11. Theory of Sparse Regularization

Chapter PDF · Complete book

11.1 Existence and Uniqueness¶

11.1.1 Existence¶

We study the penalized problem (10.10), with $\lambda>0$, and its constrained counterpart (10.11):

and the limiting problem as $\lambda\rightarrow 0$,

where $A \in \mathbb{R}^{P \times N}$, and $\mathcal{L}_y \mathrel{:=} \left\{ x \in \mathbb{R}^N \;;\; Ax=y \right\}$.

Recall that we observe noisy measurements

We seek conditions ensuring that $x_0$ solves $\mathcal{P}_0(Ax_0)$ in the noiseless case, and quantitative bounds on its distance to solutions of $\mathcal{P}_\lambda(Ax_0+w)$ when $\lambda$ is chosen according to the noise level.

The objective in ($\mathcal {P}_\lambda (y)$) is continuous and coercive because $|\!| \cdot |\!|_1$ is coercive, so a minimizer exists.

Uniqueness is not automatic: $A$ may have a nontrivial kernel, and $|\!| \cdot |\!|_1$ is not strictly convex.

If $y \in \mathop{\mathrm{Im}}(A)$, the constraint set in ($\mathcal {P}_0(y)$) is nonempty. A constrained minimizer also exists but need not be unique.

Figure 11.1 shows how the recovered coefficients change with $\lambda$, illustrating the effect of regularization on sparsity and reconstruction accuracy.

11.1.2 Polytope Projection for the Constrained Problem¶

The next proposition characterizes noiseless $\ell^1$ recovery geometrically.

For a nonzero recovered vector $x_0$, its measurement $Ax_0$ lies on the boundary of the scaled polytope $|\!| x_0 |\!|_1 A B$. When $P$ is much smaller than $N$, projection can remove faces and prevent $\ell^1$ recovery. Chapter 12 studies random projections that preserve enough of the geometry of the $\ell^1$ ball in $\mathbb{R}^N$ to allow recovery from $\mathbb{R}^P$.

Identifiability is unchanged by positive scaling: if $x_0$ is identifiable, so is $\rho x_0$ for every $\rho>0$. More generally, the relevant face of $B$, and hence condition (11.1), depends only on $\mathop{\mathrm{sign}}(x_0)$.

For this geometric discussion, assume uniqueness and let $\mathcal{A}: y \mapsto x^\star$ map each datum $y$ to the solution of ($\mathcal {P}_0(y)$).

By (11.1), $A$ and $\mathcal{A}$ are inverse maps between the cones $\mathcal{C}_s = \left\{ x \;;\; \mathop{\mathrm{sign}}(x)=s \right\}$ and $A\mathcal{C}_s$ for recoverable sign patterns $s$. Including the lower-dimensional cones and the origin, their images $A\mathcal{C}_s$ partition $\mathbb{R}^P$. Suppose that the columns $(a_j)_j$ of $A$ have unit norm. For $P=3$, intersecting the cones $A\mathcal{C}_s$ with the unit sphere in $\mathbb{R}^3$ gives a spherical Delaunay subdivision, which is a triangulation under a general-position assumption; see Figure 11.7. In higher dimensions, simplices replace triangles. The subdivision satisfies the empty-cap property: the spherical cap bounded by a triangle’s circumcircle contains no signed column $\pm a_j$ in its interior.

Figure 11.3 illustrates these conclusions in $\mathbb{R}^2$ and $\mathbb{R}^3$.

11.1.3 Optimality Conditions¶

For an index set $I \subset \{1,\ldots,N\}$, write $A=(a_i)_{i=1}^N$ for the columns of $A$ and $A_I \mathrel{:=}(a_i)_{i \in I} \in \mathbb{R}^{P \times |I|}$ for the corresponding submatrix. For a vector $x \in \mathbb{R}^N$, write $x_I \mathrel{:=}(x_i)_{i \in I} \in \mathbb{R}^{|I|}$ for its restriction.

The following proposition rephrases the first-order optimality conditions in a convenient form.

In particular, $\mathop{\mathrm{supp}}(x_\lambda) \subset \mathop{\mathrm{sat}}(\eta_\lambda)$, where $\mathop{\mathrm{sat}}(\eta)=\{i:|\eta_i|=1\}$ is the saturation set.

For the constrained case $\lambda=0$, the next proposition introduces dual certificates arising from the Lagrange multipliers associated with $\mathcal{L}_y$.

Every valid certificate satisfies $\mathop{\mathrm{supp}}(x^\star) \subset \mathop{\mathrm{sat}}(\eta)$, for $\eta\in \mathcal{D}_0(y,x^\star)$.

Writing $I=\mathop{\mathrm{supp}}(x^\star)$, one thus has

Although the notation $\mathcal{D}_0(y,x^\star)$ involves a particular solution $x^\star$, the set of certificates is the same for every primal solution. The duality argument below explains this independence.

11.1.4 Uniqueness¶

The next proposition shows that at least one minimizer uses linearly independent columns of $A$. It does not assert that every minimizer has this property.

If the columns of $A_I$ on a minimizer’s support are linearly dependent, the proof constructs another minimizer. Thus that solution of ($\mathcal {P}_\lambda (y)$) is nonunique.

For a solution $x_\lambda$ with $\ker(A_I)=\{0\}$, the optimality condition for ($\mathcal {P}_\lambda (y)$) gives the implicit formula

This expression generalizes soft thresholding (which is recovered when $A=\mathrm{Id}_N$).

Although the coefficient vector $x_\lambda$ may not be unique, its fitted data $Ax_\lambda$, also called the predictor, are unique.

Multiplying (11.4) by the restricted measurement matrix gives

Thus, up to an $O(\lambda)$ bias, this predictor is an orthogonal projection onto a low-dimensional subspace indexed by $I$.

For $\lambda>0$, a joint density for the entries of $A$ on $\mathbb{R}^{P\times N}$ implies that its columns are in general position almost surely. The Lasso solution is then unique for every $y$. Randomizing only $y$, while keeping $A$ fixed, does not in general ensure uniqueness.

11.1.5 Duality¶

The first-order conditions can also be derived from the convex duality theory in Section 18.3. This interpretation explains why certificates characterize optimal recovery.

For $\lambda>0$, the dual objective (11.6) is strongly concave. Completing the square shows that its maximizer $p_\lambda$ is an orthogonal projection:

Thus the supremum is attained at a unique point for $\lambda>0$. For $\lambda=0$ and feasible data, a dual optimum also exists, though it need not be unique.

11.2 Consistency and Sparsistency¶

11.2.1 Bregman Divergence Rates for General Regularizations¶

Consider the more general regularized problem

for a proper convex regularizer $J$ and $\lambda>0$. We state estimates for any minimizer, without assuming uniqueness.

For $x_0$ with nonempty subdifferential and a chosen $\eta \in \partial J(x_0)$, define the Bregman divergence

One has $D_\eta(x_0|x_0)=0$ and, by convexity, $D_\eta(x|x_0)\geqslant 0$; see Figure 11.5.

If $J$ is differentiable, then $\partial J(x_0)=\{\nabla J(x_0)\}$, and the divergence becomes

If $J$ is strictly convex, then $D(x|x_0)=0$ exactly when $x=x_0$. Thus $D(\cdot|\cdot)$ separates points, although it need not be symmetric or satisfy the triangle inequality.

If $J=|\!| \cdot |\!|^2$, then $D(x|x_0) = |\!| x-x_0 |\!|^2$ is the squared Euclidean distance.

For the negative-entropy functional $J(x)=\sum_i x_i(\log x_i-1)+\iota_{\mathbb{R}_+^N}(x)$, one obtains

This is the generalized Kullback–Leibler divergence for $x_0\in\mathbb{R}_{++}^N$ and $x\in\mathbb{R}_+^N$, with the convention $0\log0=0$. For probability vectors, the linear terms sum to zero.

The following estimate, associated with the work of Burger and Osher, yields a linear noise rate measured by Bregman divergence.

When $|\!| w |\!|>0$ and $p\neq0$, choose $\lambda=|\!| w |\!|/|\!| p |\!|$ in (11.10). The resulting estimate $D_\eta( x_\lambda|x_0 ) \leqslant 2 |\!| w |\!||\!| p |\!|$ is linear in the noise level.

For the simple case of a quadratic regularizer $J(x)=|\!| x |\!|^2/2$, as used in Section 9.3.2, the source condition (11.9) becomes

This is (9.12) with $\beta=1$. Under this condition, (11.10) gives the sublinear $\ell^2$ estimate

This agrees with Theorem 9.4 under the convention $x_0\in\mathop{\mathrm{Im}}((A^*A)^{\beta/2})$.

The source condition (11.9) is sufficient for $x_0$ to minimize $J$ subject to $Ax=Ax_0$. In finite dimension, it is also necessary under a subdifferential sum-rule qualification, for example when $J$ is finite and continuous on $\mathbb{R}^N$. Without such a qualification, existence of a Lagrange multiplier is an additional assumption.

11.2.2 Linear Rates in Norms for $\ell^1$ Regularization¶

For a regularizer $J$ that is not strictly convex, the Bregman bound (11.10) need not control the error in norm. We now examine how additional certificate conditions restore such control for the $\ell^1$ penalty $J=|\!| \cdot |\!|_1$.

The next lemma quantifies this control through the gap between $|\eta_i|$ and one.

We use the convention $|\!| \eta_\emptyset |\!|_\infty=0$. The margin $1-|\!| \eta_{J^c} |\!|_\infty>0$ measures how far $\eta$ lies from saturation on the complementary coordinates. A larger margin gives stronger control of the coefficient error.

This lemma yields a norm convergence rate under the certificate and injectivity assumptions of Proposition 11.7, with $x^\star=x_0$.

This theorem guarantees uniqueness of the noiseless solution $x_0$, but does not imply uniqueness of $x_\lambda$.

The source condition (11.13), combined with the injectivity condition $\ker(A_{J})=\{0\}$, gives the $\ell^2$ stability estimate (11.14).

The estimate (11.14) controls the error linearly relative to the fixed noiseless signal $x_0$. It does not by itself establish pairwise Lipschitz continuity between arbitrary noisy solutions.

Compare this linear rate with the sublinear rates for quadratic regularization in Theorem 9.4.

The source conditions for linear regularization (9.12) and nonlinear regularization (11.13) express different structural assumptions.

In the quadratic case with $\beta=1$, the source condition is $x_0\in\mathop{\mathrm{Im}}(A^*)=\ker(A)^\bot$ in finite dimension. The signal itself must have no component in $\ker(A)$ because squared-norm regularization cannot recover such a component.

The nonlinear condition instead requires a subgradient $\eta$ to belong to $\mathop{\mathrm{Im}}(A^*)$. It can therefore permit recovery of signal components in $\ker(A)$ when the prior supplies sufficient structure.

11.2.3 Sparsistency for Low Noise¶

Theorem 11.11 gives an abstract guarantee whose certificate assumptions may be difficult to verify.

We now construct a candidate certificate explicitly. When it satisfies a strict off-support bound, it guarantees both a linear rate and support stability.

For any solution $x_\lambda$ of $(\mathcal{P}_\lambda(y))$, define the dual certificate as in (11.2). Uniqueness of the fitted data makes it independent of the chosen solution:

The next proposition shows that $p_\lambda$ converges to the dual solution of minimum norm for the constrained problem.

The certificates in $\mathcal{D}_0(y,x_0)$ for $\lambda=0$ may be nonunique, but the limit $\lambda\rightarrow 0$ selects one of them. This selected certificate governs support stability when both the noise and $\lambda$ are small.

The inequality constraint $|\!| \eta_0 |\!|_\infty\leqslant 1$ makes the minimum-norm problem (11.15) difficult to solve explicitly. Dropping this $\ell^\infty$ constraint while retaining interpolation on the support gives the minimum-norm precertificate:

The notation “$\eta_F$” refers to the “Fuchs” certificate, named after J.-J. Fuchs, who first used it to study $\ell^1$ minimization.

The vector $p_F$ may fail dual feasibility because $|\!| \eta_F |\!|_\infty \leqslant 1$ need not hold. This is why the construction gives a precertificate.

If the interpolation system is consistent, $p_F$ is its minimum-norm solution: $A_I^*p=\mathop{\mathrm{sign}}(x_{0,I})$. The pseudoinverse gives $p_F=A_I^{*,+} \mathop{\mathrm{sign}}(x_{0,I})$; see Proposition 9.2. Under $\ker(A_I)=\{0\}$, this simplifies to

With the Gram matrix $C \mathrel{:=}A^* A$, the certificate becomes

The next proposition relates $\eta_F$ to $\eta_0$: whenever $\eta_F$ is feasible, it equals $\eta_0$.

The condition $|\!| \eta_F |\!|_\infty \leqslant 1$ implies that $x_0$ is a solution to ($\mathcal {P}_0(y)$).

A strict off-support bound on $\eta_F$ strengthens this recovery condition and gives both a linear error rate and support stability for small noise.

The feasibility condition $|\!| \eta_F |\!|_\infty \leqslant 1$ has the empty-cap interpretation shown in Figure 11.7. We thank Charles Dossal for pointing out this connection to spherical Delaunay triangulations.

Theorem 11.15 requires small correlated noise $|\!| A^*w |\!|_\infty$ and a small $\lambda$, with their ratio controlled, so that regularization does not eliminate the nonzero coefficients of $x_0$. Theorem 11.11, by contrast, applies at any noise level $|\!| w |\!|$, but gives neither support control nor uniqueness of $x_\lambda$, and its constants can be less favorable.

The proof provides explicit constants in terms of three parameters $K,L,S$:

• $K$ accounts for the conditioning of the operator on the support $I$;

• $L$ is the largest sum of absolute correlations between an off-support atom and all atoms on the support;

• $S$ measures the margin by which $\eta_F$ avoids saturation outside the support.

The bounds on $|\!| A^*w |\!|_\infty/\lambda$ and on $\lambda$ are given by (11.22). For fixed $\gamma>0$ and sufficiently small $\delta>0$, choosing $\lambda=(1+\gamma)R\delta/S$ gives

The strict factor $1+\gamma$ ensures the required strict inequality. Such a choice is primarily theoretical because its constants involve the unknown support.

For a class of signals $x_0$, analyzing $\ell^1$ support recovery therefore starts with testing the Fuchs precertificate $\eta_F$: feasibility requires $|\!| \eta_F |\!|_\infty \leqslant 1$, and support stability requires a strict bound outside the support.

Figure 11.9 illustrates the behavior for random $A$: the precertificate $\eta_F$ becomes nondegenerate when $P$ is sufficiently large. Section 12.2 analyzes this phenomenon.

11.2.4 Sparsistency for Arbitrary Noise¶

The small-noise theorem fixes the sign of the solution in advance. To obtain support containment for arbitrary noise levels, one can instead control all sign patterns on a prescribed support $I$. This gives the exact recovery coefficient of Tropp [32]:

The maximum over an empty complement is taken to be zero. Assume $A_I$ is injective and $\operatorname{ERC}(I)>0$. Let $P_I=A_IA_I^+$ and choose $\lambda>0$ so that

Then the Lasso has a unique solution supported within $I$; some coefficients may vanish, so equality of supports is not asserted.

Indeed, minimize the Lasso objective over vectors supported on $I$. The resulting vector $\hat x_I$ is unique, since $A_I$ is injective. Its optimality condition gives $q_I\in\partial|\!| \cdot |\!|_1(\hat x_I)$ and

For $j\notin I$, the corresponding dual coefficient is bounded by

Thus the restricted minimizer satisfies the full optimality conditions with strict inequality outside $I$, and injectivity proves uniqueness. If $y=Ax_0+w$ and $\mathop{\mathrm{supp}}(x_0)\subset I$, the orthogonal residual above is $(\mathrm{Id}-P_I)w$. The condition therefore states explicitly how large $\lambda$ must be relative to the noise.

11.3 Sparse Deconvolution Case Study¶

For random $A$, Chapter 12 gives probabilistic conditions under which $\eta_F$ is a valid certificate.

Super-resolution exposes a limitation of the Fuchs construction. The columns $(a_i)_i$ of $A$ are often samples of a smooth kernel, so nearby columns are highly correlated.

We consider $a_i=\varphi(z_i)$, where $(z_i)_i \subset \mathbb{X}$ is a sampling grid in a domain $\mathbb{X}$ and $\varphi: \mathbb{X}\rightarrow \mathcal{H}$ is a smooth map. One has

The coefficients of a sparse vector $x$ are the weights of the discrete measure $m_{x} \mathrel{:=}\sum_{i=1}^N x_i \delta_{z_i}$, whose atoms lie on the reconstruction grid.

The matrix $A$ discretizes an operator on Radon measures, $\mathcal{A}: m \in \mathcal{M}(\mathbb{X}) \mapsto y = \mathcal{A}m \in \mathcal{H}$, defined by

For a discrete measure, the two representations agree: $\mathcal{A}(m_x)=A x$.

For example, take $\mathcal{H}=L^2(\mathbb{X})$ on $\mathbb{X}=\mathbb{R}^d$ or $\mathbb{X}=\mathbb{T}^d$ and let $\varphi(z) = \tilde\varphi(\cdot-z)$. The resulting operator is convolution:

The observation space $\mathcal{H}$ is infinite-dimensional in this example. Sampling the output instead gives $\varphi(z) = (\tilde\varphi(r_j-z))_{j=1}^P$. The observation grid $r \in \mathbb{X}^P$ need not coincide with the reconstruction grid $z \in \mathbb{X}^N$.

A closely related example uses $\varphi(z)=(e^{\mathrm{i}k z})_{k=-f_c}^{f_c}$ on $\mathbb{X}=\mathbb{T}$, so that $\mathcal{A}$ corresponds to computing the $2f_c+1$ low-frequency coefficients of the Fourier transform of the measure

The operator $\mathcal{A}^*\mathcal{A}$ is a convolution against an ideal low-pass (Dirichlet) kernel. By weighting the Fourier coefficients, one can model low-pass filters on the torus.

On $\mathbb{X}=\mathbb{R}^+$, another example is the Laplace transform:

Define the continuous covariance kernel by

The function $\mathcal{C}$ is the kernel of $\mathcal{A}^* \mathcal{A}$.

The discrete covariance, defined on the computational grid, is $C=(\mathcal{C}(z_i,z_{i'}))_{i,i'} \in \mathbb{R}^{N \times N}$, while its restriction to some support set $I$ is $C_{I,I}=(\mathcal{C}(z_i,z_{i'}))_{(i,i')\in I^2}\in\mathbb{R}^{|I|\times|I|}$.

By (11.17), the vector $\eta_F$ samples the continuous precertificate $\tilde\eta_F$ on the grid:

Thus $\eta_F$ samples a linear combination of $|I|$ kernel functions $(\mathcal{C}(x,z_i))_{i \in I}$.

The question is whether $|\!| \eta_F |\!|_{\ell^\infty} \leqslant 1$. Along increasingly dense grids, continuity implies that a uniform bound on the sampled certificate requires $|\!| \tilde\eta_F |\!|_{L^\infty}\leqslant 1$. A finite grid may miss an overshoot between samples. The construction constrains $\tilde\eta_F$ only to interpolate $\mathop{\mathrm{sign}}(x_{0,i})$ at support points $z_i$; it does not prevent the function from leaving $[-1,1]$ nearby. Figure 11.11 illustrates this fact.

In one dimension, a certificate bounded in magnitude by one must satisfy $\eta'(z_i)=0$ at every interior support point $i \in I$, where it already interpolates the sign. Adding these necessary conditions gives the minimum-norm precertificate with vanishing derivatives:

Here we restrict to a one-dimensional domain and assume that the interpolation constraints are independent. The adjoint is $(\mathcal{A}^*p)(z)=\langle\varphi(z),p\rangle_\mathcal{H}$, and $\eta_V=(\tilde\eta_V(z_i))_{i=1}^N$. A derivative constraint alone does not ensure $|\tilde\eta_V|\leqslant 1$; this bound must be checked separately. As in (11.23), this precertificate is a linear combination of kernel functions, now with $2|I|$ terms:

where $\partial_2\mathcal{C}$ is the derivative of $\mathcal{C}$ with respect to the second variable, and $(b,c)$ solve a $2|I| \times 2|I|$ linear system

Evaluating $\tilde\eta_V=\mathcal{A}^*p_V$ on the reconstruction grid gives the discrete precertificate $\eta_V$.

Figure 11.11 shows better behavior of $\eta_V$ than of $\eta_F$. If $|\!| \eta_V |\!|_\infty\leqslant 1$ and $A$ is injective on its saturation set, Theorem 11.11 gives a linear rate in the grid-based $\ell^2$ norm. As the grid is refined, however, this coefficient $\ell^2$ norm does not provide an intrinsic distance between spike measures, which need not have $L^2$ densities.

When $\eta_V$ differs from $\eta_F$, one cannot directly apply Theorem 11.15: exact support stability on increasingly fine grids can fail in super-resolution problems.

Methods without a fixed grid provide a more suitable framework for such stability results. In a grid-free formulation, a valid vanishing-derivative precertificate coincides with the minimum-norm certificate under the corresponding interpolation and nondegeneracy assumptions. This leads to stability results for spike locations and amplitudes, rather than a fixed-grid coefficient norm.