The Convex Geometry of Linear Inverse Problems

In applications throughout science and engineering one is often faced with the challenge of solving an ill-posed inverse problem, where the number of available measurements is smaller than the dimension of the model to be estimated. However in many practical situations of interest, models are constrained structurally so that they only have a few degrees of freedom relative to their ambient dimension. This paper provides a general framework to convert notions of simplicity into convex penalty functions, resulting in convex optimization solutions to linear, underdetermined inverse problems. The class of simple models considered includes those formed as the sum of a few atoms from some (possibly infinite) elementary atomic set; examples include well-studied cases from many technical fields such as sparse vectors (signal processing, statistics) and low-rank matrices (control, statistics), as well as several others including sums of a few permutation matrices (ranked elections, multiobject tracking), low-rank tensors (computer vision, neuroscience), orthogonal matrices (machine learning), and atomic measures (system identification). The convex programming formulation is based on minimizing the norm induced by the convex hull of the atomic set; this norm is referred to as the atomic norm. The facial structure of the atomic norm ball carries a number of favorable properties that are useful for recovering simple models, and an analysis of the underlying convex geometry provides sharp estimates of the number of generic measurements required for exact and robust recovery of models from partial information. These estimates are based on computing the Gaussian widths of tangent cones to the atomic norm ball. When the atomic set has algebraic structure the resulting optimization problems can be solved or approximated via semidefinite programming. The quality of these approximations affects the number of measurements required for recovery, and this tradeoff is characterized via some examples. Thus this work extends the catalog of simple models (beyond sparse vectors and low-rank matrices) that can be recovered from limited linear information via tractable convex programming.

1. A spherical cap is a subset of the sphere obtained by intersecting the sphere \(\mathbb{S}^{p-1}\) with a halfspace.

2. While Proposition 3.15 follows as a consequence of the general result in Corollary 3.14, one can remove the constant factor 9 in the statement of Proposition 3.15 by carrying out a more refined analysis of the Birkhoff polytope.

This work was supported in part by AFOSR grant FA9550-08-1-0180, in part by a MURI through ARO grant W911NF-06-1-0076, in part by a MURI through AFOSR grant FA9550-06-1-0303, in part by NSF FRG 0757207, in part through ONR award N00014-11-1-0723, and NSF award CCF-1139953.

We gratefully acknowledge Holger Rauhut for several suggestions on how to improve the presentation in Sect. 3, and Amin Jalali for pointing out an error in a previous draft. We thank Santosh Vempala, Joel Tropp, Bill Helton, Martin Jaggi, and Jonathan Kelner for helpful discussions. Finally, we acknowledge the suggestions of the associate editor Emmanuel Candès as well as the comments and pointers to references made by the reviewers, all of which improved our paper.

Authors and Affiliations

1. Department of Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA, 91125, USA

   Venkat Chandrasekaran

2. Computer Sciences Department, University of Wisconsin, Madison, WI, 53706, USA

   Benjamin Recht

3. Laboratory for Information and Decision Systems, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA, 02139, USA

   Pablo A. Parrilo & Alan S. Willsky

Corresponding author

Correspondence to Venkat Chandrasekaran.

Communicated by Emmanuel Candès.

Appendix A: Proof of Proposition 3.6

Proof

First note that the Gaussian width can be upper-bounded as follows:

where \(\mathcal{B}(0,1)\) denotes the unit Euclidean ball. The expression on the right-hand side inside the expected value can be expressed as the optimal value of the following convex optimization problem for each g∈ℝ^p:

We now proceed to form the dual problem of (31) by first introducing the Lagrangian

where \(\mathbf{u}\in\mathcal{C}^{\ast}\) and γ≥0 is a scalar. To obtain the dual problem we maximize the Lagrangian with respect to z, which amounts to setting

Putting this into the Lagrangian above gives the dual problem

Solving this optimization problem with respect to γ we find that \(\gamma= \frac{1}{2} \|\mathbf{g}-\mathbf{u}\|\), which gives the dual problem to (31)

Under very mild assumptions about \(\mathcal{C}\), the optimal value of (32) is equal to that of (31) (for example as long as \(\mathcal{C}\) has a nonempty relative interior, strong duality holds). Hence we have derived

This equation combined with the bound (30) gives us the desired result. □

Appendix B: Proof of Theorem 3.9

Proof

We set \(\beta= \tfrac{1}{\varTheta}\). First note that if \(\beta\geq \tfrac{1}{4} \exp\{\tfrac{p}{9}\}\) then the width bound exceeds \(\sqrt{p}\), which is the maximal possible value for the width of \(\mathcal{C}\). Thus, we will assume throughout that \(\beta\leq \tfrac{1}{4} \exp\{\tfrac{p}{9}\}\).

Using Proposition 3.6 we need to upper-bound the expected distance to the polar cone. Let \(\mathbf{g}\sim\mathcal {N}(0,I)\) be a normally distributed random vector. Then the norm of g is independent from the angle of g. That is, ∥g∥ is independent from g/∥g∥. Moreover g/∥g∥ is distributed as a uniform sample on \(\mathbb{S}^{p-1}\), and \(\operatorname{\mathbb{E}}_{\mathbf{g}}[\|\mathbf{g}\|]\leq\sqrt{p}\). Thus we have

where u is sampled uniformly on \(\mathbb{S}^{p-1}\).

To bound the latter quantity, we will use isoperimetry. Suppose A is a subset of \(\mathbb{S}^{p-1}\) and B is a spherical cap with the same volume as A. Let N(A,r) denote the locus of all points in the sphere of Euclidean distance at most r from the set A. Let μ denote the Haar measure on \(\mathbb{S}^{p-1}\) and let μ(A;r) denote the measure of N(A,r). Then spherical isoperimetry states that μ(A;r)≥μ(B;r) for all r≥0 (see, for example, [48, 53]).

Let B now denote a spherical cap with \(\mu(B)=\mu(\mathcal{C}^{\ast}\cap\mathbb{S}^{p-1})\). Then we have

(35)

(36)

(37)

where the first equality is the integral form of the expected value and the last inequality follows by isoperimetry. Hence we can bound the expected distance to the polar cone intersecting the sphere using only knowledge of the volume of spherical caps on \(\mathbb{S}^{p-1}\).

To proceed let v(φ) denote the volume of a spherical cap subtending a solid angle φ. An explicit formula for v(φ) is

where \(z_{p} = \int_{0}^{\pi}\sin^{p-1}(\vartheta)\,\mathrm{d}\vartheta\) [45]. Let φ(β) denote the minimal solid angle of a cap such that β copies of that cap cover \(\mathbb {S}^{p-1}\). Since the geodesic distance on the sphere is always greater than or equal to the Euclidean distance, if K is a spherical cap subtending ψ radians, μ(K;t)≥v(ψ+t). Therefore

We can proceed to simplify the right-hand side integral:

(40)

(41)

(42)

(43)

(44)

(45)

(46)

(43) follows by switching the order of integration, and the rest of these equalities follow by straightforward integration and some algebra.

Using the inequalities that \(z_{p} \geq\frac{2}{\sqrt{p-1}}\) (see [48]) and sin(x)≤exp(−(x−π/2)²/2) for x∈[0,π], we can bound the last integral as

(47)

Performing the change of variables \(a = \sqrt{p-1}(\vartheta-\tfrac {\pi}{2})\), we are left with the integral

(48)

(49)

(50)

In this final bound, we bounded the first term by dropping the upper integrand, and for the second term we used the fact that

We are now left with the task of computing a lower bound for φ(β). We need to first reparameterize the problem. Let K be a spherical cap. Without loss of generality, we may assume that

for some h∈[0,1]. Here h is the height of the cap over the equator. Via elementary trigonometry, the solid angle that K subtends is given by π/2−sin⁻¹(h). Hence, if h(β) is the largest number such that β caps of height h cover \(\mathbb {S}^{p-1}\), then h(β)=sin(π/2−φ(β)).

The quantity h(β) may be estimated using the following estimate from [11]. For h∈[0,1], let γ(p,h) denote the volume of a spherical cap of \(\mathbb{S}^{p-1}\) of height h.

Lemma B.1

(See [11])

For \(1\geq h\geq\frac{2}{\sqrt{p}}\),

Note that for \(h \geq\frac{2}{\sqrt{p}}\),

So if

then h≤1 because we have assumed \(\beta\leq\tfrac{1}{4} \exp\{ \tfrac{p}{9}\}\) and p≥9. Moreover, \(h\geq\frac{2}{\sqrt{p}}\) and the volume of the cap with height h is less than or equal to 1/β. That is

Combining the estimate (50) with Proposition 3.6, and using our estimate for φ(β), we get the bound

(57)

This expression can be simplified by using the following bounds. First, sin⁻¹(x)≥x lets us upper-bound the first term by \(\sqrt{\frac{p}{p-1}}\frac{1}{8\beta}\). For the second term, using the inequality \(\sin^{-1}(x)\leq\tfrac{\pi}{2}x\) results in the upper bound

For p≥9 the upper bound can be expressed simply as \(w(\mathcal{C})\leq3\sqrt{\log(4 \beta)}\). We recall that \(\beta= \tfrac{1}{\varTheta}\), which completes the proof of the theorem.  □

Appendix C: Direct Width Calculations

We first give the proof of Proposition 3.10.

Proof

Let x ^⋆ be an s-sparse vector in ℝ^p with ℓ ₁ norm equal to 1, and let \(\mathcal{A}\) denote the set of unit-Euclidean-norm one-sparse vectors. Let Δ denote the set of coordinates where x ^⋆ is nonzero. The normal cone at x ^⋆ with respect to the ℓ ₁ ball is given by

(59)

(60)

Here Δ^c represents the zero entries of x ^⋆. The minimum squared distance to the normal cone at x ^⋆ can be formulated as a one-dimensional convex optimization problem for arbitrary z∈ℝ^p

(61)

(62)

where

is the ℓ ₁-shrinkage function. Hence, for any fixed t≥0 independent of g, we have

(64)

(65)

Now we directly integrate the second term, treating each summand individually. For a zero-mean, unit-variance normal random variable g,

(66)

(67)

(68)

(69)

(70)

The first simplification follows because the shrink function and Gaussian distributions are symmetric about the origin. The second equality follows by integrating by parts. The inequality follows by a tight bound on the Gaussian Q-function:

Using this bound, we get

Setting \(t= \sqrt{2\log(p/s)}\) gives

The last inequality follows because

whenever 0≤s≤p. □

Next we give the proof of Proposition 3.11.

Proof

Let x ^⋆ be an m ₁×m ₂ matrix of rank r with singular value decomposition UΣV ^∗, and let \(\mathcal{A}\) denote the set of rank-one unit-Euclidean-norm matrices of size m ₁×m ₂. Without loss of generality, impose the conventions m ₁≤m ₂, Σ is r×r, U is m ₁×r, and V is m ₂×r, and assume the nuclear norm of x ^⋆ is equal to 1.

Let u _k (respectively v _k ) denote the k’th column of U (respectively V). It is convenient to introduce the orthogonal decomposition \({\mathbb{R}}^{m_{1} \times m_{2}} = \Delta\oplus\Delta ^{\perp}\) where Δ is the linear space spanned by elements of the form u _k z ^T and \(\mathbf{y}\mathbf{v}_{k}^{T}\), 1≤k≤r, where z and y are arbitrary, and Δ^⊥ is the orthogonal complement of Δ. The space Δ^⊥ is the subspace of matrices spanned by the family (yz ^T), where y (respectively z) is any vector orthogonal to all the columns of U (respectively V). The normal cone of the nuclear norm ball at x ^⋆ is given by the cone generated by the subdifferential at x ^⋆:

(75)

(76)

Note that here \(\|Z\|_{\mathcal{A}}^{\ast}\) is the operator norm, equal to the maximum singular value of Z [63].

Let G be a Gaussian random matrix with i.i.d. entries, each with mean zero and unit variance. Then the matrix

is in the normal cone at x ^⋆. We can then compute

(78)

(79)

(80)

Here (79) follows because \(\mathcal{P}_{\Delta }(G)\) and \(\mathcal{P}_{\Delta^{\perp}}(G)\) are independent. The final line follows because dim(T)=r(m ₁+m ₂−r) and the Frobenius (i.e., Euclidean) norm of UV ^∗ is \(\|UV^{*}\|_{F}=\sqrt{r}\). Due to the isotropy of Gaussian random matrices, \(\mathcal{P}_{\Delta^{\perp}}(G)\) is identically distributed as an (m ₁−r)×(m ₂−r) matrix with i.i.d. Gaussian entries each with mean zero and variance one. We thus know that

(see, for example, [22]). To bound the latter expectation, we again use the integral form of the expected value. Letting \(\mu_{T^{\perp}}\) denote the quantity \(\sqrt{m_{1}-r}+\sqrt{m_{2}-r}\), we have

(82)

(83)

(84)

(85)

(86)

(87)

Using this bound in (80), we obtain

(88)

(89)

(90)

where the second inequality follows from the fact that (a+b)²≤2a ²+2b ². We conclude that 3r(m ₁+m ₂−r) random measurements are sufficient to recover a rank r, m ₁×m ₂ matrix using the nuclear norm heuristic. □

Reprints and permissions