# Covariance-Free Sparse Bayesian Learning

Alexander Lin, Andrew H. Song,
Berkin Bilgic, and Demba Ba A. Lin and D. Ba are with the School of
Engineering and Applied Sciences, Harvard University, Cambridge, MA, 02138 USA
(e-mail<EMAIL_ADDRESS>demba@seas.harvard.edu). A. H. Song is with
the Electrical Engineering and Computer Science, Massachusetts Institute of
Technology, Cambridge, MA, 02138 USA (email: andrew90@mit.edu).B. Bilgic is
with Harvard-MIT Health Sciences and Technology, Massachusetts Institute of
Technology, Cambridge, MA, USA, Athinoula A. Martinos Center for Biomedical
Imaging, Charlestown, MA, USA Department of Radiology, Harvard Medical School,
Boston, MA, USA.

###### Abstract

Sparse Bayesian learning (SBL) is a powerful framework for tackling the sparse
coding problem while also providing uncertainty quantification. The most
popular inference algorithms for SBL exhibit prohibitively large computational
costs for high-dimensional problems due to the need to maintain a large
covariance matrix. To resolve this issue, we introduce a new method for
accelerating SBL inference – named covariance-free expectation maximization
(CoFEM) – that avoids explicit computation of the covariance matrix. CoFEM
solves multiple linear systems to obtain unbiased estimates of the posterior
statistics needed by SBL. This is accomplished by exploiting innovations from
numerical linear algebra such as preconditioned conjugate gradient and a
little-known diagonal estimation rule. For a large class of compressed sensing
matrices, we provide theoretical justifications for why our method scales well
in high-dimensional settings. Through simulations, we show that CoFEM can be
up to thousands of times faster than existing baselines without sacrificing
coding accuracy. Through applications to calcium imaging deconvolution and
multi-contrast MRI reconstruction, we show that CoFEM enables SBL to tractably
tackle high-dimensional sparse coding problems of practical interest.

##  I Introduction

Sparse coding is a fundamental problem in signal processing that seeks a
sparse solution $\boldsymbol{z}$ to the equation

$\displaystyle\boldsymbol{y}=\boldsymbol{\Phi}\boldsymbol{z}+\boldsymbol{\varepsilon},$
(1)

where $\boldsymbol{z}\in\mathbb{R}^{D}$ is a latent sparse vector,
$\boldsymbol{y}\in\mathbb{R}^{N}$ is an observed measurement vector,
$\boldsymbol{\Phi}\in\mathbb{R}^{N\times D}$ is a known dictionary, and
$\boldsymbol{\varepsilon}\in\mathbb{R}^{N}$ is an unknown noise vector.

_Sparse Bayesian learning_ (SBL) is an effective methodology for sparse
coding. It has been employed in several popular models, such as sparse
Bayesian regression [1], relevance vector machines [2], and Bayesian
compressed sensing [3, 4, 5]. The practical applications of SBL are numerous,
encompassing diverse examples such as medical image reconstruction [6, 7],
direction of arrival estimation [8], human pose estimation [9], structural
health monitoring [10], seismic exploration [11], and visual tracking [12].

SBL offers several advantages compared to other common approaches to sparse
coding (e.g. $\ell_{0}$ or $\ell_{1}$ regularization). As a Bayesian method,
SBL provides uncertainty quantification and the ability to recover credible
intervals over $\boldsymbol{z}$ instead of a single point solution. Moreover,
SBL does not need to tune regularization parameters since it can learn them or
integrate them out using hyperpriors [4]. As a generative model, SBL can be
embedded as a submodule within a larger framework to enforce more complex
structure for $\boldsymbol{z}$ (e.g. group sparsity [13], block sparsity
[14]). Finally, SBL has favorable optimization properties, such as a sparser
global minimum than $\ell_{1}$ methods and fewer local minima than $\ell_{0}$
methods [15].

One often-noted limitation of SBL is the heavy computational cost of its
inference algorithm [2, 6]. On the one hand, the fact that SBL requires more
computation than non-Bayesian approaches should be unsurprising, since SBL
recovers an entire distribution instead of a single point estimate. On the
other hand, the most widely-used options for SBL inference scale poorly to
high-dimensional problems, which are becoming increasingly common in many
domains. This limitation threatens to render the SBL paradigm obsolete for
large-scale settings, as inference cannot be performed in an acceptable
timeframe for practical applications.

The principal inference procedure of SBL is the expectation-maximization (EM)
algorithm [1, 2]. Each iteration of EM is computationally expensive, requiring
$O(D^{3})$-time and $O(D^{2})$-space to invert a large covariance matrix,
where $D$ is the dimension of the sparse codes. There have been several
attempts to reduce the costs of EM, using iteratively reweighted least-squares
(IRLS) [16], approximate message passing (AMP) [17], and variational inference
(VI) [18]. A popular alternative to EM called the sequential algorithm (Seq)
[19] is often faster in practice. However, these methods lack either
_scalability_ at very high dimensions $D$ or _accuracy_ in recovering ground-
truth sparse codes when compared to EM.

###  I-A Contributions

We introduce a novel method for accelerating SBL’s EM algorithm for high-
dimensional problems. We call our method covariance-free expectation
maximization (CoFEM) because it eliminates the main bottleneck of EM – i.e.
the storage and inversion of the covariance matrix. We demonstrate that CoFEM
has significant advantages in both scalability and accuracy over other SBL
approaches (reviewed in Section II). Our contributions are categorized into
four types – methodology (Section III), theory (Section IV), experimentation
(Section V), and applications (Section VI).

_Methodology:_ We develop CoFEM, which solves multiple linear systems and
exploits a little-known diagonal estimation rule to obtain _unbiased_
estimates of the posterior moments needed by EM. The multiple systems can be
solved in parallel using an iterative solver, such as the conjugate gradient
(CG) algorithm, without constructing the covariance matrix. This simple, yet
powerful principle reduces EM’s per-iteration time complexity from $O(D^{3})$
to $O(UK\tau_{D})$, where $\tau_{D}$ is the time needed for matrix-vector
multiplication and $U,K\ll D$ are hyperparameters for the number of CG steps
and the number of linear systems. Furthermore, CoFEM reduces EM’s space
complexity from $O(D^{2})$ to $O(DK)$.

_Theory:_ We prove new theorems that further characterize CoFEM’s asymptotic
complexities. In applying this theory to matrices $\boldsymbol{\Phi}$ that
satisfy the restricted isometry property (RIP) and commonly occur in
compressed sensing applications, we show that, in the limit of a large number
of EM iterations, CoFEM’s hyperparameters $U$ and $K$ can be kept small even
as the dimensionality of the problem $D$ grows very large. We use this theory
to justify practical choices for $U$ and $K$ that need not increase with $D$,
leading to signficant computational savings compared to existing SBL
algorithms.

_Experimentation:_ We perform simulations comparing CoFEM to several existing
SBL algorithms [2], [16], [17], [18], [19], for high-dimensional sparse signal
recovery. Across all of our settings, CoFEM is able to maintain the same level
of accuracy as EM due to its unbiased estimation of posterior moments. In
addition, CoFEM’s highly-parallelized and space-saving design enables further
acceleration via graphics processing units (GPUs); most other approaches
cannot fully realize this benefit due to their heavy memory requirements. In
practice, CoFEM with GPU acceleration can reduce hours of computation for
covariance-based algorithms to a few seconds.

_Applications:_ We apply CoFEM-equipped SBL in two settings of practical
interest: (1) calcium deconvolution and (2) multi-contrast MRI reconstruction.
In these high-dimensional settings, CoFEM attains competitive computation time
with non-Bayesian approaches, while providing several benefits (e.g. superior
performance, uncertainty quantification). These applications require
extensions of the SBL model to multi-task learning and to settings with non-
negativity constraints, which we demonstrate CoFEM is flexible enough to
handle.

##  II Sparse Bayesian Learning

###  II-A Generative Model

To solve Eq. (1), SBL imposes the following model:

$\displaystyle\boldsymbol{z}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{0},\mathrm{diag}\\{\boldsymbol{\alpha}\\}^{-1}),$
$\displaystyle\boldsymbol{y}\>|\>\boldsymbol{z}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{\Phi}\boldsymbol{z},1/\beta\mathbf{I}),$
(2)

where $\beta\in\mathbb{R}$ is the inverse of the variance of the noise
(commonly called precision) and $\mathbf{I}$ is the identity matrix. Given
$\boldsymbol{y}$, SBL performs inference on this model to recover
$\boldsymbol{z}$.

The main identifying feature of SBL is the diagonal Gaussian prior with
precision parameters $\boldsymbol{\alpha}\in\mathbb{R}^{D}$ for
$\boldsymbol{z}$. The notation $\mathrm{diag}\\{\boldsymbol{\alpha}\\}$ in Eq.
(2) maps $\boldsymbol{\alpha}$ to a $D\times D$ matrix with
$\boldsymbol{\alpha}$ along its diagonal and zero elsewhere. SBL performs type
II maximum likelihood estimation [20] by integrating out $\boldsymbol{z}$ and
optimizing $\boldsymbol{\alpha}$. Thus, SBL learns a posterior distribution
with uncertainty over $\boldsymbol{z}$. The overall learning objective is:

$\displaystyle\max_{\boldsymbol{\alpha}}$ $\displaystyle\log
p(\boldsymbol{y}\>|\>\boldsymbol{\alpha})=\log\int_{\boldsymbol{z}}p(\boldsymbol{y}\>|\>\boldsymbol{z})p(\boldsymbol{z}\>|\>\boldsymbol{\alpha})d\boldsymbol{z}.$
(3)

As this objective is optimized, many of the elements of $\boldsymbol{\alpha}$
diverge to $\infty$. For these elements, the independent Gaussian priors
converge to point masses at zero, forcing their respective posteriors to
follow suit. Thus, upon convergence of $\boldsymbol{\alpha}$ to
$\boldsymbol{\hat{\alpha}}$, the recovered posterior distribution
$p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\hat{\alpha}})$ is often
highly sparse. This phenomenon is called _automatic relevance determination_
[15] because SBL learns which elements of $\boldsymbol{z}$ are “relevant”
(i.e. non-zero) from the data.

In the remainder of this section, we review several inference schemes that
have been proposed to optimize Eq. (3) and comment on their respective
shortcomings.

###  II-B Expectation-Maximization Algorithm

Most SBL inference schemes are built on the _expectation-maximization_ (EM)
algorithm, a general framework for parameter estimation in the presence of
latent variables [21, 2]. EM iteratively alternates between an expectation
step (E-Step) and a maximization step (M-Step). Let
$\boldsymbol{\alpha}^{(t)}$ be the solution at the start of the $t$-th
iteration. The E-Step computes the expectation
$Q(\boldsymbol{\alpha};\boldsymbol{\alpha}^{(t)})$ of the complete data log-
likelihood $\log p(\boldsymbol{z},\boldsymbol{y}\>|\>\boldsymbol{\alpha})$
with respect to the latent posterior
$p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\alpha}^{(t)})$:

$\displaystyle
Q(\boldsymbol{\alpha};\boldsymbol{\alpha}^{(t)})=\mathbb{E}_{p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\alpha}^{(t)})}[\log
p(\boldsymbol{z},\boldsymbol{y}\>|\>\boldsymbol{\alpha})]$
$\displaystyle=\mathbb{E}_{p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\alpha}^{(t)})}[\log
p(\boldsymbol{z}\>|\>\boldsymbol{\alpha})+\log
p(\boldsymbol{y}\>|\>\boldsymbol{z})]$
$\displaystyle\propto\sum_{j=1}^{D}\log\alpha_{j}-\alpha_{j}\cdot\mathbb{E}_{p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\alpha}^{(t)})}\left[z_{j}^{2}\right]+\text{const},$
(4)

where “const” absorbs all terms that are constant with respect to
$\boldsymbol{\alpha}$. The posterior
$p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\alpha}^{(t)})$ is a
multivariate Gaussian distribution
$\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ with mean and covariance
parameters

$\displaystyle\boldsymbol{\mu}=\beta\boldsymbol{\Sigma}\boldsymbol{\Phi}^{\top}\boldsymbol{y},$
$\displaystyle\boldsymbol{\Sigma}=(\beta\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\mathrm{diag}\\{\boldsymbol{\alpha}^{(t)}\\})^{-1}.$
(5)

The second moment of each $z_{j}$ in Eq. (4) can be decomposed into a sum over
the squared mean and variance, i.e.

$\displaystyle\mathbb{E}_{p(\boldsymbol{z}|\boldsymbol{y},\boldsymbol{\alpha}^{(t)})}[z_{j}^{2}]$
$\displaystyle=\mathbb{E}_{p(\boldsymbol{z}|\boldsymbol{y},\boldsymbol{\alpha}^{(t)})}[z_{j}]^{2}+\mathbb{V}\mathrm{ar}_{p(\boldsymbol{z}|\boldsymbol{y},\boldsymbol{\alpha}^{(t)})}[z_{j}]$
$\displaystyle=\mu_{j}^{2}+\Sigma_{j,j},$ (6)

The M-Step updates each ${\alpha}^{(t)}_{j}$ by maximizing Eq. (4) with
respect to $\alpha_{j}$. Given $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$,
this can be done in closed-form by differentiating $Q$:

$\displaystyle\frac{\partial
Q}{\partial\alpha_{j}}=\frac{1}{\alpha_{j}}-{\mu_{j}^{2}+\Sigma_{j,j}}=0\implies{\alpha}_{j}^{(t+1)}=\frac{1}{\mu_{j}^{2}+\Sigma_{j,j}},$
(7)

EM repeats Eqs. (5) and (7) for $T$ iterations until convergence, while
guaranteeing non-negative change in the log-likelihood objective of Eq. (3) at
each step.

Despite its simplicity, the EM algorithm is limited by its computational cost.
The E-Step of Eq. (5) is expensive for large $D$. Storing
$\boldsymbol{\Sigma}$ requires $O(D^{2})$-space and computing it through
matrix inversion requires $O(D^{3})$-time. This makes the standard EM
algorithm challenging at high dimensions.

###  II-C Previous Attempts to Accelerate EM

There have been several attempts to accelerate the EM algorithm for SBL
inference. One approach is based on _iteratively reweighted least squares_
(IRLS) [16], which applies the Woodbury matrix identity to
$\boldsymbol{\Sigma}$ in Eq. (5), yielding

$\displaystyle\boldsymbol{\Sigma}=\mathbb{C}-\mathbb{C}\mathbb{\Phi}^{\top}(\mathbb{I}+\mathbb{\Phi}\mathbb{C}\mathbb{\Phi}^{\top})^{-1}\mathbb{\Phi}\mathbb{C},$
(8)

where $\mathbb{C}=\mathrm{diag}\\{\boldsymbol{\alpha}^{(t)}\\}^{-1}$. By
inverting an $N\times N$ matrix (as opposed to a $D\times D$ matrix), IRLS
reduces the per-iteration time complexity of EM to $O(N^{3}+N^{2}D)$. For
problems in which $N\ll D$, IRLS is more efficient than EM. However, if
$N=O(D)$, the time complexity is still $O(D^{3})$. Furthermore, like EM, IRLS
requires storage of $\boldsymbol{\Sigma}$, which remains an expensive
$O(D^{2})$-space cost.

Another line of work is based on variants of _approximate message passing_
(AMP). AMP combines quadratic and Taylor series approximations with loopy
belief propagation to iteratively estimate the means and variances of
$\boldsymbol{z}$ in Eq. (6), circumventing matrix inversion [17, 22]. Though
faster than EM in practice, AMP is known to diverge for dictionaries
$\boldsymbol{\Phi}$ that do not satisfy zero-mean, sub-Gaussian criteria [23].

A third common strategy is to employ _variational inference_ (VI), which
approximates the true posterior
$p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\alpha})$ with a simpler
surrogate $q(\boldsymbol{z})$ [24, 25, 26]. This allows for SBL inference that
is inverse-free [18, 27]. However, VI approaches optimize a lower bound on Eq.
(3) instead of the true log-likelihood objective. Thus, they may converge to a
sub-optimal solution for $\boldsymbol{\alpha}$. Both AMP and VI are ultimately
limited by the fact that their approximations to the means and variances of
$\boldsymbol{z}$ can be _biased_ for general dictionaries $\boldsymbol{\Phi}$
[23, 27]. In Section III, we present a new method that ensures an _unbiased_
estimation of these moments, regardless of the structure of
$\boldsymbol{\Phi}$.

###  II-D Sequential Algorithm

The sequential (Seq) algorithm, is a notable alternative to EM that reduces
computation time and space in practice [19, 3]. Seq maintains a set
$\mathcal{S}\subseteq\\{1,2,\ldots,D\\}$ of “active” indices such that
$\alpha_{j}$ is finite for each $j\in\mathcal{S}$ and $\alpha_{j}=\infty$ for
all $j\neq\mathcal{S}$. Initially, $\mathcal{S}=\emptyset$. Indices are
sequentially added to or deleted from $\mathcal{S}$ if making such a change
can increase the log-likelihood objective (Eq. (3)).

At any given point, Seq only needs to store parts of the mean vector
$\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$ corresponding
to $\mathcal{S}$; all other components are assumed to be zero. Thus, for truly
sparse vectors $\boldsymbol{z}$ with $d$ non-zero components such that $d\ll
D$, Seq is more efficient than EM. Yet unlike EM, the number of iterations
needed for Seq depends on $d$, since at least $d$ steps must be taken to fully
recover $\boldsymbol{z}$. The overall time complexity is $O(d^{2}D)$ and the
space complexity is $O(d^{2}+D)$ [4].

However, Seq has several limitations. It is often the case that $d$ is a
fraction or percentage of $D$. If $d$ grows linearly with $D$, the asymptotic
time cost is $O(D^{3})$, similar to EM. This may explain why Seq can still be
up to hundreds of times slower than non-Bayesian methods for large $D$ [6].
Also, the algorithm’s sequential nature hinders its potential for speedup on
parallel machines. Lastly, for large $D$, it remains costly to store parts of
the quadratically-sized covariance matrix $\boldsymbol{\Sigma}$.

##  III Covariance-Free Expectation-Maximization

We introduce _covariance-free expectation-maximization_ (CoFEM), a new SBL
inference scheme for accelerating EM that removes the need to compute (let
alone invert) the covariance matrix $\boldsymbol{\Sigma}$. Our method is built
on several advances from numerical linear algebra [28, 29].

Our key observation is that not all elements of $\boldsymbol{\Sigma}$ are
needed for the M-Step in Eq. (7). Indeed, we only need $\boldsymbol{\mu}$
(which depends on $\boldsymbol{\Sigma}$ via Eq. (5)) and the the diagonal
elements of $\boldsymbol{\Sigma}$ to update $\boldsymbol{\hat{\alpha}}$. Thus,
we propose a simplified E-Step that can estimate the two desired quantities
$\\{\mu_{j},\Sigma_{j,j}\\}$ for all $j$ from _solutions to linear systems_ ,
thereby avoiding the need for matrix inversion. We can re-express Eq. (5) for
$\boldsymbol{\mu}$ as

$\displaystyle\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}=\beta\boldsymbol{\Phi}^{\top}\boldsymbol{y},$
(9)

where
$\boldsymbol{\Sigma}^{-1}=\beta\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\mathrm{diag}\\{\boldsymbol{\alpha}^{(t)}\\}$.
Thus, $\boldsymbol{\mu}$ is the solution $\boldsymbol{x}$ to the linear system
$\mathbf{A}\boldsymbol{x}=\boldsymbol{b}$ for
$\mathbf{A}:=\boldsymbol{\Sigma}^{-1}$ and
$\boldsymbol{b}:=\beta\boldsymbol{\Phi}^{\top}\boldsymbol{y}$.

###  III-A Diagonal Estimation

We leverage a result from [28] to obtain the diagonal of
$\boldsymbol{\Sigma}$.

######  Proposition 1 (Diagonal Estimation Rule).

Let $\mathbf{M}\in\mathbb{R}^{D\times D}$ be a square matrix. Let
$\boldsymbol{p}_{1},\boldsymbol{p}_{2},\ldots,\boldsymbol{p}_{K}\in\mathbb{R}^{D}$
be random probe vectors, where each $\boldsymbol{p}_{k}$ comprises independent
and identically distributed elements such that $\mathbb{E}[p_{k,j}]=0$ for all
$j=1,\ldots,D$. For each $\boldsymbol{p}_{k}$, let
$\boldsymbol{x}_{k}=\mathbb{M}\boldsymbol{p}_{k}$. Consider the estimator
$\boldsymbol{s}\in\mathbb{R}^{D}$, where, for each $j=1,\ldots,D$,

$\displaystyle s_{j}=\frac{\sum_{k=1}^{K}p_{k,j}\cdot
x_{k,j}}{\sum_{k=1}^{K}p_{k,j}^{2}}.$ (10)

Then, each $s_{j}$ is an unbiased estimator of $\mathrm{M}_{j,j}$.

###### Proof.

Consider $s_{j}$, the $j$-th element of $\boldsymbol{s}$. We have

$\displaystyle s_{j}$
$\displaystyle=\frac{\sum_{k=1}^{K}\left(p_{k,j}\cdot\left(\sum_{j^{\prime}=1}^{D}\mathrm{M}_{j,j^{\prime}}\cdot
p_{k,j^{\prime}}\right)\right)}{\sum_{k=1}^{K}p_{k,j}^{2}}$
$\displaystyle=\mathrm{M}_{j,j}+\sum_{j^{\prime}\neq
j}\mathrm{M}_{j,j^{\prime}}\cdot\frac{\sum_{k=1}^{K}p_{k,j}\cdot
p_{k,j^{\prime}}}{\sum_{k=1}^{K}p_{k,j}^{2}}.$ (11)

Thus, $\mathbb{E}[s_{j}]$ is equal to the following:

$\displaystyle\mathrm{M}_{j,j}+\sum_{j^{\prime}\neq
j}\mathrm{M}_{j,j^{\prime}}\cdot\left(\sum_{k=1}^{K}\underbrace{\mathbb{E}[p_{k,j^{\prime}}]}_{=0}\cdot\mathbb{E}\left[\frac{p_{k,j}}{\sum_{k=1}^{K}p_{k,j}^{2}}\right]\right),$
(12)

where we have applied the fact that the $j$ and $j^{\prime}$ components of
$\boldsymbol{p}_{k}$ are independent to arrive at a product of expectations.
Since $\mathbb{E}[p_{k,j^{\prime}}]=0$ for all $k$ and $j^{\prime}$, it
follows that $\boldsymbol{\mathbb{E}}[s_{j}]=\mathrm{M}_{j,j}$. ∎

We apply this rule to $\boldsymbol{\Sigma}$ to estimate its diagonal elements.
Following [28], we use the _Rademacher distribution_ for drawing the probe
vector $\boldsymbol{p}_{k}$, where each $p_{k,j}$ is either $-1$ or $+1$ with
equal probability. This simplifies Eq. (10) to

$\displaystyle s_{j}=\frac{1}{K}\sum_{k=1}^{K}p_{k,j}\cdot x_{k,j}.$ (13)

The diagonal estimation rule turns an _explicit_ property of a matrix (i.e.
the diagonal elements) into an _implicit_ quantity that can be estimated
without physically forming the matrix. To exploit this rule, we only need a
method for applying $\boldsymbol{\Sigma}$ to each vector $\boldsymbol{p}_{k}$
to obtain $\boldsymbol{x}_{k}$; this can be done by solving a linear system
$\mathbb{A}\boldsymbol{x}=\boldsymbol{b}$, where
$\mathbb{A}:=\boldsymbol{\Sigma}^{-1}$ and
$\boldsymbol{b}:=\boldsymbol{p}_{k}$.

In summary, the quantities $\\{\mu_{j},\Sigma_{j,j}\\}$ needed for the
simplified E-Step update can be obtained by solving $K+1$ separate linear
systems. These systems can be solved in parallel by considering the matrix
equation $\mathbb{A}\mathbb{X}=\mathbb{B}$ with inputs
$\mathbb{A}\in\mathbb{R}^{D\times D}$ and
$\mathbb{B}\in\mathbb{R}^{D\times(K+1)}$ defined as follows:

$\displaystyle\mathbb{A}$
$\displaystyle:=\beta\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\mathrm{diag}\\{\boldsymbol{\alpha}^{(t)}\\},$
$\displaystyle\mathbb{B}$
$\displaystyle:=\begin{bmatrix}\boldsymbol{p}_{1}\>|\>\boldsymbol{p}_{2}\>|\>\ldots\>|\>\boldsymbol{p}_{K}\>|\>\beta\boldsymbol{\Phi}^{\top}\boldsymbol{y}\end{bmatrix}.$
(14)

If we enumerate the columns of the solution matrix
$\mathbb{X}\in\mathbb{R}^{D\times(K+1)}$ as
$\boldsymbol{x}_{1},\boldsymbol{x}_{2},\ldots,\boldsymbol{x}_{K},\boldsymbol{\mu}$,
then our desired quantities for the simplified E-Step are $\boldsymbol{\mu}$
and $\boldsymbol{s}$, as calculated by Eq. (13). We can then perform the
M-Step update in Eq. (7) as

$\displaystyle{\alpha}^{(t+1)}_{j}=\frac{1}{\mu_{j}^{2}+s_{j}},$ (15)

completely avoiding the need to compute or invert the covariance matrix
$\boldsymbol{\Sigma}$. Algorithm 1 gives the full CoFEM algorithm.

The diagonal estimator in Eq. (13) is unbiased for any ${K\geq 1}$, yet its
variance is proportional to $1/K$, as we will show in Section IV-A. We will
provide theoretical justification that small $K$ suffices for CoFEM in
practice, and that this $K$ can be constant with respect to the dimensionality
$D$ for a large class of compressed sensing matrices $\boldsymbol{\Phi}$.

Algorithm 1 CovarianceFreeEM($\boldsymbol{y}$, $\boldsymbol{\Phi}$, $\beta$,
$T$, $K$)

1:Initialize $\alpha_{j}^{(1)}\leftarrow 1$ for $j=1,\ldots,D$.

2:for $t=1,2,\ldots,T$ do

3:  _// Simplified E-Step_

4: Define
$\mathbb{A}\leftarrow\beta\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\mathrm{diag}\\{\boldsymbol{\alpha}^{(t)}\\}$.

5: Draw $\boldsymbol{p}_{1},\boldsymbol{p}_{2},\ldots,\boldsymbol{p}_{K}\sim$
Rademacher distribution.

6: Define
$\mathbb{B}\leftarrow[\boldsymbol{p}_{1}\>|\>\boldsymbol{p}_{2}\>|\>\ldots\>|\>\boldsymbol{p}_{K}\>|\>\beta\boldsymbol{\Phi}^{\top}\boldsymbol{y}]$.

7:
$[\boldsymbol{x}_{1}\>|\>\boldsymbol{x}_{2}\>|\>\ldots\>|\>\boldsymbol{x}_{K}\>|\>\boldsymbol{\mu}]\leftarrow$
LinearSolver($\mathbb{A},\mathbb{B}$).

8: Compute $s_{j}\leftarrow 1/K\sum_{k=1}^{K}p_{k,j}\cdot x_{k,j}$ for
$j=1,\ldots,D$.

9:  _// M-Step_

10: if $t<T$ then

11: Update $\alpha_{j}^{(t+1)}\leftarrow 1/(\mu_{j}^{2}+s_{j})$ for
$j=1,\ldots,D$.

12: end if

13:end for

14:return $\boldsymbol{\alpha}^{(T)},\boldsymbol{\mu},\boldsymbol{s}$

###  III-B Parallel Conjugate Gradient

Among potential options for the linear solver in Algorithm 1, we choose the
conjugate gradient (CG) algorithm due to its efficiency and flexibility [29,
30]. CG is an iterative approach with a series of matrix-vector multiplication
steps to solve a linear system $\mathbb{A}\boldsymbol{x}=\boldsymbol{b}$. At
each step, CG does not need a physical manifestation of $\mathbb{A}$; it only
requires a way to apply $\mathbb{A}$ to an arbitrary vector $\boldsymbol{v}$
to yield $\mathbb{A}\boldsymbol{v}$. For SBL,
$\mathbb{A}:=\beta\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\mathrm{diag}\\{\boldsymbol{\alpha}^{(t)}\\}$,
which implies the complexity of CG is dominated by the time $\tau_{D}$ it
takes to apply $\boldsymbol{\Phi}$ and its transpose to $\boldsymbol{v}$. In
many applications, $\boldsymbol{\Phi}$ is a structured and matrix-free
operation. Examples include convolution, discrete cosine transform, Fourier
transform, and wavelet transform – all of which require at most
$\tau_{D}=O(D\log D)$-time [31].

CG can also easily generalize to solving _multiple_ linear systems
$\mathbb{A}\mathbb{X}=\mathbb{B}$ by simply replacing the matrix-vector
multiplications with matrix-matrix multiplications. For faster computing,
these operations can be parallelized on GPUs. In addition, CG is space-
efficient and only needs $O(D)$-space to solve the linear system; this is the
minimum requirement for any solver given that the output
$\boldsymbol{x}\in\mathbb{R}^{D}$.

###  III-C Preconditioning

The time complexity of CG depends on the number of CG steps $U$. It is
guaranteed that $U\leq D$ [29], yet $D$ can be very large for high-dimensional
problems. In general, if $\mathbb{A}$ has a small condition number
$\kappa(\mathbb{A}):=\lambda_{\text{max}}(\mathbb{A})/\lambda_{\text{min}}(\mathbb{A})$
where $\lambda_{\text{max}}(\mathbb{A})$ and
$\lambda_{\text{min}}(\mathbb{A})$ are the largest and smallest eigenvalues of
$\mathbb{A}$, then $U\ll D$ steps are needed to find an $\boldsymbol{\hat{x}}$
such that
${\lVert\mathbb{A}\boldsymbol{\hat{x}}-\boldsymbol{b}\rVert_{2}/\lVert\boldsymbol{b}\rVert_{2}\leq\epsilon_{\text{max}}}$
for small $\epsilon_{\text{max}}$. However, optimizing the SBL objective
function pushes many diagonal elements of $\mathbb{A}$ to $\infty$, resulting
in large $\kappa(\mathbb{A})$ and necessitating large $U$.

To resolve this issue, we incorporate a _preconditioner_ matrix
$\mathbb{M}\in\mathbb{R}^{D\times D}$ in CG. Instead of directly solving
$\mathbb{A}\boldsymbol{x}=\boldsymbol{b}$, we solve an equivalent system
$\mathbb{A}^{\prime}\boldsymbol{x}^{\prime}=\boldsymbol{b}^{\prime}$, where
$\mathbb{A}^{\prime}:=\mathbb{M}^{-1/2}\mathbb{A}\mathbb{M}^{-1/2}$,
$\boldsymbol{b}^{\prime}:=\mathbb{M}^{-1/2}\boldsymbol{b}$ and
$\boldsymbol{x}^{\prime}:=\mathbb{M}^{1/2}\boldsymbol{x}$. We consider two
criteria for choosing $\mathbb{M}$ in CoFEM. First, we want
$\kappa(\mathbb{A}^{\prime})\ll\kappa(\mathbb{A})$ to reduce the number of
steps $U$ that are necessary to achieve small $\epsilon_{\text{max}}$. Next,
we want to maintain the scalability of CoFEM (e.g. absence of matrix
inversion, $O(D)$-space complexity) by avoiding a dense matrix for
$\mathbb{M}$.

We thus propose to use a _diagonal preconditioner_
$\mathbb{M}=\mathrm{diag}\\{\beta\boldsymbol{\theta}+\boldsymbol{\alpha}^{(t)}\\}$,
where $\boldsymbol{\theta}\in\mathbb{R}^{D}$ is a set of customizable positive
values. This $\mathbb{M}$ is easy to invert, only requires $O(D)$-space, and
can be quickly applied to vectors. By varying $\boldsymbol{\theta}$, we can
arrive at different choices for $\mathbb{M}$. For example, setting
$\theta_{j}=\sum_{i=1}^{N}\mathrm{\Phi}_{i,j}^{2}$ leads to the popular
_Jacobi preconditioner_ satisfying $\mathrm{M}_{j,j}=\mathrm{A}_{j,j}$ [30].
In Section IV-B, we provide novel theoretical analysis for our diagonal
preconditioner within the context of SBL. We illustrate that for a large class
of compressed sensing dictionaries $\boldsymbol{\Phi}$, setting $\theta_{j}=1$
for all $j$ is a favorable choice that leads to small
$\kappa(\mathbb{A}^{\prime})$ and enables $U$ to be constant with respect to
the dimensionality $D$.

Algorithm 2 summarizes the parallel CG algorithm for inputs
$\mathbb{A}\in\mathbb{R}^{D\times D}$, $\mathbb{M}\in\mathbb{R}^{D\times D}$,
and $\mathbb{B}\in\mathbb{R}^{D\times Q}$, where $Q$ is the number of parallel
systems. For CoFEM, we have $Q=K+1$. The computation is dominated by line 5,
in which $\mathbb{A}$ is applied to vectors stored as columns of a matrix.

Algorithm 2
ParallelConjGradient($\mathbb{A},\mathbb{B},\mathbb{M},U,\epsilon_{\text{max}}$)

1:Initialize $\mathbb{X}$ as a $D\times Q$ matrix of all zeros.

2:Initialize $\mathbb{R}\leftarrow\mathbb{B}$ and
$\mathbb{P}\leftarrow\mathbb{B}$ and
$\mathbb{W}\leftarrow\mathbb{M}^{-1}\mathbb{B}$.

3:Compute
$\rho_{q}\leftarrow\sum_{j=1}^{D}\mathrm{R}_{j,q}\cdot\mathrm{W}_{j,q}$ for
$q=1,\ldots,Q$.

4:for $u=1,2,\ldots,U$ do

5: Compute $\boldsymbol{\Psi}\leftarrow\mathbb{A}\mathbb{P}$. _// Apply
matrix._

6: Compute $\pi_{q}\leftarrow\sum_{j=1}^{D}\mathrm{P}_{d,q}\cdot\Psi_{d,q}$
for $q=1,\ldots,Q$.

7: Compute $\gamma_{q}\leftarrow\rho_{q}/\pi_{q}$ for $q=1,\ldots,Q$.

8: Update $\mathbb{X}\leftarrow\mathbb{X}+\mathbb{P}\mathbb{\Gamma}$, where
$\mathbb{\Gamma}=\mathrm{diag}\\{\boldsymbol{\gamma}\\}$.

9: Update $\mathbb{R}\leftarrow\mathbb{R}-\mathbb{\Psi}\mathbb{\Gamma}$, where
$\mathbb{\Gamma}=\mathrm{diag}\\{\boldsymbol{\gamma}\\}$.

10: Let $\delta\leftarrow||\mathbb{R}||_{F}/||\mathbb{B}||_{F}$, where $F$ is
Frobenius norm.

11: if $\delta\leq\epsilon_{\text{max}}$ then

12: return $\mathbb{X}$

13: end if

14: Compute $\mathbb{W}\leftarrow\mathbb{M}^{-1}\mathbb{R}$. _// Apply
preconditioner._

15: Let $\rho^{\text{old}}_{q}\leftarrow\rho_{q}$ for $q=1,\ldots,Q$.

16: Compute
$\rho_{q}\leftarrow\sum_{j=1}^{D}\mathrm{R}_{j,q}\cdot\mathrm{W}_{j,q}$ for
$q=1,\ldots,Q$.

17: Compute $\eta_{q}\leftarrow\rho_{q}/\rho^{\text{old}}_{q}$ for
$q=1,\ldots,Q$.

18: Update $\mathbb{P}\leftarrow\mathbb{R}+\mathbb{P}\mathbb{H}$, where
$\mathbb{H}=\mathrm{diag}\\{\boldsymbol{\eta}\\}$.

19:end for

20:return $\mathbb{X}$

###  III-D Complexity Comparisons

Each of the $T$ iterations of CoFEM requires at most $U$ steps of CG – in
which we apply $\boldsymbol{\Phi}$ (and $\boldsymbol{\Phi}^{\top}$) in
$\tau_{D}$-time to $K$ vectors – giving us an overall time complexity of
$O(T\tau_{D}UK)$. CoFEM’s space complexity is dominated by CG, which requires
$O(D)$-space for each of the $(K+1)$ systems. Table I shows the complexities
of CoFEM and other SBL inference schemes. While many other methods improve
upon EM, they introduce dependencies on $N$ or $d$, which typically grow with
$D$. For example, if the size of the signal $\boldsymbol{z}$ is doubled, we
may also expect the number of measurements $N$ to be doubled (to achieve same
reconstruction error), as well as the number $d$ of non-zero values in
$\boldsymbol{z}$. Thus, increasing $D$ compounds the increase in complexities
of these algorithms. In contrast, CoFEM’s dependencies on $U$ and $K$ can be
held constant as $D$ increases, which we demonstrate in Section V.

TABLE I: Computational complexities of SBL inference schemes. Method | Time | Space
---|---|---
EM [2] | $O(TD^{3})$ | $O(D^{2})$
IRLS [16] | $O(T(DN^{2}+N^{3}))$ | $O(D^{2})$
AMP [17] | $O(TDNT_{\text{amp}})$ | $O(DN)$
VI [18] | $O(T\tau_{D})$ | $O(D)$
Seq [19] | $O(Dd^{2})$ | $O(D+d^{2})$
CoFEM (ours) | $O(T\tau_{D}UK)$ | $O(DK)$

###  III-E CoFEM for SBL Extensions

To further highlight the flexibility of CoFEM, we show how it can handle two
common extensions of the SBL model.

####  III-E1 Multi-Task Learning

In multi-task learning, there are $L$ different sparse vector recovery
problems that one wishes to solve at once. These problems may have different
observation-dictionary pairs
$(\boldsymbol{y}_{1},\boldsymbol{\Phi}_{1}),(\boldsymbol{y}_{2},\boldsymbol{\Phi}_{2}),\ldots,(\boldsymbol{y}_{L},\boldsymbol{\Phi}_{L})$,
yet the tasks are related in the sense that their corresponding vectors
$\boldsymbol{z}_{1},\boldsymbol{z}_{2},\ldots,\boldsymbol{z}_{L}$ have similar
non-zero supports. Some examples include multiple measurements vector (MMV)
problems [32], multi-task compressed sensing [4], and sparse Bayesian learning
with complex numbers [33].

A simple way to enforce joint sparsity among all tasks in SBL is to have them
share a common $\boldsymbol{\alpha}$ vector:

$\displaystyle\boldsymbol{z}_{\ell}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{0},\mathrm{diag}\\{\boldsymbol{\alpha}\\}^{-1}),$
$\displaystyle\ell=1,2,\ldots,L,$ $\displaystyle\boldsymbol{y}_{\ell}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{\Phi}_{\ell}\boldsymbol{z}_{\ell},1/\beta\mathbf{I}),$
$\displaystyle\ell=1,2,\ldots,L.$ (16)

Learning takes place through the task-separable objective:

$\displaystyle\max_{\boldsymbol{\alpha}}\log
p(\boldsymbol{y}_{1},\boldsymbol{y}_{2},\ldots,\boldsymbol{y}_{L}\>|\>\boldsymbol{\alpha})=\sum_{\ell=1}^{L}\log
p(\boldsymbol{y}_{\ell}\>|\>\boldsymbol{\alpha}).$ (17)

To optimize Eq. (17), EM runs a E-Step for each task $\ell$,

$\displaystyle\boldsymbol{\mu}_{\ell}=\beta\boldsymbol{\Sigma}_{\ell}\boldsymbol{\Phi}_{\ell}^{\top}\boldsymbol{y}_{\ell},$
$\displaystyle\boldsymbol{\Sigma}_{\ell}=(\beta\boldsymbol{\Phi}_{\ell}^{\top}\boldsymbol{\Phi}_{\ell}+\text{diag}\\{\boldsymbol{\alpha}^{(t)}\\})^{-1},$
(18)

followed by a M-Step to combine these moments,

$\displaystyle{\alpha}^{(t+1)}_{j}=\frac{L}{\sum_{\ell=1}^{L}\mu_{\ell,j}^{2}+\Sigma_{\ell,j,j}}.$
(19)

To accelerate EM, CoFEM can simply replace the E-Step for each task $\ell$
with a covariance-free version, as described in Section III-A. We can further
parallelize these $L$ E-Steps by solving their $L\cdot(K+1)$ systems all at
once through Alg. 2.

####  III-E2 Non-Negativity Constraints

In some applications, we may want to enforce non-negativity on
$\boldsymbol{z}$. In these cases, we can use an independent _rectified
Gaussian_ prior $\mathcal{N}^{R}(0,1/\alpha_{j})$ for each component $z_{j}$
of $\boldsymbol{z}$ [34], which places zero probability mass on the negative
values. Specifically, we have

$\displaystyle
p(z_{j}\>|\>\alpha_{j})=\begin{cases}\sqrt{2\alpha_{j}/\pi}\exp(-\alpha_{j}z_{j}^{2}/2),&z_{j}>0,\\\
1/2,&z_{j}=0,\\\ 0,&z_{j}<0.\end{cases}$ (20)

Due to conjugacy between the rectified Gaussian and Gaussian distributions,
the posterior $p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\hat{\alpha}})$
is also a rectified Gaussian. However, this posterior’s density function is
not analytically tractable, so we can follow [34] and approximate it with a
diagonal rectified Gaussian whose second moment is

$\displaystyle\mathbb{E}_{p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\hat{\alpha}})}[z_{j}^{2}]=\mu_{j}^{2}+\Sigma_{j,j}+\mu_{j}\cdot\sqrt{\frac{\Sigma_{j,j}}{\pi}}\cdot\frac{\exp(-\xi_{j}^{2})}{\mathrm{erfc}(-\xi_{j})},$
(21)

where $\xi_{j}=\mu_{j}/\sqrt{2\Sigma_{j,j}}$ and
$\mathrm{erfc}(x)=2/\sqrt{\pi}\int_{x}^{\infty}\exp(-t^{2})dt$ is the
complimentary error function. To compute Eq. (21), CoFEM can replace
$\Sigma_{j,j}$ with $s_{j}$ from Eq. (13). The M-Step then updates
$\alpha_{j}$ as the reciprocal of Eq. (21).

##  IV Theoretical Analysis of CoFEM

The two main hyperparameters of CoFEM are the number of probe vectors $K$ and
the number of conjugate gradient steps $U$. These values determine the per-
iteration time complexity of CoFEM as $O(UK\tau_{D})$ and its space complexity
as $O(DK)$. Thus, it is important to control $U$ and $K$. In this section, we
present new theoretical results illustrating that, for a large class of
dictionaries $\boldsymbol{\Phi}$, these hyperparameters can be kept small even
as the dimensionality of the problem $D$ grows very large.

The ultimate goal of SBL is to solve a Bayesian variant of the sparse coding
problem in which we (a) identify which $z_{j}=0$ and (b) provide uncertainty
quantification for the non-zero $z_{j}$. To achieve true sparsity for a
particular $z_{j}$, it is necessary for its prior parameter
$\alpha_{j}\to\infty$, which is only attained as the number of iterations
$t\to\infty$. Thus, our study of $U$ and $K$ is based on the following SBL
Convergence property:

######  Definition 1 (SBL Convergence).

In Alg. 1, consider the sequence of iterates $\boldsymbol{\alpha}^{(t)}$ for
$t=1,2,\ldots,$ and let
$\boldsymbol{\hat{\alpha}}:=\lim_{t\to\infty}\boldsymbol{\alpha}^{(t)}$. Then,
$(\mathcal{S},\mathcal{U},\boldsymbol{\hat{\alpha}})$-convergence is satisfied
for SBL if the indices $\mathbb{N}_{D}:=\\{1,2,\ldots,D\\}$ can be partitioned
into an “active” set $\mathcal{S}\subseteq\mathbb{N}_{D}$ and an “inactive”
set $\mathcal{U}:=\mathbb{N}_{D}\setminus\mathcal{S}$, where
$\hat{\alpha}_{j}>0$ is finite if $j\in\mathcal{S}$ and
$\hat{\alpha}_{j}=\infty$ if $j\in\mathcal{U}$.

Of course, it is impossible to run $t\to\infty$ iterations in practice and to
the best of the authors’ knowledge, there does not exist a formal proof
characterizing EM’s behavior after finitely many iterations in the context of
SBL. However, many works within the SBL literature [2, 5, 3, 4] have
illustrated (and even relied on) a well-accepted phenomenon in which the
inactive $\alpha_{j}$ grow very large and can be treated as “reaching
infinity” after $T$ finite iterations. This justifies the practical
applicability of Definition 1 and our ensuing theoretical analysis.

Notation For any positive integer $P$, let
$\mathbb{N}_{P}:=\\{1,2,\ldots,P\\}$. For any vector
$\boldsymbol{v}\in\mathbb{R}^{P}$ and set
$\mathcal{B}\subseteq\mathbb{N}_{P}$, define the sub-vector
$\boldsymbol{v}_{\mathcal{B}}:=[v_{b}\>|\>b\in\mathcal{B}]\in\mathbb{R}^{|\mathcal{B}|}$.
Similarly, for any matrix
$\mathbb{M}=[\boldsymbol{v}^{(1)},\ldots,\boldsymbol{v}^{(P)}]\in\mathbb{R}^{Q\times
P}$, define the sub-matrix
$\mathbb{M}_{\mathcal{B}}:=[\boldsymbol{v}^{(b)}\>|\>b\in\mathcal{B}]\in\mathbb{R}^{Q\times|\mathcal{B}|}$.
For any two sets $\mathcal{A}\subseteq\mathbb{N}_{Q}$ and
$\mathcal{B}\subseteq\mathbb{N}_{P}$, define the sub-matrix block
$\mathbb{M}_{\mathcal{A},\mathcal{B}}:=[\boldsymbol{v}^{(b)}_{\mathcal{A}}\>|\>b\in\mathcal{B}]\in\mathbb{R}^{|\mathcal{A}|\times|\mathcal{B}|}$.
Let $\lVert\boldsymbol{v}\rVert_{2}$ denote the Euclidean norm of vector
$\boldsymbol{v}$ and
$\lVert\boldsymbol{v}\rVert_{\mathbb{M}}:=\sqrt{\boldsymbol{v}^{\top}\mathbb{M}\boldsymbol{v}}$
denote the $\mathbb{M}$-weighted norm of $\boldsymbol{v}$. Let
$\lVert\mathbb{M}\rVert_{2}$ denote the spectral norm (i.e. largest singular
value) of matrix $\mathbb{M}$. Let
$\sigma_{\text{min}}(\mathbb{M}),\lambda_{\text{max}}(\mathbb{M})$ and
$\lambda_{\text{min}}(\mathbb{M})$ denote the smallest singular value, largest
eigenvalue, and smallest eigenvalue of $\mathbb{M}$, respectively. For a
matrix $\mathbb{M}$, let
$\kappa(\mathbb{M}):=\lVert\mathbb{M}\rVert_{2}/\sigma_{\text{min}}(\mathbb{M})$
be the condition number of $\mathbb{M}$.

###  IV-A A Theory for the Number of Probe Vectors $K$

First, we analyze the dependency of the diagonal estimator $\boldsymbol{s}$ on
the number of probe vectors $K$. We aim to characterize the standard deviation
(i.e. the square root of the variance) of each $s_{j}$, which leads to this
lemma:

######  Lemma 1.

Let $\mathbb{M}\in\mathbb{R}^{D\times D}$. Consider the estimator
$\boldsymbol{s}$ defined in Prop. 1, where
$\boldsymbol{p}_{1},\ldots,\boldsymbol{p}_{K}$ are independent Rademacher
variables. Then, the standard deviation $\nu_{j}$ of $s_{j}$ is

$\displaystyle\nu_{j}:=\sqrt{\mathbb{E}[(s_{j}-\mathbb{E}[s_{j}])^{2}]}=\sqrt{\frac{1}{K}\sum_{j^{\prime}\neq
j}\mathrm{M}_{j,j^{\prime}}^{2}}.$ (22)

###### Proof.

Within Eq. (22), we substitute the expression for $s_{j}$ from Eq. (11) and
the fact that $\mathbb{E}[s_{j}]=\mathrm{M}_{j,j}$ to yield

$\displaystyle\nu_{j}$
$\displaystyle=\sqrt{\mathbb{E}\left[\left(\sum_{j^{\prime}\neq
j}\mathrm{M}_{j,j^{\prime}}\cdot\frac{\sum_{k=1}^{K}p_{k,j}\cdot
p_{k,j^{\prime}}}{\sum_{k=1}^{K}p_{k,j}^{2}}\right)^{2}\right]}$
$\displaystyle=\sqrt{\mathbb{E}\left[\sum_{j^{\prime}\neq
j}\sum_{j^{\prime\prime}\neq
j}\mathrm{M}_{j,j^{\prime}}\cdot\mathrm{M}_{j,j^{\prime\prime}}\cdot\frac{g_{j,j^{\prime}}}{g_{j,j}}\cdot\frac{g_{j,j^{\prime\prime}}}{g_{j,j}}\right]},$
(23)

where we define $g_{j,\ell}:=\sum_{k=1}^{K}p_{k,j}\cdot p_{k,\ell}$ for all
$(j,\ell)\in\mathbb{N}_{D}\times\mathbb{N}_{D}$. In the denominator, we have
$g_{j,j}=K$ for all $j$ since $p^{2}=1$ for a Rademacher variable. In the
numerator, if $j^{\prime}=j^{\prime\prime}$, we have
$\mathbb{E}[g_{j,j^{\prime}}\cdot
g_{j,j^{\prime\prime}}]=\mathbb{E}[g_{j,j^{\prime}}^{2}]=K$. Otherwise, if
$j\neq j^{\prime\prime}$, $\mathbb{E}[g_{j,j^{\prime}}\cdot
g_{j,j^{\prime\prime}}]=0$. Thus, Eq. (23) simplifies to Eq. (22). ∎

Lemma 1 tells us that the standard deviation of our estimator decreases with
$K$ and increases with the norm of the off-diagonal entries. Analyzing Lemma 1
within the context of SBL leads to our first main theoretical result, as
stated below.

######  Theorem 1.

Let
$\boldsymbol{\Sigma}^{(t)}:=(\beta\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\mathrm{diag}\\{\boldsymbol{\alpha}^{(t)}\\})^{-1}$
be the SBL covariance matrix at the $t$-th iteration of Alg. 1. Let
$\boldsymbol{s}^{(t)}$ be the Rademacher diagonal estimator for
$\boldsymbol{\Sigma}^{(t)}$ defined in Eq. (13) with $K$ probe vectors. Let
${\nu}^{(t)}_{j}$ denote the standard deviation of ${s}^{(t)}_{j}$. Assume
that $(\mathcal{S},\mathcal{U},\boldsymbol{\hat{\alpha}})$-convergence is
satisfied. Then, for any inactive index $j\in\mathcal{U}$, we have

$\displaystyle\lim_{t\to\infty}{\nu}^{(t)}_{j}=0,$ (24)

and for any active index $j\in\mathcal{S}$, we have

$\displaystyle\lim_{t\to\infty}{\nu}^{(t)}_{j}\leq\frac{1}{\sqrt{K}}\cdot\frac{\inf_{\boldsymbol{\Theta}\in\mathcal{O}}\lVert\boldsymbol{\Theta}^{-1}\boldsymbol{\Phi}^{\top}_{\mathcal{S}}\boldsymbol{\Phi}_{\mathcal{S}}-\mathbb{I}\rVert_{2}}{\beta\cdot\sigma^{2}_{\mathrm{min}}(\boldsymbol{\Phi}_{\mathcal{S}})},$
(25)

where $\mathcal{O}$ is the set of $|\mathcal{S}|\times|\mathcal{S}|$ diagonal
matrix with positive diagonal elements and $\mathbb{I}$ is the identity
matrix.

Theorem 1 offers several insights in the limit of EM iterations: (1) the
estimator becomes deterministic with zero standard deviation for the inactive
indices, (2) $K$ only affects the estimator’s standard deviation for the
active indices, and (3) if
$\boldsymbol{\Phi}^{\top}_{\mathcal{S}}\boldsymbol{\Phi}_{\mathcal{S}}$ is
close to any diagonal matrix $\boldsymbol{\Theta}$ (i.e. the columns of
$\boldsymbol{\Phi}_{\mathcal{S}}$ are close to orthogonal), then the standard
deviation for the active indices converge to a small quantity. The proof of
Theorem 1 is given in Appendix A.

###  IV-B A Theory for the Number of Conjugate Gradient Steps $U$

Next, we analyze the number of CG steps $U$ needed for convergence. We build
on the following well-known result [30]:

######  Lemma 2 (CG Convergence).

Consider the CG algorithm for solving
$\mathbb{A}\boldsymbol{x}=\boldsymbol{b}$, where
$\mathbb{A}\in\mathbb{R}^{D\times D}$ is a positive definite matrix and
$\boldsymbol{b}\in\mathbb{R}^{D}$. Let $\boldsymbol{x}_{0}\in\mathbb{R}^{D}$
be the initial solution. Let $\boldsymbol{x}_{u}$ denote the algorithm’s
solution and
$\boldsymbol{r}_{u}:=\boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}_{u}$ denote
the algorithm’s residual at the $u$-th step of CG. Then,

$\displaystyle\lVert\boldsymbol{r}_{u}\rVert_{\mathbb{A}^{-1}}\leq
2\left(\frac{\sqrt{\kappa(\mathbb{A})}-1}{\sqrt{\kappa(\mathbb{A})}+1}\right)^{u}\lVert\boldsymbol{r}_{0}\rVert_{\mathbb{A}^{-1}},$
(26)

where
$\lVert\boldsymbol{r}\rVert_{\mathbb{A}^{-1}}:=\sqrt{\boldsymbol{r}^{\top}\mathbb{A}^{-1}\boldsymbol{r}}$
for any vector $\boldsymbol{r}\in\mathbb{R}^{D}$.

######  Corollary 1.

Let CG with matrix $\mathbb{A}$ and any $\boldsymbol{b}\in\mathbb{R}^{D}$
start with the initialization $\boldsymbol{x}_{0}=\boldsymbol{0}$. Let
$\epsilon:=\lVert\boldsymbol{r}_{U}\rVert_{\mathbb{A}^{-1}}/\lVert\boldsymbol{b}\rVert_{\mathbb{A}^{-1}}$
denote the relative residual error of CG after $U$ steps111Theoretical results
for CG often define $\epsilon$ in terms of the norm weighted by
$\mathbb{A}^{-1}$ (or $\mathbb{A}$), even though
$\epsilon_{2}:=\lVert\boldsymbol{r}_{U}\rVert_{2}/\lVert\boldsymbol{b}\rVert_{2}$
is used to determine convergence in software. We follow this convention while
noting that one can exploit the relation
$\epsilon_{2}\leq\epsilon\sqrt{\kappa(\mathbb{A})}$ to extend our results for
$\epsilon_{2}$.. Then,

$\displaystyle\epsilon\leq 2\exp(-U/\sqrt{\kappa(\mathbb{A})}).$ (27)

###### Proof.

Since $\boldsymbol{r}_{0}=\boldsymbol{b}$, we can apply Eq. (26) to obtain

$\displaystyle\epsilon\leq
2\left(\frac{\sqrt{\kappa(\mathbb{A})}-1}{\sqrt{\kappa(\mathbb{A})}+1}\right)^{U}\leq
2\left(1-\frac{1}{\sqrt{\kappa(\mathbb{A})}}\right)^{U}.$ (28)

The inequality $1-1/x\leq\exp(-1/x)$ holds for any $x\in\mathbb{R}$. Taking
$x=\sqrt{\kappa(\mathbb{A})}$ in Eq. (28) gives the result. ∎

Corollary 1 indicates that the relative residual error of CG decreases
exponentially with $U$, yet the precise exponent depends on
$\kappa(\mathbb{A})$. In the $t$-th iteration of CoFEM, we wish to solve
linear systems with
$\mathbb{A}^{(t)}:=\beta\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\mathrm{diag}\\{\boldsymbol{\alpha}^{(t)}\\}$.
However, since $\boldsymbol{\alpha}_{j}^{(t)}\rightarrow\infty$ for
$j\in\mathcal{U}$, this leads to $\kappa(\mathbb{A}^{(t)})\to\infty$ and a
vacuous bound on $\epsilon$ for any finite number of steps $U$. Therefore,
CoFEM introduces a preconditioner $\mathbb{M}$ to construct a new CG matrix
$\mathbb{A^{\prime}}^{(t)}$ with a reduced condition number. In our second
main theoretical result, we analyze $\kappa(\mathbb{A^{\prime}}^{(t)})$, the
bound that it induces on error, and the implications for $U$.

######  Theorem 2.

Let
$\mathbb{A}^{(t)}:=\beta\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\mathrm{diag}\\{\boldsymbol{\alpha}^{(t)}\\}$
denote the SBL inverse-covariance matrix at the $t$-th iteration of Alg. 1.
Let
$\mathbb{M}^{(t)}:=\mathrm{diag}\\{\beta\boldsymbol{\theta}+\boldsymbol{\alpha}^{(t)}\\}$
denote the preconditioner, where $\boldsymbol{\theta}\in\mathbb{R}^{D}$ is a
vector of positive values. Define the preconditioned matrix
$\mathbb{A^{\prime}}^{(t)}:=(\mathbb{M}^{(t)})^{-1/2}\mathbb{A}^{(t)}(\mathbb{M}^{(t)})^{-1/2}$.
Let $\boldsymbol{b}^{(t)}\in\mathbb{R}^{D}$ be any vector and
$\boldsymbol{b^{\prime}}^{(t)}:=(\mathbb{M}^{(t)})^{-1/2}\boldsymbol{b}^{(t)}$.
Let ${\epsilon}^{(t)}$ be the relative residual error after running $U$
conjugate gradient steps to solve the system
$\mathbb{A^{\prime}}^{(t)}\boldsymbol{x^{\prime}}^{(t)}=\boldsymbol{b^{\prime}}^{(t)}$
with $\boldsymbol{x^{\prime}}^{(t)}_{0}=\boldsymbol{0}$. Given
$(\mathcal{S},\mathcal{U},\boldsymbol{\hat{\alpha}})$-convergence, it follows
that

$\displaystyle\lim_{t\to\infty}{\epsilon}^{(t)}\leq
2\exp\left(-U\sqrt{\frac{1-\lVert\boldsymbol{\Theta}^{-1}\boldsymbol{\Phi}_{\mathcal{S}}^{\top}\boldsymbol{\Phi}_{\mathcal{S}}-\mathbb{I}\rVert_{2}}{1+\lVert\boldsymbol{\Theta}^{-1}\boldsymbol{\Phi}_{\mathcal{S}}^{\top}\boldsymbol{\Phi}_{\mathcal{S}}-\mathbb{I}\rVert_{2}}}\right),$
(29)

where $\mathbb{\Theta}=\mathrm{diag}\\{\boldsymbol{\theta}_{\mathcal{S}}\\}$
and $\mathbb{I}$ is the identity matrix.

Eq. (29) indicates that faster convergence is achieved for
$\boldsymbol{\theta}_{\mathcal{S}}\in\mathbb{R}^{|\mathcal{S}|}$ that
minimizes
$\lVert\mathrm{diag}\\{\boldsymbol{\theta}_{\mathcal{S}}\\}^{-1}\boldsymbol{\Phi}_{\mathcal{S}}^{\top}\boldsymbol{\Phi}_{\mathcal{S}}-\mathbb{I}\rVert_{2}$.
In practice, the support $\mathcal{S}$ is not known in advance, so we may
instead choose $\boldsymbol{\theta}\in\mathbb{R}^{D}$ to minimize
$\lVert\mathrm{diag}\\{\boldsymbol{\theta}\\}^{-1}\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}-\mathbb{I}\rVert_{2}$.
For example, if $\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}$ is a diagonal
matrix (i.e. all columns of $\boldsymbol{\Phi}$ are mutually orthogonal), then
the optimal choice is $\theta_{j}=\sum_{i=1}^{N}\mathrm{\Phi}_{i,j}^{2}$,
which corresponds to the Jacobi preconditioner. The proof of Theorem 2 is
given in Appendix B.

###  IV-C A Theory for $U$ and $K$ for Compressed Sensing Matrices

We now show that our bounds in Theorems 1 and 2 can be simplified for a large
class of matrices $\boldsymbol{\Phi}$ satisfying the restricted isometry
property (RIP) and commonly employed in compressed sensing applications.

######  Definition 2 (Restricted Isometry Property).

Let $\boldsymbol{\Phi}\in\mathbb{R}^{N\times D}$, $d\leq D$, and $\delta>0$.
Then, $\mathbb{\Phi}$ satisfies $(d,\delta)$-RIP if for every set
$\mathcal{C}\subseteq\\{1,\ldots,D\\}$ of size $|\mathcal{S}|=d$ and every
vector $\boldsymbol{v}\in\mathbb{R}^{d}$,

$\displaystyle(1-\delta)\lVert\boldsymbol{v}\rVert_{2}^{2}\leq\lVert\boldsymbol{\Phi}_{\mathcal{C}}\boldsymbol{v}\rVert_{2}^{2}\leq(1+\delta)\lVert\boldsymbol{v}\rVert_{2}^{2}.$
(30)

######  Corollary 2.

Let $\boldsymbol{\Phi}\in\mathbb{R}^{N\times D}$ satisfy $(d,\delta)$-RIP. For
any $\mathcal{C}\subseteq\mathbb{N}_{D}$ with $|\mathcal{C}|=d$, let
$\boldsymbol{\Delta}_{\mathcal{C},\mathcal{C}}:={\boldsymbol{\Phi}_{\mathcal{C}}^{\top}\boldsymbol{\Phi}_{\mathcal{C}}}-\mathbb{I}$.
Then, $\lVert\mathbb{\Delta}_{\mathcal{C},\mathcal{C}}\rVert_{2}\leq\delta$.

###### Proof.

From Eq. (30), we have for any vector $\boldsymbol{v}\in\mathbb{R}^{d}$,

$\displaystyle 1-\delta$
$\displaystyle\leq\frac{\boldsymbol{v}^{\top}\boldsymbol{\Phi}_{\mathcal{C}}^{\top}\boldsymbol{\Phi}_{\mathcal{C}}\boldsymbol{v}}{\boldsymbol{v}^{\top}\boldsymbol{v}}\leq
1+\delta\implies\left|\frac{\boldsymbol{v}^{\top}\mathbb{\Delta}_{\mathcal{C},\mathcal{C}}\boldsymbol{v}}{\boldsymbol{v}^{\top}\boldsymbol{v}}\right|\leq\delta.$
(31)

Eq. (31) bounds the Rayleigh quotient (and eigenvalues) of
$\mathbb{\Delta}_{\mathcal{C},\mathcal{C}}$ between $[-\delta,\delta]$. The
spectral norm of a symmetric matrix is equal to its largest absolute
eigenvalue. ∎

######  Corollary 3.

In Theorem 1, let $\boldsymbol{\Phi}\in\mathbb{R}^{N\times D}$ satisfy
$(d,\delta)$-RIP, where $d=|\mathcal{S}|$. Then, for $j\in\mathcal{S}$, the
standard deviation of the diagonal estimator ${\nu}^{(t)}_{j}$ satisfies

$\displaystyle\lim_{t\to\infty}{\nu}^{(t)}_{j}\leq\frac{1}{\sqrt{K}}\cdot\frac{\delta}{\beta(1-\delta)}.$
(32)

###### Proof.

In Eq. (25), take $\boldsymbol{\Theta}=\mathbb{I}$, which reduces the
numerator to
$\lVert\mathbb{\Phi}^{\top}_{\mathcal{S}}\mathbb{\Phi}_{\mathcal{S}}-\mathbb{I}\rVert_{2}$.
By Corollary 2, this quantity is at most $\delta$. Finally, by Definition 2,
$\sigma_{\text{min}}^{2}(\boldsymbol{\Phi}_{\mathcal{S}})\geq(1-\delta)$. ∎

######  Corollary 4.

In Theorem 2, let $\boldsymbol{\Phi}\in\mathbb{R}^{N\times D}$ satisfy
$(d,\delta)$-RIP, where $d=|\mathcal{S}|$. Let $\theta_{j}=1$ for all
$j\in\mathbb{N}_{D}$. Then, the relative residual error ${\epsilon}^{(t)}$ of
conjugate gradient satisfies

$\displaystyle\lim_{t\to\infty}{\epsilon}^{(t)}\leq
2\exp\left(-U\sqrt{\frac{1-\delta}{1+\delta}}\right).$ (33)

###### Proof.

This follows from applying Corollary 2 to Eq. (29). ∎

One may have expected that as $N$ and $D$ increase, CoFEM’s hyperparameters
$U$ and $K$ must also increase proportionally to ensure that $\nu_{j}$ and
$\epsilon$ remain small. However, Corollaries 3 and 4 illustrate that for RIP
matrices222RIP is one mathematical notion for the idea of “close to
orthonormality” for a set of dictionary columns. There exist other notions
(e.g. incoherence, null space property) [35]. We conjecture there are bounds
similar to Corollaries 3 and 4 that hold for matrices that satisfy these other
properties., the bounds on $U$ and $K$ only depend on $\delta$ (not $N$ or
$D$). In compressed sensing, $\delta$ can be small even as $N$ and $D$ grow
very large [35]. In Section V, we use this insight to demonstrate that even
for very large $D$, CoFEM can accurately perform sparse coding with small
constant values for $U$ and $K$.

##  V Simulated Experiments

In this section, we run a series of experiments to compare the accuracy and
scalability of CoFEM to that of other SBL inference methods across a broad
range of different settings.

Figure 1: Comparing the accuracy of different SBL inference schemes. Each
point represents the mean of 25 trials.

###  V-A Experimental Setup

####  V-A1 General Structure

We design all simulations with a structure inspired by Section V-A of [3]. We
form a ground-truth latent vector $\boldsymbol{z}^{*}\in\mathbb{R}^{D}$ of
spikes by drawing $d$ of its components from a distribution $\mathcal{P}$ and
setting the other $D-d$ components to zero. The location of the spikes are
chosen uniformly at random. Given a dictionary
$\boldsymbol{\Phi}\in\mathbb{R}^{N\times D}$, the observed data
$\boldsymbol{y}\in\mathbb{R}^{N}$ is generated as
$\boldsymbol{y}=\boldsymbol{\Phi}\boldsymbol{z}^{*}+\boldsymbol{\varepsilon}$,
where
$\boldsymbol{\varepsilon}\sim\mathcal{N}(0,\sigma^{2}\mathbf{I})\in\mathbb{R}^{N}$
with $\sigma=0.01$. The goal is to apply SBL for reconstructing
$\boldsymbol{z}^{*}$. Success is measured through minimization of normalized
root mean squared error (NRMSE),
$||\boldsymbol{\hat{z}}-\boldsymbol{z}^{*}||_{2}/||\boldsymbol{z}^{*}||_{2}\times
100$, where $\boldsymbol{\hat{z}}=\boldsymbol{\mu}$ is the mean of the
distribution $p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\hat{\alpha}})$
upon convergence.

####  V-A2 Dictionary

We consider three types of dictionaries $\mathbb{\Phi}$:

  * •

Dense Gaussian: We draw each element $\Phi_{i,j}$ independently from
$\mathcal{N}(0,1/N)$ to form a dense matrix $\boldsymbol{\Phi}$. The spike
distribution of $\boldsymbol{z}^{*}$ is $\mathcal{P}:=\mathrm{Uniform}(-2,2)$.

  * •

DCT: We let $\boldsymbol{\Phi}=\mathbb{M}\mathbb{\Omega}^{-1}$, where
$\mathbb{\Omega}\in\mathbb{R}^{D\times D}$ is the matrix corresponding to the
1D discrete cosine transform of size $D$ and $\mathbb{M}\in\mathbb{R}^{N\times
D}$ is an undersampling operator that selects $N$ out of $D$ components (where
$N\leq D$). The spike distribution is $\mathcal{P}:=\mathcal{N}(0,5)$.

  * •

Convolutional: We set $N=D$ and $\boldsymbol{\Phi}\in\mathbb{R}^{D\times D}$
to a convolution in which its columns are delayed (and truncated) repetitions
of an exponentially decaying filter $\boldsymbol{\phi}\in\mathbb{R}^{D}$. That
is, $\phi_{j}:=(1-\rho)^{j-1}$ for $0<\rho<1$. Thus, $\boldsymbol{\Phi}$ is a
lower triangular matrix in which the $j$-th column is a concatenation of $j-1$
zeros and $\\{\phi_{1},\ldots,\phi_{N-(j-1)}\\}$. The spike distribution is
$\mathcal{P}:=\mathrm{Exponential}(1.5)$.

###  V-B Accuracy Analysis

For our first analysis, we fix $D=1024$ and evaluate the accuracy of various
SBL inference schemes across different settings. We compare among EM/IRLS [2],
AMP [17], VI [18], Seq [19], and CoFEM. Note that EM and IRLS always yield the
same result; IRLS is just the Woodbury identity (Eq. (8)) applied to EM. All
EM-based methods (EM/IRLS, AMP, VI, CoFEM) are executed for $T=50$ iterations.
AMP employs $T_{\text{amp}}=10$ inner loops. For CoFEM, Corollary 3 tells us
that we can keep the number of probes $K$ very small since $\delta\approx 0$
and $\beta=1/(0.01)^{2}=10{,}000$. We use $K=20$ probe vectors, though we have
found that even smaller values for $K$ do not change the results. We employ
$U=400$ maximum CG steps with early termination if the residual error drops
below the threshold $\epsilon_{\text{max}}=10^{-4}$. Results are displayed in
Fig. 1.

####  V-B1 Performance vs. Sparsity

In Fig. 1(a), we consider the dense Gaussian dictionary with $N=\lfloor
D/4\rfloor$, which is typical in a _compressed sensing_ setting. Let the
_sparsity factor_ $f\in[0,1]$ determine the number of non-zero coefficients in
the latent signal $\boldsymbol{z}^{*}$ as $d=\lfloor f\cdot D\rfloor$. We vary
$f$ and observe its impact on NRMSE. At low $f$, all algorithms perform well
except VI. As $f$ increases, EM/IRLS, AMP, and CoFEM exhibit the same decay in
performance, while Seq decays more rapidly. In Fig. 1(b), we consider the DCT
dictionary with $N=\lfloor D/3\rfloor$ measurements, again varying $f$. We see
the same overall trend as in Fig. 1(a). In Fig. 1(c), we compare NRMSE versus
sparsity level $f$ for the convolutional dictionary. The decay rate is set to
$\rho=0.04$. EM/IRLS, Seq, and CoFEM have near-perfect NRMSE at all $f$.
However, VI fails again due to its biased objective and AMP fails due to the
convolutional dictionary not being zero-mean and sub-Gaussian.

####  V-B2 Performance vs. Undersampling

In Fig. 1(d), we revisit the Gaussian dictionary and fix $f=0.06$. We vary the
_undersampling rate_ $r>0$, where $N=\lfloor D/r\rfloor$. We see that EM/IRLS,
AMP, and CoFEM have similar performance at all $r$, while VI and Seq degrade
more rapidly with increasing $r$. A similar trend is shown in Fig. 1(e), in
which we vary $r$ for the DCT dictionary while fixing $f=0.12$.

####  V-B3 Performance vs. Decay Rate

Finally, in Fig. 1(f), we consider the convolutional dictionary with fixed
$f=0.2$ and vary the decay factor $\rho$. Smaller $\rho$ leads to slower decay
of the exponential filter, which increases the correlations between the
columns of $\boldsymbol{\Phi}$. We observe that EM/IRLS, Seq, and CoFEM are
all robust to changes in $\rho$. The performance of VI is heavily correlated
with $\rho$, while AMP diverges for all values of $\rho$.

CoFEM is the only algorithm that is as robust as EM/IRLS to changes in
sparsity level, undersampling rate, and correlation between dictionary
columns. We attribute this robustness to the fact that CoFEM performs an
unbiased estimation of the posterior variances, regardless of the dictionary
structure.

Figure 2: Comparing the scalability of different SBL inference schemes on log-
log scales. Some lines omit points for high dimensions due to memory issues.

###  V-C Scalability Analysis

We now analyze the scalability of CoFEM and other algorithms, by considering
how computation time and memory requirements change as $D$ is increased. In
all settings, CoFEM has $K=20$ probes and $U=400$ maximum CG steps (with
$\epsilon_{\text{max}}=10^{-4}$), regardless of $D$. The preconditioner
employs $\theta_{j}=1$ for all $j$. Results are in Fig. 2.

####  V-C1 CoFEM vs. EM-based Algorithms

We begin by comparing CoFEM against the other EM-based algorithms (EM, IRLS,
AMP). In each setting, we run each EM-based algorithm for $T=30$ EM
iterations, as this is sufficient to obtain low NRMSE (i.e. $<2\%$). VI is
omitted because it is unable to reach this threshold in many cases, as
illustrated by Fig. 1.

For our first setting (Fig. 2(a)), we consider the dense Gaussian
$\boldsymbol{\Phi}$ with $D=2^{p}$ for $p=\\{9,10,\ldots,15\\}$. For each $D$,
we let $N=\lfloor D/4\rfloor$ and $d=\lfloor 0.06D\rfloor$. As $D$ increases,
we observe that CoFEM becomes faster than EM by an order of magnitude. CoFEM
is slower than AMP and IRLS, yet the gap decreases for large $D$. This is
because it takes $\tau_{D}=O(DN)$-time to apply a dense $\boldsymbol{\Phi}$ to
a vector, so CoFEM has the same asymptotic complexity as AMP for the dense
case (see Table I). Furthermore, due to its space-saving and parallelization-
friendly design, CoFEM can exploit a GPU333We use a 16-GB Nvidia T4 GPU and a
32-GB, 2.3 GHz Intel Xeon CPU. to be up to $7\times$ faster than IRLS/AMP and
$100\times$ faster than EM.

Next, in Fig. 2(b) and Fig. 2(c), we repeat the experiment of increasing $D$
for the two structured $\boldsymbol{\Phi}$ (DCT and convolutional). In both
cases, there are fast algorithms for applying $\boldsymbol{\Phi}$ to an
arbitrary vector in $\tau_{D}=O(D\log D)$-time. As a result, CoFEM is faster
at high $D$ than all other algorithms. In the DCT case, we let $d=\lfloor
0.12D\rfloor$ and $N=\lfloor D/3\rfloor$ for all $D$. We observe that CoFEM
can be faster than EM by up to $360\times$ on the CPU and $3800\times$ on the
GPU, reducing over two hours of computation for EM at $D=2^{15}$ to two
seconds. In the convolutional case, we let $d=\lfloor 0.2D\rfloor$. For
$D=2^{15}$, EM and IRLS have memory issues due to the large $N=D$. AMP is
unable to yield a sensible solution (see Fig. 1(c)). In contrast, CoFEM (on
both CPU and GPU) is accurate while being faster than EM/IRLS and not
experiencing memory issues at high $D$.

####  V-C2 CoFEM vs. Seq

Finally, we compare CoFEM to the sequential algorithm. Both of these
algorithms have low space complexity (see Table I), preventing memory issues
at very high $D$. Since CoFEM and Seq have different optimization procedures,
we run both algorithms until they reach 2% NRMSE. We repeat the experimental
settings from Fig. 2(a)-(c) at higher dimensions $D=2^{p}$ for
$p=\\{12,13,\ldots,18\\}$.

For the dense dictionary (Fig. 2(d))444We could not run $D>2^{16}$ for the
dense case because there is not enough memory on our devices to store all the
entries of $\boldsymbol{\Phi}$., we observe that Seq is faster than CoFEM on
the CPU for low $D$, but the gap decreases for high $D$. On the GPU, CoFEM is
faster than Seq at moderate $D$ due to its greater ability to exploit
parallelized hardware. For the two structured dictionaries (Fig. 2(e) and
2(f)), CoFEM is faster across all settings. We further observe that for many
of the larger dimensions, Seq inevitably encounters memory issues as a
covariance-based method. In contrast, CoFEM has no such issue due to its low
space complexity. For example, at $D=2^{18}$, CoFEM can leverage the GPU to be
$5000\times$ faster than Seq. In summary, CoFEM’s ability to obviate
covariance computation provides many advantages in terms of scalability over
existing SBL inference schemes.

We emphasize that across all of our experiments, $U$ and $K$ are kept at small
values despite substantially increasing $D$, demonstrating CoFEM’s scalability
at very high dimensions.

Figure 3: Empirical insights into $U$ and $K$ based on theory of CoFEM. Figure
4: Results of running SBL-CoFEM for calcium imaging and comparison with
FOOPSI.

###  V-D Empirical Insights from the Theory of CoFEM

Finally, we give some insights into our theory for CoFEM’s hyperparameters $U$
and $K$ (Section IV). We use the DCT dictionary setup of Section V-A with
$N=\lfloor D/3\rfloor$ measurements.

In Fig. 3(a), we let $D=1024$ and plot the relationship between (1) the number
of probes $K$ and (2) the maximum standard deviation
$\max_{j=1}^{D}{\nu}^{(T)}_{j}$ over all $D$ coordinates of the diagonal
estimator at iteration $T=50$ of CoFEM. We present a “theoretical” curve
calculated using Theorem 1 with $\boldsymbol{\Theta}=\mathbb{I}$ and an
“empirical” curve in which each ${\nu}^{(T)}_{j}$ is estimated through the
empirical standard deviation of $1000$ Monte Carlo trials. Fig. 3(a) shows
that our theoretical bound (1) captures the true decay of the standard
deviation as a function of $K$ and (2) accurately upperbounds the empirical
curve after $T$ _finite_ iterations despite Theorem 1’s condition of
$T\to\infty$.

In Fig. 3(b), we consider different orders of magnitude for $D$ and plot the
number of CG steps $U$ needed to achieve small error
$\epsilon_{\text{max}}=10^{-4}$ at the $T=50$-th and $T=100$-th iterations of
CoFEM. We explore the impact of our diagonal preconditioner (Section III-C) on
$U$. The figure reveals that our theoretical insights from Section IV-B and
Theorem 2 hold in practice: (1) when there is no preconditioner and
$\epsilon_{\text{max}}$ is fixed, $U$ needs to grow substantially with
increasing iterations $T$ and/or dimensionality $D$, yet (2) when a
preconditioner is used, $U$ can be small and constant for large $T$ and large
$D$.

##  VI Real-Data Experiments

We now demonstrate the utility of CoFEM for two real data applications –
calcium deconvolution and MRI reconstruction.

###  VI-A Calcium Deconvolution

Calcium imaging is a widely used tool in neuroscience for monitoring the
electrical activity of neurons [36]. It is a method for indirectly observing
the spiking activity of a neuron through a fluorescence trace
$\boldsymbol{y}\in\mathbb{R}^{D}$, approximated as the convolution of an
intrinsically sparse spiking pattern $\boldsymbol{z}^{*}\in\mathbb{R}^{D}$
with a decaying calcium response $\boldsymbol{\phi}\in\mathbb{R}^{D}$. In
calcium deconvolution, we aim to recover $\boldsymbol{z}^{*}$ from
$\boldsymbol{y}$ and $\boldsymbol{\phi}$.

####  VI-A1 SBL Model

We can cast calcium deconvolution as a SBL problem with a dictionary
$\boldsymbol{\Phi}\in\mathbb{R}^{D\times D}$ consisting of delayed (and
truncated) versions of $\boldsymbol{\phi}$ as its columns, similar to the
setting in Section V-A2. The data is then assumed to be generated as
$\boldsymbol{y}=\boldsymbol{\Phi}\boldsymbol{z}^{*}+\boldsymbol{\varepsilon}$,
where each $\varepsilon_{i}\sim\mathcal{N}(0,1/\beta)$. Since
$\boldsymbol{z}^{*}$ for calcium deconvolution is a non-negative vector, we
employ SBL with _non-negativity constraints_ (Section III-E2).

####  VI-A2 Spike Inference

We use the CoFEM inference algorithm adapted for non-negativity constraints.
Since $\boldsymbol{\Phi}$ represents a discrete-time convolution, we can
efficiently apply $\boldsymbol{\Phi}$ to any vector $\boldsymbol{v}$ through
fast Fourier transforms and an element-wise product. Non-negative SBL yields a
rectified Gaussian posterior
$p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\hat{\alpha}})$ over the
latent spikes $\boldsymbol{z}$. To obtain a point estimate
$\boldsymbol{\hat{z}}$, we find a _filtered mode_.Specifically, we first
filter $\boldsymbol{z}$ by selecting components $z_{j}$ that are highly likely
to be non-zero, i.e. $z_{j}$ such that
$p(z_{j}=0\>|\>\boldsymbol{y},\boldsymbol{\hat{\alpha}})<q$, where $q$ is some
small percentile (e.g. $0.05$, $0.01$). This query is possible only because
SBL models uncertainty in $\boldsymbol{z}$. Setting the unselected components
of $\boldsymbol{z}$ to zero, we then find the most likely values for all
selected $z_{j}$’s according to the posterior, resulting in
$\boldsymbol{\hat{z}}$. More details can be found in Appendix C-1. This is
analogous to thresholding heuristics commonly used by $\ell_{1}$-based sparse
coding algorithms [37]. However, unlike those value-based strategies, the
percentile filtering for SBL is value-agnostic and instead operates on the
learned posterior.

####  VI-A3 Data and Hyperparameters

We apply SBL to five fluoresence traces from the GENIE dataset [38]. Each
$\boldsymbol{y}$ contains data at a sampling rate of $\nu=60$ Hz for a total
of $D=14{,}400$ time points, which is a high-dimensional problem. The response
$\boldsymbol{\phi}$ has $\phi_{i}=(1-\psi)^{i-1}$ for $\psi=1/(\nu\times
0.7)=0.0238$, a widely-used value for the calcium indicator GCaMP6f. CoFEM
employs $T=20$, $K=20$ and $U=400$. The noise precision $\beta$ is estimated
from $\boldsymbol{y}$ through a Fourier domain procedure, as described in
[39]. To obtain $\boldsymbol{\hat{z}}$, we use a filter percentile of
$q=0.05$.

####  VI-A4 Results

To evaluate $\boldsymbol{\hat{z}}$, we employ the following standard practice
[39]: The GENIE dataset has ground-truth times for neural spikes. Let
$\boldsymbol{z}^{*}\in\mathbb{R}^{D}$ be a zero-one vector indicating when
true spiking occurred. For bin length $b$, we reduce $\boldsymbol{z}^{*}$ and
$\boldsymbol{\hat{z}}$ to vectors $\boldsymbol{c}^{*}$ and
$\boldsymbol{\hat{c}}$ of length $\lceil D/b\rceil$ by summing across windows
of $b$ consecutive components. We then compute the Pearson correlation
coefficient $\rho_{b}$ between $\boldsymbol{c}^{*}$ and
$\boldsymbol{\hat{c}}$. A high value for $\rho_{b}$ indicates agreement
between $\boldsymbol{\hat{z}}$ and $\boldsymbol{z}^{*}$.

Fig. 4(a) plots sample SBL-CoFEM outputs and compares its inferred spike times
with the ground truth. Fig. 4(b) shows an averaged curve over the five traces
of $\rho_{b}$ versus $b$ at various bin lengths $b\in\\{10,20,30,40,50,60\\}$
for SBL-CoFEM and a popular $\ell_{1}$-based method called FOOPSI [40]. The
figure shows that SBL-CoFEM outperforms FOOPSI, with the gap growing for
larger bin sizes $b$.

###  VI-B Multi-Contrast MRI Reconstruction

Magnetic resonance imaging (MRI) is one of the dominant modalities for imaging
the human body [41]. The standard data acquisition practice samples a set of
points (called “$k$-space”) $\boldsymbol{k}\in\mathbb{C}^{N}$ from the two-
dimensional Fourier transform (2DFT) of the image
$\boldsymbol{x}\in\mathbb{R}^{D}$.555Though MRIs are generally complex-valued,
the data we use here are real-valued. See [6] for how to generalize SBL to the
complex case. In practice, one may aim to collect $N<D$ points to reduce the
amount of time a patient needs to remain in the scanner. However, doing so
leads to an ill-posed inverse problem
$\mathbb{M}\mathbb{F}\boldsymbol{x}=\boldsymbol{k}$ for $\boldsymbol{x}$,
where $\mathbb{F}\in\mathbb{C}^{D\times D}$ is the 2DFT and
$\mathbb{M}\in\mathbb{R}^{N\times D}$ is an undersampling operator. Thus,
compressed sensing strategies often exploit the sparsity of $\boldsymbol{x}$
with respect to some transform for accurate reconstruction.

In _multi-contrast_ MRI reconstruction, there are $L$ images
$\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{L}$ of an object that one wishes to
recover from corresponding undersampled $k$-space measurements
$\boldsymbol{k}_{1},\ldots,\boldsymbol{k}_{L}$. Bilgic et al. [6] showed that
SBL with _multi-task learning_ (Section III-E1) achieves successful joint
recovery of the multiple images. They outperformed $\ell_{1}$-based methods by
exploiting common sparsity patterns among the horizontal/vertical image
gradients (i.e. row-wise/column-wise finite differences) of
$\boldsymbol{x}_{1},\boldsymbol{x}_{2},\ldots,\boldsymbol{x}_{L}$. However,
the main drawback is the computation time for reconstruction, which is many
times slower than $\ell_{1}$ methods. We demonstrate how CoFEM can accelerate
this method while maintaining its superior performance.

Figure 5: Undersampling $k$-space masks and SBL-CoFEM reconstrutions for the
SRI24 atlas. Error maps are scaled by 15$\times$ to aid visualization.

####  VI-B1 SBL Model and Inference

For each contrast $\ell$, let
$\boldsymbol{\Phi}_{\ell}=\mathbb{M}_{\ell}\mathbb{F}$ denote the $\ell$-th
undersampled 2DFT operator with undersampling mask $\mathbb{M}_{\ell}$. The
task is to infer the sparse latent vector
$\boldsymbol{z}^{\text{horz}}_{\ell}=\partial^{\text{horz}}\boldsymbol{x}_{\ell}\in\mathbb{R}^{D}$,
where $\partial^{\text{horz}}$ denotes the horizontal image gradient operator.
Note that $\partial^{\text{horz}}$ is a convolution, implying that it
corresponds to a diagonal matrix $\boldsymbol{\Delta}^{\text{horz}}$ in the
Fourier domain. Let the observed data be
$\boldsymbol{y}^{\text{horz}}_{\ell}=\mathbb{M}_{\ell}\boldsymbol{\Delta}^{\text{horz}}\mathbb{M}_{\ell}^{\top}\boldsymbol{k}_{\ell}$
for all $\ell$. We impose a multi-task SBL model on
$(\boldsymbol{z}^{\text{horz}}_{\ell},\boldsymbol{y}^{\text{horz}}_{\ell},\boldsymbol{\Phi}_{\ell})$,
as per Eq. (16). We use a shared $\boldsymbol{\alpha}^{\text{horz}}$ to ensure
that we will learn grouped sparsity patterns among
$\\{\boldsymbol{z}_{\ell}^{\text{horz}}\\}_{\ell=1}^{L}$. The same procedure
is repeated for $\boldsymbol{z}_{\ell}^{\text{vert}}$ based on the operator
$\partial^{\text{vert}}$. Further details are in [6].

We employ CoFEM for multi-task SBL (Section III-E1) to recover
$p(\boldsymbol{z}_{\ell}^{\text{horz}}\>|\>\boldsymbol{y}^{\text{horz}}_{\ell},\boldsymbol{\hat{\alpha}}^{\text{horz}})$
and
$p(\boldsymbol{z}_{\ell}^{\text{vert}}\>|\>\boldsymbol{y}_{\ell}^{\text{vert}},\boldsymbol{\hat{\alpha}}^{\text{vert}})$
for all $\ell$. The reconstructed image $\hat{\boldsymbol{x}}_{\ell}$ are
computed from these posterior distributions. More details are given in
Appendix D-1.

TABLE II: Results on Multi-Contrast MRI Reconstruction Algorithm | NRMSE | Computation Time
---|---|---
SparseMRI | 5.4% | 22.3 min
SBL-Seq | 3.4% | 89.9 min
SBL-CoFEM | 2.9% | 2.5 min (CPU), 0.2 min (GPU)

####  VI-B2 Data and Hyperparameters

We consider the SRI24 atlas [42], a set of $L=3$ MRI contrasts with dimensions
$200\times 200$ for a total of $D=40{,}000$ pixels. For each image
$\boldsymbol{x}_{\ell}^{*}$, we undersample its 2DFT by a factor of four in
the horizontal dimension, observing $N=10{,}000$ points to form
$\boldsymbol{k}_{\ell}$. The mask $\mathbb{M}_{\ell}$ is randomly determined
according to a power rule favoring the center of $k$-space [43]. For CoFEM, we
have $T=15$, $\beta=10^{6}$, $K=8$, $\epsilon_{\text{max}}=10^{-5}$ and
$U=200$. The algorithm is not sensitive to variations in these values.

####  VI-B3 Results

Fig. 5 provides images of the masks $\mathbb{M}_{\ell}$ and SBL-CoFEM’s
reconstructions $\boldsymbol{\hat{x}}_{\ell}$. Success is measured through
NRMSE between the vectorized forms of
$\boldsymbol{\hat{x}}\in\mathbb{R}^{D\cdot L}$ and
$\boldsymbol{x}^{*}\in\mathbb{R}^{D\cdot L}$. Table II compares SBL-CoFEM
against SparseMRI ($\ell_{1}$-based compressed sensing [43]) and SBL-Seq (SBL
with sequential algorithm [6]). Although SBL-Seq has lower NRMSE than
SparseMRI, it requires high computation time. In contrast, SBL-CoFEM attains
the lowest error and can exploit the GPU to be 450$\times$ faster than SBL-
Seq.

Finally, the bottom half of Fig. 5 displays error maps of absolute differences
between $\boldsymbol{\hat{x}}_{\ell}$ and $\boldsymbol{x}^{*}_{\ell}$, along
with _variance maps_ for each image. Each variance map captures the model’s
confidence over different areas of its reconstruction; pixels with high
variance indicate more potential to deviate from the point estimate
$\boldsymbol{\hat{x}}_{\ell}$. SBL can create variance maps because it models
uncertainty; non-Bayesian methods that do not model uncertainty (e.g.
$\ell_{1}$ methods) cannot generate these maps. Appendix D-2 explains how SBL-
CoFEM can generate these variance maps using the diagonal estimation rule. The
variance maps bear similarity to the ground-truth error maps, suggesting that
SBL-CoFEM can predict its errors in reconstruction.

##  VII Conclusion

We developed covariance-free expectation-maximization (CoFEM) to accelerate
sparse Bayesian learning (SBL). By solving linear systems to circumvent matrix
inversion, CoFEM exhibits superior time-efficiency and space-efficiency over
existing SBL inference schemes, especially when the dictionary
$\boldsymbol{\Phi}$ admits fast matrix-vector multiplication. We theoretically
analyzed CoFEM’s hyperparameters, such as the number of linear systems and
number of solver steps, showing that they can remain small even at high
dimensions. Coupled with GPU acceleration, CoFEM can be up to thousands of
times faster than other SBL methods without sacrificing sparse coding
accuracy. Finally, we used CoFEM for real-data applications, showing that it
can adapt to multi-task learning and non-negativity constraints, while
enabling SBL to be competitive with non-Bayesian methods in accuracy and
scalability.

##  Appendix A Proof of Theorem 1

###### Proof.

We begin by characterizing
$\boldsymbol{\hat{\Sigma}}:=\lim_{t\to\infty}\boldsymbol{\Sigma}^{(t)}$. By
the block matrix inversion formula [44], any symmetric positive definite
matrix with diagonal blocks $\mathbb{X},\mathbb{Z}$ and off-diagonal block
$\mathbb{Y}$ and can be inverted as

$\displaystyle\begin{bmatrix}\mathbb{X}&\mathbb{Y}\\\
\mathbb{Y}^{\top}&\mathbb{Z}\end{bmatrix}^{-1}=\begin{bmatrix}\mathbb{X}^{-1}+\mathbb{X}^{-1}\mathbb{Y}\mathbb{W}\mathbb{Y}^{\top}\mathbb{X}^{-1}&-\mathbb{X}^{-1}\mathbb{Y}\mathbb{W}\\\
-\mathbb{W}\mathbb{Y}^{\top}\mathbb{X}^{-1}&\mathbb{W}\end{bmatrix},$ (34)

where
$\mathbb{W}:=(\mathbb{Z}-\mathbb{Y}^{\top}\mathbb{X}^{-1}\mathbb{Y})^{-1}$ is
the inverse Schur complement. We apply Eq. (34) to
$\boldsymbol{\Sigma}^{(t)}$, with
$\mathbb{X}:=\beta\mathbb{\Phi}^{\top}_{\mathcal{S}}\boldsymbol{\Phi}_{\mathcal{S}}+\mathrm{diag}\\{\boldsymbol{\alpha}^{(t)}_{\mathcal{S}}\\}$,
$\mathbb{Y}:=\beta\boldsymbol{\Phi}^{\top}_{\mathcal{S}}\boldsymbol{\Phi}_{\mathcal{U}}$,
and
$\mathbb{Z}:=\beta\mathbb{\Phi}^{\top}_{\mathcal{U}}\boldsymbol{\Phi}_{\mathcal{U}}+\mathrm{diag}\\{\boldsymbol{\alpha}^{(t)}_{\mathcal{U}}\\}$.
As $t\to\infty$, we have ${\alpha}^{(t)}_{j}\to\infty$ for $j\in\mathcal{U}$,
which forces $\mathbb{W}\to\boldsymbol{0}$, where $\mathbb{0}$ is the zero-
matrix. Then, by Eq. (34),

$\displaystyle\boldsymbol{\hat{\Sigma}}=\lim_{t\to\infty}\boldsymbol{\Sigma}^{(t)}=\begin{bmatrix}(\beta\mathbb{\Phi}^{\top}_{\mathcal{S}}\boldsymbol{\Phi}_{\mathcal{S}}+\mathrm{diag}\\{\boldsymbol{\hat{\alpha}}_{\mathcal{S}}\\})^{-1}&\mathbb{0}\\\
\mathbb{0}&\mathbb{0}\end{bmatrix}.$ (35)

Eq. (24) follows from applying Eq. (35) to Lemma 1 for $j\in\mathcal{U}$;
since all rows of $\boldsymbol{\hat{\Sigma}}$ that correspond to $\mathcal{U}$
are zero, the estimator’s standard deviation must also converge to zero.

We now prove Eq. (25) for an active index $j\in\mathcal{S}$. Let $\psi_{j}>0$
be any positive real number. Using Lemma 1 and Eq. (35), we can bound
$\hat{\nu}_{j}:=\lim_{t\to\infty}{\nu}^{(t)}_{j}$ with

$\displaystyle\hat{\nu}_{j}=\sqrt{\frac{1}{K}\sum_{\begin{subarray}{c}j^{\prime}\in\mathcal{T}_{j}\end{subarray}}\mathrm{\hat{\Sigma}}_{j,j^{\prime}}^{2}}\leq\frac{1}{\sqrt{K}}\sqrt{({\hat{\Sigma}}_{j,j}-\psi_{j})^{2}+\sum_{\begin{subarray}{c}j^{\prime}\in\mathcal{T}_{j}\end{subarray}}\mathrm{\hat{\Sigma}}_{j,j^{\prime}}^{2}},$
(36)

where $\mathcal{T}_{j}:=\mathcal{S}\setminus\\{j\\}$. Let
$\mathbb{\Psi}\in\mathbb{R}^{|\mathcal{S}|\times|\mathcal{S}|}$ be a diagonal
matrix with $\psi_{j}$ for all $j\in\mathcal{S}$ along its diagonal. Let
$\boldsymbol{e}_{j}$ be the standard unit vector in
$\mathbb{R}^{|\mathcal{S}|}$ corresponding to $j$. From Eq. (36), we have

$\displaystyle\hat{\nu}_{j}$
$\displaystyle\leq\tfrac{1}{\sqrt{K}}\lVert(\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}-\mathbb{\Psi})\boldsymbol{e}_{j}\rVert_{2}$
$\displaystyle\leq\tfrac{1}{\sqrt{K}}\lVert\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}-\mathbb{\Psi}\rVert_{2}=\tfrac{1}{\sqrt{K}}\lVert\mathbb{\Psi}(\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}^{-1}-\mathbb{\Psi}^{-1})\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}\rVert_{2}$
$\displaystyle\leq\tfrac{1}{\sqrt{K}}\lVert\mathbb{\Psi}(\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}^{-1}-\mathbb{\Psi}^{-1})\rVert_{2}\lVert\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}\rVert_{2}$
$\displaystyle=\tfrac{1}{\sqrt{K}}\lVert\mathbb{\Psi}(\beta\boldsymbol{\Phi}^{\top}_{\mathcal{S}}\boldsymbol{\Phi}_{\mathcal{S}}+\mathrm{diag}\\{\boldsymbol{\hat{\alpha}}_{\mathcal{S}}\\}-\mathbb{\Psi}^{-1})\rVert_{2}\lVert\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}\rVert_{2},$
(37)

where the last step uses Eq. (35) to expand
$\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}^{-1}$.

We now perform a change-of-variables: Let
$\mathbb{\Theta}\in\mathbb{R}^{|\mathcal{S}|\times|\mathcal{S}|}$ be any
diagonal matrix of positive diagonal values. We define
$\boldsymbol{\Psi}:=(\beta\boldsymbol{\Theta}+\mathrm{diag}\\{\boldsymbol{\hat{\alpha}}_{\mathcal{S}}\\})^{-1}$
and re-write Eq. (37) as

$\displaystyle\hat{\nu}_{j}$
$\displaystyle\leq\tfrac{1}{\sqrt{K}}\lVert(\beta\boldsymbol{\Theta}+\mathrm{diag}\\{\boldsymbol{\hat{\alpha}}_{\mathcal{S}}\\})^{-1}(\beta\boldsymbol{\Phi}^{\top}_{\mathcal{S}}\boldsymbol{\Phi}_{\mathcal{S}}-\beta\boldsymbol{\Theta})\rVert_{2}\lVert\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}\rVert_{2}$
$\displaystyle\leq\tfrac{1}{\sqrt{K}}\lVert(\beta\boldsymbol{\Theta})^{-1}(\beta\boldsymbol{\Phi}^{\top}_{\mathcal{S}}\boldsymbol{\Phi}_{\mathcal{S}}-\beta\boldsymbol{\Theta})\rVert_{2}\lVert\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}\rVert_{2}$
$\displaystyle=\tfrac{1}{\sqrt{K}}\lVert\boldsymbol{\Theta}^{-1}\boldsymbol{\Phi}^{\top}_{\mathcal{S}}\boldsymbol{\Phi}_{\mathcal{S}}-\mathbb{I}\rVert_{2}\lVert\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}\rVert_{2}.$
(38)

Finally, we can bound the last term of Eq. (38) with

$\displaystyle\lVert\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}\rVert_{2}=\frac{1}{\lambda_{\text{min}}(\boldsymbol{\hat{\Sigma}}_{\mathcal{S},\mathcal{S}}^{-1})}\leq\frac{1}{\beta\cdot\lambda_{\text{min}}(\boldsymbol{\Phi}_{\mathcal{S}}^{\top}\boldsymbol{\Phi}_{\mathcal{S}})},$
(39)

where we use the fact that the singular values coincide with the eigenvalues
for a symmetric positive definite matrix. ∎

##  Appendix B Proof of Theorem 2

###### Proof.

The statement follows from Eq. (27) and showing that
$\hat{\kappa}:=\lim_{t\to\infty}\kappa(\mathbb{A^{\prime}}^{(t)})\leq(1+\xi)/(1-\xi)$,
where
$\xi:=\lVert\boldsymbol{\Theta}^{-1}\boldsymbol{\Phi}_{\mathcal{S}}^{\top}\boldsymbol{\Phi}_{\mathcal{S}}-\mathbb{I}\rVert$.
Let
$\mathbb{\Delta}:=\beta\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}-\mathrm{diag}\\{\beta\boldsymbol{\theta}\\}$.
Then,

$\displaystyle\mathbb{A}^{(t)}=\mathbb{M}^{(t)}+\mathbb{\Delta}\implies\mathbb{A^{\prime}}^{(t)}=\mathbb{I}+(\mathbb{M}^{(t)})^{-\frac{1}{2}}\mathbb{\Delta}(\mathbb{M}^{(t)})^{-\frac{1}{2}}$
$\displaystyle\implies\hat{\mathbb{A}}^{\prime}:=\lim_{t\to\infty}\mathbb{A^{\prime}}^{(t)}=\mathbb{I}+\begin{bmatrix}\hat{\mathbb{M}}_{\mathcal{S},\mathcal{S}}^{-\frac{1}{2}}\mathbb{\Delta}_{\mathcal{S},\mathcal{S}}\hat{\mathbb{M}}_{\mathcal{S},\mathcal{S}}^{-\frac{1}{2}}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{0}\end{bmatrix},$ (40)

where
$\hat{\mathbb{M}}:=\mathrm{diag}\\{\beta\boldsymbol{\theta}+\boldsymbol{\hat{\alpha}}\\}$.

Our goal is to bound
$\hat{\kappa}=\lambda_{\text{max}}(\hat{\mathbb{A}}^{\prime})/\lambda_{\text{min}}(\hat{\mathbb{A}}^{\prime})$.
Eq. (40) shows that if $\eta$ is an eigenvalue of
$\hat{\mathbb{M}}_{\mathcal{S},\mathcal{S}}^{-1/2}\mathbb{\Delta}_{\mathcal{S},\mathcal{S}}\hat{\mathbb{M}}_{\mathcal{S},\mathcal{S}}^{-1/2}$,
then $1+\eta$ is an eigenvalue of $\hat{\mathbb{A}}^{\prime}$. This reduces
our task to bounding $\eta$. A matrix $\mathbb{X}\in\mathbb{R}^{D\times D}$ is
similar to another matrix $\mathbb{Y}\in\mathbb{R}^{D\times D}$ if there
exists an invertible matrix $\mathbb{Z}\in\mathbb{R}^{D\times D}$ such that
$\mathbb{Y}=\mathbb{Z}^{-1}\mathbb{X}\mathbb{Z}$, and this further implies
that $\mathbb{X}$ and $\mathbb{Y}$ have the same eigenvalues [44]. Since
$\mathbb{X}:=\hat{\mathbb{M}}_{\mathcal{S},\mathcal{S}}^{-1/2}\mathbb{\Delta}_{\mathcal{S},\mathcal{S}}\hat{\mathbb{M}}_{\mathcal{S},\mathcal{S}}^{-1/2}$
is similar to
$\mathbb{Y}:=\hat{\mathbb{M}}_{\mathcal{S},\mathcal{S}}^{-1}\mathbb{\Delta}_{\mathcal{S},\mathcal{S}}$
for $\mathbb{Z}:=\hat{\mathbb{M}}_{\mathcal{S},\mathcal{S}}^{1/2}$, it follows
that $\eta$ is also an eigenvalue of $\mathbb{Y}$. The absolute eigenvalues of
a matrix cannot exceed its spectral norm, implying

$\displaystyle|\eta|\leq\lVert\hat{\mathbb{M}}_{\mathcal{S},\mathcal{S}}^{-1}\mathbb{\Delta}_{\mathcal{S},\mathcal{S}}\rVert_{2}=\lVert\hat{\mathbb{M}}_{\mathcal{S},\mathcal{S}}^{-1}(\beta\boldsymbol{\Phi}_{\mathcal{S}}^{\top}\boldsymbol{\Phi}_{\mathcal{S}}-\beta\mathbb{\Theta})\rVert_{2}$
$\displaystyle\leq\lVert(\beta\boldsymbol{\Theta})^{-1}(\beta\boldsymbol{\Phi}_{\mathcal{S}}^{\top}\boldsymbol{\Phi}_{\mathcal{S}}-\beta\mathbb{\Theta})\rVert_{2}=\lVert\boldsymbol{\Theta}^{-1}\boldsymbol{\Phi}_{\mathcal{S}}^{\top}\boldsymbol{\Phi}_{\mathcal{S}}-\mathbb{I}\rVert_{2}.$
(41)

It follows that $\lambda_{\text{max}}(\hat{\mathbb{A}}^{\prime})\leq
1+\lVert\boldsymbol{\Theta}^{-1}\boldsymbol{\Phi}_{\mathcal{S}}^{\top}\boldsymbol{\Phi}_{\mathcal{S}}-\mathbb{I}\rVert_{2}$
and $\lambda_{\text{min}}(\hat{\mathbb{A}}^{\prime})\geq
1-\lVert\boldsymbol{\Theta}^{-1}\boldsymbol{\Phi}_{\mathcal{S}}^{\top}\boldsymbol{\Phi}_{\mathcal{S}}-\mathbb{I}\rVert_{2}$,
which bounds $\kappa(\hat{\mathbb{A}}^{\prime})$. ∎

##  Appendix C Details of SBL for Calcium Deconvolution

####  C-1 Filtered Mode

Let $\mathcal{S}\subseteq\\{1,2,\ldots,D\\}$ be the set of selected indices
after percentile filtering of
$p(\boldsymbol{z}\>|\>\boldsymbol{y},\boldsymbol{\hat{\alpha}})$ recovered by
CoFEM. Let
$\boldsymbol{\Phi}_{\mathcal{S}}\in\mathbb{R}^{N\times|\mathcal{S}|}$ denote
the sub-matrix of $\boldsymbol{\Phi}$ composed of the columns corresponding to
$\mathcal{S}$. Then, the _filtered mode_ is the solution to the following
problem:

$\displaystyle\boldsymbol{\hat{u}}=\arg\min_{\boldsymbol{u}\geq\boldsymbol{0}\in\mathbb{R}^{|\mathcal{S}|}}||\boldsymbol{y}-\boldsymbol{\Phi}_{\mathcal{S}}\boldsymbol{u}||_{2}^{2}+\sum_{j\in\mathcal{S}}(\hat{\alpha}_{j}/\beta)u_{j}^{2},$
(42)

which can be obtained using non-negative least-squares solvers. Solving Eq.
(42) is fast in practice, because it is a low-dimensional problem (i.e.
$|\mathcal{S}|\ll D$). Our point estimate solution is $\boldsymbol{\hat{z}}$,
where $\hat{z}_{j}=u_{j}$ for $j\in\mathcal{S}$ and $\hat{z}_{j}=0$ for
$j\not\in\mathcal{S}$.

##  Appendix D Details of SBL for MRI Reconstruction

####  D-1 Final Reconstruction

Let $\boldsymbol{\mu}_{\ell}^{\text{horz}}$ and
$\boldsymbol{\mu}_{\ell}^{\text{vert}}$ denote the respective means of
$p(\boldsymbol{z}_{\ell}^{\text{horz}}\>|\>\boldsymbol{y}^{\text{horz}}_{\ell},\boldsymbol{\hat{\alpha}}^{\text{horz}})$
and
$p(\boldsymbol{z}_{\ell}^{\text{vert}}\>|\>\boldsymbol{y}_{\ell}^{\text{vert}},\boldsymbol{\hat{\alpha}}^{\text{vert}})$.
These quantities are combined through solving a constrained least-squares
problem to yield a final reconstruction $\boldsymbol{\hat{x}}_{\ell}$:

$\displaystyle\boldsymbol{\hat{x}}_{\ell}=\arg$
$\displaystyle\min_{\boldsymbol{x}_{\ell}}\lVert\partial^{\text{horz}}\boldsymbol{x}_{\ell}-\boldsymbol{\mu}_{\ell}^{\text{horz}}\rVert_{2}^{2}+\lVert\partial^{\text{vert}}\boldsymbol{x}_{\ell}-\boldsymbol{\mu}_{\ell}^{\text{vert}}\rVert_{2}^{2},$
$\displaystyle\text{ s.t.
}\mathbb{M}_{\ell}\mathbb{F}\boldsymbol{x}_{\ell}=\boldsymbol{k}_{\ell}.$ (43)

We use Parseval’s Theorem [31] to cast Eq. (43) to the Fourier domain. This
converts $\partial^{\text{horz}}$ and $\partial^{\text{vert}}$ into diagonal
matrices $\mathbb{\Delta}^{\text{horz}}$ and $\mathbb{\Delta}^{\text{vert}}$
in the Fourier domain, giving an element-wise separable problem with a closed-
form solution [6].

####  D-2 Variance Map

For each MRI contrast $\ell$, SBL learns posterior distributions
$\mathcal{N}(\boldsymbol{\mu}^{\text{horz}}_{\ell},\boldsymbol{\Sigma}^{\text{horz}}_{\ell})$
and
$\mathcal{N}(\boldsymbol{\mu}^{\text{vert}}_{\ell},\boldsymbol{\Sigma}^{\text{vert}}_{\ell})$
for the image gradients. We use the fact that for
$\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ and a
matrix $\mathbb{E}$, we have
$\mathbb{E}\boldsymbol{z}\sim\mathcal{N}(\mathbb{E}\boldsymbol{\mu},\mathbb{E}\mathbb{\Sigma}\mathbb{E}^{\top})$.
The solution to Eq. (43) has the form
$\boldsymbol{\hat{x}}_{\ell}=\mathbb{E}_{1}\boldsymbol{\mu}_{\ell}^{\text{horz}}+\mathbb{E}_{2}\boldsymbol{\mu}_{\ell}^{\text{horz}}+\boldsymbol{e}$
for some matrices $\mathbb{E}_{1},\mathbb{E}_{2}\in\mathbb{R}^{D\times D}$ and
some vector $\boldsymbol{e}\in\mathbb{R}^{D}$. To obtain variance maps for
$\boldsymbol{\hat{x}}_{\ell}$, we find the diagonal elements of its covariance
matrix
$\boldsymbol{\Psi}:=\mathbb{E}_{1}\boldsymbol{\Sigma}^{\text{horz}}_{\ell}\mathbb{E}_{1}^{\top}+\mathbb{E}_{2}\boldsymbol{\Sigma}^{\text{vert}}_{\ell}\mathbb{E}_{2}^{\top}$.
We can do this by drawing probes and applying the diagonal estimation rule
(Section III-A). In doing so, we need to apply
$\boldsymbol{\Sigma}^{\text{horz}}_{\ell}=(\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\boldsymbol{\hat{\alpha}}^{\text{horz}})^{-1}$
and
$\boldsymbol{\Sigma}^{\text{vert}}_{\ell}=(\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\boldsymbol{\hat{\alpha}}^{\text{vert}})^{-1}$
to arbitrary vectors, which can be done via parallel CG (Section III-B).

## References

  * [1] D. J. MacKay, “Bayesian methods for backpropagation networks,” in _Models of neural networks III_. Springer, 1996.
  * [2] M. E. Tipping, “Sparse bayesian learning and the relevance vector machine,” _Journal of machine learning research_ , vol. 1, no. Jun, pp. 211–244, 2001\.
  * [3] S. Ji, Y. Xue, and L. Carin, “Bayesian compressive sensing,” _IEEE Transactions on signal processing_ , vol. 56, no. 6, pp. 2346–2356, 2008.
  * [4] S. Ji, D. Dunson, and L. Carin, “Multitask compressive sensing,” _IEEE Transactions on Signal Processing_ , vol. 57, no. 1, pp. 92–106, 2008.
  * [5] D. P. Wipf and B. D. Rao, “Sparse bayesian learning for basis selection,” _IEEE Transactions on Signal processing_ , vol. 52, no. 8, pp. 2153–2164, 2004.
  * [6] B. Bilgic, V. K. Goyal, and E. Adalsteinsson, “Multi-contrast reconstruction with bayesian compressed sensing,” _Magnetic resonance in medicine_ , vol. 66, no. 6, pp. 1601–1615, 2011.
  * [7] S. Liu, J. Jia, Y. D. Zhang, and Y. Yang, “Image reconstruction in electrical impedance tomography based on structure-aware sparse bayesian learning,” _IEEE transactions on medical imaging_ , vol. 37, no. 9, pp. 2090–2102, 2018\.
  * [8] P. Chen, Z. Cao, Z. Chen, and X. Wang, “Off-grid doa estimation using sparse bayesian learning in mimo radar with unknown mutual coupling,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 1, pp. 208–220, 2018.
  * [9] Y. Ma, Y. Konishi, K. Kinoshita, S. Lao, and M. Kawade, “Sparse bayesian regression for head pose estimation,” in _18th International Conference on Pattern Recognition (ICPR’06)_ , vol. 3. IEEE, 2006, pp. 507–510.
  * [10] Z. Zhang, T.-P. Jung, S. Makeig, and B. D. Rao, “Compressed sensing for energy-efficient wireless telemonitoring of noninvasive fetal ecg via block sparse bayesian learning,” _IEEE Transactions on Biomedical Engineering_ , vol. 60, no. 2, pp. 300–309, 2012.
  * [11] S. Yuan, S. Wang, M. Ma, Y. Ji, and L. Deng, “Sparse bayesian learning-based time-variant deconvolution,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 55, no. 11, pp. 6182–6194, 2017.
  * [12] O. Williams, A. Blake, and R. Cipolla, “Sparse bayesian learning for efficient visual tracking,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 27, no. 8, pp. 1292–1304, 2005.
  * [13] Q. Wu, Y. D. Zhang, M. G. Amin, and B. Himed, “Space–time adaptive processing and motion parameter estimation in multistatic passive radar using sparse bayesian learning,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 54, no. 2, pp. 944–957, 2015.
  * [14] J. Fang, Y. Shen, H. Li, and P. Wang, “Pattern-coupled sparse bayesian learning for recovery of block-sparse signals,” _IEEE Transactions on Signal Processing_ , vol. 63, no. 2, pp. 360–372, 2014.
  * [15] D. P. Wipf, S. S. Nagarajan, J. Platt, D. Koller, and Y. Singer, “A new view of automatic relevance determination.” in _NIPS_ , 2007, pp. 1625–1632.
  * [16] D. Wipf and S. Nagarajan, “Iterative reweighted $\ell_{1}$and $\ell_{2}$ methods for finding sparse solutions,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 4, no. 2, pp. 317–329, 2010.
  * [17] J. Fang, L. Zhang, and H. Li, “Two-dimensional pattern-coupled sparse Bayesian learning via generalized approximate message passing,” _IEEE Transactions on Image Processing_ , vol. 25, no. 6, pp. 2920–2930, 2016.
  * [18] H. Duan, L. Yang, J. Fang, and H. Li, “Fast inverse-free sparse bayesian learning via relaxed evidence lower bound maximization,” _IEEE Signal Processing Letters_ , vol. 24, no. 6, pp. 774–778, 2017.
  * [19] M. E. Tipping, A. C. Faul _et al._ , “Fast marginal likelihood maximisation for sparse bayesian models.” in _AISTATS_ , 2003.
  * [20] D. P. Wipf, B. D. Rao, and S. Nagarajan, “Latent variable bayesian models for promoting sparsity,” _IEEE Transactions on Information Theory_ , vol. 57, no. 9, pp. 6236–6255, 2011.
  * [21] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the em algorithm,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 39, no. 1, pp. 1–22, 1977\.
  * [22] M. Al-Shoukairi and B. Rao, “Sparse bayesian learning using approximate message passing,” in _2014 48th Asilomar Conference on Signals, Systems and Computers_. IEEE, 2014, pp. 1957–1961.
  * [23] M. Al-Shoukairi, P. Schniter, and B. D. Rao, “A gamp-based low complexity sparse bayesian learning algorithm,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 2, pp. 294–308, 2017.
  * [24] C. M. Bishop and M. Tipping, “Variational relevance vector machines,” _Uncertainty in Artificial Intelligence_ , 2000.
  * [25] S. D. Babacan, M. Luessi, R. Molina, and A. K. Katsaggelos, “Low-rank matrix completion by variational sparse bayesian learning,” in _2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2011, pp. 2188–2191.
  * [26] D. Shutin, T. Buchgraber, S. R. Kulkarni, and H. V. Poor, “Fast variational sparse bayesian learning with automatic relevance determination for superimposed signals,” _IEEE Transactions on Signal Processing_ , vol. 59, no. 12, pp. 6257–6261, 2011.
  * [27] B. Worley, “Scalable mean-field sparse bayesian learning,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 24, pp. 6314–6326, 2019.
  * [28] C. Bekas, E. Kokiopoulou, and Y. Saad, “An estimator for the diagonal of a matrix,” _Applied numerical mathematics_ , vol. 57, no. 11-12, pp. 1214–1229, 2007.
  * [29] M. R. Hestenes, E. Stiefel _et al._ , _Methods of conjugate gradients for solving linear systems_. NBS Washington, DC, 1952, vol. 49, no. 1.
  * [30] J. R. Shewchuk _et al._ , “An introduction to the conjugate gradient method without the agonizing pain,” 1994.
  * [31] A. V. Oppenheim, J. R. Buck, and R. W. Schafer, _Discrete-time signal processing. Vol. 2_. Upper Saddle River, NJ: Prentice Hall, 2001.
  * [32] D. P. Wipf and B. D. Rao, “An empirical bayesian strategy for solving the simultaneous sparse approximation problem,” _IEEE Transactions on Signal Processing_ , vol. 55, no. 7, pp. 3704–3716, 2007.
  * [33] Q. Wu, Y. D. Zhang, M. G. Amin, and B. Himed, “Complex multitask bayesian compressive sensing,” in _2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2014, pp. 3375–3379.
  * [34] A. Nalci, I. Fedorov, M. Al-Shoukairi, T. T. Liu, and B. D. Rao, “Rectified gaussian scale mixtures and the sparse non-negative least squares problem,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 12, pp. 3124–3139, 2018.
  * [35] S. Foucart and H. Rauhut, _A Mathematical Introduction to Compressive Sensing._ , 2013.
  * [36] C. Grienberger and A. Konnerth, “Imaging calcium in neurons,” _Neuron_ , vol. 73, no. 5, pp. 862–885, 2012.
  * [37] J. Friedrich, P. Zhou, and L. Paninski, “Fast online deconvolution of calcium imaging data,” _PLoS computational biology_ , vol. 13, no. 3, p. e1005423, 2017.
  * [38] J. Akerboom, T.-W. Chen, T. J. Wardill, L. Tian, J. S. Marvin, S. Mutlu, N. C. Calderón, F. Esposti, B. G. Borghuis, X. R. Sun _et al._ , “Optimization of a gcamp calcium indicator for neural activity imaging,” _Journal of neuroscience_ , vol. 32, no. 40, pp. 13 819–13 840, 2012.
  * [39] E. A. Pnevmatikakis, D. Soudry, Y. Gao, T. A. Machado, J. Merel, D. Pfau, T. Reardon, Y. Mu, C. Lacefield, W. Yang _et al._ , “Simultaneous denoising, deconvolution, and demixing of calcium imaging data,” _Neuron_ , vol. 89, no. 2, pp. 285–299, 2016.
  * [40] J. T. Vogelstein, A. M. Packer, T. A. Machado, T. Sippy, B. Babadi, R. Yuste, and L. Paninski, “Fast nonnegative deconvolution for spike train inference from population calcium imaging,” _Journal of neurophysiology_ , vol. 104, no. 6, pp. 3691–3704, 2010.
  * [41] D. G. Nishimura, _Principles of magnetic resonance imaging_. Standford Univ., 2010.
  * [42] T. Rohlfing, N. M. Zahr, E. V. Sullivan, and A. Pfefferbaum, “The sri24 multichannel atlas of normal adult human brain structure,” _Human brain mapping_ , vol. 31, no. 5, pp. 798–819, 2010.
  * [43] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid mr imaging,” _Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine_ , vol. 58, no. 6, pp. 1182–1195, 2007.
  * [44] G. Strang, _Linear algebra and its applications._ Belmont, CA: Thomson, Brooks/Cole, 2006.