# Dynamic Gap: Formal Guarantees for Safe Gap-based Navigation in Dynamic
Environments

Max Asselmeier1, Ye Zhao1 and Patricio A. Vela2 1M. Asselmeier, and Y. Zhao
are with the School of Mechanical Engineering, Georgia Institute of
Technology, Atlanta, GA 30308, USA<EMAIL_ADDRESS>P.A. Vela is with the
School of Electrical and Computer Engineering and the Institute for Robotics
and Intelligent Machines, Georgia Institute of Technology, Atlanta, GA 30308,
USA.*This work supported in part by NSF Award #1849333, #2144309 and Georgia
Tech IRIM/IPaT Aware Home Seed Grant.

###### Abstract

This paper extends the family of gap-based local planners to unknown dynamic
environments through generating provable collision-free properties for
hierarchical navigation systems. Existing perception-informed local planners
that operate in dynamic environments rely on emergent or empirical robustness
for collision avoidance as opposed to providing formal guarantees for safety.
In addition to this, the obstacle tracking that is performed in these existent
planners is often achieved with respect to a global inertial frame, subjecting
such tracking estimates to transformation errors from odometry drift. The
proposed local planner, called dynamic gap, shifts the tracking paradigm to
modeling how the free space, represented as gaps, evolves over time. Gap
crossing and closing conditions are developed to aid in determining the
feasibility of passage through gaps, and Bézier curves are used to define a
safe navigable gap that encapsulates both local environment dynamics and ego-
robot reachability. Artificial Harmonic Potential Field (AHPF) methods that
guarantee collision-free convergence to the goal are then leveraged to
generate safe local trajectories. Monte Carlo benchmarking experiments are run
in structured simulation worlds with dynamic agents to showcase the benefits
that such formal safety guarantees provide.

##  I Introduction

Collision avoidance is of utmost importance for safe robot navigation. This
task is typically handled by a local planner which utilizes sensory
information to evade obstacles. One family of local planners is the gap-based
planner [1], which identifies passable regions, or “gaps”, and synthesizes
motion commands through them. With this emphasis on free space, gap-based
planners are an approach based on the affordances of the environment [2], and
they have shown great promise with capabilities of respecting dynamic, visual,
and geometric constraints [3, 4, 5, 6, 7] as well as generating provably
collision-free trajectories [8].

Despite this success, gap-based planners have yet to be explicitly extended to
handling dynamic obstacle avoidance, a challenge that accurately reflects the
unknown, changing environments of the real world. Reactive planners may
provide collision avoidance in certain circumstances, but without formal
safety guarantees, the performance of such planners is limited. This research
deficit suggests that gap-based planners stand to benefit greatly from the
field of motion planning in dynamic environments.

Figure 1: Visualization of gaps and trajectories generated by dynamic gap. The
central blue dot represents the ego-robot while the red dots represent dynamic
agents. Bolder dots denote the current agent poses while the more transparent
dots denote future agent poses. Magenta arcs represent the instantaneous gaps
obtained during Gap Detection (Section III-A1), and magenta arrows represent
the gap state estimates (Section III-A3). Orange points represent the local
gap goals, and green lines represent the Bézier curves that define the
navigable gaps - labeled A-G in the figure - used during trajectory synthesis
(Section III-C2). Dashed blue lines running through the navigable gaps labeled
A-F represent the candidate trajectories generated at the current time step
(Section III-C4). During the Gap Feasibility Analysis step (Section III-B),
Gap G is deemed infeasible, resulting in the gap being discarded prior to
trajectory synthesis.

This paper details an extension to prior work on the potential gap planner
[8], involving augmentations to the planning framework to provide formal
guarantees for safe navigation in dynamic environments. This extension is
referred to hereafter as the dynamic gap local planner and can be visualized
in Figure 1. Egocentric free space tracking is integrated to develop
predictions of how gaps will evolve over time. These gap predictions along
with known kinematic constraints are used to ensure that the defined workspace
will remain as free space throughout the entire local time horizon. The
employment of this safe workspace in junction with collision-free trajectory
generation methods provides formal safety guarantees in ideal settings. For
this work, ideal settings involve isolated gaps (no overlap between gaps
during the local time horizon), and a first-order, point-mass, holonomic
robot.

##  II Related Work

###  II-A Motion Planning in Dynamic Environments

Methods for tracking and predicting dynamic obstacles can be categorized into
two approaches: model-based and motion-based. Model-based methods involve
classifying obstacles from a predefined set of categories such as pedestrians
or automobiles. This usually involves visual inferencing from a learned
classifier [9]. Model-based methods are advantageous in that prediction models
can be catered towards specific behaviors or classes [10], and they can
identify objects in the scene that may not currently be moving but could move
in the future. However, these approaches usually require a high level of
confidence in the object classifier.

Motion-based methods involve detecting dynamic obstacles from sensor-derived
motion cues. Gap-based planners fall into this category. These methods are
advantageous in that they will capture any moving artifact in the environment,
but they are disadvantageous in that they model all artifacts identically. To
make obstacle predictions, motion-based methods perform data clustering, data
association, and data estimation. Data clustering involves grouping input data
into individual obstacles or subsections of obstacles. Data association
involves finding pairings of clusters between consecutive time steps. The most
common method is to treat this step as a Rectangular Assignment problem and
solve via the Hungarian Algorithm [11]. Other methods may be used that are
better suited for the problem domain [12]. Data estimation involves updating
prediction models of the dynamic obstacles. Obstacles can be geometrically
modeled in various ways: polygons or lines [13], circles [14, 15], or sets of
points [16]. Prediction models can be factored into planning decisions through
capturing them in cost functions or constraints in MPC formulations [17, 13],
or embedding them in occupancy maps that can be searched over for paths [18,
19]. The methods cited here assign prediction models to individual obstacles
and represent these prediction models with respect to a fixed frame.
Performing obstacle tracking in the egocentric frame allows the local planner
to leverage the benefits of operating in the perception space.

###  II-B Perception Space and Gap-based Navigation

A perception-space approach to planning involves keeping sensory input in its
raw egocentric form to take advantage of the computational benefits that come
with foregoing intensive data processing. All local planning steps downstream
are then performed through an egocentric point of view. Prior work from the
authors has been done to develop a hierarchical perception-space navigation
stack. Specifically, efforts have been made towards efficient collision
checking [20], egocentric environment representations [21], and local path
planning [22, 8]. This hierarchical egocentric navigation stack, including the
contributions highlighted in this work, can be seen in Figure 2.

Figure 2: Information flow for egocentric hierarchical navigation stack,
including dynamic gap as the local planner. Red boxes represent input
perception data, blue boxes represent steps maintained from potential gap, and
green boxes represent novel steps that were either added or augmented.

Gap-based planners typically function by operating on 1D laser scans and
generating egocentric gaps: regions of collision-free space defined by either
leading or trailing edges of obstacles. These planners either synthesize local
trajectories or reactive control inputs through these gaps [3, 1, 7, 5, 4].
Potential Gap explored how to simplify and manipulate gap representations to
uphold collision-free gap passage guarantees for manually-crafted artificial
potential fields. To the authors’ knowledge, there are no gap-based planners
in the literature that are explicitly designed to operate within dynamic
environments.

###  II-C Artificial Harmonic Potential Fields

Artificial Potential Fields (APFs) are motion planning methods that obtain
control inputs from the gradient of a potential function that returns higher
values closer to obstacles and lower values closer to the goal [23, 24]. APFs
encapsulate both safety and goal convergence through fast and simple means,
but a major drawback of APFs is their tendency to have local minima due to
complex world geometries. Many innovations have been made to alleviate such
issues, one being Artificial Harmonic Potential Fields (AHPFs) which by design
exhibit no local minima within a specified workspace [25, 26]. This is due to
the fact that the critical points of a harmonic function strictly occur at
saddle points, meaning that the robot may momentarily exist at these points
but it will not settle at them. Recent work has explored algorithms that
provide AHPFs for safe navigation in unknown workspaces by constraining the
resultant vector field to point into the workspace at discrete points along
the workspace boundary. [27]. It has been proven that assuming sufficient
discretization, the resultant vector field will point inwards along the entire
boundary, guaranteeing collision-free convergence to the goal position [28].

##  III Dynamic Gap Local Planning Module

This section details the approach that the proposed work takes to planning
through feasible dynamic gaps. By virtue of the reachability-informed
workspace, the resultant vector field guarantees safe gap passage under ideal
conditions.

###  III-A Gap Detection, Association, and Estimation

####  III-A1 Gap Detection and Simplification

Potential gap defines a gap detection policy that operates on incoming
egocircle [22] data to detect the instantaneous static gaps in the local
environment. Gaps can either be classified as radial, in which there is a
significant instantaneous change in range between consecutive scan points, or
swept, in which there is a significant span of scan indices with maximum range
values. Gap detection involves merging radial gaps into swept gaps, the latter
of which possess favorable properties for planning such as high line-of-sight
visibility. Dynamic gap adopts these detection and simplification policies,
ultimately yielding a set of simplified gaps, $\mathcal{G}_{\rm simp}$.

####  III-A2 Gap Association

In static environments, solely gap detection is required for safe navigation,
but association is now needed to understand how gaps evolve over time.

Each gap $G\in\mathcal{G}_{\rm simp}$ is comprised of a left and right point
$p_{l},p_{r}\in\mathbb{R}^{2}$. Association is performed on the left and right
points of the simplified gaps, this set can be referred to as $P_{\rm simp}$.
The association step is represented as a rectangular assignment problem, where
the cost is equivalent to the distance between points across consecutive
steps, $P_{\rm simp}^{t-1}$ and $P_{\rm simp}^{t}$. This assignment problem is
then solved with the Hungarian Algorithm, producing a minimum distance one-to-
one mapping between $P_{\rm simp}^{t-1}$ and $P_{\rm simp}^{t}$. This
association determines which points to transition models from and to across
consecutive time steps. If the distance between two associated points exceeds
a threshold $\tau_{assoc}$, then that point-to-point association is discarded.

####  III-A3 Gap Estimation

Given the desire to perform gap estimation in the perception space, the
prediction models must represent the gap points with respect to the local
robot frame which is constantly changing. Therefore, the constant velocity
dynamics model taken with respect to the inertial frame is augmented to allow
for translations and rotations of the robot. The state vector is defined as:

$\textbf{X}=\begin{bmatrix}p_{s/b}\\\ v_{s/b}\\\
\end{bmatrix}=\begin{bmatrix}p_{s}-p_{b}\\\ v_{s}-v_{b}\end{bmatrix},$ (1)

where $p_{s/b}$ and $v_{s/b}$ represent the two-dimensional position and
velocity of the gap side $s$ (left or right) relative to the robot body $b$,
respectively. The system dynamics are then:

$\dot{\textbf{X}}=\begin{bmatrix}\dot{p}_{s/b}\\\ \dot{v}_{s/b}\\\
\end{bmatrix}=\begin{bmatrix}v_{s/b}-\omega_{b}\times p_{s/b}\\\
-a_{b}-\omega_{b}\times v_{s/b}\\\ \end{bmatrix},$ (2)

where $\omega_{b}$ represents the angular velocity of the robot and $a_{b}$
represents the linear acceleration of the robot. All values used for the
estimator are represented in the body frame.

###  III-B Gap Feasibility Analysis

Within static environments, gap feasibility solely depends on whether or not
the gap is wide enough to fit the robot. For dynamic environments, additional
considerations must be made to determine whether or not the robot can pass
through the gap before the gap shuts. By pruning away infeasible gaps, this
step not only conserves energy of the robot, but it also avoids potentially
dangerous gaps through which it would be difficult, if not impossible, for the
robot to pass.

In order to remove the ego-robot motion from the gap state and solely analyze
the motion of the gap itself, the gap-only dynamics are recovered from the
prediction models by adding the robot’s ego-velocity $v_{b}$ to the relative
velocity estimate $v_{s/b}$ to obtain the gap-only velocity $v_{s}$.

####  III-B1 Gap Categorization

Gap points can be categorized as static, expanding, or shrinking depending on
their side, left or right, and their bearing rate of change, $\dot{\beta}$:

$\text{Category}=\begin{cases}\text{Static}&\text{if $\dot{\beta}_{l}=0$ (or
$\dot{\beta}_{r}=0$)}\\\ \text{Expanding}&\text{if $\dot{\beta}_{l}>0$ (or
$\dot{\beta}_{r}<0$)},\\\ \text{Shrinking}&\text{if $\dot{\beta}_{l}<0$ (or
$\dot{\beta}_{r}>0$)}\end{cases}$ (3)

where $\dot{\beta}$ can be calculated as follows:

$\dot{\beta}_{s}=\frac{p_{s,x}v_{s,y}-p_{s,y}v_{s,x}}{p_{s,x}^{2}+p_{s,y}^{2}}.$
(4)

Gaps themselves can also be placed into the same categories depending on the
overall rate of change of the gap angle. Define the angular size of gap
$G\in\mathcal{G}_{\rm simp}$ as $\alpha$. Then:

$\text{Category}=\begin{cases}\text{Static}&\text{if $\dot{\alpha}=0$}\\\
\text{Expanding}&\text{if $\dot{\alpha}>0$}\\\ \text{Shrinking}&\text{if
$\dot{\alpha}<0$}\end{cases},$ (5)

where $\dot{\alpha}=\dot{\beta_{l}}-\dot{\beta_{r}}$. Shrinking gaps are the
only category that can be deemed infeasible given that they are the only
category in which the gap can shut, signified by the left and right bearings
crossing one another. Determining at what time this crossing occurs dictates
whether or not the gap is kinematically feasible for the robot.

####  III-B2 Gap Propagation

Gap models are integrated forward for a local time horizon $T$ to determine
the terminal gap points. Gap points can terminate prior to the end of the time
horizon in two ways: if the gap closes, or if an expanding gap point becomes
reachable. If neither of these conditions apply nor occur, the terminal points
are set to the propagated points at the end of the time horizon.

At each time step $t$ during integration, a crossing condition is checked. A
bearing crossing is interpreted as the clockwise angle of the gap switching
from convex to non-convex when both the left and right bearings pass across
the center of the gap. The clockwise angle condition is evaluated as:

$(\alpha^{t-1}_{cw}>0)\wedge(\alpha^{t}_{cw}<0),$ (6)

where

$\alpha^{t}_{cw}=\arctan(\frac{\eta^{t}_{l}\cdot\eta^{t}_{r}}{\eta^{t}_{l}\times\eta^{t}_{r}}),$
(7)

and the unit bearing vector for the gap sides at time $t$, $\eta_{s}^{t}$, is
defined as:

$\eta_{s}^{t}=\frac{\langle p_{s,x}^{t},p_{s,y}^{t}\rangle}{\lVert
p_{s}^{t}\rVert}.$ (8)

The condition for both bearings passing across the gap center is evaluated as:

$(\eta^{t}_{l}\cdot\eta^{t-1}_{c}>0)\wedge(\eta^{t}_{l}\cdot\eta^{t-1}_{c}>0),$
(9)

where $\eta_{c}^{t}$ refers to the unit bearing vector of the gap center. The
crossing condition is visualized in Figure 3.

Figure 3: Crossing condition of a gap. At time step $t-1$, the left and right
sides of the gap have yet to cross. At time step $t$, the gap sides have
crossed, as shown by the clockwise angle and gap center conditions.

In the case that the gap bearings have crossed, a further evaluation is
required to check if the gap has closed. For this, a range proximity condition
is checked:

$\lVert p_{l}^{t}-p_{r}^{t}\rVert<2r_{\rm infl},$ (10)

where $r_{\rm infl}$ is the inflated radius of the robot. If the left and
right points are sufficiently close together at the crossing, then the gap is
closed. If a gap is deemed closed, the last time step at which the robot can
fit between the gap points is used to set the terminal gap points.

The second terminal condition, for when expanding gap points first become
reachable, is used to preserve the collision-free property of the gap. By
terminating expanding gap points when they become reachable, it is ensured
that no part of the gap point’s future trajectory that is reachable will
overlap with the gap, therefore guaranteeing that the gap is a collision-free
space. This condition is evaluated as:

$\lVert p_{s}^{t}\rVert<v_{\rm max}t,$ (11)

where $v_{\rm max}$ is the maximum linear velocity of the robot. This is a
strict notion of reachability, assuming that the robot will travel in a
straight line at its maximum velocity to intersect with the gap point’s
trajectory.

####  III-B3 Gap Feasibility Testing

Once the initial and terminal gap points are set, a trajectory can be
approximated to determine if it is feasible for the robot to pass through the
gap during the local time horizon. This trajectory is estimated using a cubic
spline that is constructed from the initial and terminal conditions of the
robot and gap. If the maximum velocity along the spline does not exceed the
linear velocity constraints of the robot, then the gap is deemed feasible.
This feasibility analysis returns a set of gaps, $\mathcal{G}_{\rm feas}$, in
which candidate trajectories will be synthesized.

###  III-C Gradient Field Construction

####  III-C1 Gap Manipulation

Potential gap applies manipulation steps to the gaps in $\mathcal{G}_{\rm
simp}$ to ensure required gap properties for guaranteeing collision-free
passage such as convexity. However, the trajectory synthesis method used in
this work [27] assumes an unknown workspace and only requires that the
boundary can be defined by a set of continuously differentiable functions.
Therefore, the only manipulation step required is to inflate the gap boundary.

In static environments, it is only required to inflate a gap angularly to
account for the radius of the robot. However, in dynamic environments, it is
also required to inflate gaps outwards radially. Gaps may be attached to
moving agents, so inflating the points outwards ensures that gap passage is
synonymous with the robot moving completely beyond any agents potentially
attached to the gap. This manipulation step yields a set of gaps
$\mathcal{G}_{\rm manip}$ that are suitable for use in trajectory synthesis.

####  III-C2 Navigable Gap Construction

When planning in static environments, all free space within each gap
$G\in\mathcal{G}_{\rm manip}$ is navigable given that $G$ persists without
change for all time. However, we must now consider what subset of the gap is
reachable at each future time step. To factor in this aspect, the left and
right sides of the navigable gaps are modeled as weighted Bézier curves:

$B_{s}(u)=(1-u^{2})O+2u(1-u)w_{s}^{0}p_{s}^{0}+u^{2}p_{s}^{1},u\in[0,1],$ (12)

where $O\in\mathbb{R}^{2}$ is the origin of the Bézier curve,
$p_{s}^{0},p_{s}^{1}\in\mathbb{R}^{2}$ are the initial and terminal gap points
of $G$, and $u$ is the sliding variable used to parameterize the Bézier curve.
A weight $w^{0}_{s}$ is applied to the initial gap point and set to
$w_{s}^{0}=\lVert v_{s}\rVert/\lVert v_{\rm max}\rVert$. This weighting term
works to lessen the curvature of the Bézier curve, making the gap easier to
pass through while still respecting the gap dynamics.

####  III-C3 Gap Goal Placement

Navigable gap goals are placed within the angular span of the terminal gap
bearings. For closing gaps, goals are placed beyond the point at which the gap
points cross each other. For gaps that cross but do not close, goals are
placed between the crossing points. Otherwise, goal placement is biased
towards a navigation waypoint taken from the global plan.

####  III-C4 Artificial Harmonic Potential Field Construction

Dynamic environments can produce non-convex and unfavorable gap geometries in
which it is difficult to guarantee robustness for a manually designed gradient
field. To address this issue, we opt to leverage methods of automatically
synthesizing artificial harmonic potential fields online which provide
collision-free convergence guarantees. The involved trajectory synthesis
method will be briefly explained below, and more implementation details can be
found at [27].

Harmonic function terms are centered at $K$ discrete points along the weighted
Bézier curves with one additional term placed at the reachable gap goal, which
will be denoted as $p_{0}$. These harmonic terms take the form of

$v_{k}(p;p_{k})=\ln(\lVert p-p_{k}\rVert),$ (13)

where $p\in\mathbb{R}^{2}$ is the point within the workspace at which the term
is evaluated and $p_{k}\in\mathbb{R}^{2}$ is the point at which the term is
centered. The potential function $\Phi(p;p_{C}):\mathbb{R}^{2\times
K+1}\rightarrow\mathbb{R}^{2}$ is constructed by taking a linear weighted
combination of these harmonic terms:

$\Phi(p;p_{C})=w^{T}v(p;p_{C}),$ (14)

where $p_{C}\in\mathbb{R}^{2\times K+1}$ represents the vector of harmonic
term centers and $w\in\mathbb{R}^{K+1}$ is a vector of weights applied to the
harmonic terms. The weights are obtained by solving a Quadratic Program that
minimizes $\lVert w\rVert$ while enforcing an inward-pointing vector field at
the discrete points along the boundary. The linear velocity commands obtained
from the vector field are then

$u(p)=-||p-p_{0}||^{2}\nabla\Phi(p;p_{C}).$ (15)

An example navigable gap along with its corresponding AHPF is shown in Figure
4.

Figure 4: Example navigable gap and AHPF. The blue points correspond to the
trajectory of the left and right gap points, with the fading of the points
corresponding to the gap evolving over time.

Through this parameterization of the gap, collision-free convergence to the
gap goal is guaranteed.

Figure 5: Collision percentage of dynamic gap (DG) and dynamic TEB (DT) vs.
field of view. Single gap benchmarks are shown in subplot (a). Navigation
benchmark performance along with observed failure modes for dynamic gap and
dynamic TEB are shown in subplots (b) and (c) respectively.

####  III-C5 Proof of Collision-Free Passage

Assumptions made include an isolated, feasible gap characterized by constant
velocity points, an ideal robot (first-order, point-mass, holonomic), and
sufficient discretization of the gap boundary.

First, it will be shown that the navigable gap is comprised solely of free
space over the local time horizon. This can be done by proving each gap side
category separately. For static gap sides, this is trivial given that the gap
side does not evolve over time. For closing gap sides, the Bézier curve
divides the inside of the curve from the line segment between the two control
points that represent the gap point trajectory. Therefore, the entire gap
point trajectory lies outside of the navigable gap. For opening gap sides, if
the terminal gap point is not reachable by the robot, then no part of the gap
point trajectory is at risk of collision. If the terminal gap point is set to
the first reachable point along the gap point’s trajectory, gap inflation
prevents this point from leading to a collision. Therefore, the entire
navigable gap is collision-free during the local time horizon.

Next, it will be shown that the resultant vector field generates collision-
free trajectories. Prior work demonstrates that the employed AHPF guarantees
safety over the entire boundary of the workspace as long as the workspace
boundary is discretized with a sufficient number of points [28]. The minimum
number of points required, $N_{0}$ can be found from an inequality derived
from a Taylor Series expansion of the curves that define the workspace
boundary. This inequality is detailed in the Theorem 1 of [28]. We implemented
this inequality on the gaps used in this work and found that $N_{0}$ typically
falls on the range of $[10,15]$. In implementation, the discrete number of
points is set to $N=25$ to avoid the computation of repeatedly evaluating this
inequality, and this is assumed to be sufficient discretization.

Finally, it must also be shown for closing gaps that the robot reaches the
goal prior to the gap shutting. Since the gap is assumed to be feasible, there
exists a feasible timescale that can be applied to the synthesized trajectory
to ensure that the robot passes the gap before it closes.

###  III-D Trajectory Synthesis, Scoring, and Switching

####  III-D1 Trajectory Synthesis

In practice, the raw signal obtained from the AHPF was found to be small in
magnitude, so the resulant velocities are scaled to provide a reasonable
control signal. Forward integration is then performed on the resulting AHPFs
for the local time horizon $T$ to obtain a set of candidate trajectories.

####  III-D2 Trajectory Scoring

Trajectories are scored in a similar manner as potential gap: based on their
proximity to obstacles and their distance from the desired direction to the
goal. A trajectory cost is comprised of a cumulative egocentric pose-wise cost
based on the proximity to the current egocircle and a terminal pose cost based
on proximity to a waypoint obtained from the global plan. In order to account
for the evolution of the egocircle during the execution of the candidate
trajectory, the odometries of dynamic agents are forward propagated and the
egocircle is augmented accordingly. The highest-scoring candidate trajectory
is compared to the currently executing trajectory to determine if a trajectory
change should occur, as detailed in the next section.

####  III-D3 Trajectory Switching

The core idea behind safe hierarchical planners involves chaining together
multiple safe local trajectories. Therefore, a method of triggering a switch
to a newly synthesized local trajectory must be defined. Local trajectory
switching is performed if the currently executing trajectory is deemed unsafe
(either due to a future collision or due to the current gap being deemed
infeasible), if the currently executing trajectory has ended, or if the
prediction models that characterize the currently executing gap have been
discarded during the gap point association step. If any of these situations
arise, a switch to the highest-scoring candidate trajectory is triggered. If
there are no valid candidate trajectories to switch to, then an obstacle
avoidance policy is enacted until a candidate trajectory is synthesized.

##  IV Experiments and Outcomes

This section details the experiments performed between the proposed dynamic
gap planner and the state-of-the-art Timed-Elastic-Bands (TEB) planner [13],
which includes an extension for handling dynamic obstacles [29]. This
extension performs blob detection on the fixed frame local costmap to identify
obstacles and employs constant velocity Cartesian Kalman filters for tracking.

First, single gap benchmarks were performed in the Simple Two-Dimensional
Robot (STDR) simulator [30] to confirm the collision-free properties that
dynamic gap provides. Pairs of agents are randomly assigned a gap category
(static, shrinking, or expanding) along with a speed (0.15, 0.30, or 0.45
m/s). Agent pairs are then spawned in an empty world with the ego-robot, and
the ego-robot is to pass through the gap. This set up is meant to closely
recreate the ideal settings under which the proof of collision-free gap
passage is given. 100 trials were run for both dynamic gap and dynamic TEB,
and the results are shown in subplot (a) of Figure 5. As expected, zero
collisions were observed for dynamic gap whereas 19 were observed for dynamic
TEB.

Second, navigation benchmarks are also run in STDR. Random unknown obstacles,
both static and dynamic, are spawned throughout the world. Dynamic agents are
given random start and goal poses along with a simple global path tracking
controller. Four benchmarking worlds are utilized: the Dense, Campus, Office,
and Sector worlds. The Campus and Office worlds can be seen in Figure 6. An
agent density of 0.035 $\mathrm{agents/m^{2}}$ of free space is used in an
attempt to normalize the difficulty of navigation across the worlds. Depending
on the map size, this density yielded $16-74$ agents in total. The navigation
planners are also tested at three fields of view: $180\degree,270\degree,$ and
$360\degree$. Each setting is run with 50 unique random initial seeds in each
world. Experiments are run on an Intel Core i7 processor.

Subplots (b) and (c) in Figure 5 provide the results of dynamic gap and
dynamic TEB respectively along with a breakdown of the failure modes observed.
These modes include aborts when the planner fails to find a realizable path to
the goal, collisions with static artifacts in the environment such as walls or
scattered obstacles, and collisions with the dynamic agents. Despite these
benchmarks violating the assumptions for collision-free gap passage, the
underlying framework of dynamic gap still leads the planner to outperform TEB
by roughly 30% across the three fields of view tested. TEB leverages a soft-
constraint time-optimal approach to local path planning. Therefore, no
guarantees for collision avoidance can be made. This factor contributes
significantly to the high collision rates observed for TEB. A majority of the
collisions observed across both planners were due to collisions with the
dynamic agents. A small number of the dynamic gap runs resulted in collisions
with the static environment, but no runs result in an abort. Dynamic TEB shows
the opposite trend, with some runs ending in an abort, but no runs resulting
in static collisions.

Figure 6: Campus (left) and Office (right) World in STDR simulator, with
dynamic agents (red) and ego-robot (blue).

##  V Conclusion

This paper proposes an extension to an existent egocentric hierarchical
navigation stack which involves adapting local gap-based planners to dynamic
environments. This is accomplished through tracking the free space evolution
relative to the egocentric frame, constructing collision-free navigable gaps
based on gap dynamics and robot reachability, and synthesizing collision-free
trajectories through Artificial Harmonic Potential Field methods. Navigation
benchmarks showcase that the proposed work outperforms a state-of-the-art
local planner that is also extended to dynamic environments. This code is
open-sourced as a ROS package [31].

This work investigates methods of ensuring formal guarantees for safe
navigation regarding trajectory generation for gap-based planners in dynamic
environments. Future work will involve extending this framework to
additionally provide formal guarantees for online trajectory execution as well
as non-ideal robot dynamics such as nonholonomic or legged systems through the
utilization of barrier functions.

## References

  * [1] M. Mujahad, D. Fischer, B. Mertsching, and H. Jaddu, “Closest gap based (cg) reactive obstacle avoidance navigation for highly cluttered environments,” in _2010 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2010, pp. 1805–1812.
  * [2] J. J. Gibson, _The Ecological Approach to Visual Perception_. Boston: Houghton Mifflin, 1979.
  * [3] V. Sezer and M. Gokasan, “A novel obstacle avoidance algorithm: “follow the gap method”,” _Robotics and Autonomous Systems_ , vol. 60, no. 9, pp. 1123–1134, 2012.
  * [4] M. Mujahed and B. Mertsching, “The admissible gap (ag) method for reactive collision avoidance,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , 2017, pp. 1916–1921.
  * [5] M. Mujahed, D. Fischer, and B. Mertsching, “Safe gap based (sg) reactive navigation for mobile robots,” in _2013 European Conference on Mobile Robots_ , 2013, pp. 325–330.
  * [6] M. Mujahed and B. Mertsching, “A new gap-based collision avoidance method for mobile robots,” in _2016 IEEE International Symposium on Safety, Security, and Rescue Robotics (SSRR)_ , 2016, pp. 220–226.
  * [7] M. Mujahed, H. Jaddu, D. Fischer, and B. Mertsching, “Tangential closest gap based (tcg) reactive obstacle avoidance navigation for cluttered environments,” in _2013 IEEE International Symposium on Safety, Security, and Rescue Robotics (SSRR)_ , 2013, pp. 1–6.
  * [8] R. Xu, S. Feng, and P. A. Vela, “Potential gap: A gap-informed reactive policy for safe hierarchical navigation,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 8325–8332, 2021.
  * [9] T. Salzmann, B. Ivanovic, P. Chakravarty, and M. Pavone, _Trajectron++: Dynamically-Feasible Trajectory Forecasting with Heterogeneous Data_ , 12 2020, pp. 683–700.
  * [10] L. Zhao and C. Thorpe, “Qualitative and quantitative car tracking from a range image sequence,” in _Proceedings. 1998 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (Cat. No.98CB36231)_ , 1998, pp. 496–501.
  * [11] H. W. Kuhn, “The hungarian method for the assignment problem,” _Naval Research Logistics Quarterly_ , vol. 2, no. 1-2, pp. 83–97, 1955.
  * [12] D. Z. Wang, I. Posner, and P. Newman, “Model-free detection and tracking of dynamic objects with 2d lidar,” _The International Journal of Robotics Research_ , vol. 34, no. 7, pp. 1039–1063, 2015.
  * [13] C. Rösmann, F. Hoffmann, and T. Bertram, “Timed-elastic-bands for time-optimal point-to-point nonlinear model predictive control,” in _2015 European Control Conference (ECC)_ , 2015, pp. 3352–3357.
  * [14] P. Fiorini and Z. Shiller, “Motion planning in dynamic environments using the relative velocity paradigm,” in _[1993] Proceedings IEEE International Conference on Robotics and Automation_ , 1993, pp. 560–565 vol.1.
  * [15] J. van den Berg, M. Lin, and D. Manocha, “Reciprocal velocity obstacles for real-time multi-agent navigation,” in _2008 IEEE International Conference on Robotics and Automation_ , 2008, pp. 1928–1935.
  * [16] S. Vaskov, S. Kousik, H. Larson, F. Bu, J. Ward, S. Worrall, M. Johnson-Roberson, and R. Vasudevan, “Towards provably not-at-fault control of autonomous robots in arbitrary dynamic environments,” 02 2019.
  * [17] M. Gaertner, M. Bjelonic, F. Farshidian, and M. Hutter, “Collision-free mpc for legged robots in static and dynamic scenes,” 2021-04, p. 1550.
  * [18] R. Senanayake, L. Ott, S. O' Callaghan, and F. T. Ramos, “Spatio-temporal hilbert maps for continuous occupancy representation in dynamic environments,” in _Advances in Neural Information Processing Systems_ , D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, Eds., vol. 29. Curran Associates, Inc., 2016\.
  * [19] V. Guizilini, R. Senanayake, and F. Ramos, “Dynamic hilbert maps: Real-time occupancy predictions in changing environments,” in _2019 IEEE International Conference on Robotics and Automation (ICRA)_ , 2019, pp. 4091–4097.
  * [20] J. S. Smith and P. Vela, “Pips: Planning in perception space,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , 2017, pp. 6204–6209.
  * [21] J. S. Smith, S. Feng, F. Lyu, and P. A. Vela, _Real-Time Egocentric Navigation Using 3D Sensing_. Cham: Springer International Publishing, 2020, pp. 431–484.
  * [22] J. S. Smith, R. Xu, and P. Vela, “egoteb: Egocentric, perception space navigation using timed-elastic-bands,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 2703–2709.
  * [23] D. Koditschek, “Exact robot navigation by means of potential functions: Some topological considerations,” in _Proceedings. 1987 IEEE International Conference on Robotics and Automation_ , vol. 4, 1987, pp. 1–6.
  * [24] O. Khatib, “Real-time obstacle avoidance for manipulators and mobile robots,” in _Proceedings. 1985 IEEE International Conference on Robotics and Automation_ , vol. 2, 1985, pp. 500–505.
  * [25] J.-O. Kim and P. Khosla, “Real-time obstacle avoidance using harmonic potential functions,” _IEEE Transactions on Robotics and Automation_ , vol. 8, no. 3, pp. 338–349, 1992.
  * [26] S. G. Loizou, “Closed form navigation functions based on harmonic potentials,” in _2011 50th IEEE Conference on Decision and Control and European Control Conference_ , 2011, pp. 6361–6366.
  * [27] P. Rousseas, C. P. Bechlioulis, and K. J. Kyriakopoulos, “Trajectory planning in unknown 2d workspaces: A smooth, reactive, harmonics-based approach,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 1992–1999, 2022\.
  * [28] P. Rousseas, C. Bechlioulis, and K. J. Kyriakopoulos, “Harmonic-based optimal motion planning in constrained workspaces using reinforcement learning,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 2005–2011, 2021\.
  * [29] C. Rösmann, F. Hoffmann, and T. Bertram, “Track and include dynamic obstacles via costmap_converter,” 2018. [Online]. Available: http://wiki.ros.org/teb_local_planner
  * [30] A. T. M. Tsardoulias and C. Zalidis, “stdr_simulator - ros_wiki,” 2014. [Online]. Available: http://wiki.ros.org/stdr_simulator
  * [31] M. Asselmeier, Y. Zhao, and P. A. Vela, “Dynamic gap.” [Online]. Available: https://github.com/ivaROS/DynamicGap


# A canonical polyadic tensor basis for fast Bayesian estimation of multi-
subject fMRI activation patterns

Michelle F. Miranda
Department of Mathematics and Statistics
University of Victoria
Victoria, BC, Canada
<EMAIL_ADDRESS>

###### Abstract

Task-evoked functional magnetic resonance imaging studies, such as the Human
Connectome Project (HCP), are a powerful tool for exploring how brain activity
is influenced by cognitive tasks like memory retention, decision-making, and
language processing. A fast Bayesian function-on-scalar model is proposed for
estimating population-level activation maps linked to the working memory task.
The model is based on the canonical polyadic (CP) tensor decomposition of
coefficient maps obtained for each subject. This decomposition effectively
yields a tensor basis capable of extracting both common features and subject-
specific features from the coefficient maps. These subject-specific features,
in turn, are modeled as a function of covariates of interest using a Bayesian
model that accounts for the correlation of the CP-extracted features. The
dimensionality reduction achieved with the tensor basis allows for a fast MCMC
estimation of population-level activation maps. This model is applied to one
hundred unrelated subjects from the HCP dataset, yielding significant insights
into brain signatures associated with working memory.

_K_ eywords tensor decomposition $\cdot$ neuroimaging $\cdot$ fMRI $\cdot$ CP
decomposition $\cdot$ Bayesian modeling $\cdot$ functional regression model

##  1 Introduction

Task-evoked functional magnetic resonance imaging (fMRI) studies, such as the
Human Connectome Project (HCP), are a powerful tool for exploring how brain
activity is influenced by cognitive tasks such as memory retention, decision-
making, and language processing. In the HCP, a core set of tasks was devised
aiming to delineate the functionally unique nodes within the human brain to
establish connections between tasks and activation signatures in critical
network nodes [2]. These data have motivated the work presented here.

Brain imaging data is often represented in the form of a multi-dimensional
array, called a tensor. These arrays can be seen as functions of space and/or
time. In multi-subject studies, it is of interest to relate these functional
response data to a set of scalar predictors through a function-on-scalar
model. Due to the high-dimensionality of the functional response, most models
in the literature will employ a form of regularization. Regularization
involves borrowing strength across observations within a function. Although
regularization can be achieved through various methods such as penalized
regression methods or/and sparse-inducing priors, one very important approach
is by the use of basis functions. Each basis function defines a linear
combination among the locations within the function, establishing a framework
for borrowing strength between observations. Commonly used basis functions are
splines, Fourier series, wavelets, and principal components (PCs)[18].

Some efforts have been devoted to analyzing neuroimaging as tensor data. In
one line of research, the CP decomposition is used as a descriptive tool to
represent multisubject and multisession fMRI data. The resulting spatial
patterns or modes of the tensor representation are, in turn, validated by
visual inspection. The extracted components in each spatial, temporal, runs,
among other dimensions, can be further inspected to describe variability in
each of these directions [1]. Although this descriptive decomposition is
useful, it is often of interest to relate the tensor neuroimaging data with
other variables such as age, gender, and sex, disease status, through a
regression analysis framework.

Previous research that tackles tensor regression are models in which the
response is a scalar and the predictor is a tensor [27, 17, 12, 7]. Commonly
in these approaches, the goal is to understand the changes of the response
outcome as the tensor variable changes. Often, a tensor decomposition is
implemented as tool for dimensionality reduction of tensor covariates. For
example, the authors in [17] use this approach to predict Alzheimer’s Disease
status as a function of magnetic resonance imaging data through a CP
decomposition of the imaging covariate.

Meanwhile, tensor response regression aims to study the change of the image as
the predictors such as the disease status and age vary. Research in this
direction often propose a low-rank tensor decomposition of the coefficients
associated with covariates [25, 23]. In a different application context, [14]
propose a likelihood-based latent variable tensor response regression, where
the tensor data corresponds to latent variables informed by additional
covariates. The authors apply their model to a 2D facial image database with
facial descriptor.

Although aiming to address the same type of problem, our proposal is
methodologically distinct from these previous approaches. Here, we use the CP
tensor decomposition to construct a tensor basis function that not only
determines the mechanisms for borrowing information across voxels but also
effectively reduces the dimensionality of the brain data, which is projected
onto the basis space. The projected data is in turn modeled as functions of
covariates of interest using a Bayesian MCMC algorithm. This approach is
inspired by the wavelet functional mixed models proposed in [20] and further
developed in [19, 26, 28, 29]. The novelty of the methodology here is that the
CP decomposition captures the 3D spatial brain intricacies across different
subjects and concatenation of the spatial directions does not happen as in
wavelet basis or principal components. We apply this approach to a second-
level analyses of fMRI data from the one hundred unrelated subjects of the HCP
data. The method, however, extends to a vast range of imaging modalities such
as structural MRI images, and electroencephalography data.

The rest of the paper is organized as follows. In Section 2 we introduce the
CP decomposition and estimation method, including a brief review of tensor
concepts and properties (Section 2.1) and a review of the method used to
obtain subject-level brain maps for fMRI task data (Section 2.2, and the
design the CP tensor basis functions with the Bayesian modeling approach
(Section 2.3). In Section 3 we apply the proposed approach to the Working
Memory data from HCP to find population-level maps of a particular Working
Memory contrast. In Section 4 we provide some concluding remarks and
directions for future work.

##  2 Methods

The proposed approach involves a two-stage procedure. In the first stage, the
4D data for each subject are modeled through the Bayesian composite-hybrid
basis model proposed in [16] and described in Section 2.2. In the second
stage, posterior means of contrast brain maps are modeled by the CP tensor
basis model. This modeling strategy is flexible, as subject-level maps can be
estimated by any known single-subject estimation technique and combined into
coefficients maps for any contrasts of interest. Furthermore, these maps are
assembled into a 4D tensor (3D volume and subjects) to be associated with
population-level covariates, e.g., age and sex.

###  2.1 Background

We review a few important concepts and operations necessary for the
comprehension of this manuscript. For an exhaustive review on tensor
decomposition and applications, see [10].

####  2.1.1 Working with tensors

A tensor is a multidimensional array with order given by the number of its
dimensions, also known as modes. For example, a matrix is a tensor of order 2.
A tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\ldots\times I_{K}}$ of
order $K$ is said to be a rank-one tensor if it can be written as an outer
product of $K$ vectors, i.e.,

$\mathcal{X}={\mbox{\boldmath${a}$}}^{(1)}\circ{\mbox{\boldmath${a}$}}^{(2)}\circ\ldots\circ{\mbox{\boldmath${a}$}}^{(K)},$
(1)

where $\circ$ is the outer product of vectors. This is equivalent to writing
each element of $x_{i_{1}i_{2}\ldots
i_{K}}=a^{(1)}_{i_{1}}a^{(2)}_{i_{2}}\ldots a^{(K)}_{i_{K}}$, for $1\leq
i_{k}\leq i_{K}$.

Matricization or unfolding of a tensor is the process of reordering the tensor
into a matrix. We will denote the mode-$n$ matricization of the tensor
$\mathcal{X}$ above as ${\mbox{\boldmath${X}$}}_{(n)}$ and arrange the
mode-$n$ dimension to be the columns of ${\mbox{\boldmath${X}$}}_{(n)}$.
Please see [10] for how the observations are arranged for a general tensor
order. As an example, consider the matricization of a tensor
$\mathcal{X}\in\mathbb{R}^{2\times 3\times 4}$. The mode-$3$ matricization
${\mbox{\boldmath${X}$}}_{(3)}$ is a matrix of size $6\times 4$ where, for a
fixed column $k$ in the matrix, we have elements
$x_{11k},x_{21k},x_{x12k},x_{22k},x_{13k},x_{23k}$ in the $k$th row.

Tensor decompositions are higher-order generalizations of the singular value
decomposition and principal component analysis. Although there are various
tensor decomposition methods, we focus here on the Canonical Polyadic (CP).
The CP decomposition factorizes a tensor as a sum of rank-one tensors, as
illustrated in Figure 1 for a tensor of order four.

Figure 1: CP decomposition of a 4D array into the sum of rank-one tensors.

Mathematically, we write the CP decomposition above as

$\mathcal{X}\approx\llbracket{\mbox{\boldmath${\lambda}$}};{\mbox{\boldmath${A}$}},{\mbox{\boldmath${B}$}},{\mbox{\boldmath${C}$}},{\mbox{\boldmath${D}$}}\rrbracket\equiv\sum_{r=1}^{R}\lambda_{r}{\mbox{\boldmath${a}$}}_{r}\circ{\mbox{\boldmath${b}$}}_{r}\circ{\mbox{\boldmath${c}$}}_{r}\circ{\mbox{\boldmath${d}$}}_{r}$
(2)

The factorization means that every element of the four-dimensional array
$\mathcal{Y}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{3}\times J_{4}}$ can be
written as
$y_{j_{1}j_{2}j_{3}j_{4}}=\sum_{r=1}^{R}\lambda_{r}a_{j_{1}r}b_{j_{2}r}c_{j_{3}r}d_{j_{4}r}$,
where the vectors ${\mbox{\boldmath${a}$}}_{r}\in\mathbb{R}^{J_{1}}$,
${\mbox{\boldmath${b}$}}_{r}\in\mathbb{R}^{J_{2}}$,
${\mbox{\boldmath${c}$}}_{r}\in\mathbb{R}^{J_{3}}$, and
${\mbox{\boldmath${d}$}}_{r}\in\mathbb{R}^{J_{4}}$ form the columns of the
factor matrices
${\mbox{\boldmath${A}$}},{\mbox{\boldmath${B}$}},{\mbox{\boldmath${C}$}},$ and
${D}$.

The factor matrices refer to the combination of the vectors from the rank-one
components, i.e.,
${\mbox{\boldmath${A}$}}=[{\mbox{\boldmath${a}$}}_{1}{\mbox{\boldmath${a}$}}_{2}\ldots{\mbox{\boldmath${a}$}}_{R}]$
and likewise for B, C, and D. We assume that the columns of A, B, C, and D are
normalized to length one, with the weights absorbed into the vector
${\mbox{\boldmath${\lambda}$}}\in\mathbb{R}^{R}$. The approximation refers to
the estimation of the factor matrices that minimize the Frobenius norm of the
difference between $\mathcal{X}$ and
$\llbracket{\mbox{\boldmath${\lambda}$}};{\mbox{\boldmath${A}$}},{\mbox{\boldmath${B}$}},{\mbox{\boldmath${C}$}},{\mbox{\boldmath${D}$}}\rrbracket$
[4].

###  2.2 Subject-level Modeling

For each voxel $v=1,\ldots,N_{v}$ in a volumetric image, the task fMRI
response function is a time series of $T$ time points based on the Blood
Oxygen Level Dependent (BOLD) contrast [21]. These measurements are an
indirect and non-invasive measure of brain activity and can be modeled as a
function of the external stimuli:

$y_{t}(v)=\sum_{j=1}^{p}b_{j}(v)x_{j}(t)+e_{t}(v),$ (3)

where $x_{j}(t)$ is the covariate for stimulus $j$, representing the indicator
function of stimulus $j$ at time $t$ convolved with the hemodynamic response
of a neural event at time $t$ [13], and $e_{t}(v)$ is a measurement error. The
coefficients $b_{j}(v)$ characterize the relationship of brain activity during
stimulus $j$. The primary goal of activation studies is to identify which
voxels are differentially activated by specified stimuli.

There are many ways to estimate the coefficient maps $b_{j}(v)$, including the
popular generalized linear model available through the software SPM [6]. The
model in (3) can be written in matrix form as:

$\mathbb{Y}=\mathbb{X}{\mbox{\boldmath${B}$}}+\mathbb{E}$ (4)

To estimate the brain maps ${\mbox{\boldmath${B}$}}\in\mathbb{R}^{p\times
N_{v}}$, we will use the Bayesian approach proposed in [16], which consists of
fitting a basis-space model:

$\mathbb{Y^{*}}=\mathbb{X^{*}}{\mbox{\boldmath${B}$}}^{*}+\mathbb{E}^{*}$ (5)

This can be written in vectorized form as:

$\text{vec}(\mathbb{Y}^{*})=({\mbox{\boldmath${I}$}}_{S}\otimes\mathbb{X}^{*})({\mbox{\boldmath${I}$}}_{S}\otimes{\mbox{\boldmath${B}$}}^{*})+\text{vec}(\mathbb{E}^{*}),$
(6)

with
$\text{vec}(\mathbb{Y}^{*})=({\mbox{\boldmath${\Psi}$}}^{T}{\mbox{\boldmath${\Phi}$}}^{T}\otimes{\mbox{\boldmath${W}$}})\text{vec}({\mathbb{Y}})$
and
$\text{vec}(\mathbb{E}^{*})=({\mbox{\boldmath${\Psi}$}}^{T}{\mbox{\boldmath${\Phi}$}}^{T}\otimes{\mbox{\boldmath${W}$}})\text{vec}(\mathbb{E})$.
In summary, the authors’ tensor basis approach consists of choosing the set of
bases
$\\{{\mbox{\boldmath${\Psi}$}},{\mbox{\boldmath${\Phi}$}},{\mbox{\boldmath${W}$}}\\}$
to account for global and local spatial correlations and long-memory temporal
correlations. A Bayesian MCMC is proposed and samples from the coefficient
maps are obtained in the basis space, and the inverse transform is applied to
obtain the MCMC samples in the voxel space. For more details, refer to [16].

###  2.3 Multi-subject Modeling

Let $\mathcal{Y}$ be the 4D tensor representing each coefficient map or
contrast map for all subjects. Then $\mathcal{Y}\in\mathbb{R}^{p_{1}\times
p_{2}\times p_{3}\times N}$. We can model $\mathcal{Y}$ as a function of $p$
subject-level covariates ${\mbox{\boldmath${Z}$}}\in\mathbb{R}^{p}$ as
follows:

${\mbox{\boldmath${Y}$}}={\mbox{\boldmath${Z}$}}{\mbox{\boldmath${\gamma}$}}+{\mbox{\boldmath${E}$}},$
(7)

where ${Y}$ is a $N\times N_{v}$ matrix representing the mode four unfolding
of the tensor $\mathcal{Y}$; ${Z}$ is a $N\times p$ design matrix of subject-
level covariates; ${\gamma}$ is a $p\times N_{v}$ matrix of coefficients; and
${E}$ is a $N\times N_{v}$ matrix of random errors.

Properly modeling the data as in (7) is challenging due to its very large
number of voxels, typically around a million, and the spatial dependency
inherited from the brain volume data. Most Bayesian models rely on either
downsampling the volume data to reduce the number of voxels or rely on
simplified assumptions on the prior of ${\gamma}$ and/or ${E}$. The solution
proposed here involves projecting the data into a basis space through a tensor
basis obtained from the CP decomposition of the tensor $\mathcal{Y}$. This
decomposition is low rank and effectively extracts brain volume features that
are common to all subjects. By extracting the common features, the remaining
information is subject-specific and can be then related to the covariates of
interest.

The CP decomposition of $\mathcal{Y}$ is given by:

$\mathcal{Y}\approx\sum_{r=1}^{R}\lambda_{r}{\mbox{\boldmath${a}$}}^{(1)}_{r}\circ{\mbox{\boldmath${a}$}}^{(2)}_{r}\circ{\mbox{\boldmath${a}$}}^{(3)}_{r}\circ{\mbox{\boldmath${g}$}}_{r}\equiv\llbracket{\mbox{\boldmath${\lambda}$}};{\mbox{\boldmath${A}$}}^{(1)},{\mbox{\boldmath${A}$}}^{(2)},{\mbox{\boldmath${A}$}}^{(3)},{\mbox{\boldmath${G}$}}\rrbracket,$
(8)

where the vectors ${\mbox{\boldmath${a}$}}^{(1)}_{r}\in\mathbb{R}^{p_{1}}$,
${\mbox{\boldmath${a}$}}^{(2)}_{r}\in\mathbb{R}^{p_{2}}$,
${\mbox{\boldmath${a}$}}^{(3)}_{r}\in\mathbb{R}^{p_{3}}$, and
${\mbox{\boldmath${g}$}}_{r}\in\mathbb{R}^{N}$. The matrices
${\mbox{\boldmath${A}$}}^{(1)}$, ${\mbox{\boldmath${A}$}}^{(2)}$,
${\mbox{\boldmath${A}$}}^{(3)}$, and ${G}$ are called factor matrices and
their columns are associated with these vectors. From [10], we can write ${Y}$
from model (7) as:

${\mbox{\boldmath${Y}$}}={\mbox{\boldmath${G}$}}{\mbox{\boldmath${\Lambda}$}}({\mbox{\boldmath${A}$}}^{(3)}\odot{\mbox{\boldmath${A}$}}^{(2)}\odot{\mbox{\boldmath${A}$}}^{(1)})^{T}+{\mbox{\boldmath${\epsilon}$}},$
(9)

where ${G}$, ${\mbox{\boldmath${A}$}}^{(1)}$, ${\mbox{\boldmath${A}$}}^{(2)}$,
and ${\mbox{\boldmath${A}$}}^{(3)}$ are estimated to minimize the Frobenius
norm of ${\epsilon}$ [10]. Let
${\mbox{\boldmath${L}$}}={\mbox{\boldmath${\Lambda}$}}({\mbox{\boldmath${A}$}}^{(3)}\odot{\mbox{\boldmath${A}$}}^{(2)}\odot{\mbox{\boldmath${A}$}}^{(1)})^{T}\in\mathbb{R}^{N_{v}\times
R}$, then ${\mbox{\boldmath${P}$}}={\mbox{\boldmath${L}$}}^{\dagger}$, the
Moore-Penrose pseudoinverse of ${L}$, provides a projection matrix formed by
spatial basis that can be used to effectively reduce data from $N_{v}$ to $R$
spatial locations.

Based on the projection matrix, we multiply both sides of equation (7) by
${P}$ to obtain the basis space model:

${\mbox{\boldmath${G}$}}={\mbox{\boldmath${Z}$}}{\mbox{\boldmath${\gamma}$}}^{*}+{\mbox{\boldmath${E}$}}^{*},$
(10)

where
${\mbox{\boldmath${\gamma}$}}^{*}={\mbox{\boldmath${\gamma}$}}{\mbox{\boldmath${P}$}}$
is the $p\times R$ matrix of coefficients in the basis space, and
${\mbox{\boldmath${E}$}}^{*}={\mbox{\boldmath${E}$}}{\mbox{\boldmath${P}$}}$
is a $N\times R$ matrix of random errors.

####  2.3.1 Basis-space Bayesian Modeling

The columns of ${G}$ are correlated due to the nature of the CP decomposition,
and model assumptions should reflect that. For each row of
${\mbox{\boldmath${E}$}}^{*}$, we assume ${\mbox{\boldmath${E}$}}_{i}^{*}\sim
N(0,{\mbox{\boldmath${\Sigma}$}}_{\epsilon})$ and priors:

${\mbox{\boldmath${p}$}}({\mbox{\boldmath${\gamma}$}}^{*},{\mbox{\boldmath${\Sigma}$}}_{\epsilon})={\mbox{\boldmath${p}$}}({\mbox{\boldmath${\Sigma}$}}_{\epsilon}){\mbox{\boldmath${p}$}}({\mbox{\boldmath${\gamma}$}}^{*}/{\mbox{\boldmath${\Sigma}$}}_{\epsilon}),$
(11)

with
${\mbox{\boldmath${\gamma}$}}^{*}\sim\text{MN}({\mbox{\boldmath${g}$}}_{0},{\mbox{\boldmath${L}$}}_{0}^{-1},{\mbox{\boldmath${\Sigma}$}}_{\epsilon})$
and
${\mbox{\boldmath${\Sigma}$}}_{\epsilon}\sim\text{W}^{-1}({\mbox{\boldmath${V}$}}_{0},\nu_{0})$,
where MN indicates the matrix normal distribution and $\text{W}^{-1}$
indicates the Inverse Wishart distribution. The joint posterior distribution
of
${\mbox{\boldmath${\gamma}$}}^{*},{\mbox{\boldmath${\Sigma}$}}_{\epsilon}/{\mbox{\boldmath${G}$}},{\mbox{\boldmath${Z}$}}$
factorizes as a product of an Inverse Wishart and a Matrix Normal
distributions as follows:

${\mbox{\boldmath${\Sigma}$}}_{\epsilon}/{\mbox{\boldmath${G}$}},{\mbox{\boldmath${Z}$}}\sim\text{W}^{-1}({\mbox{\boldmath${V}$}}_{n},\nu_{n}),$
${\mbox{\boldmath${\gamma}$}}^{*}/{\mbox{\boldmath${G}$}},{\mbox{\boldmath${Z}$}},{\mbox{\boldmath${\Sigma}$}}_{\epsilon}\sim\text{MN}({\mbox{\boldmath${g}$}}_{n},{\mbox{\boldmath${L}$}}_{n}^{-1},{\mbox{\boldmath${\Sigma}$}}_{\epsilon}),$
(12)

with:

${\mbox{\boldmath${V}$}}_{n}={\mbox{\boldmath${V}$}}_{0}+({\mbox{\boldmath${Y}$}}-{\mbox{\boldmath${Z}$}}{\mbox{\boldmath${g}$}}_{n})^{T}({\mbox{\boldmath${Y}$}}-{\mbox{\boldmath${Z}$}}{\mbox{\boldmath${g}$}}_{n})+({\mbox{\boldmath${g}$}}_{n}-{\mbox{\boldmath${g}$}}_{0})^{T}{\mbox{\boldmath${L}$}}_{0}({\mbox{\boldmath${g}$}}_{n}-{\mbox{\boldmath${g}$}}_{0}),$
$\nu_{n}=\nu_{0}+N,$
${\mbox{\boldmath${g}$}}_{n}=({\mbox{\boldmath${Z}$}}^{T}{\mbox{\boldmath${Z}$}}+{\mbox{\boldmath${L}$}}_{0})^{-1}({\mbox{\boldmath${Z}$}}^{T}{\mbox{\boldmath${G}$}}+{\mbox{\boldmath${L}$}}_{0}{\mbox{\boldmath${g}$}}_{0}),$
${\mbox{\boldmath${L}$}}_{n}={\mbox{\boldmath${Z}$}}^{T}{\mbox{\boldmath${Z}$}}+{\mbox{\boldmath${L}$}}_{0}.$

Sampling from these distributions is straightforward and done through Gibbs
sampling.

####  2.3.2 Rank Selection

To choose the appropriate rank for the CP decomposition, we propose a 10-fold
cross-validation procedure with the following steps for each of the
training/test set combination:

  1. 1.

Obtain the factor matrices of the CP decomposition for the training data,

  2. 2.

Proceed to estimate ${\mbox{\boldmath${\gamma}$}}^{*}$ from model (10)
following the MCMC algorithm detailed in Section 2.3.1,

  3. 3.

Using $\widehat{{\mbox{\boldmath${\gamma}$}}^{*}}$ from the previous step and
${\mbox{\boldmath${Z}$}}_{\text{Test}}$, estimate
${\mbox{\boldmath${G}$}}_{\text{Test}}={\mbox{\boldmath${Z}$}}_{\text{Test}}\widehat{{\mbox{\boldmath${\gamma}$}}^{*}}$,

  4. 4.

Obtain
$\widehat{\mathcal{Y}}_{\text{Test}}=\llbracket{\mbox{\boldmath${\lambda}$}};{\mbox{\boldmath${A}$}}^{(1)},{\mbox{\boldmath${A}$}}^{(2)},{\mbox{\boldmath${A}$}}^{(3)},{\mbox{\boldmath${G}$}}_{\text{Test}}\rrbracket$
as in (8),

  5. 5.

Choose the rank that minimizes the Frobenius norm of the average of the
$\mathcal{Y}_{\text{Test}}-\widehat{\mathcal{Y}}_{\text{Test}}$, obtained
across the training/test sets.

####  2.3.3 Inference on the Brain Space

For every MCMC sample of ${\mbox{\boldmath${\gamma}$}}^{*}$, we obtain
${\mbox{\boldmath${\gamma}$}}={\mbox{\boldmath${\gamma}$}}^{*}{\mbox{\boldmath${L}$}}$
in the original voxel space. To control for the experimentwise error rate in
the voxel space, joint credible bands as in [24] can be calculated. For each
voxel $v$, a joint $100(1-\alpha)\%$ credible interval for any group-level
contrast ${C}$ is given by

$I_{\alpha}(v)=\hat{C}(v)\pm
q_{(1-\alpha)}[\widehat{\mbox{std}}(\hat{C}(v))],$ (13)

where the quantile $q_{(1-\alpha)}$ is obtained from $z^{(1)},\ldots,z^{(M)}$
with
$z^{(m)}=\max_{v\in\mathcal{V}}\left|\frac{C^{(m)}(v)-\hat{C}(v)}{\widehat{\mbox{std}}(\hat{C}(v))}\right|$,
for $m=1,\ldots,M$ and $M$ is the total number of MCMC samples after burn-in
and thinning [16].

These joint credible bands account for multiple testing across the voxel space
and can be inverted to flag activation signatures based on voxels for which
the $100(1-\alpha)\%$ joint credible bands exclude zero [15, 16]. Denote by
$P_{SimBas}$ the minimum $\alpha$ at which each interval excludes zero
($P_{SimBas}=\min\\{\alpha:0\notin I_{\alpha}(v)\\}$). $P_{SimBas}$ can be
directly computed by

$P_{SimBas}(v)=\frac{1}{M}\sum_{m=1}^{M}1\left\\{\left|\frac{\hat{C}(v)}{\widehat{\mbox{std}}(\hat{C}(v))}\right|\leq
z^{(m)}\right\\}.$ (14)

##  3 Application to the HCP data

The proposed model is applied to the 100 unrelated subjects data from the
Working Memory task of the Human Connectome Project. The volumes considered
were collected from the right-left phase. The experiment consists of a total
of 8 blocks, half corresponding to the 2-back task and the other half the
0-back task. In each block, the participants were shown a series of 10 images.
Images corresponding to 4 stimuli (tools, places, body parts, faces) were
embedded within the memory task. Each image was shown for 2.5s, followed by a
15s inter-block interval. fMRI volumes were obtained every 720 ms (TR). Each
volume consisted of images of size $91\times 109\times 91$ for a total of 405
time frames.

Subject-level analysis. We follow the analysis described in Section 2.2. To
obtain the local spatial basis, we partitioned the brain into ROIs using a
digitalized version of the original Talairach structural labelling that was
registered into the MNI152 space. The atlas can be obtained from the FSL atlas
library with information found at
https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Atlases [8]. We considered ROIs that
had at least 125 voxels and obtained a total of 298 regions. Within each
region, the local basis ${\Phi}$ was formed by the eigenvectors of the SVD
decomposition for the matrix data in each ROI as described in [16]. Next, the
global spatial basis ${\Psi}$ was obtained by first projecting the brain data
using the local basis and then calculating the SVD decomposition of the
project data, exactly as in [16]. The temporal basis ${W}$ was chosen to be
the Haar wavelet. The estimation was performed at the cluster Cedar managed by
the Digital Research Alliance of Canada (former Compute Canada). The jobs were
submitted in parallel for the 100 subjects and results took an hour to finish.

For each subject, after estimation, we constructed a contrast map of the 2
back versus 0 back working memory task. The maps were obtained from the
posterior means of the coefficients, and assembled into a 4D array
$\mathcal{Y}\in\mathbb{R}^{91\times 109\times 91\times 405}$.

Multi-subject analysis. We take the tensor $\mathcal{Y}$ containing the
posterior means of the contrast 2back versus 0back for the maps obtained in
the single-level analysis. Next, we follow the method proposed in Section 2.3
to estimate the coefficients maps ${\gamma}$ for the mean of 2back versus
0back and the covariate Sex. We are interested in finding common activation
maps for all subjects and aim to investigate if there are differences in
working memory task based on sex. Next, we obtain the CP decomposition of the
tensor $\mathcal{Y}$, taking $R=30$ selected based on the cross-validation
procedure described in Section 2.3.2. Figure 2 shows the Frobenius norm as a
function of the rank, for rank values $R=10,15,20,25,\ldots,80$. Following the
decomposition, we aim to identify significant locations in the brain for the
Working Memory task 2back versus 0 back, we proceed as follows.

  1. (i)

Obtain MCMC samples for ${\mbox{\boldmath${\gamma}$}}^{*}$, as detailed in
Section 2.3.1.

  2. (ii)

Use the inverse transform of ${P}$,
${\mbox{\boldmath${L}$}}={\mbox{\boldmath${\Lambda}$}}({\mbox{\boldmath${A}$}}^{(1)}\odot{\mbox{\boldmath${A}$}}^{(2)}\odot{\mbox{\boldmath${A}$}}^{(3)})^{T}\in\mathbb{R}^{N_{v}\times
R}$, to obtain MCMC samples for
${\mbox{\boldmath${\gamma}$}}={\mbox{\boldmath${\gamma}$}}^{*}{\mbox{\boldmath${L}$}}$
in the voxel space.

  3. (iii)

For each subject-level covariate, use the MCMC samples of ${\gamma}$ to obtain
$P_{SimBas}$ as described in (14) to control for the experimentwise error rate
on multiple testing.

  4. (iv)

For $P_{SimBas}<\alpha=0.01$ flag the voxel as different than zero, indicating
activation.

  5. (v)

Obtain clusters of significant voxels that are at least 125 voxels large.

The computational time to run this analysis is 1.62 minutes on a MacBook Pro
2.8 GHz Quad-Core Intel Core i7 with 16 GB 2133 MHz LPDDR3 memory. For the
cross-validation procedure, it took 38 minutes for a grid of 15 rank values.
For faster computational time on rank selection, one could run a five-fold
cross validation instead and/or a less fine grid of ranks.

Figure 2: Cross-validation to choose the rank R. We want R to minimize the
Frobenius norm.

We consider two maps of coefficients, the first is the intercept, which is the
average activation for the contrast 2 back versus 0 back of the Working Memory
task. The second is the map associated with the variable sex. We did not find
any significant differences in the working memory contrast between females and
males. Figure 3 shows the results for the intercept revealing 4445 active
voxels. The larger cluster in red shows 1268 voxels encompassing 24 distinct
brain regions as subdivided by the Tailarach Atlas, (labels can be found here
http://www.talairach.org/about.html#Labels). Within this cluster, larger
active regions are in the right parietal lobe, including inferior, superior,
precuneous parietal regions. The cluster in green is formed by 1182 voxels,
with active regions on the left parietal lobe, including the inferior,
superior, precuneous, and cuneous regions. The parietal lobe has been
previously identified as playing a key role in working memory [3, 11, 9].

The next cluster in blue shows the right frontal lobe, including the superior
and middle frontal gyrus. The cluster in yellow shows 345 voxels in the right
frontal lobe, including the superior and middle frontal gyrus. The cluster in
light orange shows 276 voxels on the left frontal lobe in the middle frontal
gyrus. These regions have also been identified as playing a critical role in
working memory [22]. The cluster in burgundy indicates 313 voxels on the right
cerebellum at the posterior lobe region. A smaller cluster of 139 voxels was
also identified on the posterior lobe at the left cerebellum. The authors in
[5] found that the cerebellum can influence the accurate maintenance of
working memory.

Figure 3: Working Memory contrast of 2back versus 0back.

Next, we compare the results obtained from the Bayesian CP tensor basis
approach with the results obtained from the Statistical Parametric Map (SPM)
software [6]. We ran both the first and the second-level analysis in SPM with
parameter settings to reflect no spatial normalization and an AR(1) noise
component. The results for the Average of 2 back versus 0 back contrast are
displayed in Figure 4. The two largest clusters were sizes 64 and 58 voxels,
respectively. Overlaid onto the template is a mask with all significant voxels
based on a voxel-level threshold that controls for the family-wise error rate
of 0.05. Inspecting Figure 4 reveals that this SPM voxelwise model only
capture small parts of the parietal and frontal lobe and has significant less
power than the CP tensor basis approach.

Figure 4: Working Memory contrast of 2back versus 0back obtained from SPM.

##  4 Final Remarks

In this paper, we propose a fast Bayesian model that takes tensor responses as
a function of population-level covariates. The model is based on using the CP
tensor decomposition approach to construct a basis function with the goal of
modeling the spatial structure of the tensor data. This approach allows us to
borrow information across voxels, increasing the power for the statistical
tests performed. Because the estimation of the modeling components is
performed in the basis space, we have a very fast Bayesian procedure. The
Bayesian MCMC estimation method for the second-level analysis, including
cross-validation for rank selection, estimation, and inference, takes 40
minutes on a regular laptop. This time can be much faster if using a cluster
with more powerful computing capabilities.

We apply this approach to a second-level analysis of fMRI data from one
hundred unrelated subjects of the HCP dataset. The results obtained are
consistent with findings in the literature showing the role of the parietal
lobe, the middle frontal gyrus, and the posterior lobe of the cerebellum in
working memory.

One limitation of the proposed approach is that CP decomposition applies only
to volumetric data. For more modern imaging files, such as the CIFTI files
available in the Human Connectome Project, the CP decomposition would need to
aggregate spatial information as one large dimension. In this case, the
benefits of the CP decomposition are unclear as the CP decomposition might be
similar to other basis approaches, including principal components and
wavelets. The authors plan to conduct an investigation on that front.

## Funding

This research was funded by Natural Sciences and Engineering Research Council
of Canada grant number RGPIN-2020-06941.

## References

  * [1] Anders H Andersen and William S Rayens.  Structure-seeking multilinear methods for the analysis of fmri data.  NeuroImage, 22(2):728–739, 2004.
  * [2] Deanna M. Barch, Gregory C. Burgess, Michael P. Harms, Steven E. Petersen, Bradley L. Schlaggar, Maurizio Corbetta, Matthew F. Glasser, Sandra Curtiss, Sachin Dixit, Cindy Feldt, Dan Nolan, Edward Bryant, Tucker Hartley, Owen Footer, James M. Bjork, Russ Poldrack, Steve Smith, Heidi Johansen-Berg, Abraham Z. Snyder, and David C. Van Essen.  Function in the human connectome: Task-fmri and individual differences in behavior.  NeuroImage, 80:169–189, 2013.  Mapping the Connectome.
  * [3] Marian E. Berryhill and Ingrid R. Olson.  The right parietal lobe is critical for visual working memory.  Neuropsychologia, 46(7):1767–1774, 2008.  Part Special Issue: What is the Parietal Lobe Contribution to Human Memory?
  * [4] J Douglas Carroll and Jae-On Chang.  Analysis of individual differences in multidimensional scaling via an n-way generalization of “eckart-young” decomposition.  Psychometrika, 35(3):283–319, 1970.
  * [5] Ben Deverett, Mikhail Kislin, David W. Tank, et al.  Cerebellar disruption impairs working memory during evidence accumulation.  Nature Communications, 10:3128, 2019.
  * [6] K. J. Friston, William D. Penny, John Ashburner, Stefan J. Kiebel, and Thomas E. Nichols.  Statistical Parametric Mapping: The Analysis of Functional Brain Images.  Academic Press, London, 2006.
  * [7] Rajarshi Guhaniyogi, Shaan Qamar, and David B. Dunson.  Bayesian tensor regression.  Journal of Machine Learning Research, 18(79):1–31, 2017.
  * [8] Mark Jenkinson, Christian F. Beckmann, Timothy E.J. Behrens, Mark W. Woolrich, and Stephen M. Smith.  Fsl.  NeuroImage, 62(2):782–790, 2012.  20 YEARS OF fMRI.
  * [9] Michael Koenigs, Aron K. Barbey, Bradley R. Postle, and Jordan Grafman.  Superior parietal cortex is critical for the manipulation of information in working memory.  Journal of Neuroscience, 29(47):14980–14986, 2009.
  * [10] Tamara G. Kolda and Brett W. Bader.  Tensor decompositions and applications.  SIAM Review, 51(3):455–500, 2009.
  * [11] Bohee Lee, June Sic Kim, and Chun Kee Chung.  Parietal and medial temporal lobe interactions in working memory goal-directed behavior.  Cortex, 150:126–136, 2022.
  * [12] Xiaoshan Li, Da Xu, Hua Zhou, and Lexin Li.  Tucker tensor regression and neuroimaging analysis.  Statistics in Biosciences, 10(3):520–545, December 2018.
  * [13] Martin A. Lindquist.  The statistical analysis of fmri data.  Statist. Sci., 23(4):439–464, 11 2008.
  * [14] Eric F. Lock and Gen Li.  Supervised multiway factorization.  Electronic Journal of Statistics, 12(1):1150 – 1180, 2018.
  * [15] Mark J. Meyer, Brent A. Coull, Francesco Versace, Paul Cinciripini, and Jeffrey S. Morris.  Bayesian function-on-function regression for multilevel functional data.  Biometrics, 71(3):563–574, 2015.
  * [16] Michelle F Miranda and Jeffrey S Morris.  Novel bayesian method for simultaneous detection of activation signatures and background connectivity for task fmri data.  arXiv preprint arXiv:2109.00160, 2021.
  * [17] Michelle F. Miranda, Hongtu Zhu, and Joseph G. Ibrahim.  TPRM: Tensor partition regression models with applications in imaging biomarker detection.  The Annals of Applied Statistics, 12(3):1422 – 1450, 2018.
  * [18] Jeffrey S. Morris.  Functional regression.  Annual Review of Statistics and Its Application, 2(1):321–359, 2015.
  * [19] Jeffrey S. Morris, Veerabhadran Baladandayuthapani, Richard C. Herrick, Pietro Sanna, and Howard Gutstein.  Automated analysis of quantitative image data using isomorphic functional mixed models, with application to proteomics data.  Annals of Applied Statistics, 5:605–1125, 2011.
  * [20] Jeffrey S. Morris and Raymond J. Carroll.  Wavelet-based functional mixed models.  Journal of the Royal Statistical Society, Series B, 68:179–199, 2006.
  * [21] S Ogawa, T M Lee, A R Kay, and D W Tank.  Brain magnetic resonance imaging with contrast dependent on blood oxygenation.  Proceedings of the National Academy of Sciences, 87(24):9868–9872, 1990.
  * [22] Adrian M. Owen, Chantal E. Stern, Rodney B. Look, Irene Tracey, Bruce R. Rosen, and Michael Petrides.  Functional organization of spatial and nonspatial working memory processing within the human lateral frontal cortex.  Proceedings of the National Academy of Sciences, 95(13):7721–7726, 1998.
  * [23] Guillaume Rabusseau and Hachem Kadri.  Low-rank regression with tensor responses.  In Advances in Neural Information Processing Systems, volume 29, pages 1867–1875, 2016.
  * [24] David Ruppert, M. P. Wand, and R. J. Carroll.  Semiparametric Regression.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, 2003.
  * [25] Daniel Spencer, Rajarshi Guhaniyogi, and Raquel Prado.  Joint bayesian estimation of voxel activation and inter-regional connectivity in fmri experiments.  Psychometrika, 85(4):845–869, December 2020.
  * [26] Lin Zhang, Veerabhadran Baladandayuthapani, Hongxiao Zhu, Keith A. Baggerly, Tadeusz Majewski, Bogdan A. Czerniak, and Jeffrey S. Morris.  Functional car models for large spatially correlated functional datasets.  Journal of the American Statistical Association, 111(514):772–786, 2016.
  * [27] Hua Zhou, Lexin Li, and Hongtu Zhu.  Tensor regression with applications in neuroimaging data analysis.  Journal of the American Statistical Association, 108(502):540–552, 2013.  PMID: 24791032.
  * [28] Hongxiao Zhu, Philip J. Brown, and Jeffrey S. Morris.  Robust, adaptive functional regression in functional mixed model framework.  Journal of the American Statistical Association, 106(495):1167–1179, 2011.
  * [29] Hongxiao Zhu, Francesco Versace, Paul M. Cinciripini, Philip Rausch, and Jeffrey S. Morris.  Robust and gaussian spatial functional regression models for analysis of event-related potentials.  NeuroImage, 181:501–512, 2018.