Let J(x)J(x) be a convex sparsity-promoting cost function (i.e., J(x)J(x) is small for sparse xx.) To recover a sparse signal representation x^x^ from measurements y=Φx,Φ∈RM×Ny=Φx,Φ∈RM×N, we may either solve

when there is no noise, or solve

when there is noise in the measurements. Here, HH is a cost function that penalizes the distance between the vectors ΦxΦx and yy. For an appropriate penalty parameter μμ, Equation 2 is equivalent to the unconstrained formulation:

for some μ>0μ>0. The parameter μμ may be chosen by trial-and-error, or by statistical techniques such as cross-validation [2].

For convex programming algorithms, the most common choices of JJ and HH are usually chosen as follows: J(x)=∥x∥1J(x)=∥x∥1, the ℓ1ℓ1-norm of xx, and H(Φx,y)=12∥Φx-y∥22H(Φx,y)=12∥Φx-y∥22, the ℓ2ℓ2-norm of the error between the observed measurements and the linear projections of the target vector xx. In statistics, minimizing this HH subject to ∥x∥1≤δ∥x∥1≤δ is known as the Lasso problem. More generally, J(·)J(·) acts as a regularization term and can be replaced by other, more complex, functions; for example, the desired signal may be piecewise constant, and simultaneously have a sparse representation under a known basis transform ΨΨ. In this case, we may use a mixed regularization term:

It might be tempting to use conventional convex optimization packages for the above formulations (Equation 1, Equation 2, and Equation 3). Nevertheless, the above problems pose two key challenges which are specific to practical problems encountered in CS: (i) real-world applications are invariably large-scale (an image of a resolution of 1024×10241024×1024 pixels leads to optimization over a million variables, well beyond the reach of any standard optimization software package); (ii) the objective function is nonsmooth, and standard smoothing techniques do not yield very good results. Hence, for these problems, conventional algorithms (typically involving matrix factorizations) are not effective or even applicable. These unique challenges encountered in the context of CS have led to considerable interest in developing improved sparse recovery algorithms in the optimization community.

In the noiseless case, the ℓ1ℓ1-minimization problem (obtained by substituting J(x)=∥x∥1J(x)=∥x∥1 in Equation 1) can be recast as a linear program (LP) with equality constraints. These can be solved in polynomial time (O(N3)O(N3)) using standard interior-point methods [3]. This was the first feasible reconstruction algorithm used for CS recovery and has strong theoretical guarantees, as shown earlier in this course. In the noisy case, the problem can be recast as a second-order cone program (SOCP) with quadratic constraints. Solving LPs and SOCPs is a principal thrust in optimization research; nevertheless, their application in practical CS problems is limited due to the fact that both the signal dimension NN, and the number of constraints MM, can be very large in many scenarios. Note that both LPs and SOCPs correspond to the constrained formulations in Equation 1 and Equation 2 and are solved using first order interior-point methods.

A newer algorithm called “l1_ls" [12] is based on an interior-point algorithm that uses a preconditioned conjugate gradient (PCG) method to approximately solve linear systems in a truncated-Newton framework. The algorithm exploits the structure of the Hessian to construct their preconditioner; thus, this is a second order method. Computational results show that about a hundred PCG steps are sufficient for obtaining accurate reconstruction. This method has been typically shown to be slower than first-order methods, but could be faster in cases where the true target signal is highly sparse.

As opposed to solving the constrained formulation, an alternate approach is to solve the unconstrained formulation in Equation 3. A widely used method for solving ℓ1ℓ1-minimization problems of the form

for a convex and differentiable HH, is an iterative procedure based on shrinkage (also called soft thresholding; see Equation 6 below). In the context of solving Equation 5 with a quadratic HH, this method was independently proposed and analyzed in [1], [9], [13], [14], and then further studied or extended in [4], [5], [6], [7], [11], [16]. Shrinkage is a classic method used in wavelet-based image denoising. The shrinkage operator on any scalar component can be defined as follows:

This concept can be used effectively to solve Equation 5. In particular, the basic algorithm can be written as following the fixed-point iteration: for i=1,...,Ni=1,...,N, the i th i th coefficient of xx at the (k+1) th (k+1) th time step is given by

where τ>0τ>0 serves as a step-length for gradient descent (which may vary with kk) and μμ is as specified by the user. It is easy to see that the larger μμ is, the larger the allowable distance between xk+1xk+1 and xkxk. For a quadratic penalty term H(·)H(·), the gradient ▽H▽H can be easily computed as a linear function of xkxk; thus each iteration of Equation 7 essentially boils down to a small number of matrix-vector multiplications.

The simplicity of the iterative approach is quite appealing, both from a computational, as well as a code-design standpoint. Various modifications, enhancements, and generalizations to this approach have been proposed, both to improve the efficiency of the basic iteration in Equation 7, and to extend its applicability to various kinds of JJ [8], [10], [16]. In principle, the basic iteration in Equation 7 would not be practically effective without a continuation (or path-following) strategy [11], [16] in which we choose a gradually decreasing sequence of values for the parameter μμ to guide the intermediate iterates towards the final optimal solution.

This procedure is known as continuation; in [11], the performance of an algorithm known as Fixed-Point Continuation (FPC) has been compared favorably with another similar method known as Gradient Projection for Sparse Reconstruction (GPSR) [10] and “l1_ls” [12]. A key aspect to solving the unconstrained optimization problem is the choice of the parameter μμ. As discussed above, for CS recovery, μμ may be chosen by trial and error; for the noiseless constrained formulation, we may solve the corresponding unconstrained minimization by choosing a large value for μμ.

In the case of recovery from noisy compressive measurements, a commonly used choice for the convex cost function H(x)H(x) is the square of the norm of the residual. Thus we have:

For this particular choice of penalty function, Equation 7 reduces to the following iteration:

which is run until convergence to a fixed point. The algorithm is detailed in pseudocode form below.

Inputs: CS matrix ΦΦ, signal measurements yy, parameter sequence μnμn
Outputs: Signal estimate x^x^
initialize: x^0=0x^0=0, r=yr=y, k=0k=0.
while ħalting criterion false do
    1. k←k+1k←k+1
    2. x←x^-τΦTrx←x^-τΦTr {take a gradient step}
    3. x^← shrink (x,μkτ)x^← shrink (x,μkτ) {perform soft thresholding}
    4. r←y-Φx^r←y-Φx^ {update measurement residual}
end while
return x^←x^x^←x^

It turns out that an efficient method to obtain the solution to the constrained optimization problem in Equation 1 can be devised by solving a small number of the unconstrained problems in the form of Equation 3. These subproblems are commonly referred to as Bregman iterations. A simple version can be written as follows:

The problem in the second step can be solved by the algorithms reviewed above. Bregman iterations were introduced in [15] for constrained total variation minimization problems, and was proved to converge for closed, convex functions J(x)J(x). In [17], it is applied to Equation 1 for J(x)=∥x∥1J(x)=∥x∥1 and shown to converge in a finite number of steps for any μ>0μ>0. For moderate μμ, the number of iterations needed is typically lesser than 5. Compared to the alternate approach that solves Equation 1 through directly solving the unconstrained problem in Equation 3 with a very large μμ, Bregman iterations are often more stable and sometimes much faster.

All the methods discussed in this section optimize a convex function (usually the ℓ1ℓ1-norm) over a convex (possibly unbounded) set. This implies guaranteed convergence to the global optimum. In other words, given that the sampling matrix ΦΦ satisfies the conditions specified in "Signal recovery via ℓ1ℓ1 minimization", convex optimization methods will recover the underlying signal xx. In addition, convex relaxation methods also guarantee stable recovery by reformulating the recovery problem as the SOCP, or the unconstrained formulation.