tablesection algorithmsection

# Hierarchical Bayesian Inverse Problems:
A High-Dimensional Statistics Viewpoint

D. Sanz-Alonso and N. Waniorek

(University of Chicago)

###### Abstract

This paper analyzes hierarchical Bayesian inverse problems using techniques
from high-dimensional statistics. Our analysis leverages a property of
hierarchical Bayesian regularizers that we call _approximate decomposability_
to obtain non-asymptotic bounds on the reconstruction error attained by
_maximum a posteriori_ estimators. The new theory explains how hierarchical
Bayesian models that exploit sparsity, group sparsity, and sparse
representations of the unknown parameter can achieve accurate reconstructions
in high-dimensional settings.

Keywords: hiearchical Bayesian inverse problems, high dimensional statistics,
MAP estimation, non-asymptotic error bounds

AMS classification: 65M32, 62C10, 65F22

##  1 Introduction

The problem of estimating a high-dimensional parameter from a few noisy and
indirect measurements is of central importance in science and engineering [29,
18, 51] and arises, for instance, in genomics [6], spatial statistics [31],
imaging [27, 40], and data assimilation [38]. High-dimensionality —which in
many inverse problems results from the unknown parameter being a fine
discretization of a function [12, 54]— presents both statistical and
computational challenges. From a statistical perspective, consistent
estimation of a $d$-dimensional unknown from $n\ll d$ measurements is only
possible when the problem exhibits some underlying low-dimensional structure
[11, 59]. From a computational perspective, the scalability of optimization
and sampling algorithms hinges on adequately exploiting such structure [26,
52, 56, 1].

In the Bayesian approach to inverse problems, low-dimensional structure can be
imposed on the reconstruction through the choice of prior distribution [51,
36]. Hierarchical modelling [32] provides an expressive approach to prior
design. The idea is to model the parameters in the prior as random variables
distributed according to a hyperprior, thus specifying a joint prior
distribution for the unknown of interest and the prior parameters.
Conditioning on the observed data yields the posterior distribution, which
represents the Bayesian solution to the inverse problem [54]. We consider a
popular family of sparsity-promoting hierarchical models with a mean-zero
Gaussian prior and gamma hyperpriors on the prior variances. This model was
originally proposed for inverse problems in medical imaging, such as cerebral
source localization [13] and positron emission tomography [4]. Heuristically,
sparsity is promoted in the following way: the gamma hyperprior imposes the
belief that most of the variances are small except for a few outliers;
consequently, the entries of prior draws will be nearly zero except for those
with outlier variances. Beyond the original motivation in imaging, similar
hierarchical models have been applied to deconvolution problems [17],
dictionary learning [46, 61, 9], multiple measurement inverse problems [33],
adaptive meshing [8], and identification of dynamical systems [2].

Hierarchical Bayesian models with gamma hyperpriors have been extensively
studied from a computational viewpoint. For linear inverse problems, the
posterior mode —known as the _maximum a posteriori_ (MAP) estimator— can be
efficiently computed using an iterative alternating scheme that converges
linearly [14]. Path-following methods for MAP estimation were developed in
[53]. For nonlinear inverse problems, the MAP estimator can be approximated
using ensemble Kalman methods [37]. Recently, algorithms for sampling [19] and
variational inference [2] have been proposed. While this body of work
compellingly showcases the computational benefits afforded by these
hierarchical models, reconstruction error bounds are not yet available.

In this paper, we present a new statistical analysis of a family of
hierarchical Bayesian models with gamma hyperpriors, obtaining the first error
bounds on the reconstruction error of MAP estimators. Our bounds are _non-
asymptotic_ : they are applicable for any dimension and number of
measurements. We view MAP estimators as particular instances of regularized
M-estimators, that is, as being defined by minimizing an objective comprising
a loss term and a regularization term. M-estimators play a central role in the
statistical literature [30]. A canonical example is the Lasso estimator [55],
which minimizes the sum of a least-squares loss and an $\ell_{1}$ penalty that
enforces sparsity in the reconstruction of the unknown [7]. Other widely used
M-estimators include the group Lasso [3], graphical Lasso [35], sparse
generalized linear models [58], and nuclear norm penalization for matrix
completion [21, 49]. The similarities between regularized M-estimation and MAP
estimation are immediate, and in fact many M-estimators can be derived as MAP
estimators of certain posterior distributions. For instance, the Lasso
corresponds to the MAP estimator of a posterior derived from a linear forward
model with additive Gaussian noise and a Laplace prior [45]. Our new theory
for hierarchical Bayesian models builds on the general framework for
M-estimation developed in [43], which identifies two key properties of the
loss and the penalty under which accurate reconstructions are possible in
high-dimensional settings: _restricted strong convexity_ of the loss and
_decomposability_ of the regularizer. As we shall see, decomposability does
not hold for the penalties that stem from the hierarchical Bayesian models for
inverse problems that we consider. To tackle this challenge, we will introduce
a weaker notion of _approximate decomposability_ under which we obtain non-
asymptotic bounds on the deviation between the MAP estimator and the true
unknown generating the data.

### Outline and Main Contributions

Section 2 formalizes the inverse problem of interest and introduces three
guiding examples of hierarchical Bayesian models that promote sparsity, group
sparsity, and sparse representations. These examples will showcase the broad
applicability of our general theory, and will also contribute to making the
presentation accessible to a wide audience. Section 3 sets forth our unified
theoretical framework to establish reconstruction bounds leveraging
approximate decomposability of hierarchical Bayesian regularizers. Section 4
applies the general theory to our three guiding examples. To streamline the
presentation, we defer proofs to Section 5 and several appendices. We close in
Section 6 with conclusions and questions for future research.

This paper contributes to bridge the fields of inverse problems and high-
dimensional statistics. In addition to providing a unified and accessible
perspective of recent literature on inverse problems and statistical
M-estimation theory, this paper makes the following original contributions:

  * •

Theorem 2.6 proves linear convergence of an iterative alternating scheme to
compute the MAP estimator for a hierarchical model that promotes group
sparsity. Example 2.7 proposes a novel hierarchical model to promote sparse
representations in an overcomplete dictionary. We derive an iterative
alternating scheme to compute the MAP estimator, and we prove its linear
convergence in Theorem 2.8.

  * •

Definition 3.1 introduces the notion of approximate decomposability, and Lemma
3.8 highlights an important consequence of this property. Our main result,
Theorem 3.11, establishes reconstruction error bounds for M-estimators with
approximately decomposable regularizers. Our theory generalizes the one in
[43], which only holds for exactly decomposable regularizers and is not
applicable to the hierarchical models we consider.

  * •

Theorems 4.1 and 4.3 provide the first known reconstruction error bounds for
MAP estimation using sparsity-promoting conditionally Gaussian models with
gamma hyperpriors under both well-specified and misspecified settings. For
Gaussian noise models and forward maps with Gaussian columns, we sharply
characterize the reconstruction errors under well-specified and misspecified
settings in Corollaries 4.2 and 4.4. Our results show that with an appropriate
choice of hyperparameters, the MAP estimator achieves the minimax rate for
high-dimensional sparse linear regression. Analogous results for group
sparsity are established in Theorem 4.5 and Corollary 4.6, and for sparse
representations in Theorem 4.7 and Corollary 4.8.

##  2 Problem Formulation

Consider the inverse problem of reconstructing an unknown $u\in\mathbb{R}^{d}$
from data $y\in\mathbb{R}^{n}$ related by

$\displaystyle y=Au+\varepsilon,\quad\varepsilon\sim\mathcal{N}(0,I).$ (2.1)

Here, $A\in\mathbb{R}^{n\times d}$ is a given forward map and $\varepsilon$
denotes measurement noise. We assume white noise without loss of generality
—our theory can be readily extended to correlated Gaussian noise using a
whitening transformation. We are interested in high-dimensional settings where
$d\gg n.$ In this underdetermined regime, accurate reconstruction hinges on
low-dimensional structure that arises, for instance, when the unknown is
sparse or admits a sparse representation.

In the Bayesian approach to inverse problems, inference is based on the
posterior distribution, which combines the Gaussian likelihood
$y|u\sim\mathcal{N}(Au,I)$ implied by (2.1) with a prior distribution that
encodes structural beliefs on the unknown. We focus on sparsity-promoting
hierarchical priors, where the level of sparsity is controlled by an auxiliary
parameter $\theta\in\mathbb{R}^{k}.$ Specifically, we consider priors of the
form

$\pi(u,\theta)=\pi(u|\theta)\pi_{\eta}(\theta),$

where the hyperprior $\pi_{\eta}$ on $\theta$ depends on a parameter $\eta$
that will play an important role in the theory. Combining the likelihood and
the hierarchical prior, we obtain the posterior distribution

$\displaystyle\begin{split}\pi(u,\theta|y)\propto\pi(y|u)\pi(u|\theta)\pi_{\eta}(\theta)\propto\exp\bigl{(}-\mathsf{J}(u,\theta)\bigr{)},\end{split}$
(2.2)

where $\mathsf{J}$ denotes (up to an additive constant) the negative log-
density of the posterior. The _maximum a posteriori_ (MAP) estimator for
$(u,\theta)$ is then given by

$\displaystyle(\hat{u},\hat{\theta})=\operatorname*{arg\,max}_{u,\theta}\pi(u,\theta|y)=\operatorname*{arg\,min}_{u,\theta}\mathsf{J}(u,\theta).$
(2.3)

This paper conducts a high-dimensional analysis of MAP estimators for
hierarchical Bayesian models under the frequentist assumption that the data
are generated from a fixed true value $u^{\star}$ of the unknown, so that

$\displaystyle y=Au^{\star}+\varepsilon,\quad\varepsilon\sim\mathcal{N}(0,I).$
(2.4)

Our goal is to derive bounds on the reconstruction error
$\hat{\Delta}:=\hat{u}-u^{\star}$ between the MAP estimator $\hat{u}$ and the
true unknown parameter $u^{\star}$. In a unified framework, we will analyze
several hierarchical models that satisfy the following two key properties:

######  Property 2.1.

The MAP estimator $(\hat{u},\hat{\theta})$ can be efficiently computed via an
alternating optimization scheme of the form:

  1. 1.

Initialize $\theta^{0},$ $\ell=0.$

  2. 2.

Iterate until convergence

    1. (i)

$u$-update:
$u^{\ell+1}:=\operatorname*{arg\,min}_{u}{\mathsf{J}}(u,\theta^{\ell}).$

    2. (ii)

$\theta$-update:
$\theta^{\ell+1}:=\operatorname*{arg\,min}_{\theta}{\mathsf{J}}(u^{\ell+1},\theta).$

    3. (iii)

$\ell\to\ell+1.$ ∎

######  Property 2.2.

If $(\hat{u},\hat{\theta})$ is the MAP estimator, then $\hat{u}$ is the
minimizer of an objective $\mathsf{F}:\mathbb{R}^{d}\to\mathbb{R}$ of the form

$\mathsf{F}(u):=\frac{1}{2n}\|y-Au\|_{2}^{2}+\lambda\mathsf{R}_{\eta}(u),$
(2.5)

where $\lambda\mathsf{R}_{\eta}(\cdot)$ is a convex regularization term
determined by the hierarchical prior. ∎

Previous work has largely focused on Property 2.1, motivating hierarchical
Bayesian models by the simplicity and effectiveness of alternating
optimization algorithms for MAP estimation [20]. On the other hand, in this
paper we focus on Property 2.2 and show that many regularizers
$\mathsf{R}_{\eta}$ arising from hierarchical Bayesian models satisfy an
approximate decomposability property which enables accurate reconstructions
under appropriate model assumptions. We next provide three important examples
of hierarchical models that satisfy Properties 2.1 and 2.2. In Section 4, we
will apply to these three examples the unified theory introduced in Section 3.

######  Example 2.3 (Hierarchical Models for Sparsity and the IAS Algorithm).

First, we consider a popular hierarchical Bayesian model with gamma
hyperpriors to reconstruct sparse parameters, given by

$\displaystyle\begin{split}u|\theta&\sim\mathcal{N}\Bigl{(}0,\frac{\sqrt{2}}{\lambda
n}D_{\theta}\Bigr{)},\qquad D_{\theta}=\text{diag}(\theta),\\\
\theta_{j}&\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}\text{Gamma}\mathopen{}\mathclose{{}\left(\frac{\sqrt{2}\theta_{j}^{*}}{\lambda
n},\frac{3}{2}+\frac{\sqrt{2}}{2}\lambda n\eta}\right),\quad 1\leq j\leq
d,\end{split}$ (2.6)

with scale parameter $\frac{\sqrt{2}\theta_{j}^{*}}{\lambda n}>0$ and shape
parameter $\frac{3}{2}+\frac{\sqrt{2}}{2}\lambda n\eta>0.$ Here, the auxiliary
parameter $\theta\in\mathbb{R}^{k},$ with $k:=d,$ controls the variance of the
coordinates of $u,$ and thereby the level of sparsity imposed on the
reconstruction. The hierarchical Bayesian model with gamma hyperpriors in
(2.6) was proposed and investigated in [14, 20]. From (2.6), we obtain the
joint prior

$\displaystyle\pi(u,\theta)\propto\exp\mathopen{}\mathclose{{}\left(-\frac{\lambda
n}{2\sqrt{2}}\sum_{j=1}^{d}\frac{u_{j}^{2}}{\theta_{j}}-\frac{\lambda
n}{\sqrt{2}}\sum_{j=1}^{d}\mathopen{}\mathclose{{}\left[\frac{\theta_{j}}{\theta_{j}^{*}}-\eta\log\frac{\theta_{j}}{\theta_{j}^{*}}}\right]}\right).$

We assume for the simplicity of our subsequent analysis that the scale
parameters are chosen such that $\theta_{j}^{*}=\frac{n}{\|A_{j}\|_{2}^{2}},$
where $A_{j}\in\mathbb{R}^{n}$ denotes the $j^{th}$ column of the forward map.
We then make the transformation
$\mathopen{}\mathclose{{}\left(\theta_{j},u_{j},A}\right)\mapsto\mathopen{}\mathclose{{}\left(\frac{\theta_{j}}{\theta_{j}^{*}},\frac{u_{j}}{\sqrt{\theta_{j}^{*}}},AD_{\theta^{*}}^{1/2}}\right)$.
Notice that our assumption on the scale parameters is therefore tantamount to
assuming that the forward map is column-normalized so that
$\frac{\|A_{j}\|_{2}}{\sqrt{n}}=1$ and that $\theta_{j}^{*}=1$ for each $1\leq
j\leq d.$ As such, we consider a posterior with negative log-density given, up
to an additive constant, by

$\mathsf{J}(u,\theta):=\mathrlap{\underbrace{\phantom{\frac{1}{2n}\|y-Au\|_{2}^{2}+\frac{\lambda}{2\sqrt{2}}\sum_{j=1}^{d}\frac{u_{j}^{2}}{\theta_{j}}}}_{(a)}}\frac{1}{2n}\|y-Au\|_{2}^{2}+\overbrace{\frac{\lambda}{2\sqrt{2}}\sum_{j=1}^{d}\frac{u_{j}^{2}}{\theta_{j}}+\frac{\lambda}{\sqrt{2}}\sum_{j=1}^{d}\mathopen{}\mathclose{{}\left(\theta_{j}-\eta\log\theta_{j}}\right)}^{(b)}\,\,.$
(2.7)

The iterative alternating scheme (IAS) computes the MAP estimator for this
objective by minimizing the least-squares objective (a) for the $u$-update and
minimizing $(b)$ for the $\theta$-update. Given $(u^{\ell},\theta^{\ell})$
both updates have a closed form solution:

$\displaystyle\begin{split}u^{\ell+1}&:=\mathopen{}\mathclose{{}\left(\frac{A^{T}A}{n}+\frac{\lambda}{\sqrt{2}}D_{\theta^{\ell}}^{-1}}\right)^{-1}\frac{A^{T}y}{n},\\\
\theta_{j}^{\ell+1}&:=\frac{\eta}{2}+\sqrt{\frac{\eta^{2}}{4}+\frac{(u_{j}^{\ell+1})^{2}}{2}},\qquad
1\leq j\leq d.\end{split}$ (2.8)

For large problems, the $u$-update can be efficiently computed using conjugate
gradient with early stopping based on Morozov’s discrepancy principle [14, 15,
20]. The closed form expression for the $\theta$-update was derived in [20].
The IAS algorithm (2.8) was introduced in [14, 20] and has been widely
studied, see e.g. [14, 15, 20, 16, 17]. The following known result verifies
Property 2.2 and summarizes other key computational properties.

######  Theorem 2.4.

For $\eta>0,$ the objective function ${\mathsf{J}}$ in (2.7) is strictly
convex over $\mathbb{R}^{d}\times\mathbb{R}^{d}_{+}$ and thus admits a unique
MAP estimator $(\hat{u},\hat{\theta}).$ Moreover, the $u$-iterates $u^{\ell}$
of the IAS algorithm (2.8) converge at least linearly to $\hat{u}$ in the
Mahalanobis norm $\|\cdot\|_{D_{\hat{\theta}}}.$ Furthermore, $\hat{u}$ is the
global minimizer of $\mathsf{F}$ in (2.5) with

$\displaystyle\mathsf{R}_{\eta}(u):=\frac{1}{\sqrt{2}}\sum_{j=1}^{d}\biggl{[}\frac{u_{j}^{2}}{2f_{j}(u)}+f_{j}(u)-\eta\log
f_{j}(u)\biggr{]},\quad\quad
f_{j}(u):=\frac{\eta}{2}+\sqrt{\frac{\eta^{2}}{4}+\frac{u_{j}^{2}}{2}}\,.$
(2.9)

As illustrated by Theorem 2.4, previous works focus on the convexity of
${\mathsf{J}}$ and the convergence of IAS iterates $u^{\ell}$ to the MAP
estimator $\hat{u}.$ In this work we focus instead on understanding the
convergence of the MAP estimator $\hat{u}$ to the true parameter $u^{\star}$
generating the data. ∎

######  Example 2.5 (Hierarchical Models for Group Sparsity and the GS-IAS
Algorithm).

The second hierarchical Bayesian model we consider promotes group sparsity. A
particular case of the model where each group consists of precisely three
indices was proposed in [14], and the general setting of arbitrary group sizes
was later formulated in [9]. We assume the a priori knowledge that components
of the unknown $u$ are partitioned into groups $\\{g_{1},\ldots,g_{k}\\}$
where $|g_{j}|=p_{j}$. The model also postulates a conditionally Gaussian
prior on each group and the variance of each group is modeled via a gamma
hyperprior,

$\displaystyle\begin{split}u_{g_{j}}|\theta_{j}&\sim\mathcal{N}\Bigl{(}0,\frac{\sqrt{2}\theta_{j}}{\lambda
n}C_{j}\Bigr{)},\qquad C_{j}\succ 0,\\\
\theta_{j}&\sim\text{Gamma}\mathopen{}\mathclose{{}\left(\frac{\sqrt{2}\theta_{j}^{*}}{\lambda
n},\frac{p_{j}+2}{2}+\frac{\sqrt{2}}{2}\lambda n\eta}\right),\quad 1\leq j\leq
k,\end{split}$ (2.10)

with scale parameter $\frac{\sqrt{2}\theta_{j}^{*}}{\lambda n}>0$ and shape
parameter $\frac{p_{j}+2}{2}+\frac{\sqrt{2}}{2}\lambda n\eta$. Here, the
auxiliary parameter $\theta\in\mathbb{R}^{k}$ controls the level of group
sparsity in the reconstructed solution. From (2.10), we have the joint prior
distribution

$\displaystyle\pi(u,\theta)\propto\exp\mathopen{}\mathclose{{}\left(-\frac{\lambda
n}{2\sqrt{2}}\sum_{j=1}^{k}\frac{\|u_{g_{j}}\|_{C_{j}}^{2}}{\theta_{j}}-\frac{\lambda
n}{\sqrt{2}}\sum_{j=1}^{k}\mathopen{}\mathclose{{}\left[\frac{\theta_{j}}{\theta_{j}^{*}}-\eta\log\frac{\theta_{j}}{\theta_{j}^{*}}}\right]}\right).$

To simplify our analysis, we assume that the scale parameters are chosen such
that $\theta_{j}^{*}=\frac{n}{\|A_{g_{j}}\|_{C_{j},2}^{2}}$, where
$A_{g_{j}}\in\mathbb{R}^{n\times p_{j}}$ denotes the matrix obtained by
keeping the columns from $A$ in group $g_{j},$ and
$\|A_{g_{j}}\|_{C_{j},2}=\max_{\|v\|_{C_{j}}=1}\|A_{g_{j}}\|_{2}.$ Then, via
the transformation
$(\theta_{j},u_{g_{j}},A_{g_{j}})\mapsto\mathopen{}\mathclose{{}\left(\frac{\theta_{j}}{\theta_{j}^{*}},\frac{u_{j}}{\sqrt{\theta_{j}^{*}}},A_{g_{j}}\sqrt{\theta_{j}^{*}}}\right)$
we see that our assumption on the scale parameters is equivalent to assuming
the forward map is block-normalized so that
$\frac{\|A_{g_{j}}\|_{C_{j},2}}{\sqrt{n}}=1$ with $\theta_{j}^{*}=1.$ Applying
Bayes rule, we consider a posterior with negative log-density given, up to an
additive constant, by

$\mathsf{J}(u,\theta):=\mathrlap{\underbrace{\phantom{\frac{1}{2n}\|y-Au\|_{2}^{2}+\frac{\lambda}{\sqrt{2}}\|u\|_{D_{\theta}}^{2}}}_{(a)}}\frac{1}{2n}\|y-Au\|_{2}^{2}+\overbrace{\frac{\lambda}{2\sqrt{2}}\|u\|_{D_{\theta}}^{2}+\frac{\lambda}{\sqrt{2}}\sum_{j=1}^{k}\mathopen{}\mathclose{{}\left(\theta_{j}-\eta\log\theta_{j}}\right)}^{(b)}\,\,,$
(2.11)

where $D_{\theta}$ is the block-diagonal matrix given by

$\displaystyle D_{\theta}:=\begin{bmatrix}\theta_{1}C_{1}&&\\\ &\ddots&\\\
&&\theta_{k}C_{k}\end{bmatrix}.$

To compute the MAP estimator for this objective, we again use an alternating
optimization scheme, minimizing the least-squares objective $(a)$ for the
$u$-update and minimizing $(b)$ for the $\theta$-update. Given
$(u^{\ell},\theta^{\ell})$, these updates have closed form solutions:

$\displaystyle\begin{split}u^{\ell+1}&:=\mathopen{}\mathclose{{}\left(\frac{A^{T}A}{n}+\frac{\lambda}{\sqrt{2}}D_{\theta^{\ell}}^{-1}}\right)\frac{A^{T}y}{n},\\\
\theta_{j}^{\ell+1}&:=\frac{\eta}{2}+\sqrt{\frac{\eta^{2}}{4}+\frac{\|u_{g_{j}}\|_{C_{j}}^{2}}{2}},\qquad
1\leq j\leq k.\end{split}$ (2.12)

The group-sparse iterative alternating scheme (GS-IAS) described by (2.12) was
derived in [14] and it was shown that the iterates converge to the minimizer
of 2.11, although no rate of convergence was derived. This objective and the
GS-IAS algorithm admit a number of computationally convenient properties,
summarized in the following theorem.

######  Theorem 2.6.

For $\eta>0,$ the objective function $\mathsf{J}$ in $\eqref{eq:GS
functional}$ is strictly convex over $\mathbb{R}^{d}\times\mathbb{R}^{k}_{+},$
and thus admits a unique MAP estimator $(\hat{u},\hat{\theta}).$ Moreover, the
$u$-iterates $u^{\ell}$ of the GS-IAS algorithm (2.12) converge at least
linearly to $\hat{u}$ in the Mahalanobis norm $\|\cdot\|_{D_{\hat{\theta}}}$.
Furthermore, $\hat{u}$ is the global minimizer of $\mathsf{F}$ in (2.5) with

$\displaystyle\mathsf{R}_{\eta}(u):=\frac{1}{\sqrt{2}}\sum_{j=1}^{k}\biggl{[}\frac{\|u_{g_{j}}\|^{2}_{C_{j}}}{2f_{j}(u)}+f_{j}(u)-\eta\log
f_{j}(u)\biggr{]},\qquad
f_{j}(u):=\frac{\eta}{2}+\sqrt{\frac{\eta^{2}}{4}+\frac{\|u_{g_{j}}\|_{C_{j}}^{2}}{2}}.$
(2.13)

For completeness, we include a proof of Theorem 2.6 in Appendix A. ∎

######  Example 2.7 (Hierarchical Models for Sparse Representations and the
O-IAS Algorithm).

Our final motivating example is a novel hierarchical Bayesian model to promote
the sparse representation of the unknown in some redundant dictionary,
$W\in\mathbb{R}^{d\times k}$. In this model, we assume the dictionary to be
overcomplete, meaning that $k\gg d$. For simplicity, we assume that $W$ is a
tight frame [60] normalized so that $WW^{T}=I$. Relevant examples of such
dictionaries include Gabor frames [41], curvelet frames [23], and wavelet
frames [50]. To encode our a priori information that the coefficients of $u$
in $W$ are sparse, we place a Gaussian prior on $W^{T}u$ and gamma hyperpriors
on the variances,

$\displaystyle\begin{split}W^{T}u|\theta&\sim\mathcal{N}\Bigl{(}0,\frac{\sqrt{2}}{\lambda
n}D_{\theta}\Bigr{)},\qquad D_{\theta}=\text{diag}(\theta),\\\
\theta_{j}&\sim\text{Gamma}\mathopen{}\mathclose{{}\left(\frac{\sqrt{2}\theta_{j}^{*}}{\lambda
n},\frac{3}{2}+\frac{\sqrt{2}}{2}\lambda n\eta}\right),\quad 1\leq j\leq
k,\end{split}$ (2.14)

with scale parameters $\frac{\sqrt{2}\theta_{j}^{*}}{\lambda n}$ and shape
parameters $\frac{3}{2}+\frac{\sqrt{2}}{2}\lambda n\eta$ for $\eta,\lambda>0$.
Here, the auxiliary parameter $\theta\in\mathbb{R}^{k}$ controls the variance
of the coefficients of $u$ in the dictionary $W$, and consequently the level
of sparsity imposed upon the coefficients. From (2.14) we have the joint prior
distribution

$\displaystyle\pi(W^{T}u,\theta)\propto\exp\mathopen{}\mathclose{{}\left(-\frac{\lambda
n}{2\sqrt{2}}\sum_{j=1}^{k}\frac{(W^{T}u)_{j}^{2}}{\theta_{j}}-\frac{\lambda
n}{\sqrt{2}}\sum_{j=1}^{k}\mathopen{}\mathclose{{}\left[\frac{\theta_{j}}{\theta_{j}^{*}}-\eta\log\frac{\theta_{j}}{\theta_{j}^{*}}}\right]}\right).$

For the convenience of our analysis in the sequel, we assume the scale
parameters are chosen such that
$\theta_{j}^{*}=\frac{n}{\|(AW)_{j}\|_{2}^{2}}$. Making
$\mathopen{}\mathclose{{}\left(\theta_{j},(W^{T}u)_{j},(AW)_{j}}\right)\mapsto\mathopen{}\mathclose{{}\left(\frac{\theta_{j}}{\theta_{j}^{*}},\frac{(W^{T}u)_{j}}{\sqrt{\theta_{j}^{*}}},(AW)_{j}\sqrt{\theta_{j}^{*}}}\right),$
we see that our assumption on the scale parameters is equivalent to assuming
that the forward map is normalized so that
$\frac{\|(AW)_{j}\|_{2}}{\sqrt{n}}=1$ with $\theta_{j}^{*}=1$. Consequently,
we consider a posterior with negative log-density given, up to an additive
constant, by

$\mathsf{J}(u,\theta):=\mathrlap{\underbrace{\phantom{\frac{1}{2n}\|y-Au\|_{2}^{2}+\frac{\lambda}{2\sqrt{2}}\sum_{j=1}^{k}\frac{(W^{T}u)_{j}^{2}}{\theta_{j}}}}_{(a)}}\frac{1}{2n}\|y-Au\|_{2}^{2}+\overbrace{\frac{\lambda}{2\sqrt{2}}\sum_{j=1}^{k}\frac{(W^{T}u)_{j}^{2}}{\theta_{j}}+\frac{\lambda}{\sqrt{2}}\sum_{j=1}^{k}\mathopen{}\mathclose{{}\left(\theta_{j}-\eta\log\theta_{j}}\right)}^{(b)}\,\,.$
(2.15)

To compute the MAP estimator, we use the following overcomplete iterative
alternating scheme (O-IAS), which minimizes $(a)$ for the $u$-update and $(b)$
for the $\theta$-update:

$\displaystyle\begin{split}u^{\ell+1}&:=\mathopen{}\mathclose{{}\left(\frac{A^{T}A}{n}+\frac{\lambda}{\sqrt{2}}WD_{\theta^{\ell}}^{-1}W^{T}}\right)\frac{A^{T}y}{n},\\\
\theta_{j}^{\ell+1}&:=\frac{\eta}{2}+\sqrt{\frac{\eta^{2}}{4}+\frac{(W^{T}u)_{j}^{2}}{2}},\qquad
1\leq j\leq k.\end{split}$ (2.16)

The objective (2.15) and the accompanying O-IAS algorithm (2.16) admit a
number of key computationally convenient properties.

######  Theorem 2.8.

For $\eta>0$, the objective function ${\mathsf{J}}$ in (2.15) is strictly
convex over $\mathbb{R}^{d}\times\mathbb{R}^{k}_{+}$, and thus admits a unique
MAP estimator $(\hat{u},\hat{\theta})$. Moreover, the $u$-iterates $u^{\ell}$
of the O-IAS algorithm converge at least linearly to $\hat{u}$ in the
Mahalanobis norm $\|\cdot\|_{(WD_{\hat{\theta}}^{-1}W^{T})^{-1}}.$
Furthermore, $\hat{u}$ is the global minimizer of $\mathsf{F}$ in (2.5) with

$\displaystyle\mathsf{R}_{\eta}(u)=\frac{1}{\sqrt{2}}\sum_{j=1}^{k}\mathopen{}\mathclose{{}\left[\frac{(W^{T}u)_{j}}{2f_{j}(u)}+f_{j}(u)-\eta\log
f_{j}(u)}\right],\quad\quad
f_{j}(u)=\frac{\eta}{2}+\sqrt{\frac{\eta^{2}}{4}+\frac{(W^{T}u)_{j}^{2}}{2}}.$
(2.17)

For completeness, we include a proof of Theorem 2.8 in Appendix A. ∎

##  3 Unified Theoretical Framework

This section develops a framework to obtain reconstruction bounds for
hierarchical Bayesian models under appropriate model assumptions. In
Subsection 3.1 we introduce the notion of approximate decomposabilility, and
show that this property is satisfied by the regularizers in the examples in
Section 2. Then, in Subsection 3.2 we present the main result, Theorem 3.11.

###  3.1 Approximate Decomposability

We introduce a model subspace $\mathcal{M}\subset\mathbb{R}^{d}$ to capture
the sparsity constraints of the problem. We refer to the orthogonal complement
$\mathcal{M}^{\perp}$ of the model subspace as the perturbation subspace.
Following [43, 59], we say that a regularizer
$\mathsf{R}:\mathbb{R}^{d}\to\mathbb{R}^{+}$ is decomposable with respect to
$\mathcal{M}$ if

$\displaystyle\mathsf{R}(u+v)=\mathsf{R}(u)+\mathsf{R}(v),\qquad\forall
u\in\mathcal{M},\forall v\in\mathcal{M}^{\perp}.$ (3.1)

For norm-based regularizers, it necessarily holds that
$\mathsf{R}(u+v)\leq\mathsf{R}(u)+\mathsf{R}(v).$ Decomposability requires
that equality holds for perturbations $v\in\mathcal{M}^{\perp}$ away from
vectors $u\in\mathcal{M},$ thus ensuring that deviations from the model
subspace along the perturbation subspace are heavily penalized. The
regularizers $\mathsf{R}_{\eta}$ discussed in Section 2 determined by the
hierarchical priors are not decomposable over a natural model subspace.
However, in the limit as $\eta\to 0$, each of these regularizers approach a
decomposable regularizer. Motivated by this observation, we introduce the
following notion of approximate decomposability.

######  Definition 3.1.

A family of regularizers $\mathsf{R}_{\eta}:\mathbb{R}^{d}\to\mathbb{R}^{+}$
parameterized by $\eta>0$ is approximately decomposable with respect to
$\mathcal{M}$ if there exists a norm-based regularizer $\mathsf{R}$
decomposable with respect to $\mathcal{M}$ and non-negative functions
$c_{i}^{L}(\eta)$ and $c_{i}^{U}(\eta),$ $i=1,2,$ such that, for any
$u\in\mathbb{R}^{d},$

$\displaystyle\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)}\right)\mathsf{R}(u)-c_{2}^{L}(\eta)\leq\mathsf{R}_{\eta}(u)\leq\mathopen{}\mathclose{{}\left(1+c_{1}^{U}(\eta)}\right)\mathsf{R}(u)+c_{2}^{U}(\eta),$
(3.2)

where $c_{i}^{L}(\eta)\to 0$ and $c_{i}^{U}(\eta)\to 0$ as $\eta\to 0$.

The functions $c_{i}^{L}$ and $c_{i}^{U}$ quantify how fast the regularizer
$\mathsf{R}_{\eta}(u)$ converges to the decomposable regularizer
$\mathsf{R}(u)$ as $\eta$ is taken to be close to zero. We remark that other
generalizations of decomposability, such as weakly decomposable norms [57] and
atomic norms [25, 22], have previously been studied. Our notion of approximate
decomposability differs most significantly from these other generalizations in
that it does not require the regularizer itself be a norm, merely that the
regularizer be close to a norm for small $\eta$. We now discuss the choice of
model subspace and establish approximate decomposability for the three example
models in Section 2.

######  Example 3.2 (Approximate Decomposability: Hierarchical Models for
Sparsity).

The natural modelling assumption for the IAS algorithm is that of sparsity.
Let $S\subset\\{1,\ldots,d\\}$ be a subset of indices with cardinality
$|S|=s.$ We define the model subspace

$\displaystyle\mathcal{M}(S):=\mathopen{}\mathclose{{}\left\\{u\in\mathbb{R}^{d}:u_{j}=0\text{
for all }j\notin S}\right\\}.$ (3.3)

Any vector in the model subspace $\mathcal{M}(S)$ has at most $s$ nonzero
entries. In this setting, the perturbation subspace is the space of vectors
with support contained entirely outside of $S$:

$\displaystyle\mathcal{M}^{\perp}(S):=\mathopen{}\mathclose{{}\left\\{u\in\mathbb{R}^{d}:u_{j}=0\text{
for all }j\in S}\right\\}.$

One can verify that, for any subset of indices $S$, the $\ell_{1}$ norm
$\mathsf{R}(u)=\|u\|_{1}$ is decomposable with respect to $\mathcal{M}(S)$.
The following lemma shows that the regularizer $\mathsf{R}_{\eta}$ defined in
(2.9) is approximately decomposable with respect to $\mathcal{M}(S).$

######  Lemma 3.3.

Consider, for $0<\eta<\frac{1}{2},$ the family of regularizers
$\mathsf{R}_{\eta}$ in (2.9). For any $u\in\mathbb{R}^{d},$

$\Bigl{(}1-\frac{\eta}{2}\Bigr{)}\|u\|_{1}-\frac{d\eta}{\sqrt{2}}\leq\mathsf{R}_{\eta}(u)\leq\|u\|_{1}+\frac{d\eta}{\sqrt{2}}(1-\log\eta).$
(3.4)

The proof can be found in Appendix B. In the notation of Definition 3.1,
approximate decomposability holds with $c_{1}^{L}(\eta)=\frac{\eta}{2}$,
$c_{2}^{L}(\eta)=\frac{d\eta}{\sqrt{2}}$, $c_{1}^{U}(\eta)=0,$ and
$c_{2}^{U}(\eta)=\frac{d\eta}{\sqrt{2}}(1-\log\eta).$ ∎

######  Example 3.4 (Approximate Decomposability: Hierarchical Models for
Group Sparsity).

We partition the index set into groups $\mathcal{G}=\\{g_{1},\ldots g_{k}\\}$
where each $g_{j}\subset\\{1,\ldots,d\\}$, $g_{i}\cap g_{j}=\emptyset$ for
$i\neq j$, and $\cup_{j=1}^{k}g_{j}=\\{1,\ldots,d\\}.$ We denote the size of
each group as $|g_{j}|=p_{j}.$ For any set of group indices
$S_{\mathcal{G}}\subset\mathcal{G}$, we define the model subspace

$\displaystyle\mathcal{M}(S_{\mathcal{G}}):=\mathopen{}\mathclose{{}\left\\{u\in\mathbb{R}^{d}:u_{g}=0\text{
for all }g\notin S_{\mathcal{G}}}\right\\}.$

This is precisely the subspace of vectors supported only on groups in the
index set $S_{\mathcal{G}}$. Conversely, the perturbation subspace is the
space of vectors with support contained in groups not in $S_{\mathcal{G}}$:

$\displaystyle\mathcal{M}^{\perp}(S_{\mathcal{G}}):=\mathopen{}\mathclose{{}\left\\{u\in\mathbb{R}^{d}:u_{g}=0\text{
for all }g\in S_{\mathcal{G}}}\right\\}.$

A decomposable regularizer for this model subspace is the group sparse norm
[62], given by

$\displaystyle\mathsf{R}(u)=\sum_{j=1}^{k}\|u_{g_{j}}\|_{C_{j}},$

where $C_{j}\in\mathbb{R}^{p_{j}\times p_{j}}$ is a positive definite matrix.
Indeed, we see that, for any $u\in\mathcal{M}(S_{\mathcal{G}})$ and
$v\in\mathcal{M}^{\perp}(S_{\mathcal{G}}),$ it holds that

$\displaystyle\mathsf{R}(u+v)=\sum_{g_{j}\in\mathcal{G}}\|u_{g_{j}}\|_{C_{j}}+\sum_{g_{j}\notin\mathcal{G}}\|v_{g_{j}}\|_{C_{j}}=\mathsf{R}(u)+\mathsf{R}(v).$

The following lemma shows that the regularizer $\mathsf{R}_{\eta}$ defined in
(2.13) is approximately decomposable with respect to
$\mathcal{M}(S_{\mathcal{G}})$.

######  Lemma 3.5.

Consider, for $0<\eta<\frac{1}{2},$ the family of regularizers
$\mathsf{R}_{\eta}$ in (2.13). For any $u\in\mathbb{R}^{d},$

$\Bigl{(}1-\frac{\eta}{2}\Bigr{)}\sum_{j=1}^{k}\|u_{g_{j}}\|_{C_{j}}-\frac{k\eta}{\sqrt{2}}\leq\mathsf{R}_{\eta}(u)\leq\sum_{j=1}^{k}\|u_{g_{j}}\|_{C_{j}}+\frac{k\eta}{\sqrt{2}}(1-\log\eta).$

The proof can be found in Appendix B. In the notation of Definition 3.1,
approximate decomposability holds with $c_{1}^{L}(\eta)=\frac{\eta}{2}$,
$c_{2}^{L}(\eta)=\frac{k\eta}{\sqrt{2}}$, $c_{1}^{U}(\eta)=0$, and
$c_{2}^{U}(\eta)=\frac{k\eta}{\sqrt{2}}(1-\log\eta).$ ∎

######  Example 3.6 (Approximate Decomposability: Hierarchical Models for
Sparse Representations).

We consider the situation where the data generating parameter itself is not
sparse, but admits a sparse representation in an overcomplete dictionary,
$W\in\mathbb{R}^{d\times k},$ where $k\gg d$. As is common in the compressed
sensing literature [24], we assume that $W$ is a tight frame such that
$WW^{T}=I.$ We let $S_{W}\subset\\{1,\ldots,k\\}$ be a subset of indices of
cardinality $|S_{W}|=s$. We then define the model subspace

$\displaystyle\mathcal{M}(S_{W}):=\mathopen{}\mathclose{{}\left\\{u\in\mathbb{R}^{d}:(W^{T}u)_{j}=0\text{
for all }j\notin S_{W}}\right\\}.$ (3.5)

This is precisely the subspace of vectors whose representation in $W$ has
support contained in $S_{W}$. The perturbation subspace is then given by

$\displaystyle\mathcal{M}^{\perp}(S_{W}):=\mathopen{}\mathclose{{}\left\\{u\in\mathbb{R}^{d}:(W^{T}u)_{j}=0\text{
for all }j\in S_{W}}\right\\}.$

A decomposable regularizer for this model subspace is given by

$\displaystyle\mathsf{R}(u)=\|W^{T}u\|_{1},$

since, for any $u\in\mathcal{M}(S_{W})$ and $v\in\mathcal{M}^{\perp}(S_{W}),$
it holds that

$\displaystyle\mathsf{R}(u+v)=\|W^{T}u+W^{T}v\|_{1}=\|W^{T}u\|_{1}+\|W^{T}v\|_{1}=\mathsf{R}(u)+\mathsf{R}(v).$

The following lemma shows that the regularizer $\mathsf{R}_{\eta}$ defined in
(2.17) is approximately decomposable with respect to $\mathcal{M}(S_{W}).$

######  Lemma 3.7.

Consider, for $0<\eta<\frac{1}{2},$ the family of regularizers
$\mathsf{R}_{\eta}$ in (2.17). For any $u\in\mathbb{R}^{d}$,

$\Bigl{(}1-\frac{\eta}{2}\Bigr{)}\|W^{T}u\|_{1}-\frac{k\eta}{\sqrt{2}}\leq\mathsf{R}_{\eta}(u)\leq\|W^{T}u\|_{1}+\frac{k\eta}{\sqrt{2}}(1-\log\eta).$

The proof can be found in Appendix B. In the notation of Definition 3.1,
approximate decomposability holds with $c_{1}^{L}(\eta)=\frac{\eta}{2}$,
$c_{2}^{L}(\eta)=\frac{k\eta}{\sqrt{2}},$ $c_{1}^{U}(\eta)=0$, and
$c_{2}^{U}(\eta)=\frac{k\eta}{\sqrt{2}}(1-\log\eta).$ ∎

We now show how the notion of approximate decomposability can be leveraged to
control the error vector $\hat{\Delta}=\hat{u}-u^{\star}$. To do so, we recall
that the dual norm of $\mathsf{R}$ with respect to the Euclidean inner product
is defined by

$\mathsf{R}^{*}(u):=\sup_{v\in\mathbb{R}^{d}\setminus\\{0\\}}\frac{\langle
u,v\rangle}{\mathsf{R}(v)}=\sup_{\mathsf{R}(v)\leq 1}\langle u,v\rangle.$

In the following lemma, we show that if the regularization parameter $\lambda$
is taken to be large enough relative to the dual norm of the noise vector,
then the reconstruction error lies in a specific set. For vector
$u\in\mathbb{R}^{d}$ and subspace $\mathcal{M}\subset\mathbb{R}^{d},$ we
denote by $u_{\mathcal{M}}$ the Euclidean projection of $u$ onto
$\mathcal{M}.$

######  Lemma 3.8.

Let $\hat{u}$ be a minimizer of the objective

$\displaystyle\mathsf{F}(u)=\frac{1}{2n}\|y-Au\|_{2}^{2}+\lambda\mathsf{R}_{\eta}(u),$
(3.6)

where $\mathsf{R}_{\eta}$ is an approximately decomposable regularizer with
respect to $\mathcal{M}$. If
$\lambda\geq\frac{2}{n}\mathsf{R}^{*}\mathopen{}\mathclose{{}\left(A^{T}\varepsilon}\right)$
and $\eta$ is small enough so that $c_{1}^{L}(\eta)<\frac{1}{4}$, then the
error $\hat{\Delta}=\hat{u}-u^{\star}$ lies in the set

$\displaystyle\mathbb{C}_{\eta}(\mathcal{M}):=\mathopen{}\mathclose{{}\left\\{\Delta\in\mathbb{R}^{d}:\mathsf{R}(\Delta_{\mathcal{M}^{\perp}})\leq
7\mathsf{R}(\Delta_{\mathcal{M}})+8\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})+4c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})+4c_{2}(\eta)}\right\\},$

where $c_{i}(\eta):=c_{i}^{L}(\eta)+c_{i}^{U}(\eta).$

The proof can be found in Appendix C. The fact that the error lies in the set
$\mathbb{C}_{\eta}(\mathcal{M})$ allows us to control the components of the
error outside the model subspace by the components of the error within the
model subspace along with terms that get smaller as the true vector is assumed
to lie within the model subspace and as $\eta$ is taken to be smaller. In the
next subsection we leverage this lemma to get a bound on the reconstruction
error.

###  3.2 General Bound on the Reconstruction Error

To bound the reconstruction error, we will need to quantify the cost of
translating between bounds on a vector in $\mathcal{M}$ given by
$\|\cdot\|_{2}$ and bounds given by $\mathsf{R}(\cdot).$ To that end, we will
rely on the notion of subspace Lipschitz constant.

######  Definition 3.9.

For a subspace $\mathcal{M}\in\mathbb{R}^{d}$, the subspace Lipschitz constant
with respect to $(\mathsf{R},\|\cdot\|_{2})$ is given by

$\displaystyle\Psi(\mathcal{M}):=\sup_{u\in\mathcal{M}\backslash\\{0\\}}\frac{\mathsf{R}(u)}{\|u\|_{2}}.$

We will also need a condition on the forward map. In the high-dimensional
setting, the likelihood term $\|y-Au\|_{2}^{2}$, necessarily has $d-n$
directions in which it has no curvature up to second order. Consequently, we
need an assumption to guarantee that in the relevant directions, chiefly the
directions given by vectors in $\mathbb{C}_{\eta}(\mathcal{M})$, the
likelihood has significant curvature at $u^{\star}.$ One such assumption is
that of restricted strong convexity.

######  Definition 3.10.

For a given regularizer $\mathsf{R}$, a linear map
$A:\mathbb{R}^{d}\to\mathbb{R}^{n}$ satisfies the restricted strong convexity
(RSC) condition with curvature $\kappa$ and tolerance $\tau^{2}$ if

$\displaystyle\frac{\|A\Delta\|_{2}^{2}}{2n}\geq\frac{\kappa}{2}\|\Delta\|_{2}^{2}-\tau^{2}\mathsf{R}^{2}(\Delta),\qquad\forall\Delta\in\mathbb{R}^{d}.$
(3.7)

If the matrix $A$ were not rank degenerate, then the RSC condition would hold
with $\tau=0,$ which recovers the usual definition of strong convexity. In the
high-dimensional setting, this cannot be the case, and a nonzero tolerance
term is required. In general, we expect the tolerance to be decreasing in the
number $n$ of observations.

Given these definitions, we now state our main result, which generalizes the
bounds in [59, Chapter 9] to hold for approximately decomposable regularizers.
The proof can be found in Appendix C.

######  Theorem 3.11.

Let $u^{\star}$ be the true parameter generating the data via (2.4) and let
$\hat{u}$ be the minimizer of the objective (3.6). Assume that
$\mathsf{R}_{\eta}$ is approximately decomposable with respect to
$\mathcal{M}$ and that $\eta$ is small enough such that
$c_{1}^{L}(\eta)<\frac{1}{4}$ and that the forward map $A$ satisfies the RSC
condition with curvature $\kappa$ and tolerance $\tau^{2}$. Then, if
$\tau\Psi(\mathcal{M})\leq\frac{\sqrt{\kappa}}{32}$ and $\lambda\geq
2\mathsf{R}^{*}(\frac{1}{n}A^{T}\varepsilon),$ we have the bound

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim\frac{\lambda^{2}}{\kappa^{2}}\Psi^{2}(\mathcal{M})+\frac{\tau^{2}}{\kappa}\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})^{2}+\frac{\lambda}{\kappa}\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})+\mathcal{E}_{\eta},$
(3.8)

where

$\mathcal{E}_{\eta}:=\frac{\tau^{2}}{\kappa}\Bigl{(}c_{1}(\eta)^{2}\mathsf{R}(u^{\star}_{\mathcal{M}})^{2}+c_{2}(\eta)^{2}\Bigr{)}+\frac{\lambda}{\kappa}\Bigl{(}c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})+c_{2}(\eta)\Bigr{)},$

and $c_{i}(\eta):=c_{i}^{L}(\eta)+c_{i}^{U}(\eta),\,i=1,2.$

We emphasize that this bound is deterministic, and makes no specific
distributional assumptions on the noise or forward map. The first three terms
are precisely the terms that appear in the bound proven in [59] for exactly
decomposable regularizers. The first term,

$\displaystyle\mathcal{E}_{\text{est}}:=\frac{\lambda^{2}}{\kappa^{2}}\Psi^{2}(\mathcal{M}),$

can be thought of as the statistical estimation error which captures the
inherent difficulty of estimating a parameter in $\mathcal{M}$. The second and
third terms together,

$\displaystyle\mathcal{E}_{\text{app}}:=\frac{\tau^{2}}{\kappa}\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})^{2}+\frac{\lambda}{\kappa}\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}}),$

represent the approximation error. These terms capture the modelling error
incurred when the data-generating parameter does not lie exactly in
$\mathcal{M}$, thus quantifying the amount of model misspecification. Notice
that if $u^{\star}\in\mathcal{M}$, then these terms vanish. The final term,

$\mathcal{E}_{\eta}:=\frac{\tau^{2}}{\kappa}\Bigl{(}c_{1}(\eta)^{2}\mathsf{R}(u^{\star}_{\mathcal{M}})^{2}+c_{2}(\eta)^{2}\Bigr{)}+\frac{\lambda}{\kappa}\Bigl{(}c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})+c_{2}(\eta)\Bigr{)},$

is novel to our bound, and captures the error incurred by the approximate
decomposability of the regularizer. This term can be made arbitrarily small
with $\eta$, which is user-specified.

##  4 Reconstruction Error for Sparsity-Promoting Models

In this section we apply the unified theoretical framework set forth in
Theorem 3.11 to each of the examples of hierarchical models discussed in
Section 2. The results in this section provide the first known reconstruction
error bounds for these models. Proofs will be deferred to Section 5.

###  4.1 Hierarchical Models for Sparsity

Here we apply our general theory to the sparsity-promoting hierarchical model
presented in Example 2.3 under both a hard-sparsity assumption and a weaker
$\ell^{q}$-sparsity assumption. For both of these types of sparsity, we first
prove a deterministic result guaranteeing that under the RSC condition, the
MAP estimator will be close to the true data generating parameter as long as
we regularize enough relative to the noise level. We emphasize that these
results do not make specific assumptions on the noise distribution or the
forward map, besides the RSC condition. Applying the deterministic results to
the probabilistic setting of Gaussian noise and Gaussian forward maps, we
derive sharp rates in terms of only the dimension, sample size, and the
parameter $\eta$. These results hold with high probability.

The first assumption we consider is that of hard sparsity. We assume that the
true parameter generating the data has support of cardinality $s$, where $s\ll
d$. The following result holds for any $s$-sparse $u^{\star}$ so long as the
forward map satisfies the RSC condition and the regularization parameter is
taken to be large enough.

######  Theorem 4.1.

Let $u^{\star}$ be the true parameter generating data via (2.1) and suppose
that $\text{supp}(u^{\star})\subset S$ where $|S|=s.$ Let
$(\hat{u},\hat{\theta})$ be the unique minimizer of the objective function
(2.7) with $0<\eta<1/2$. Assume that
$\lambda\geq\frac{2\|A^{T}\varepsilon\|_{\infty}}{n}$ and that $A$ satisfies
the RSC condition (3.7) with curvature $\kappa$ and tolerance $\tau^{2}$.
Then, if $\tau\sqrt{s}\leq\frac{\sqrt{\kappa}}{32}$ it holds that

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim
s\frac{\lambda^{2}}{\kappa^{2}}+\mathcal{E}_{\eta},$ (4.1)

where

$\mathcal{E}_{\eta}:=\frac{\tau^{2}}{\kappa}\mathopen{}\mathclose{{}\left(\eta^{2}\|u^{\star}_{S}\|_{1}^{2}+\frac{d^{2}\eta^{2}}{2}(2-\log\eta)^{2}}\right)+\frac{\lambda}{\kappa}\mathopen{}\mathclose{{}\left(\eta\|u^{\star}_{S}\|_{1}+\frac{d\eta}{\sqrt{2}}(2-\log\eta)}\right).$

We emphasize that this theorem is deterministic. It makes no specific
distributional assumptions on the noise or forward map, only requiring that
the regularization parameter is large enough relative to the noise and that
the forward map satisfies the RSC condition. The bound illustrates that, as
long as the condition $\lambda\geq\frac{2\|A^{T}\varepsilon\|_{\infty}}{n}$ is
satisfied, we would like to take the regularization parameter $\lambda,$ the
RSC tolerance $\tau$, and the parameter $\eta$ to be as small as possible. We
remark that $\eta$ is user supplied and can be made arbitrarily small, but
doing so may slow down the convergence of the IAS algorithm. The other terms
in the bound, $\lambda$ and $\tau$, inherently depend on the noise and the
forward map. As an illustrative example, we explicitly characterize $\lambda$
and $\tau$ under the assumptions of Gaussian noise and a Gaussian forward map
in the following corollary:

######  Corollary 4.2.

Let the forward map $A\in\mathbb{R}^{n\times d}$ have rows distributed i.i.d.
according to $\mathcal{N}(0,\Sigma)$, and assume the columns of $A$ are
normalized such that $\frac{\|A_{j}\|_{2}}{\sqrt{n}}=1$ for all
$j=1,\ldots,d.$ Assume the noise is given by $\varepsilon\sim\mathcal{N}(0,I)$
and that $(\hat{u},\hat{\theta})$ is the unique minimizer of the IAS objective
functional (2.7) with $\lambda=4\sqrt{\frac{\log d}{n}}$ and $0<\eta<1/2$. If
$u^{\star}$ is $s$-sparse, then there exists a constant $c_{\Sigma}$ depending
only on $\Sigma$, and universal constants $c$ and $c^{\prime}$ such that if
$n>c_{\Sigma}s\log d$ it holds with probability at least
$1-c\exp(-c^{\prime}n\lambda^{2})$ that

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim c_{\Sigma}\frac{s\log
d}{n}+\min\mathopen{}\mathclose{{}\left\\{\varphi(\eta),\varphi(\eta)^{2}}\right\\},$

where

$\displaystyle\varphi(\eta):=\sqrt{\frac{\log
d}{n}}\mathopen{}\mathclose{{}\left(\eta\|u^{\star}_{S}\|_{1}+\frac{d\eta}{\sqrt{2}}(2-\log\eta)}\right).$

Notice that the term $\varphi(\eta)$ can be made arbitrarily small by choosing
a small enough $\eta$. Taking $\eta$ small enough to match the first term in
the bound, the MAP estimator for the hierarchical model in Example 2.3
satisfies the non-asymptotic bound

$\displaystyle\|\hat{u}-u^{\star}\|_{2}\lesssim\sqrt{c_{\Sigma}\frac{s\log
d}{n}},$

which agrees with known non-asymptotic bounds for the Lasso estimator [7, 42].
The assumptions of Gaussianity in this corollary are largely for theoretical
convenience, and are not strictly necessary. Verifying the RSC condition for
the forwards maps of interest in the inverse problems community remains a
direction of further research.

Instead of assuming the true parameter is exactly $s$-sparse, one can make the
more general assumption of $\ell^{q}$ sparsity. That is, we assume

$\displaystyle
u^{\star}\in\mathbb{B}_{q}(R_{q}):=\mathopen{}\mathclose{{}\left\\{u\in\mathbb{R}^{d}:\sum_{i=1}^{d}|u_{i}|^{q}\leq
R_{q}}\right\\},$

where $R_{q}>0$ controls the level of sparsity. Notice that if $q=0$ and
$R_{q}=s$, we recover the hard-sparsity assumption. For $q>0$ this assumption
enforces a decay rate on the entries of $u^{\star}.$ Under $\ell^{q}$
sparsity, we choose our model subspace to be vectors with support
$S_{\delta}$, where

$\displaystyle
S_{\delta}:=\mathopen{}\mathclose{{}\left\\{j\in\\{1,\ldots,d\\}:|u_{j}^{\star}|>\delta}\right\\}.$
(4.2)

This model subspace is well defined for any positive $\delta$. We will see
that our bound on the reconstruction error depends on $\delta$. Note that we
define $\delta$ only as a theoretical tool, since the model subspace is not
involved methodologically when computing the MAP estimator, $\hat{u}$. We
emphasize that in the setting of $\ell^{q}$ sparsity, we have that
$u^{\star}_{S^{C}_{\delta}}\neq 0$ in general, and consequently the model
error terms in the general bounds are non-negligible. We now present a
deterministic result analogous to Theorem 4.1 under the $\ell^{q}$-sparsity
assumption.

######  Theorem 4.3.

Let $u^{\star}$ be the true parameter generating data via (2.1) and suppose
that $u^{\star}\in\mathbb{B}_{q}(R_{q}).$ Let $(\hat{u},\hat{\theta})$ be the
unique minimizer of the objective function (2.7) with $0<\eta<1/2$. Assume
that $\lambda\geq\frac{2\|A^{T}\varepsilon\|_{\infty}}{n}$ and that $A$
satisfies the RSC condition (3.7) with curvature $\kappa$ and tolerance
$\tau^{2}$. Then, for any $\delta>0$, if
$\tau\sqrt{R_{q}\delta^{-q}}\leq\frac{\sqrt{\kappa}}{32}$ it holds that

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim\frac{\lambda^{2}}{\kappa^{2}}R_{q}\delta^{-q}+\frac{\tau^{2}}{\kappa}R_{q}\delta^{1-q}+\frac{\lambda}{\kappa}R_{q}\delta^{1-q}+\mathcal{E}_{\eta},$

where

$\mathcal{E}_{\eta}:=\frac{\tau^{2}}{\kappa}\Bigl{(}\eta^{2}R_{q}^{2/q}+d^{2}\eta^{2}(2-\log\eta)^{2}\Bigr{)}+\frac{\lambda}{\kappa}\Bigl{(}\eta
R_{q}^{1/q}+d\eta(2-\log\eta)\Bigr{)}.$

Since this deterministic error bound holds for any positive $\delta,$ Theorem
4.3 provides a family of bounds parameterized by $\delta.$ Under noise and
forward map assumptions that yield particular values for $\lambda$, $\kappa$,
and $\tau$, one can optimize $\delta$ to yield the sharpest bound. We do so in
the following corollary under the assumptions of Gaussian noise and a Gaussian
forward map.

######  Corollary 4.4.

Let the forward map $A\in\mathbb{R}^{n\times d}$ have rows distributed i.i.d.
according to $\mathcal{N}(0,\Sigma)$, and assume the columns of $A$ are
normalized such that $\frac{\|A_{j}\|_{2}}{\sqrt{n}}=1$ for all
$j=1,\ldots,d.$ Assume the noise is given by $\varepsilon\sim\mathcal{N}(0,I)$
and that $(\hat{u},\hat{\theta})$ is the unique minimizer of the IAS objective
functional (2.7) with $\lambda=4\sqrt{\frac{\log d}{n}}$ and $0<\eta<1/2$. If
$u^{\star}\in\mathbb{B}_{q}(R_{q})$, then there exists a constant $c_{\Sigma}$
depending only on $\Sigma$, and universal constants $c$ and $c^{\prime}$ such
that if $n>c_{\Sigma}\log dR_{q}^{1-q/2}$ it holds with probability at least
$1-c\exp(-c^{\prime}n\lambda^{2})$ that

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim
c_{\Sigma}R_{q}\mathopen{}\mathclose{{}\left(\frac{\log
d}{n}}\right)^{1-\frac{q}{2}}+c_{\Sigma}\min\mathopen{}\mathclose{{}\left\\{\varphi(\eta),\varphi(\eta)^{2}}\right\\},$

where

$\displaystyle\varphi(\eta):=\sqrt{\frac{\log
d}{n}}\mathopen{}\mathclose{{}\left(\eta
R_{q}^{\frac{1}{q}}+\frac{d\eta}{\sqrt{2}}(2-\log\eta)}\right).$

This corollary suggests that $\eta$ should be taken small enough so that
$\varphi(\eta)$ matches the rate of
$R_{q}\mathopen{}\mathclose{{}\left(\frac{\log d}{n}}\right)^{1-\frac{q}{2}}$.
In particular, for any $\eta$ such that

$\displaystyle\eta\log\frac{1}{\eta}\leq\min\mathopen{}\mathclose{{}\left\\{\frac{R_{q}}{d}\mathopen{}\mathclose{{}\left(\frac{\log
d}{n}}\right)^{\frac{1-q}{2}},R_{q}^{1-\frac{1}{q}}\mathopen{}\mathclose{{}\left(\frac{\log
d}{n}}\right)^{\frac{1-q}{2}}}\right\\},$

the MAP estimator of the hierarchical Bayesian model in Example 2.3 achieves
up to constant factors the minimax rate for estimation under $\ell^{q}$
sparsity [48], given by

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim
R_{q}\mathopen{}\mathclose{{}\left(\frac{\log d}{n}}\right)^{1-\frac{q}{2}}.$

That is, with this choice of $\eta$, the MAP estimator achieves the best
possible reconstruction error.

###  4.2 Hierarchical Models for Group Sparsity

In this subsection we derive non-asymptotic bounds on the reconstruction error
for the group-sparsity-promoting hierarchical model considered in Example 2.5.
We prove reconstruction bounds under a weak $\ell^{q}$-group-sparsity
assumption. We suppose that

$\displaystyle
u^{\star}\in\mathbb{B}_{q}(R_{q},\mathcal{G}):=\mathopen{}\mathclose{{}\left\\{u\in\mathbb{R}^{d}:\sum_{j=1}^{k}\|u_{g_{j}}\|_{C_{j}}^{q}\leq
R_{q}}\right\\}.$

We emphasize that in the case where $q=0$ and $R_{q}=s$, we recover the hard-
group-sparsity assumption that $u_{g}=0$ for all $g\in S_{\mathcal{G}}$ with
$|S_{\mathcal{G}}|=s.$ We consider the set of group indices with norm larger
than $\delta>0$, defined by

$\displaystyle
S_{\mathcal{G},\delta}:=\mathopen{}\mathclose{{}\left\\{g_{j}\in\mathcal{G}:\|u_{g_{j}}^{\star}\|_{C_{j}}>\delta}\right\\}.$
(4.3)

We now present a deterministic result that holds for any positive $\delta$ so
long as the regularization parameter is taken to be large enough relative to
the noise and the forward map satisfies an RSC condition.

######  Theorem 4.5.

Let $u^{\star}$ be the true parameter generating the data via (2.4) and
suppose that $u^{\star}\in\mathbb{B}_{q}(R_{q},\mathcal{G})$. Let
$(\hat{u},\hat{\theta})$ be the minimizer of the objective function (2.11)
with $0<\eta<1/2.$ Assume that
$\lambda\geq\max_{g_{j}\in\mathcal{G}}\frac{2}{n}\|(A^{T}\varepsilon)_{g_{j}}\|_{C_{j}^{-1}},$
and that $A$ satisfies the RSC condition (3.7) with curvature $\kappa$ and
tolerance $\tau^{2}$. Then, for any $\delta>0$, if
$\tau\sqrt{R_{q}\delta^{-q}}\leq\frac{\sqrt{\kappa}}{32}$ it holds that

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim\frac{\lambda^{2}}{\kappa^{2}}R_{q}\delta^{-q}+\frac{\tau^{2}}{\kappa}R_{q}\delta^{1-q}+\frac{\lambda}{\kappa}R_{q}\delta^{1-q}+\mathcal{E}_{\eta},$

where

$\displaystyle\mathcal{E}_{\eta}:=\frac{\tau^{2}}{\kappa}\mathopen{}\mathclose{{}\left(\eta^{2}R_{q}^{\frac{2}{q}}+k^{2}\eta^{2}(2-\log\eta)^{2}}\right)+\frac{\lambda}{\kappa}\mathopen{}\mathclose{{}\left(\eta
R_{q}^{\frac{1}{q}}+k\eta(2-\log\eta)}\right).$

We remark that if $q=0$ and $R_{q}=s$, this theorem recovers bounds for the
hard-sparsity assumption that $u^{\star}\in\mathcal{M}(S_{\mathcal{G}})$:

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim
s\frac{\lambda^{2}}{\kappa^{2}}+\mathcal{E}_{\eta},$

where

$\displaystyle\mathcal{E}_{\eta}:=\frac{\tau^{2}}{\kappa}\mathopen{}\mathclose{{}\left(\eta^{2}\mathopen{}\mathclose{{}\left(\sum_{j\in
S_{\mathcal{G}}}\|u^{\star}_{g_{j}}\|_{C_{j}}}\right)^{2}+k^{2}\eta^{2}(2-\log\eta)^{2}}\right)+\frac{\lambda}{\kappa}\mathopen{}\mathclose{{}\left(\eta\sum_{j\in
S_{\mathcal{G}}}\|u^{\star}_{g_{j}}\|_{C_{j}}+k\eta(2-\log\eta)}\right).$

As was the case for the models discussed in the previous section, under
distributional assumptions on the noise and forward map, we can use this
deterministic bound to sharply characterize non-asymptotic bounds on the
reconstruction error with high probability. For simplicity of presentation, we
assume that $C_{j}=I$ for $1\leq j\leq k$. This assumption is straightforward
to remove, but doing so is notationally cumbersome. We require a block-
normalization condition on $A$ that generalizes the column-normalization
assumption used in the previous subsection. Given a group
$g_{j}\in\mathcal{G}$ of size $p_{j}$, we consider the submatrix
$A_{g_{j}}:\mathbb{R}^{n\times p_{j}}$, and require that

$\displaystyle\frac{\|A_{g_{j}}\|_{op}}{\sqrt{n}}=1,\quad
j\in\\{1,\ldots,k\\}.$ (4.4)

As remarked in Example 2.5, this condition can be guaranteed by appropriately
choosing the model hyperparameters. For the following corollary, we introduce
the notation $p_{\max}:=\max_{j\in\\{1,\ldots,k\\}}p_{j}.$

######  Corollary 4.6.

Let the forward map $A\in\mathbb{R}^{n\times d}$ have rows distributed i.i.d.
according to $\mathcal{N}(0,\Sigma)$ and assume that the block normalization
condition (4.4) holds. Assume that the noise is given by
$\varepsilon\sim\mathcal{N}(0,I)$ and that $(\hat{u},\hat{\theta})$ is the
unique minimizer of the GS-IAS objective functional $\eqref{eq:GS functional}$
with
$\lambda=\mathopen{}\mathclose{{}\left(\sqrt{\frac{p_{\max}}{n}}+\sqrt{\frac{\log
k}{n}}}\right)$ and $0<\eta<1/2.$ If
$u^{\star}\in\mathbb{B}_{q}(R_{q},\mathcal{G}),$ then there exists a constant
$c_{\Sigma}$ depending only on $\Sigma$ and such that if
$\sqrt{n}>c_{\Sigma}R_{q}^{1-q/2}\mathopen{}\mathclose{{}\left(\sqrt{p_{\max}}+\sqrt{\log
k}}\right)$ it holds with probability at least $1-2/k^{2}$ that

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim
c_{\Sigma}R_{q}\mathopen{}\mathclose{{}\left(\frac{p_{\max}}{n}+\frac{\log
k}{n}}\right)^{1-\frac{q}{2}}+\min\mathopen{}\mathclose{{}\left\\{\varphi(\eta),\varphi(\eta)^{2}}\right\\},$

where

$\displaystyle\varphi(\eta):=\sqrt{\frac{p_{\max}+\log
k}{n}}\mathopen{}\mathclose{{}\left(\eta
R_{q}^{\frac{1}{q}}+k\eta(2-\log\eta)}\right).$

If $q=0$ and $R_{q}=s$, we recover a bound under a hard-group-sparsity
assumption:

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim\frac{sp_{\max}}{n}+\frac{s\log
k}{n}+\max\mathopen{}\mathclose{{}\left\\{\varphi(\eta),\varphi(\eta)^{2}}\right\\},$

where

$\displaystyle\varphi(\eta):=\sqrt{\frac{p_{\max}+\log
k}{n}}\mathopen{}\mathclose{{}\left(\eta\sum_{j=1}^{k}\|u^{\star}_{g_{j}}\|_{C_{j}}+k\eta(2-\log\eta)}\right).$

In both the $\ell^{q}$-group-sparsity and hard-group-sparsity settings, by
taking $\eta$ to be small enough to match the estimation rate in the first
half of the bound, we see that the hierarchical model in Example 2.5 and the
GS-IAS algorithm can achieve the known bounds for reconstructing group sparse
vectors [43, 39, 34].

###  4.3 Hierarchical Models for Sparse Representations

Finally, we apply our theory to derive non-asymptotic reconstruction error
bounds for the hierarchical model promoting sparse representations described
in Example 2.7. We assume that the data-generating unknown $u^{\star}$ has a
representation in a tight from $W\in\mathbb{R}^{d\times k}$, where $k$ may be
much larger than $d$. We assume this representation lies in an $\ell^{q}$
ball:

$\displaystyle
u^{\star}\in\mathbb{B}_{q}(R_{q},W):=\mathopen{}\mathclose{{}\left\\{u\in\mathbb{R}^{d}:\sum_{j=1}^{k}|(W^{T}u)_{j}|^{q}\leq
R_{q}}\right\\}.$

As in the previous examples, we recover the hard-sparsity assumption that the
support of $W^{T}u^{\star}$ has cardinality $s$ in the case where $q=0$ and
$R_{q}=s$. We define our model subspace as in (3.5), taking the index set to
be

$\displaystyle
S_{W,\delta}:=\mathopen{}\mathclose{{}\left\\{j\in\\{1,\ldots,k\\}:(W^{T}u^{\star})_{j}>\delta}\right\\}.$
(4.5)

We show a deterministic bound on the reconstruction error that holds for any
positive $\delta$ provided that the regularization parameter is large enough
relative to the noise and that the forward map satisfies the RSC condition.

######  Theorem 4.7.

Let $u^{\star}$ be the true parameter generating the data via (2.4) and
suppose that $u^{\star}\in\mathbb{B}_{q}(R_{q},W).$ Let
$(\hat{u},\hat{\theta})$ be the minimizer of the objective function (2.15)
with $0<\eta<1/2.$ Assume that
$\lambda\geq\frac{2\|W^{T}A^{T}\varepsilon\|_{\infty}}{n}$ and that $A$
satisfies the RSC condition $\eqref{eq:RSC}$ with curvature $\kappa$ and
tolerance $\tau^{2}$. Then, for any $\delta>0$, if
$\tau\sqrt{R_{q}\delta^{-q}}\leq\frac{\kappa}{32}$ it holds that

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim\frac{\lambda^{2}}{\kappa^{2}}R_{q}\delta^{-q}+\frac{\tau^{2}}{\kappa}R_{q}\delta^{1-q}+\frac{\lambda}{\kappa}R_{q}\delta^{1-q}+\mathcal{E}_{\eta},$

where

$\displaystyle\mathcal{E}_{\eta}:=\frac{\tau^{2}}{\kappa}\mathopen{}\mathclose{{}\left(\eta^{2}R_{q}^{2/q}+k^{2}\eta^{2}(2-\log(\eta)^{2}}\right)+\frac{\lambda}{\kappa}\mathopen{}\mathclose{{}\left(\eta
R_{q}^{1/q}+k\eta(2-\log\eta)}\right).$

When $q=0$ and $R_{q}=s$ this result recovers a bound under the hard-sparsity
assumption that $\text{supp}(W^{T}u^{\star})=S$ with $|S|=s:$

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim
s\frac{\lambda^{2}}{\kappa^{2}}+\mathcal{E}_{\eta},$

where

$\displaystyle\mathcal{E}_{\eta}:=\frac{\tau^{2}}{\kappa}\mathopen{}\mathclose{{}\left(\eta^{2}\|(W^{T}u^{\star})_{S}\|_{1}^{2}+k^{2}\eta^{2}(2-\log\eta)^{2}}\right)+\frac{\lambda}{\kappa}\mathopen{}\mathclose{{}\left(\eta\|(W^{T}u^{\star})_{S}\|_{1}+k\eta(2-\log\eta)}\right).$

Once again, by making assumptions on the distribution of the noise and forward
map, we can derive sharp non-asymptotic bounds on the reconstruction error
with high probability.

######  Corollary 4.8.

Let the forward map $A\in\mathbb{R}^{n\times d}$ have rows distributed i.i.d.
according to $\mathcal{N}(0,\Sigma)$ and assume that $A$ is normalized such
that $\frac{\|(AW)_{j}\|_{2}}{\sqrt{n}}=1$ for all $j=1,\ldots,k$. Assume the
noise is given by $\varepsilon\sim\mathcal{N}(0,I)$ and that
$(\hat{u},\hat{\theta})$ is the minimizer of the O-IAS objective functional
(2.15) with $\lambda=4\sqrt{\frac{\log k}{n}}$ and $0<\eta<1/2.$ If
$u^{\star}\in\mathbb{B}_{Q}(R_{q},Q)$, then there exists a constant
$c_{\Sigma}$ depending only on $\Sigma$ and universal constants $c$ and
$c^{\prime}$ such that if $n>c_{\Sigma}R_{q}^{1-q/2}\log k$ it holds with
probability at least $1-c\exp(c^{\prime}n\lambda^{2})$ that

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim
R_{q}\mathopen{}\mathclose{{}\left(\frac{\log
k}{n}}\right)^{1-\frac{q}{2}}+\min\mathopen{}\mathclose{{}\left\\{\varphi(\eta),\varphi(\eta)^{2}}\right\\},$

where

$\displaystyle\varphi(\eta):=\sqrt{\frac{\log
k}{n}}\mathopen{}\mathclose{{}\left(\eta
R_{q}^{1/q}+k\eta(2-\log\eta)}\right).$

The corresponding hard-sparsity bound under the assumption that $q=0$ and
$R_{q}=s$ is

$\displaystyle\|\hat{u}-u^{\star}\|_{2}^{2}\lesssim\frac{s\log
k}{n}+\max\mathopen{}\mathclose{{}\left\\{\varphi(\eta),\varphi(\eta)^{2}}\right\\},$

where

$\displaystyle\varphi(\eta):=\sqrt{\frac{\log
k}{n}}\mathopen{}\mathclose{{}\left(\eta\|W^{T}u^{\star}\|_{1}+k\eta(2-\log\eta)}\right).$

##  5 Proof of Reconstruction Error for Sparsity-Promoting Models

This section contains the proofs of all the results in Section 4.

###  5.1 Hierarchical Models for Sparsity

###### Proof of Theorem 4.1.

The result follows from applying Theorem 3.11. From Lemma 3.3 we have that
approximate decomposability holds with $c_{1}(\eta)=\frac{\eta}{2}$ and
$c_{2}(\eta)=\frac{d\eta}{\sqrt{2}}(2-\log\eta)$. For the model subspace
$\mathcal{M}(S)$ defined in (3.3), the subspace Lipschitz constant is given by

$\displaystyle\Psi\bigl{(}\mathcal{M}(S)\bigr{)}=\sup_{u\in\mathcal{M}(S)}\frac{\|u\|_{1}}{\|u\|_{2}}=\sqrt{s}.$

The dual of the one norm is the infinity norm, so the assumption that
$\lambda\geq\frac{2\|A^{T}\varepsilon\|_{\infty}}{n}$ satisfies the condition
of Theorem 3.11. Since we have assumed $u^{\star}\in\mathcal{M}(S)$, we have
that $u^{\star}_{\mathcal{M}^{\perp}}=0,$ and thus
$\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})=0$. Plugging each of these pieces
into the bound (3.8) gives the desired result. ∎

###### Proof of Corollary 4.2.

First we verify that the given choice of $\lambda$ satisfies the assumption on
Theorem 4.1 with high probability. It is shown in [43, Corollary 2] that under
the column normalization assumption, for any $\varepsilon$ that is
(sub-)Gaussian with identity covariance, it holds that

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left(\frac{\|A^{T}\varepsilon\|_{\infty}}{n}\geq
t}\right)\leq 2\exp\mathopen{}\mathclose{{}\left(-\frac{nt^{2}}{2}+\log
d}\right).$

Taking $t^{2}=\frac{4\log d}{n},$ we conclude that the choice of
$\lambda=4\sqrt{\frac{\log d}{n}}$ satisfies the assumption of Theorem 4.1
with probability $1-c_{1}\exp(-c_{2}n\lambda^{2}).$ Next we characterize the
curvature and tolerance in the RSC condition for a Gaussian forward map. By
[59, Theorem 7.16], for any matrix $A\in\mathbb{R}^{n\times d}$ with rows
i.i.d. from $\mathcal{N}(0,\Sigma)$ there exist positive constants
$c_{1}<1<c_{2}$ such that, with probability at least
$1-\frac{e^{-n/32}}{1-e^{-n/32}},$

$\displaystyle\frac{\|A\Delta\|_{2}^{2}}{2n}\geq
c_{1}\lambda_{\min}(\Sigma)\|\Delta\|_{2}^{2}-c_{2}\rho(\Sigma)^{2}\frac{\log
d}{n}\|\Delta\|_{1}^{2},\qquad\forall\Delta\in\mathbb{R}^{d},$

where $\rho(\Sigma)^{2}:=\max_{1\leq j\leq d}\Sigma_{jj}$. In other words, the
RSC condition holds with curvature $\kappa=c_{1}\lambda_{\min}(\Sigma)$ and
tolerance $\tau^{2}=c_{2}\rho(\Sigma)^{2}\frac{\log d}{n}$. Consequently, to
guarantee that $\tau\sqrt{s}\leq\frac{\sqrt{\kappa}}{32}$ as assumed in
Theorem 4.1, we require that

$\displaystyle\sqrt{n}\geq
32\sqrt{\frac{\frac{c_{1}}{c_{2}}\rho(\Sigma)^{2}s\log
d}{\lambda_{\min}(\Sigma)}}.$

Given this requirement on $n$ we have satisfied the hypotheses of Theorem 4.1.
Consequently, we plug in our characterizations of $\lambda$, $\kappa$, and
$\tau$ into (4.1) and the desired result follows. ∎

###### Proof of Theorem 4.3.

From Lemma 3.3 we have that approximate decomposability holds with
$c_{1}(\eta)=\frac{\eta}{2}$ and
$c_{2}(\eta)=\frac{d\eta}{\sqrt{2}}(2-\log\eta)$. For the model subspace
$\mathcal{M}(S_{\delta})$ where $S_{\delta}$ is given by (4.2), the subspace
Lipschitz constant is then given by

$\displaystyle\Psi\bigl{(}\mathcal{M}(S_{\delta})\bigr{)}=|S_{\delta}|^{1/2}.$

As such, we need to characterize the cardinality of $S_{\delta}$. We remark
that

$\displaystyle R_{q}\geq\sum_{j=1}^{d}|u_{j}^{\star}|^{q}\geq\sum_{j\in
S_{\delta}}|u_{j}^{\star}|^{q}\geq\delta^{q}|S_{\delta}|.$

It follows that

$\displaystyle|S_{\delta}|\leq R_{q}\delta^{-q}.$

Under $\ell^{q}$ sparsity, we do not have that $u^{\star}_{S_{\delta}^{C}}=0$.
However, we have the bound

$\displaystyle\|u^{\star}_{S_{\delta}^{C}}\|_{1}=\sum_{j\in
S^{c}_{\delta}}|u_{j}^{\star}|^{q}|u_{j}^{\star}|^{1-q}\leq
R_{q}\delta^{1-q},$

since $|u_{j}^{\star}|\leq\delta$ for all $j\in S_{\delta}^{C}$. Finally,
since $q<1$, we have that

$\displaystyle\|u^{\star}_{S_{\delta}}\|_{1}\leq R_{q}^{1/q}.$

Hence, we apply Theorem 3.11 and the desired result follows. ∎

###### Proof of Corollary 4.4.

Exactly as in the proof of Corollary 4.2, the choice of
$\lambda=4\sqrt{\frac{\log d}{n}}$ satisfies the assumption that
$\lambda\geq\frac{2\|A^{T}\varepsilon\|_{\infty}}{n}$ with probability
$1-c_{1}\exp(-c_{2}n\lambda^{2}).$ Again, exactly as in Corollary 4.2, the RSC
condition holds with curvature $\kappa=c_{1}\lambda_{\min}(\Sigma)$ and
tolerance $\tau^{2}=c_{2}\rho(\Sigma)^{2}\frac{\log d}{n}$ with probability at
least $1-\frac{e^{-n/32}}{1-e^{-n/32}}$. We pick $\delta$ in the definition of
$S_{\delta}$ to be $\delta=\frac{\lambda}{\lambda_{\min}(\Sigma)}$ and
plugging each of the above equations into Theorem 4.3 yields the desired
result. ∎

###  5.2 Hierarchical Models for Group Sparsity

###### Proof of Theorem 4.5.

From Lemma 3.5 we have that approximate decomposability holds with
$c_{1}(\eta)=\frac{\eta}{2}$ and
$c_{2}(\eta)=\frac{k\eta}{\sqrt{2}}(2-\log\eta).$ For the model subspace
$\mathcal{M}(S_{\mathcal{G},\delta})$, where $S_{\mathcal{G},\delta}$ is given
in (4.3), the subspace Lipschitz constant is given by

$\displaystyle\Psi\bigl{(}\mathcal{M}(S_{\mathcal{G},\delta})\bigr{)}=\sup_{u\in\mathcal{M}(S_{\mathcal{G},\delta})\backslash\\{0\\}}\frac{\sum_{j=1}^{k}\|u_{g_{j}}\|_{C_{j}}}{\|u\|_{2}}=\frac{1}{\sqrt{|S_{\mathcal{G},\delta}|}}\sum_{j\in
S_{\mathcal{G},j}}\sqrt{\frac{1}{\lambda_{\min}(C_{j})}}.$ (5.1)

Note that if $C_{j}=I$ for all $j$, this reduces to
$\Psi\bigl{(}\mathcal{M}(S_{\mathcal{G},\delta})\bigr{)}=\sqrt{|S_{\mathcal{G},\delta}|}.$
We can upper bound the subspace Lipschitz constant as

$\displaystyle\Psi\bigl{(}\mathcal{M}(S_{\mathcal{G},\delta})\bigr{)}\leq\sqrt{\frac{1}{\lambda_{\min}(C_{j^{*}})}\mathopen{}\mathclose{{}\left|S_{\mathcal{G},\delta}}\right|},$
(5.2)

where $j^{*}=\operatorname*{arg\,max}_{j\in
S_{\mathcal{G},\delta}}\frac{1}{\lambda_{\min}(C_{j})}.$ If the matrices
$C_{j}$ are all equal, then the upper bound (5.2) is an equality. We will use
this upper bound in what follows, but we remark that if the spectrum of the
matrices $C_{j}$ varies drastically between groups of variables, one may be
able to derive sharper bounds by working with (5.1) directly. The dual norm of
our group sparse norm is given by

$\displaystyle\mathsf{R}^{*}(u):=\sup_{v\neq 0}\frac{\langle
u,v\rangle}{\sum_{j=1}^{k}\|v_{g_{j}}\|_{C_{j}}}=\max_{g_{j}\in\mathcal{G}}\frac{2}{n}\|u_{g_{j}}\|_{C_{j}^{-1}},$

and we see that our assumption on $\lambda$ satisfies the condition of Theorem
3.11. As in the proof of Theorem 4.3, one can show that

$\displaystyle\mathopen{}\mathclose{{}\left|S_{\mathcal{G},\delta}}\right|$
$\displaystyle\leq R_{q}\delta^{-q},$
$\displaystyle\mathsf{R}\bigl{(}u^{\star}_{\mathcal{M}(S_{\mathcal{G},\delta}^{C})}\bigr{)}$
$\displaystyle=\sum_{j\in
S_{\mathcal{G},\delta}^{C}}\|u^{\star}_{g_{j}}\|_{C_{j}}\leq
R_{q}\delta^{1-q},$
$\displaystyle\mathsf{R}\bigl{(}u^{\star}_{\mathcal{M}(S_{\mathcal{G},\delta})}\bigr{)}$
$\displaystyle=\sum_{j\in\mathcal{M}(S_{\mathcal{G},\delta})}\|u^{\star}_{g_{j}}\|_{C_{j}}\leq
R_{q}^{\frac{1}{q}}.$

These bounds along with Theorem 3.11 give the desired result. ∎

###### Proof of Corollary 4.6.

We first want to verify that the choice of
$\lambda=2\mathopen{}\mathclose{{}\left(\sqrt{\frac{p_{\max}}{n}}+\sqrt{\frac{\log
k}{n}}}\right)$ satisfies

$\displaystyle\lambda\geq\max_{g_{j}\in\mathcal{G}}\frac{2}{n}\|(A^{T}\varepsilon)_{g_{j}}\|_{2}$
(5.3)

with high probability. We claim that if $A$ satisfies the block-normalization
condition (4.4) and $\varepsilon\sim\mathcal{N}(0,I)$, then we have

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\frac{1}{n}\max_{1\leq
j\leq k}\|X_{g_{j}}^{T}\varepsilon\|_{2}\geq
2\mathopen{}\mathclose{{}\left(\sqrt{\frac{p_{\max}}{n}}+\sqrt{\frac{\log
k}{n}}}\right)}\right]\leq 2\exp\mathopen{}\mathclose{{}\left(2-\log
k}\right).$ (5.4)

For a fixed group $g_{j}\in\mathcal{G}$ of size $p_{j}$, we consider the
submatrix $A_{g_{j}}\in\mathbb{R}^{n\times p_{j}}$. We characterize the
deviation above the mean using the same argument as [44] to get that, for all
$t>0,$

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\frac{1}{n}\|X^{T}_{g_{j}}\varepsilon\|_{2}\geq\mathbb{E}\mathopen{}\mathclose{{}\left[\frac{1}{n}\|X^{T}_{g_{j}}\varepsilon\|_{2}}\right]+t}\right]\leq
2\exp\mathopen{}\mathclose{{}\left(-\frac{nt^{2}}{2}}\right).$ (5.5)

The expectation in this bound can be bounded as follows

$\displaystyle\mathbb{E}\mathopen{}\mathclose{{}\left[\frac{1}{n}\|X^{T}_{g_{j}}\varepsilon\|_{2}}\right]\overset{\text{(i)}}{\leq}\sqrt{\frac{2}{n}}\mathbb{E}\mathopen{}\mathclose{{}\left[\sqrt{\varepsilon^{T}\varepsilon}}\right]\overset{\text{(ii)}}{\leq}\sqrt{\frac{2p_{\max}}{n}},$
(5.6)

where for (i) we use the same Gaussian comparison principle argument as [44]
and for (ii) we use Jensen’s inequality and that
$\varepsilon\sim\mathcal{N}(0,I).$ Combining (5.5) and (5.6) gives that, for
all $t>0,$

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\frac{1}{n}\|X^{T}_{g_{j}}\varepsilon\|_{2}\geq\sqrt{\frac{2p_{\max}}{n}}+t}\right]\leq
2\exp\mathopen{}\mathclose{{}\left(-\frac{nt^{2}}{2}}\right).$

Applying a union bound over all $j\in\\{1,\ldots,k\\}$ and setting
$t^{2}=\frac{4\log k}{n}$, we deduce that (5.4) holds as claimed. Hence, our
choice of $\lambda$ satisfies the bound (5.3) with probability at least
$1-2/k^{2}$. From [44, Proposition 1], we conclude that the RSC condition
holds with curvature $\kappa=\frac{1}{4}\lambda_{\min}(\Sigma)$ and tolerance
$\tau^{2}=9\max_{g_{j}\in\mathcal{G}}\|(\Sigma^{1/2}_{g_{j}}\|_{op}\mathopen{}\mathclose{{}\left(\sqrt{\frac{p_{\max}}{n}}+\sqrt{\frac{\log
k}{n}}}\right)^{2}$. Applying Theorem 4.5 with $\delta\asymp\lambda$ gives the
desired result. ∎

###  5.3 Hierarchical Models for Sparse Representations

###### Proof of Theorem 4.7.

From Lemma 3.7 we have that approximate decomposability holds with
$c_{1}(\eta)=\frac{\eta}{2}$ and
$c_{2}(\eta)=\frac{k\eta}{\sqrt{2}}(2-\log\eta).$ For the model subspace
$\mathcal{M}(S_{W,\delta})$, where $S_{\mathcal{G},\delta}$ is given in (4.5),
we have that the subspace Lipschitz constant is given by

$\displaystyle\Psi\bigl{(}\mathcal{M}(S_{W,\delta})\bigr{)}=\sup_{u\in\mathcal{M}(S_{W,\delta})\backslash\\{0\\}}\frac{\|W^{T}u\|_{1}}{\|u\|_{2}}=\sqrt{|S_{W,\delta}|}.$

The dual of $\|W^{T}\cdot\|_{1}$ is given by $\|W^{T}\cdot\|_{\infty}$, so our
assumption on $\lambda$ satisfies the condition of Theorem 3.11. By the same
arguments as in the previous settings, we get the bounds

$\displaystyle|S_{W,\delta}|$ $\displaystyle\leq R_{q}\delta^{-q},$
$\displaystyle R(u^{\star}_{\mathcal{M}(S_{W,\delta}^{C})})$
$\displaystyle=\|W^{T}u^{\star}_{\mathcal{M}(S_{W,\delta}^{C})}\|_{1}\leq
R_{q}\delta^{1-q},$ $\displaystyle R(u^{\star}_{\mathcal{M}(S_{W,\delta})})$
$\displaystyle=\|u^{\star}_{\mathcal{M}(S_{W,\delta})}\|_{1}\leq R_{q}^{1/q}.$

These bounds along with Theorem 3.11 give the desired result. ∎

###### Proof of Corollary 4.8.

By the normalization assumption and the fact that $\varepsilon$ is Gaussian,
we have that, for each $j=1,\ldots,k$,

$\displaystyle\mathbb{P}\Bigl{[}\bigl{|}\langle(AW)_{j},\varepsilon\rangle\bigr{|}\geq
t\Bigr{]}\leq 2\exp\mathopen{}\mathclose{{}\left(-\frac{nt^{2}}{2}}\right).$

A union bound then gives that

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left[\frac{\|W^{T}A^{T}\varepsilon\|_{\infty}}{n}\geq
t}\right]\leq 2\exp\mathopen{}\mathclose{{}\left(-\frac{nt^{2}}{2}}\right),$

and taking $t^{2}=\frac{4\log k}{n}$ we see that the choice of
$\lambda=4\sqrt{\frac{\log k}{n}}$ satisfies the assumption of Theorem 4.7.
The proof of [47, Theorem 1] can be adapted to show that, with probability at
least $1-\frac{e^{-n/32}}{1-e^{-n/32}},$

$\displaystyle\frac{\|A\Delta\|_{2}}{\sqrt{n}}\geq\frac{1}{4}\lambda_{\min}(\Sigma)-9\rho(W^{T}\Sigma
W)^{2}\frac{\log
k}{n}\|W^{T}\Delta\|_{1}^{2},\qquad\forall\Delta\in\mathbb{R}^{d},$

where $\rho(W^{T}\Sigma W)^{2}:=\max_{1\leq j\leq k}(W^{T}\Sigma W)_{jj}$.
Hence, the RSC condition holds with curvature
$\kappa=\frac{1}{4}\lambda_{\min}(\Sigma)$ and tolerance
$\tau^{2}=9\rho(W^{T}\Sigma W)^{2}\frac{\log k}{n}.$ Picking
$\delta\asymp\lambda$ and applying Theorem 4.7 yields the desired result. ∎

##  6 Conclusions

This paper has put forth a unified perspective on sparsity-promoting
conditionally Gaussian hypermodels as M-estimators, bridging the fields of
inverse problems and high-dimensional statistics. We have introduced the
notion of approximate decomposability, which quantifies how close the
regularization term resulting from a Bayesian hierarchical model is to a
decomposable regularizer. We have then generalized the theory in [43] to the
setting of approximately decomposable regularizers, and proved the first known
reconstruction error bounds for MAP estimators of conditionally Gaussian
hypermodels. In particular, we have obtained bounds for widely-used
hierarchical models to reconstruct sparse and group-sparse parameters, as well
as for a novel model which promotes sparse representations. Our bounds suggest
that the minimax rate can be achieved provided that the hyperpameter $\eta$ is
taken sufficiently small.

This work opens a number of important questions for future research. As noted
above, from a statistical perspective, the hyperparameter $\eta$ should be
taken as small as possible to promote sparsity. However, from a computational
perspective, the positivity of $\eta$ enforces the strict convexity of the
objective function, and taking $\eta$ to be too small may degrade the
convergence of the IAS algorithm used to compute the MAP estimators. As such,
an important avenue for research is to investigate the statistical-
computational tradeoff in the choice of $\eta$, as well as computational
comparisons between the Bayesian hypermodels and frequentist M-estimators that
also admit efficient and well-studied optimization algorithms such as ISTA
[28], FISTA [5], and ADMM [10].

The connection between conditionally Gaussian hypermodels and M-estimators can
be leveraged to motivate the design of new hierarchical Bayesian models and
algorithms inspired by M-estimator theory. We conjecture that hierarchical
models can be developed for tasks such as low-rank matrix estimation, sparse
generalized linear models, and sparsity-promoting dictionary learning, and
that these models will inherit the same convenient computational and
statistical properties as models considered in this paper. Likewise, our
theory could be extended to the broader family of generalized gamma
hyperpriors, for which greedy algorithms for efficient MAP estimation have
been developed in [17, 16].

Finally, the analysis in this work relies heavily on the forward map
satisfying the RSC condition (3.7). This condition holds for Gaussian forward
maps, which may not be a realistic assumption for many inverse problems of
interest. As such, verifying the RSC condition (possibly after
preconditioning) for important forward maps in inverse problems, such as
deconvolution or tomography [36], is an important direction for further
research.

## Acknowledgments

DSA is grateful for the support of the NSF CAREER award DMS-2237628, the DOE
grant DE-SC0022232, and the BBVA Foundation. The authors are thankful to
Jiajun Bao and Jiaheng Chen for helpful feedback on a previous version of this
manuscript, and to Omar Al-Ghattas for helpful discussions.

## References

  * [1] S. Agapiou, O. Papaspiliopoulos, D. Sanz-Alonso, and A. M. Stuart, Importance sampling: Intrinsic dimension and computational cost, Statistical Science, (2017), pp. 405–431.
  * [2] S. Agrawal, H. Kim, D. Sanz-Alonso, and A. Strang, A variational inference approach to inverse problems with gamma hyperpriors, SIAM/ASA Journal on Uncertainty Quantification, 10 (2022), pp. 1533–1559.
  * [3] F. R. Bach, Consistency of the group Lasso and multiple kernel learning, Journal of Machine Learning Research, 9 (2008), pp. 1179–1225.
  * [4] J. M. Bardsley, D. Calvetti, and E. Somersalo, Hierarchical regularization for edge-preserving reconstruction of PET images, Inverse Problems, 26 (2010), p. 035010.
  * [5] A. Beck and M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM journal on imaging sciences, 2 (2009), pp. 183–202.
  * [6] P. J. Bickel, J. B. Brown, H. Huang, and Q. Li, An overview of recent developments in genomics and associated statistical methods, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 367 (2009), pp. 4313–4337.
  * [7] P. J. Bickel, Y. Ritov, and A. B. Tsybakov, Simultaneous analysis of Lasso and Dantzig selector, The Annals of Statistics, 37 (2009), pp. 1705–1732.
  * [8] A. Bocchinfuso, D. Calvetti, and E. Somersalo, Adaptive anisotropic Bayesian meshing for inverse problems, arXiv preprint arXiv:2310.03855, (2023).
  * [9] A. Bocchinfuso, D. Calvetti, and E. Somersalo, Bayesian sparsity and class sparsity priors for dictionary learning and coding, arXiv preprint arXiv:2309.00999, (2023).
  * [10] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, et al., Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends® in Machine learning, 3 (2011), pp. 1–122.
  * [11] P. Bühlmann and S. A. van De Geer, Statistics for High-Dimensional Data: Methods, Theory and Applications, Springer Science & Business Media, 2011.
  * [12] T. Bui-Thanh, O. Ghattas, J. Martin, and G. Stadler, A computational framework for infinite-dimensional Bayesian inverse problems Part I: The linearized case, with application to global seismic inversion, SIAM Journal on Scientific Computing, 35 (2013), pp. A2494–A2523.
  * [13] D. Calvetti, H. Hakula, S. Pursiainen, and E. Somersalo, Conditionally Gaussian hypermodels for cerebral source localization, SIAM Journal on Imaging Sciences, 2 (2009), pp. 879–909.
  * [14] D. Calvetti, A. Pascarella, F. Pitolli, E. Somersalo, and B. Vantaggi, A hierarchical Krylov–Bayes iterative inverse solver for MEG with physiological preconditioning, Inverse Problems, 31 (2015), p. 125005.
  * [15] D. Calvetti, F. Pitolli, E. Somersalo, and B. Vantaggi, Bayes meets Krylov: Statistically inspired preconditioners for CGLS, SIAM Review, 60 (2018), pp. 429–461.
  * [16] D. Calvetti, M. Pragliola, and E. Somersalo, Sparsity promoting hybrid solvers for hierarchical Bayesian inverse problems, SIAM Journal on Scientific Computing, 42 (2020), pp. A3761–A3784.
  * [17] D. Calvetti, M. Pragliola, E. Somersalo, and A. Strang, Sparse reconstructions from few noisy data: analysis of hierarchical Bayesian models with generalized gamma hyperpriors, Inverse Problems, 36 (2020), p. 025010.
  * [18] D. Calvetti and E. Somersalo, Bayesian Scientific Computing, vol. 215, Springer Nature, 2023.
  * [19] D. Calvetti and E. Somersalo, Computationally efficient sampling methods for sparsity promoting hierarchical Bayesian models, arXiv preprint arXiv:2303.16988, (2023).
  * [20] D. Calvetti, E. Somersalo, and A. Strang, Hierachical Bayesian models and sparsity: $L^{2}$-magic, Inverse Problems, 35 (2019), p. 035003.
  * [21] E. Candes and B. Recht, Exact matrix completion via convex optimization, Communications of the ACM, 55 (2012), pp. 111–119.
  * [22] E. Candes and B. Recht, Simple bounds for recovering low-complexity models, Mathematical Programming, 141 (2013), pp. 577–589.
  * [23] E. J. Candès and D. L. Donoho, New tight frames of curvelets and optimal representations of objects with piecewise $C^{2}$ singularities, Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 57 (2004), pp. 219–266.
  * [24] E. J. Candes, Y. C. Eldar, D. Needell, and P. Randall, Compressed sensing with coherent and redundant dictionaries, Applied and Computational Harmonic Analysis, 31 (2011), pp. 59–73.
  * [25] V. Chandrasekaran, B. Recht, P. A. Parrilo, and A. S. Willsky, The convex geometry of linear inverse problems, Foundations of Computational Mathematics, 12 (2012), pp. 805–849.
  * [26] S. Cotter, G. Roberts, A. Stuart, and D. White, MCMC methods for functions: modifying old algorithms to make them faster, Statistical Science, 28 (2013), pp. 424–446.
  * [27] L. David, Hyperspectral image data analysis as a high dimensional signal processing problem, IEEE Signal Processing Magazine, 19 (2002), pp. 17–28.
  * [28] D. L. Donoho, De-noising by soft-thresholding, IEEE Transactions on Information Theory, 41 (1995), pp. 613–627.
  * [29] D. L. Donoho et al., High-dimensional data analysis: The curses and blessings of dimensionality, AMS Math Challenges Lecture, 1 (2000), p. 32.
  * [30] S. A. Geer, Empirical Processes in M-Estimation, vol. 6, Cambridge University Press, 2000.
  * [31] A. E. Gelfand, P. Diggle, P. Guttorp, and M. Fuentes, Handbook of Spatial Statistics, CRC Press, 2010.
  * [32] A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin, Bayesian Data Analysis, Chapman and Hall/CRC, 1995.
  * [33] J. Glaubitz and A. Gelb, Leveraging joint sparsity in hierarchical Bayesian learning, arXiv preprint arXiv:2303.16954, (2023).
  * [34] J. Huang and T. Zhang, The benefit of group sparsity, The Annals of Statistics, 38 (2010), pp. 1978–2004.
  * [35] L. Jacob, G. Obozinski, and J.-P. Vert, Group Lasso with overlap and graph Lasso, in Proceedings of the 26th Annual International Conference on Machine Learning, 2009, pp. 433–440.
  * [36] J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, vol. 160, Springer Science & Business Media, 2006.
  * [37] H. Kim, D. Sanz-Alonso, and A. Strang, Hierarchical ensemble Kalman methods with sparsity-promoting generalized gamma hyperpriors, Foundations of Data Science, 5 (2023), pp. 366–388.
  * [38] K. Law, A. Stuart, and K. Zygalakis, Data Assimilation: A Mathematical Introduction, vol. 62, Springer, 2015.
  * [39] K. Lounici, M. Pontil, A. Tsybakov, and S. Van De Geer, Taking advantage of sparsity in multi-task learning, in COLT 2009-The 22nd Conference on Learning Theory, 2009.
  * [40] M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly, Compressed sensing MRI, IEEE Signal Processing Magazine, 25 (2008), pp. 72–82.
  * [41] S. Mallat, A Wavelet Tour of Signal Processing, Elsevier, 1999.
  * [42] N. Meinshausen and B. Yu, Lasso-type recovery of sparse representations for high-dimensional data, The Annals of Statistics, 37 (2009), pp. 246–270.
  * [43] S. N. Negahban, P. Ravikumar, M. J. Wainwright, and B. Yu, A unified framework for high-dimensional analysis of M-estimators with decomposable regularizers, Statistical Science, 27 (2012), pp. 538–557.
  * [44] S. N. Negahban, P. Ravikumar, M. J. Wainwright, and B. Yu, Supplement to ‘A unified framework for high-dimensional analysis of M-estimators with decomposable regularizers’, Statistical Science, (2012).
  * [45] T. Park and G. Casella, The Bayesian Lasso, Journal of the American Statistical Association, 103 (2008), pp. 681–686.
  * [46] M. Pragliola, D. Calvetti, and E. Somersalo, Overcomplete representation in a hierarchical Bayesian framework, Inverse Problems and Imaging, 16 (2022), pp. 19–38.
  * [47] G. Raskutti, M. J. Wainwright, and B. Yu, Restricted eigenvalue properties for correlated Gaussian designs, The Journal of Machine Learning Research, 11 (2010), pp. 2241–2259.
  * [48] , Minimax rates of estimation for high-dimensional linear regression over $\ell_{q}$-balls, IEEE Transactions on Information Theory, 57 (2011), pp. 6976–6994.
  * [49] B. Recht, M. Fazel, and P. A. Parrilo, Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization, SIAM Review, 52 (2010), pp. 471–501.
  * [50] A. Ron and Z. Shen, Affine systems in $L^{2}(\mathbb{R}^{d})$: the analysis of the analysis operator, Journal of Functional Analysis, 148 (1997), pp. 408–447.
  * [51] D. Sanz-Alonso, A. Stuart, and A. Taeb, Inverse Problems and Data Assimilation, vol. 107, Cambridge University Press, 2023.
  * [52] D. Sanz-Alonso and N. Waniorek, Analysis of a computational framework for Bayesian inverse problems: ensemble Kalman updates and MAP estimators under mesh refinement, arXiv preprint arXiv:2304.09933, (2023).
  * [53] Z. Si, Y. Liu, and A. Strang, Path-following methods for Maximum a Posteriori estimators in Bayesian hierarchical models: How estimates depend on hyperparameters, arXiv preprint arXiv:2211.07113, (2022).
  * [54] A. M. Stuart, Inverse problems: a Bayesian perspective, Acta Numerica, 19 (2010), pp. 451–559.
  * [55] R. Tibshirani, Regression shrinkage and selection via the Lasso, Journal of the Royal Statistical Society Series B: Statistical Methodology, 58 (1996), pp. 267–288.
  * [56] N. G. Trillos, Z. Kaplan, T. Samakhoana, and D. Sanz-Alonso, On the consistency of graph-based Bayesian semi-supervised learning and the scalability of sampling algorithms, The Journal of Machine Learning Research, 21 (2020), pp. 1022–1068.
  * [57] S. van de Geer, Weakly decomposable regularization penalties and structured sparsity, Scandinavian Journal of Statistics, 41 (2014), pp. 72–86.
  * [58] S. A. van de Geer, High-dimensional generalized linear models and the Lasso, The Annals of Statistics, 36 (2008), pp. 614–645.
  * [59] M. J. Wainwright, High-Dimensional Statistics: A Non-Asymptotic Viewpoint, vol. 48, Cambridge University Press, 2019.
  * [60] S. F. Waldron, An Introduction to Finite Tight Frames, Springer, 2018\.
  * [61] N. Waniorek, D. Calvetti, and E. Somersalo, Bayesian hierarchical dictionary learning, Inverse Problems, 39 (2023), p. 024006.
  * [62] M. Yuan and Y. Lin, Model selection and estimation in regression with grouped variables, Journal of the Royal Statistical Society Series B: Statistical Methodology, 68 (2006), pp. 49–67.

##  Appendix A Computational Properties of GS-IAS and O-IAS

###### Proof of Theorem 2.6.

Lemma 2.2 in [14] shows that the objective function (2.11) is strictly convex
over $\mathbb{R}^{d}\times\mathbb{R}^{k}_{+}$ in the special case where each
group is of size $p_{j}=3$. The proof holds as written for non-overlapping
groups of arbitrary sizes. We denote the Hessian of the objective function
(2.11) as

$\displaystyle
H(u,\theta)=\begin{bmatrix}\nabla^{2}_{uu}{\mathsf{J}}(u,\theta)&\nabla^{2}_{u\theta}{\mathsf{J}}(u,\theta)\\\
\nabla^{2}_{\theta
u}{\mathsf{J}}(u,\theta)&\nabla^{2}_{\theta\theta}{\mathsf{J}}(u,\theta)\\\
\end{bmatrix}.$

Let $e^{\ell}:=u^{\ell}-\hat{u}.$ To show linear convergence, we apply Lemma
4.1 in [20], which gives

$\displaystyle e^{\ell+1}$
$\displaystyle=\nabla^{2}_{uu}{\mathsf{J}}(\hat{u},\hat{\theta})^{-1}\nabla^{2}_{u\theta}{\mathsf{J}}(\hat{u},\hat{\theta})\nabla^{2}_{\theta\theta}{\mathsf{J}}(\hat{u},\hat{\theta})^{-1}\nabla^{2}_{\theta
u}{\mathsf{J}}(\hat{u},\hat{\theta})e^{\ell}+o(\|e^{\ell}\|^{2}),$

where each of the blocks of the Hessian are evaluated at the MAP estimator
$(\hat{u},\hat{\theta}).$ Denoting by $\text{diag}(B_{j})$ the block-diagonal
matrix with blocks $B_{j},1\leq j\leq k,$ we thus obtain that

$\displaystyle e^{\ell+1}$
$\displaystyle=\mathopen{}\mathclose{{}\left(\frac{A^{T}A}{n}+\lambda
D_{\hat{\theta}}^{-1}}\right)^{-1}\text{diag}(\lambda\hat{\theta}_{j}C_{j}^{-1}\hat{u}_{g_{j}})\times$
$\displaystyle\hskip
56.9055pt\times\text{diag}(\lambda\hat{\theta}_{j}^{-3}\|\hat{u}_{g_{j}}\|_{C_{j}}^{2}+\hat{\theta}_{j}^{-2}\eta)^{-1}\text{diag}(\lambda\hat{\theta}_{j}C_{j}^{-1}\hat{u}_{g_{j}})e^{\ell}+o(\|e^{\ell}\|^{2})$
$\displaystyle=D_{\hat{\theta}}^{1/2}\mathopen{}\mathclose{{}\left(\frac{D_{\hat{\theta}}^{1/2}A^{T}AD_{\hat{\theta}}^{1/2}}{n\lambda}+I}\right)^{-1}\text{diag}(\hat{\theta}_{j}^{-3/2}C_{j}^{-1/2}\hat{u}_{g_{j}})\times$
$\displaystyle\hskip
56.9055pt\times\text{diag}(\hat{\theta}_{j}^{-3}\|\hat{u}_{g_{j}}\|_{C_{j}}^{2}+\hat{\theta}_{j}^{-2}\eta)^{-1}\text{diag}(\hat{\theta}_{j}C_{j}^{-1}\hat{u}_{g_{j}})e^{\ell}+o(\|e^{\ell}\|^{2}).$

Multiplying both sides by $D_{\theta}^{-1/2},$ taking norms, and neglecting
higher order terms, we get that

$\displaystyle\|e^{\ell+1}\|_{D_{\hat{\theta}}}\leq\mathopen{}\mathclose{{}\left\|\text{diag}\biggl{(}\frac{C_{j}^{1/2}\hat{u}_{g_{j}}\hat{u}_{g_{j}}^{T}C_{j}^{1/2}}{\|\hat{u}_{g_{j}}\|_{C_{j}}^{2}+\hat{\theta}_{j}\eta}\biggr{)}}\right\|_{op}\|e^{\ell}\|_{D_{\hat{\theta}}}.$

Since for every $1\leq j\leq k$ it holds that
$\|C_{j}^{1/2}\hat{u}_{g_{j}}\hat{u}_{g_{j}}^{T}C_{j}^{1/2}\|_{op}\leq\|\hat{u}_{g_{j}}\|_{C_{j}}^{2}$
and that $\hat{\theta}_{j}\eta>0$, we have linear convergence in the
Mahalanobis norm.

To show that $\hat{u}$ is the global minimizer of (2.13), we refer to Theorem
2.1 in [14], which proves the result for the special case of each group having
size $p_{j}=3$. Their proof still holds in our more general setting. ∎

###### Proof of Theorem 2.8.

We first show that the O-IAS objective function $\mathsf{J}(z)$, where
$z=(u,\theta)\in\mathbb{R}^{d}\times\mathbb{R}^{k}_{+},$ given in (2.15) is
strictly convex over $\mathbb{R}^{d}\times\mathbb{R}^{k}_{+}$. We let
$\tilde{\mathsf{J}}(\tilde{z})$ denote the IAS objective given in (2.7) over
$\tilde{z}=(\tilde{u},\theta)\in\mathbb{R}^{k}\times\mathbb{R}^{k}_{+}$, with
the forward map $\tilde{A}=AW$. From [20], we know that
$\tilde{\mathsf{J}}(\tilde{z})$ is strictly convex over
$\mathbb{R}^{k}\times\mathbb{R}^{k}_{+}$. We then observe that
$\mathsf{J}(z)=\tilde{\mathsf{J}}(C\tilde{z}),$ where

$\displaystyle C=\begin{bmatrix}W^{T}&0\\\ 0&I\\\
\end{bmatrix}\in\mathbb{R}^{2k\times d+k}.$

Since this matrix has trivial nullspace and strict convexity is preserved
under such a transformation, we conclude that $\mathsf{J}(z)$ is strictly
convex. Since $\mathsf{J}(z)$ tends to infinity as $\theta_{j}\to 0$,
$\theta_{j}\to\infty$, and $\|u\|\to\infty$, a minimizer
$(\hat{u},\hat{\theta})$ exists and is unique. As in the proof of Theorem 2.6,
we let $e^{\ell}=u^{\ell}-\hat{u}$ and apply Lemma 4.1 in [20], which gives

$\displaystyle e^{\ell+1}$
$\displaystyle=\nabla^{2}_{uu}{\mathsf{J}}(\hat{u},\hat{\theta})^{-1}\nabla^{2}_{u\theta}{\mathsf{J}}(\hat{u},\hat{\theta})\nabla^{2}_{\theta\theta}{\mathsf{J}}(\hat{u},\hat{\theta})^{-1}\nabla^{2}_{\theta
u}{\mathsf{J}}(\hat{u},\hat{\theta})e^{\ell}+o(\|e^{\ell}\|^{2}),$

where each of the blocks of the Hessian are evaluated at the MAP estimator
$(\hat{u},\hat{\theta}).$ Plugging in the expressions of the Hessian blocks,
we get that

$\displaystyle e^{\ell+1}$
$\displaystyle=\mathopen{}\mathclose{{}\left(\frac{A^{T}A}{n}+\lambda
WD_{\hat{\theta}}^{-1}W^{T}}\right)^{-1}W\times$ $\displaystyle\hskip
28.45274pt\times\text{diag}(\lambda
W^{T}\hat{u}/\hat{\theta}^{2})\mathopen{}\mathclose{{}\left(\text{diag}(\lambda(W^{T}\hat{u})^{2}/\hat{\theta}^{3}+\eta/\hat{\theta}^{2})}\right)^{-1}\text{diag}(\lambda
W^{T}\hat{u}/\hat{\theta}^{2})W^{T}e^{\ell}+o(\|e^{\ell}\|^{2}).$

Multiplying both sides by $\text{diag}(1/\hat{\theta}^{1/2})W^{T},$ we get
that

$\displaystyle\text{diag}(\hat{\theta}^{-1/2})W^{T}e^{\ell+1}$
$\displaystyle=\text{diag}(\hat{\theta}^{-1/2})W^{T}\mathopen{}\mathclose{{}\left(\frac{A^{T}A}{n}+\lambda
WD_{\hat{\theta}}^{-1}W^{T}}\right)^{-1}W\times$ $\displaystyle\hskip
56.9055pt\times\text{diag}\mathopen{}\mathclose{{}\left(\frac{\lambda(W^{T}\hat{u})^{2}}{(W^{T}\hat{u})^{2}\hat{\theta}^{1/2}+\eta\hat{\theta}^{3/2}}}\right)\text{diag}(\hat{\theta}^{-1/2})W^{T}e^{\ell}+o(\|e^{\ell}\|^{2}).$

Taking norms of both sides, we see that if

$\displaystyle\mu=\mathopen{}\mathclose{{}\left\|\text{diag}(\hat{\theta}^{-1/2})W^{T}\mathopen{}\mathclose{{}\left(\frac{A^{T}A}{n}+\lambda
WD_{\hat{\theta}}^{-1}W^{T}}\right)^{-1}W\text{diag}\mathopen{}\mathclose{{}\left(\frac{\lambda(W^{T}\hat{u})^{2}}{(W^{T}\hat{u})^{2}\hat{\theta}^{1/2}+\eta\hat{\theta}^{3/2}}}\right)}\right\|<1,$

then we have linear convergence in the Mahalanobis norm. Letting
$W\text{diag}(1/\hat{\theta})W^{T}=ZZ$, where $Z\in\mathbb{R}^{d\times d}$ is
the unique matrix square root of $W\text{diag}(1/\hat{\theta})W^{T}$, we have
that

$\displaystyle\mu$
$\displaystyle=\mathopen{}\mathclose{{}\left\|\text{diag}(1/\hat{\theta}^{1/2})W^{T}Z^{-1}\mathopen{}\mathclose{{}\left(\frac{Z^{-1}A^{T}AZ^{-1}}{\lambda
n}+I}\right)^{-1}Z^{-1}W\text{diag}\mathopen{}\mathclose{{}\left(\frac{(W^{T}\hat{u})^{2}}{(W^{T}\hat{u})^{2}\hat{\theta}^{1/2}+\eta\hat{\theta}^{3/2}}}\right)}\right\|$
$\displaystyle\leq\mathopen{}\mathclose{{}\left\|\text{diag}(1/\hat{\theta}^{1/2})W^{T}Z^{-1}Z^{-1}W\text{diag}\mathopen{}\mathclose{{}\left(\frac{(W^{T}\hat{u})^{2}}{(W^{T}\hat{u})^{2}\hat{\theta}^{1/2}+\eta\hat{\theta}^{3/2}}}\right)}\right\|$
$\displaystyle=\mathopen{}\mathclose{{}\left\|\text{diag}(1/\hat{\theta}^{1/2})W^{T}\mathopen{}\mathclose{{}\left(W\text{diag}(1/\hat{\theta})W^{T}}\right)^{-1}W\text{diag}\mathopen{}\mathclose{{}\left(\frac{(W^{T}\hat{u})^{2}}{(W^{T}\hat{u})^{2}\hat{\theta}^{1/2}+\eta\hat{\theta}^{3/2}}}\right)}\right\|.$

Since $WW^{T}=I$, the SVD of $W$ must have the form

$\displaystyle W=U\begin{bmatrix}I&0\end{bmatrix}V^{T}=U\tilde{V}^{T},$

where $\tilde{V}\in\mathbb{R}^{k\times d}$ is the matrix containing the first
$d$ columns of $V$. Plugging this expression into our upper bound for $\mu$,
we get that

$\displaystyle\mu$
$\displaystyle\leq\mathopen{}\mathclose{{}\left\|\text{diag}(1/\hat{\theta}^{1/2})\tilde{V}U^{T}\mathopen{}\mathclose{{}\left(U\tilde{V}^{T}\text{diag}(1/\hat{\theta})\tilde{V}U^{T}}\right)^{-1}U\tilde{V}^{T}\text{diag}\mathopen{}\mathclose{{}\left(\frac{(W^{T}\hat{u})^{2}}{(W^{T}\hat{u})^{2}\hat{\theta}^{1/2}+\eta\hat{\theta}^{3/2}}}\right)}\right\|$
$\displaystyle=\mathopen{}\mathclose{{}\left\|\text{diag}(1/\hat{\theta}^{1/2})\tilde{V}\mathopen{}\mathclose{{}\left(\tilde{V}^{T}\text{diag}(1/\hat{\theta})\tilde{V}}\right)^{-1}\tilde{V}^{T}\text{diag}\mathopen{}\mathclose{{}\left(\frac{(W^{T}\hat{u})^{2}}{(W^{T}\hat{u})^{2}\hat{\theta}^{1/2}+\eta\hat{\theta}^{3/2}}}\right)}\right\|.$

Making the substitution $V^{\prime}=\text{diag}(1/\hat{\theta}^{1/2})$, we get
that

$\displaystyle\mu\leq\mathopen{}\mathclose{{}\left\|V^{\prime}\mathopen{}\mathclose{{}\left(V^{\prime
T}V^{\prime}}\right)^{-1}V^{\prime
T}\text{diag}\mathopen{}\mathclose{{}\left(\frac{(W^{T}\hat{u})^{2}}{(W^{T}\hat{u})^{2}+\eta\hat{\theta}}}\right)}\right\|.$

Observing that $V^{\prime}\mathopen{}\mathclose{{}\left(V^{\prime
T}V^{\prime}}\right)^{-1}V^{\prime T}$ is an orthogonal projector onto the
range of $V^{\prime}$, we have that

$\displaystyle\mu\leq\mathopen{}\mathclose{{}\left\|\text{diag}\mathopen{}\mathclose{{}\left(\frac{(W^{T}\hat{u})^{2}}{(W^{T}\hat{u})^{2}+\eta\hat{\theta}}}\right)}\right\|<1.$

Consequently, we have shown linear convergence of the O-IAS algorithm in the
Mahalanobis norm.

We will now show that the minimizer of (2.15) over $u$ and $\theta$ can be
written as the minimizer of $\mathsf{J}\bigl{(}u,f(u)\bigr{)},$ where

$\displaystyle
f_{j}(u)=\mathopen{}\mathclose{{}\left(\frac{\eta}{2}+\sqrt{\frac{\eta^{2}}{4}+\frac{(W^{T}u)_{j}}{2}}}\right).$

Since the minimizer $(\hat{u},\hat{\theta})$ is a critical point, it follows
that

$\displaystyle\nabla_{u}\mathsf{J}(\hat{u},\hat{\theta})$
$\displaystyle=\mathopen{}\mathclose{{}\left(\frac{A^{T}A}{n}+\lambda
W\text{diag}(1/\hat{\theta})W^{T}}\right)\hat{u}-\frac{A^{T}y}{n}=0,$
$\displaystyle\nabla_{\theta}\mathsf{J}(\hat{u},\hat{\theta})$
$\displaystyle=\sum_{j=1}^{k}\mathopen{}\mathclose{{}\left(\frac{(W^{T}\hat{u})_{j}^{2}}{2\hat{\theta}_{j}^{2}}-\frac{\eta}{\hat{\theta}_{j}}-1}\right)e_{j}=0.$
(A.1)

Given $\hat{u}$, we solve (A.1) for $\hat{\theta}$ in terms of $\hat{u}$. Some
algebra shows that the solution is given by $\hat{\theta}=f(\hat{u})$, and
thus $\hat{u}$ is a critical point of

$\displaystyle\mathsf{F}(u)=\mathsf{J}\bigl{(}u,f(u)\bigr{)}.$ (A.2)

If $\tilde{u}$ is a critical point of (A.2), letting
$\tilde{\theta}=f(\tilde{u})$, we get that $(\tilde{u},\tilde{\theta})$ is a
critical point of $\mathsf{J}(u,\theta)$ and by uniqueness of the minima we
have that $\hat{u}=\tilde{u}$, completing the proof. ∎

##  Appendix B Approximate Decomposability: Examples

###### Proof of Lemma 3.3.

By the triangle inequality, it holds that

$\displaystyle\frac{|u_{j}|}{\sqrt{2}}\leq
f_{j}(u)\leq\eta+\frac{|u_{j}|}{\sqrt{2}},\quad\text{and}\quad\eta\leq
f_{j}(u).$ (B.1)

Using these inequalities, we deduce that

$\displaystyle\mathsf{R}_{\eta}(u)$
$\displaystyle=\frac{1}{2\sqrt{2}}\sum_{j=1}^{d}\mathopen{}\mathclose{{}\left[\frac{u_{j}^{2}}{f_{j}(u)}+2f_{j}(u)-2\eta\log
f_{j}(u)}\right]$
$\displaystyle\leq\frac{1}{2\sqrt{2}}\sum_{j=1}^{d}\mathopen{}\mathclose{{}\left[\sqrt{2}|u_{j}|+2\eta+\sqrt{2}|u_{j}|-2\eta\log\mathopen{}\mathclose{{}\left(\eta}\right)}\right]\leq\|u\|_{1}+\frac{d\eta}{\sqrt{2}}(1-\log\eta),$

proving the upper bound in (3.4). For the lower bound, we have

$\displaystyle\mathsf{R}_{\eta}(u)$
$\displaystyle\overset{\text{(i)}}{\geq}\frac{1}{2\sqrt{2}}\sum_{j=1}^{d}\mathopen{}\mathclose{{}\left[\frac{u_{j}^{2}}{\eta+\frac{|u_{j}|}{\sqrt{2}}}+\sqrt{2}|u_{j}|-2\eta\log\mathopen{}\mathclose{{}\left(\eta+\frac{|u_{j}|}{\sqrt{2}}}\right)}\right]$
$\displaystyle\overset{\text{(ii)}}{\geq}\frac{1}{2\sqrt{2}}\sum_{j=1}^{d}\mathopen{}\mathclose{{}\left[\frac{u_{j}^{2}}{\eta+\frac{|u_{j}|}{\sqrt{2}}}+\sqrt{2}|u_{j}|-2\eta\frac{|u_{j}|}{\sqrt{2}}}\right]$
$\displaystyle\overset{\text{(iii)}}{\geq}\frac{1}{2\sqrt{2}}\sum_{j=1}^{d}\mathopen{}\mathclose{{}\left[\sqrt{2}|u_{j}|-2\eta+\sqrt{2}|u_{j}|-\sqrt{2}\eta|u_{j}|}\right]=\Bigl{(}1-\frac{\eta}{2}\Bigr{)}\|u\|_{1}-\frac{d\eta}{\sqrt{2}},$

where (i) follows by the inequalities in (B.1), (ii) follows using that for
$0<\eta<1$ it holds that
$\log\Bigl{(}\eta+\frac{|u_{j}|}{\sqrt{2}}\Bigr{)}\leq\frac{|u_{j}|}{\sqrt{2}},$
and (iii) follows using that
$\frac{u_{j}^{2}}{\eta+\frac{|u_{j}|}{\sqrt{2}}}\geq\sqrt{2}|u_{j}|-2\eta.$ ∎

###### Proof of Lemma 3.5.

Similar to the proof of Lemma 3.3, replacing $|u_{j}|$ with
$\|u_{g_{j}}\|_{C_{j}}.$ ∎

###### Proof of Lemma 3.7.

The statement follows from Lemma 3.3 applied to $W^{T}u\in\mathbb{R}^{k}$. ∎

##  Appendix C Approximate Decomposability and General Bounds

In this appendix we prove the key property of approximately decomposable
regularizers in Lemma 3.8 and the general reconstruction bound in Theorem
3.11. To control the reconstruction error $\hat{\Delta}=\hat{u}-u^{\star}$
between the minimizer of the objective $\mathsf{F}$ in (3.6) and the true
parameter $u^{\star}$ generating the data according to (2.4), we introduce a
function $\mathcal{F}:\mathbb{R}^{d}\to\mathbb{R}$ given by

$\displaystyle\begin{split}\mathcal{F}(\Delta):&=\mathsf{F}(u^{\star}+\Delta)-\mathsf{F}(u^{\star})\\\
&=\frac{1}{2n}\Bigl{[}\|y-A(u^{\star}+\Delta)\|_{2}^{2}-\|y-Au^{\star}\|_{2}^{2}\Bigl{]}+\lambda\Bigl{[}\mathsf{R}(u^{\star}+\Delta)-\mathsf{R}(u^{\star})\Bigr{]}.\end{split}$
(C.1)

By construction, $\mathcal{F}(0)=0$. Moreover, for
$\hat{\Delta}=\hat{u}-u^{\star}$ we have the basic inequality

$\mathcal{F}(\hat{\Delta})=\mathsf{F}(\hat{u})-\mathsf{F}(u^{\star})\leq 0,$
(C.2)

since $\hat{u}$ is the minimizer of $\mathsf{F}.$ We will leverage this
inequality to control the reconstruction error.

###### Proof of Lemma 3.8.

We first prove a lower bound on the difference of the regularizer evaluated at
the minimizer and at the true parameter. Using approximate decomposability of
$\mathsf{R}_{\eta}$, reverse triangle inequality, and decomposability of
$\mathsf{R}$ over $\mathcal{M},$ it follows that

$\displaystyle\begin{split}\mathsf{R}_{\eta}(\hat{u})&\geq\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)}\right)\mathsf{R}(u^{\star}+\hat{\Delta})-c_{2}^{L}(\eta)\\\
&\geq\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)}\right)\mathopen{}\mathclose{{}\left(\mathsf{R}(u^{\star}_{\mathcal{M}}+\hat{\Delta}_{\mathcal{M}^{\perp}})-\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})-\mathsf{R}(\hat{\Delta}_{\mathcal{M}})}\right)-c_{2}^{L}(\eta)\\\
&\geq\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)}\right)\mathopen{}\mathclose{{}\left(\mathsf{R}(u^{\star}_{\mathcal{M}})+\mathsf{R}(\hat{\Delta}_{\mathcal{M}^{\perp}})-\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})-\mathsf{R}(\hat{\Delta}_{\mathcal{M}})}\right)-c_{2}^{L}(\eta).\end{split}$
(C.3)

Similarly, using approximate decomposability of $\mathsf{R}_{\eta}$ and
decomposability of $\mathsf{R},$ we obtain

$\mathsf{R}_{\eta}(u^{\star})\leq\mathopen{}\mathclose{{}\left(1+c_{1}^{U}(\eta)}\right)\Bigl{(}\mathsf{R}(u_{\mathcal{M}}^{\star})+\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})\Bigr{)}+c_{2}^{U}(\eta).$
(C.4)

Combining (C.3) and (C.4) gives

$\displaystyle\begin{split}\mathsf{R}_{\eta}(\hat{u})-\mathsf{R}_{\eta}(u^{\star})&\geq\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)}\right)\mathopen{}\mathclose{{}\left(\mathsf{R}(\hat{\Delta}_{\mathcal{M}^{\perp}})-\mathsf{R}(\hat{\Delta}_{\mathcal{M}})}\right)-c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})-2\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})-c_{2}(\eta),\end{split}$
(C.5)

where $c_{i}(\eta)=c_{i}^{L}(\eta)+c_{i}^{U}(\eta)$ and we used that
$c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})\geq 0.$ We next show a
lower bound on the difference of the likelihood evaluated at the minimizer and
at the true parameter. By direct calculation, it holds that

$\displaystyle\frac{1}{2n}\Bigl{[}\|y-A\hat{u}\|_{2}^{2}-\|y-Au^{\star}\|_{2}^{2}\Bigr{]}=\frac{1}{2n}\|A\hat{\Delta}\|_{2}^{2}-\frac{1}{n}\langle
A^{T}\varepsilon,\hat{\Delta}\rangle\geq-\frac{1}{n}\mathopen{}\mathclose{{}\left|\langle
A^{T}\varepsilon,\hat{\Delta}\rangle}\right|.$

Then, using Hölder’s inequality with the regularizer and its dual, we get that

$\displaystyle\begin{split}\frac{1}{2n}\Bigl{[}\|y-A\hat{u}\|_{2}^{2}-\|y-Au^{\star}\|_{2}^{2}\Bigr{]}&\geq-\frac{1}{n}\mathsf{R}^{*}(A^{T}\varepsilon)\mathsf{R}(\hat{\Delta})=-\frac{1}{n}\mathsf{R}^{*}(A^{T}\varepsilon)\mathopen{}\mathclose{{}\left(\mathsf{R}(\hat{\Delta}_{\mathcal{M}})+\mathsf{R}(\hat{\Delta}_{\mathcal{M}^{\perp}})}\right)\\\
&\geq-\frac{\lambda}{2}\mathopen{}\mathclose{{}\left(\mathsf{R}(\hat{\Delta}_{\mathcal{M}})+\mathsf{R}(\hat{\Delta}_{\mathcal{M}^{\perp}})}\right),\end{split}$
(C.6)

where the last inequality follows by our assumption on $\lambda.$ Finally,
(C.5) and (C.6) combined with the basic inequality (C.2) give

$\displaystyle 0\geq\mathcal{F}(\hat{\Delta})$
$\displaystyle\geq\lambda\biggl{(}\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)}\right)\mathopen{}\mathclose{{}\left(\mathsf{R}(\hat{\Delta}_{\mathcal{M}^{\perp}})-\mathsf{R}(\hat{\Delta}_{\mathcal{M}})}\right)-c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})-2\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})-c_{2}(\eta)\biggl{)}$
$\displaystyle-\frac{\lambda}{2}\mathopen{}\mathclose{{}\left(\mathsf{R}(\hat{\Delta}_{\mathcal{M}})+\mathsf{R}(\hat{\Delta}_{\mathcal{M}^{\perp}})}\right),$

which after rearranging implies that

$\displaystyle\mathopen{}\mathclose{{}\left(\frac{1}{2}-c_{1}^{L}(\eta)}\right)\mathsf{R}(\hat{\Delta}_{\mathcal{M}^{\perp}})\leq$
$\displaystyle\mathopen{}\mathclose{{}\left(\frac{3}{2}+c_{1}^{L}(\eta)}\right)\mathsf{R}(\hat{\Delta}_{\mathcal{M}})+2\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})+c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})+c_{2}(\eta).$

The desired result follows from the assumption that
$c_{1}^{L}(\eta)<\frac{1}{4}.$ ∎

Letting $\delta>0$ be an error radius that we will specify later, we define
the set

$\displaystyle\mathbb{k}(\delta):=\mathbb{C}_{\eta}(\mathcal{M})\cap\\{\|\Delta\|_{2}=\delta\\}.$

The following result from [44, Lemma 4] demonstrates that it is sufficient to
show a positive lower bound for the function $\mathcal{F}$ given in (C.1) over
the set $\mathbb{K}(\delta).$

######  Lemma C.1.

If $\mathcal{F}(\Delta)>0$ for all vectors $\Delta\in\mathbb{K}(\delta)$, then
$\|\hat{\Delta}\|_{2}\leq\delta.$

###### Proof.

The proof in [44] still holds essentially as written in our slightly modified
setting. In particular, the proof does not require that the regularizers in
$\mathcal{F}(\Delta)$ be decomposable. The regularizers in our definition of
$\mathbb{C}_{\eta}$ are decomposable with respect to $\mathcal{M}$, and thus
the argument this set is “star shaped” (see [44] or [59, Chapter 9] for a
definition and discussion of this property) holds as written. ∎

We are now ready to prove Theorem 3.11.

###### Proof of Theorem 3.11.

Fix a radius $\delta>0$, which we will specify at the conclusion of the proof.
By Lemma C.1, we proceed by deriving a positive lower bound for
$\mathcal{F}(\Delta)$ over $\mathbb{K}(\delta),$ which will yield a bound on
$\|\hat{\Delta}\|_{2}^{2}.$ By definition of $\mathcal{F}$, the RSC condition,
and the inequality (C.5)

$\displaystyle\mathcal{F}(\Delta)$
$\displaystyle=\frac{1}{2n}\Bigl{[}\|y-A(u^{\star}+\Delta)\|_{2}^{2}-\|y-Au^{\star}\|_{2}^{2}\Bigr{]}+\lambda\mathopen{}\mathclose{{}\left[\mathsf{R}(u^{\star}+\Delta)-\mathsf{R}(u^{\star})}\right]$
$\displaystyle\geq\langle
A^{T}\varepsilon,\Delta\rangle+\frac{\kappa}{2}\|\Delta\|_{2}^{2}-\tau^{2}\mathsf{R}^{2}(\Delta)+\lambda\mathopen{}\mathclose{{}\left[\mathsf{R}(u^{\star}+\Delta)-\mathsf{R}(u^{\star})}\right]$
$\displaystyle\geq\langle
A^{T}\varepsilon,\Delta\rangle+\frac{\kappa}{2}\|\Delta\|_{2}^{2}-\tau^{2}\mathsf{R}^{2}(\Delta)$
$\displaystyle+\lambda\Big{[}\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)}\right)\mathopen{}\mathclose{{}\left(\mathsf{R}(\hat{\Delta}_{\mathcal{M}^{\perp}})-\mathsf{R}(\hat{\Delta}_{\mathcal{M}})}\right)-c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})-2\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})-c_{2}(\eta)\Big{]}.$

Applying Hölder’s inequality and using the assumption on $\lambda$, we have
that $|\langle
A^{T}\varepsilon,\Delta\rangle|\leq\frac{\lambda}{2}\mathsf{R}(\Delta),$ which
gives

$\displaystyle\mathcal{F}(\Delta)\geq\frac{\kappa}{2}\|\Delta\|_{2}^{2}-\tau^{2}\mathsf{R}^{2}(\Delta)-\frac{\lambda}{2}\mathsf{R}(\Delta)$
$\displaystyle+\lambda\Big{[}\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)}\right)\mathopen{}\mathclose{{}\left(\mathsf{R}(\hat{\Delta}_{\mathcal{M}^{\perp}})-\mathsf{R}(\hat{\Delta}_{\mathcal{M}})}\right)-c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})-2\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})-c_{2}(\eta)\Big{]}$
$\displaystyle=\frac{\kappa}{2}\|\Delta\|_{2}^{2}-\tau^{2}\mathsf{R}^{2}(\Delta)$
$\displaystyle+\lambda\biggl{[}\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)-\frac{1}{2}}\right)\mathsf{R}(\hat{\Delta}_{\mathcal{M}^{\perp}})-\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)+\frac{1}{2}}\right)\mathsf{R}(\hat{\Delta}_{\mathcal{M}})-c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})-2\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})-c_{2}(\eta)\biggr{]},$

where the equality follows from the decomposability of $\mathsf{R}.$ Notice
that, for any $\Delta\in\mathbb{C}_{\eta}(\mathcal{M}),$

$\displaystyle\mathsf{R}^{2}(\Delta)=\bigl{(}\mathsf{R}(\Delta_{\mathcal{M}})+\mathsf{R}(\Delta_{\mathcal{M}^{\perp}})\bigr{)}^{2}$
$\displaystyle\leq
16\bigl{(}2\Psi(\mathcal{M})\|\Delta\|_{2}+2\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})+c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})+c_{2}(\eta)\bigr{)}^{2}$
$\displaystyle\leq
64\bigl{[}4\Psi^{2}(\mathcal{M})\|\Delta\|_{2}^{2}+4\mathsf{R}^{2}(u^{\star}_{\mathcal{M}^{\perp}})+c_{1}^{2}(\eta)\mathsf{R}^{2}(u^{\star}_{\mathcal{M}})+c_{2}^{2}(\eta)\bigr{]}.$

Plugging this bound into the previous inequality, we get

$\displaystyle\mathcal{F}(\Delta)$
$\displaystyle\geq\frac{\kappa}{2}\|\Delta\|_{2}^{2}-64\tau^{2}\bigl{[}4\Psi^{2}(\mathcal{M})\|\Delta\|_{2}^{2}+4\mathsf{R}^{2}(u^{\star}_{\mathcal{M}^{\perp}})+c_{1}^{2}(\eta)\mathsf{R}^{2}(u^{\star}_{\mathcal{M}})+c_{2}^{2}(\eta)\bigr{]}$
$\displaystyle-\lambda\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)+\frac{1}{2}}\right)\mathsf{R}(\hat{\Delta}_{\mathcal{M}})+c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})+2\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})+c_{2}(\eta)}\right]$
$\displaystyle\geq\frac{\kappa}{2}\|\Delta\|_{2}^{2}-64\tau^{2}\bigl{[}4\Psi^{2}(\mathcal{M})\|\Delta\|_{2}^{2}+4\mathsf{R}^{2}(u^{\star}_{\mathcal{M}^{\perp}})+c_{1}^{2}(\eta)\mathsf{R}^{2}(u^{\star}_{\mathcal{M}})+c_{2}^{2}(\eta)\bigr{]}$
$\displaystyle-\lambda\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)+\frac{1}{2}}\right)\Psi(\mathcal{M})\|\Delta\|_{2}+c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})+2\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})+c_{2}(\eta)}\right].$

By our assumption that $\tau\Psi(\mathcal{M})\leq\frac{\sqrt{\kappa}}{32}$, it
follows that

$\displaystyle\mathcal{F}(\Delta)$
$\displaystyle\geq\frac{\kappa}{4}\|\Delta\|_{2}^{2}-64\tau^{2}\mathopen{}\mathclose{{}\left[4\mathsf{R}^{2}(u^{\star}_{\mathcal{M}^{\perp}})+c_{1}^{2}(\eta)\mathsf{R}^{2}(u^{\star}_{\mathcal{M}})+c_{2}^{2}(\eta)}\right]$
$\displaystyle-\lambda\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(1-c_{1}^{L}(\eta)+\frac{1}{2}}\right)\Psi(\mathcal{M})\|\Delta\|_{2}+c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})+2\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})+c_{2}(\eta)}\right].$

The right-hand side of this inequality is a quadratic polynomial in
$\|\Delta\|_{2}$ that is positive for $\|\Delta\|$ sufficiently large. In
particular, the quadratic formula shows that taking

$\displaystyle\delta$
$\displaystyle=72\frac{\lambda^{2}}{\kappa^{2}}\Psi^{2}(\mathcal{M})$
$\displaystyle+\frac{8}{\kappa}\biggl{(}64\tau^{2}\mathopen{}\mathclose{{}\left[4\mathsf{R}^{2}(u^{\star}_{\mathcal{M}^{\perp}})+c_{1}^{2}(\eta)\mathsf{R}^{2}(u^{\star}_{\mathcal{M}})+c_{2}^{2}(\eta)}\right]+\lambda\mathopen{}\mathclose{{}\left[c_{1}(\eta)\mathsf{R}(u^{\star}_{\mathcal{M}})+2\mathsf{R}(u^{\star}_{\mathcal{M}^{\perp}})+c_{2}(\eta)}\right]\biggr{)}$

guarantees a positive lower bound on $\mathcal{F}(\Delta)$ over
$\mathbb{K}(\delta)$, completing the proof. ∎