Augmented L1 and Nuclear-Norm Models with a Globally Linearly Convergent Algorithm

Introduction

Sparse vector recovery and low-rank matrix recovery problems have drawn lots of attention from researchers in different fields in the past several years. They have wide applications in compressive sensing, signal/image processing, machine learning, etc. The fundamental problem of sparse vector recovery is to find the vector with (nearly) fewest nonzero entries from an underdetermined linear system ${\mathbf{A}}{\mathbf{x}}={\mathbf{b}}$, and that of low-rank matrix recovery is to find a matrix of (nearly) lowest rank from an underdetermined ${\mathcal{A}}({\mathbf{X}})={\mathbf{b}}$, where ${\mathcal{A}}$ is a linear operator.

To recover a sparse vector ${\mathbf{x}}^{0}$, a well-known model is the basis pursuit problem :

For vector ${\mathbf{b}}$ with noise or generated by an approximately sparse vector, a variant of (1) is

where $\|{\mathbf{X}}\|_{*}$ equals the summation of the singular values of ${\mathbf{X}}$. Similar to (2), a useful variant of (3) is

The nonsmooth objective functions in problems (1)–(4) pose numerical challenges. We augment or “smooth” them by adding $\frac{1}{2\alpha}\|{\mathbf{x}}\|_{2}^{2}$ or $\frac{1}{2\alpha}\|{\mathbf{X}}\|_{F}^{2}$, where $\alpha$ is a positive scalar. We argue that minimizing the augmented objective $\|{\mathbf{x}}\|_{1}+\frac{1}{2\alpha}\|{\mathbf{x}}\|_{2}^{2}$, as well as $\|{\mathbf{X}}\|_{*}+\frac{1}{2\alpha}\|{\mathbf{X}}\|_{F}^{2}$, leads to fast numerical algorithms because not only accurate solutions can be obtained by using a sufficiently large, yet not excessive large, value of $\alpha$, but the Lagrange dual problems are also continuously differentiable and subject to gradient-based acceleration techniques such as line search.

Next, we briefly review the related works and summarize the contributions of this paper. The augmented model for (1) is

which can be solved by the linearized Bregman algorithm (LBreg) , which is analyzed in . (Note that LBreg is different from the Bregman algorithm , which solves problem (1) instead of (5).)

The exact regularization property of (5) is proved in : the solution to (5) is also a solution to (1) as long as $\alpha$ is sufficiently large. The property can also be obtained from . However, neither paper tells how to select $\alpha$, whereas the size of $\alpha$ affects the numerical performance. It has been observed by several groups of researchers that a larger $\alpha$ tends to cause slower convergence. Hence, one would like to choose a moderate $\alpha$ that is just large enough for (5) to return a solution to (1). For recovering a sparse vector ${\mathbf{x}}^{0}$ and a low-rank matrix ${\mathbf{X}}^{0}$, this paper gives the simple formulae

respectively, where the operator norm $\|{\mathbf{X}}^{0}\|_{2}$ equals the maximum singular value of ${\mathbf{X}}^{0}$. Although ${\mathbf{x}}^{0}$ and ${\mathbf{X}}^{0}$ are not known when $\alpha$ must be set, $\|{\mathbf{x}}^{0}\|_{\infty}$ and $\|{\mathbf{X}}^{0}\|_{2}$ are often easy to estimate. For example, in compressive sensing, $\|{\mathbf{x}}^{0}\|_{\infty}$ is the maximum intensity of the underlying signal or the maximum sensor reading. When the total energy $\|{\mathbf{x}}^{0}\|_{2}$ is roughly known, one can apply the more conservative formula: $\alpha\geq 10\|{\mathbf{x}}^{0}\|_{2}$ since $\|{\mathbf{x}}^{0}\|_{2}\geq\|{\mathbf{x}}^{0}\|_{\infty}$. Similarly, a more conservative formula is $\alpha\geq 10\|{\mathbf{X}}^{0}\|_{F}$ for the matrix case.

This paper also shows that the Lagrange dual problem of (5) is unconstrained and differentiable, and its objective is uniformly strongly convex when restricted to certain pairs of points. Consequently, algorithm LBreg, as well as two faster variants, enjoys global linear convergence; specifically, both the objective error and solution error are bounded by $O(\mu^{k})$, where $k$ is the iteration number and $\mu$ is a constant strictly less than $1$. The value of $\mu$ depends on $\alpha$, the dynamic range of the solution’s nonzero entries, as well as some properties of ${\mathbf{A}}$. Although several first-order algorithms for (1) have been shown to have asymptotic linear convergence, this is the first global linear convergence result that comes with an explicit rate.

We shall discuss strong convexity. Many of the algorithms for recovering sparse solutions from under-determined systems of equations are observed to have a linearly converging behavior, at least on problems that are not severely “ill-conditioned”; however, their underlying objective functions do not have strong convexity – a property commonly used to ensure global linear convergence – when the linear operator ${\mathbf{A}}$ has fewer rows than columns. Specifically, the loss function in the form of $g({\mathbf{A}}{\mathbf{x}}-{\mathbf{b}})$, even for strongly convex function $g$, is “flat” along many directions. Flatness or near flatness along a direction means a small directional gradient, which can generally cause slow decrease in the objective value. However, in problems with certain types of matrix ${\mathbf{A}}$, moving along these directions will significantly change the regularization function. In the recent paper , the definition of strong convexity is extended to include a relaxation term involving the regularizer function. The paper argues that, with high probability for problems with ${\mathbf{A}}$ that is random or satisfies restricted eigenvalue or other suitable properties, their “restricted strong convexity” definition is satisfied by the sum of the regularization and loss functions, and as a result, the prox-linear or gradient projection iteration applied to minimizing the sum has a (nearly-)linear convergence behavior, specifically,

where $c<1$, ${\mathbf{x}}^{*}$ and ${\mathbf{x}}^{0}$ are the minimizer and underlying true signal, respectively, and ${\mathbf{x}}^{(k)}$ stands for the $k$th iterate. This paper presents a different approach. Due to smoothing, unmodified linear convergence to the exact solution is achievable without a probabilistic argument. The Lagrange dual of (5) is strongly convex, not in the global sense, but restricted between the current point and its projection to the solution set. This property turns out to be sufficient for global linear convergence without a modification.

Numerically, LBreg without acceleration is not very efficient because it is equivalent to the dual gradient ascent with a fixed step size, as shown in . Nonetheless, the step size can be relaxed. Since the augmentation term $\frac{1}{2\alpha}\|{\mathbf{x}}\|_{2}^{2}$ makes the dual problem unconstrained and differentiable, the dual is subject to advanced gradient-descent techniques such as Barzilai-Borwein (BB) step sizes , non-monotone line search, Nesterov’s technique , as well as semi-smooth Newton methods. Indeed, LBreg has been improved in several recent works: applies a kicking trick; considers applying BB step sizes and non-monotone line search, as well as the limited memory BFGS method ; applies the alternating direction method to the Lagrange dual of (5); applies Nesterov’s technique and obtains the convergence rate $O(1/k^{2})$. Based on the restricted strong convexity of the dual objective and some existing proofs, we theoretically show and numerically demonstrate that LBreg with BB step sizes with non-monotone line search also enjoys global linear convergence.

LBreg has also been extended to recovering simply structured matrices. The algorithms SVT for matrix completion and IT for robust principal components are of the LBreg type, namely, they are gradient iterations that solve

respectively, where $\Omega$ is the set of the observed matrix entries and $\|{\mathbf{S}}\|_{1}=\sum_{i,j}|S_{i,j}|$. shows that the exact regularization property for the vector case also holds for (6) and (7). Although this paper does not analyze (6) and (7) specifically, it gives recovery guarantees for models

assuming $\alpha\geq 10\|{\mathbf{X}}^{0}\|_{2}$.

The matlab codes and demos of LBreg, including the original, line search, and Nesterov acceleration versions, can be found from the second author’s homepage.

Eliminating ${\mathbf{z}}$ from the last equation gives the following dual problem.

It is interesting to compare (10) with the Lagrange dual of (1):

Instead of confining each component of ${\mathbf{A}}^{\top}{\mathbf{y}}$ to $$, (10) applies quadratic penalty to the violation. This leads to its advantage of being unconstrained and differentiable (despite the presence of projection).

The gradient of the last term in (10) is $\alpha{\mathbf{A}}\shrink({\mathbf{A}}^{\top}{\mathbf{y}})$. Furthermore, given a solution ${\mathbf{y}}^{*}$ to (10), one can recover the solution ${\mathbf{x}}^{*}=\alpha\shrink({\mathbf{A}}^{\top}{\mathbf{y}}^{*})$ to (5) (since (10) has a vanishing gradient ${\mathbf{A}}{\mathbf{x}}^{*}-{\mathbf{b}}=\mathbf{0}$, and ${\mathbf{x}}^{*}$ and ${\mathbf{y}}^{*}$ lead to 0-gap primal and dual objectives, respectively). Therefore, solving (10) solves (5), and it is easier than solving (1). In particular, (10) enjoys a rich set of classical techniques such as line search, Barzilai-Borwein steps , semi-smooth Newton methods, Nesterov’s acceleration , which do not directly apply to problems (1) or (11).

The objective of (13) is differentiable except at ${\mathbf{y}}=\mathbf{0}$. However, this is not an issue since ${\mathbf{y}}=\mathbf{0}$ is a solution to (13) only if ${\mathbf{x}}=\mathbf{0}$ is the solution to (12). In other words, (13) is practically differentiable and thus also amenable to classical gradient-based acceleration techniques.

Equality-constrained augmented $\|\cdot\|_{*}$: The primal and dual of the augmented model of (3) are (8) and

respectively, where ${\mathcal{A}}^{*}{\mathbf{y}}:=\sum_{i=1}^{m}y_{i}{\mathbf{A}}_{i}$ and $\{{\mathbf{X}}:\|{\mathbf{X}}\|_{2}\leq 1\}$ is the set of $n_{1}$-by-$n_{2}$ matrices with spectral norms no more than 1. In (14), inside the Frobenius norm is the singular value soft-thresholding of ${\mathcal{A}}^{*}{\mathbf{y}}$.

The primal and dual of the augmented model (4) are (9) and

respectively. Like the augmented models for vectors, problems (14) and (15) are practically differentiable and thus also amenable to advanced optimization techniques for unconstrained differentiable problems.

Recovery Guarantees

First, we present some numerical simulations to motivate the subsequent analysis.

We are interested in comparing model (5) to model (1), whose the performance on recovering sparse solutions have been widely studied. To this end, we conducted three sets of simulations. Without loss of generality, we fixed $\|{\mathbf{x}}^{0}\|_{\infty}=1$ and solved (1) and then (5) with $\alpha=1,10$, and $25$ to reconstruct signals of $n=400$ dimensions. We set the signal sparsity $k=1,2,\ldots,80$ and the number of measurements $m=40,41,\ldots,200$. The entries of ${\mathbf{A}}$ were sampled from the standard Gaussian distribution.

It turns out that the recovery performance of (5) depends on the decay speed of the nonzero entries of the signal ${\mathbf{x}}^{0}$. So, we tested three decay speeds: (i) flat magnitude — no decay, (ii) independent Gaussian — moderate decay, and (iii) power law — fast decay. In the power-law decay, the $i$th largest entry had magnitude $i^{-2}$ and a random sign.

For each $(m,k)$, 100 independent tests were run, and the average of

was recorded, where ${\mathbf{x}}^{*}$ stands for a solution of either (1) or (5). The slightly smoothed cut-off curves at two different levels of relative errors are depicted in Figure 1. Above each curve is the region where a model fails to recover the signals to the specified average relative error. Hence, a higher curve means fewer fails and thus better recovery performance.

In all tests, the best curve is from BP or model (1). Closely following it are those of $\alpha=25$ and $\alpha=10$ of model (5). As long as $\alpha\geq 10$, model (5) is as good as model (1) up to a negligible difference.

The curve of $\alpha=1$ is noticeably lower than others when the signal is flat or decays slowly. For this reason, we do not recommend using $\alpha=\|{\mathbf{x}}^{0}\|_{\infty}$ for model (5) unless when the underlying signals decay very fast.

The differences of the fours curves are very similar across the two levels $10^{-3}$ and $10^{-5}$ of relative errors. We tested other levels and found the same. Therefore, the performance differences are independent of the error level chosen to plot the curves.

2 Null space property

Matrix ${\mathbf{A}}$ satisfies the NSP if

holds for all ${\mathbf{h}}\in\Null({\mathbf{A}})$ and coordinate sets ${\mathcal{S}}\subset\{1,2,\cdots,n\}$ of cardinality $|{\mathcal{S}}|\leq k$. If so, problem (1) recovers all $k$-sparse vectors ${\mathbf{x}}^{0}$ from measurements ${\mathbf{b}}={\mathbf{A}}{\mathbf{x}}^{0}$. The NSP is also necessary for exact recovery of all $k$-sparse vectors uniformly. The wide use of NSP can be found in, e.g., . Note that it holds regardless the value of $\|{\mathbf{x}}^{0}\|_{\infty}$. We now give a necessary and sufficient condition for problem (5).

Assume $\|{\mathbf{x}}^{0}\|_{\infty}$ is fixed. Problem (5) uniquely recovers all $k$-sparse vectors ${\mathbf{x}}^{0}$ with the fixed $\|{\mathbf{x}}^{0}\|_{\infty}$ from measurements ${\mathbf{b}}={\mathbf{A}}{\mathbf{x}}^{0}$ if and only if

holds for all vectors ${\mathbf{h}}\in\Null({\mathbf{A}})$ and coordinate sets ${\mathcal{S}}$ of cardinality $|{\mathcal{S}}|\leq k$.

where the first inequality follows from the triangle inequality, and the second follows from $\|{\mathbf{h}}_{\mathcal{S}}\|_{2}^{2}+\|{\mathbf{h}}_{\mathcal{Z}}\|_{2}^{2}=\|{\mathbf{h}}\|_{2}^{2}$ and $\langle{\mathbf{x}}^{0}_{\mathcal{S}},{\mathbf{h}}_{\mathcal{S}}\rangle\geq-\|{\mathbf{x}}^{0}_{\mathcal{S}}\|_{\infty}\|{\mathbf{h}}_{\mathcal{S}}\|_{1}=-\|{\mathbf{x}}^{0}\|_{\infty}\|{\mathbf{h}}_{\mathcal{S}}\|_{1}$.

Since $\|{\mathbf{h}}\|_{2}^{2}>0$, $\|{\mathbf{x}}^{0}+{\mathbf{h}}\|_{1}+\frac{1}{2\alpha}\|{\mathbf{x}}^{0}+{\mathbf{h}}\|_{2}$ is strictly larger than $\|{\mathbf{x}}^{0}\|_{1}+\frac{1}{2\alpha}\|{\mathbf{x}}^{0}\|_{2}$ provided that the second block of (19) is nonnegative. Hence, condition (18) is sufficient for ${\mathbf{x}}^{0}$ to be the unique minimizer of (5) .

for all $0<\tau\leq 1$, which in turn requires (18) to hold. ∎

For any finite $\alpha>0$, (18) is stronger than (17) due to the extra term $\frac{\|{\mathbf{x}}^{0}_{\mathcal{S}}\|_{\infty}}{\alpha}$. Since various uniform recovery results establish conditions that guarantee (17), one can tighten these conditions so that they guarantee (18) and thus the uniform recovery by problem (5). How much tighter these conditions have to be depends on the value $\frac{\|{\mathbf{x}}^{0}_{\mathcal{S}}\|_{\infty}}{\alpha}$.

3 Restricted isometry property

In this subsection, we first review the RIP-based sparse recovery guarantees and then show that given certain RIP conditions, any $\alpha\geq 10\|{\mathbf{x}}^{0}\|_{2}$ guarantees exact and stable recovery by (5) and (12), respectively.

The RIP constant $\delta_{k}$ of matrix ${\mathbf{A}}$ is the smallest value such that

For (1) to recover any $k$-sparse vector uniformly, shows the sufficiency of $\delta_{2k}<0.4142$, which is later improved to $\delta_{2k}<0.4531$ , $\delta_{2k}<0.4652$ , $\delta_{2k}<0.4721$ , as well as $\delta_{2k}<0.4931$ . The bound is still being improved. Adapting results in , we give the uniform recovery conditions for (5) below.

For $\delta_{2k}=0.4404$, we obtain $\left(\theta_{2k}^{-1}-1\right)^{-1}\|{\mathbf{x}}^{0}\|_{\infty}\approx 9.9849\|{\mathbf{x}}^{0}\|_{\infty}\leq\alpha$, which proves the theorem. ∎

Different values of $\delta_{2k}$ are associated with different conditions on $\alpha$. Following (23), if $\delta_{2k}\leq 0.4715$, $\alpha\geq 25\|{\mathbf{x}}^{0}\|_{\infty}$ guarantees exact recovery. If $\delta_{2k}\leq 0.1273$, $\alpha\geq\|{\mathbf{x}}^{0}\|_{\infty}$ guarantees exact recovery. In general, a smaller $\delta_{2k}$ allows a smaller $\alpha$.

Next we study the case where ${\mathbf{b}}$ is noisy or ${\mathbf{x}}^{0}$ is not exactly sparse, or both. For comparison, we present two inequalities next to each other for problems (2) and (5) each, where the first one is easy to obtain; see for example.

where $\|{\mathbf{x}}^{0}_{\mathcal{Z}}\|_{1}$ is the best $k$-term approximation error of ${\mathbf{x}}^{0}$ and

We only show (25). Since ${\mathbf{x}}^{*}={\mathbf{x}}^{0}+{\mathbf{h}}$ is the minimizer of (12), we have

where the first inequality follows from the triangle inequality, and the second from $\langle{\mathbf{a}},{\mathbf{b}}\rangle\leq\|{\mathbf{a}}\|_{\infty}\|{\mathbf{b}}\|_{1}$. Combining (27) and (28), we obtain

and thus (25) after dropping the nonnegative term $\frac{1}{2\alpha}\|{\mathbf{h}}\|^{2}$. ∎

We now present the stable recovery guarantee.

Assume the setting of Lemma 1. Let ${\mathbf{b}}:=A{\mathbf{x}}^{0}+{\mathbf{n}}$, where ${\mathbf{n}}$ is an arbitrary noisy vector with $\|{\mathbf{n}}\|_{2}\leq\sigma$. If ${\mathbf{A}}$ satisfies RIP with $\delta_{2k}\leq 0.3814$, then the solution ${\mathbf{x}}^{*}$ of (12) with any $\alpha\geq 10\|{\mathbf{x}}^{0}\|_{\infty}$ satisfies

where $C_{1}$, $C_{2}$, $\bar{C}_{1}$, and $\bar{C}_{2}$ are given in (33a)–(34b) as functions of only $\delta_{2k}$, $C_{3}$, and $C_{4}$ in (26).

We follow an argument similar to that in . According to Lemma 4.3 of , from $\|{\mathbf{A}}{\mathbf{h}}\|_{2}=\|{\mathbf{A}}{\mathbf{x}}^{*}-{\mathbf{A}}{\mathbf{x}}^{0}\|_{2}=\|{\mathbf{A}}{\mathbf{x}}^{*}-{\mathbf{b}}+{\mathbf{n}}\|_{2}\leq\|A{\mathbf{x}}^{*}-{\mathbf{b}}\|_{2}+\|{\mathbf{n}}\|_{2}\leq 2\|{\mathbf{n}}\|_{2}$ and $\delta_{2k}<2/3$, we obtain

where $\theta_{2k}$ is defined in (22) as a function of $\delta_{2k}$. It is easy to verify that with the choice of $\delta_{2k}\leq 0.3814$ and $\alpha$ in the theorem, $C_{3}\theta_{2k}<1$ holds for all nonzero ${\mathbf{x}}^{0}$. Hence, combining (25) of Lemma 1 and (31) yield the bound of $\|{\mathbf{h}}_{\mathcal{Z}}\|_{1}$:

To prove (30), we apply (32) to the inequality (Page 7 of )

A key inequality in the proof above is $C_{3}\theta_{2k}<1$, where $C_{3}$ (cf. (26)) depends on $\alpha$, $\|{\mathbf{x}}^{0}_{\mathcal{S}}\|_{\infty}$, and $\|{\mathbf{x}}^{0}_{\mathcal{Z}}\|_{\infty}$, and $\theta_{2k}$ (cf. (22)) depends on $\delta_{2k}$. If the nonzeros of ${\mathbf{x}}^{0}$ decay faster in magnitude, $C_{3}$ becomes smaller and thus the condition $C_{3}\theta_{2k}<1$ is easier to hold. Therefore, a faster decaying ${\mathbf{x}}^{0}$ is easier to recover. This is consistent with the numerical simulation in subsection 3.1. In Theorem 3, the condition on $\delta_{2k}$ and bound on $\alpha$ are given for the worst case corresponding to no decay, namely, $\|{\mathbf{x}}^{0}_{\mathcal{S}}\|_{\infty}=\|{\mathbf{x}}^{0}_{\mathcal{Z}}\|_{\infty}$. If $\|{\mathbf{x}}^{0}_{\mathcal{S}}\|_{\infty}>\|{\mathbf{x}}^{0}_{\mathcal{Z}}\|_{\infty}$, one can allow a larger $\delta_{2k}$ for each fixed $\alpha$ or, equivalently, a smaller $\alpha$ for each fixed $\delta_{2k}$. For example, if $\|{\mathbf{x}}^{0}_{\mathcal{S}}\|_{\infty}\geq 10\|{\mathbf{x}}^{0}_{\mathcal{Z}}\|_{\infty}$, one only needs $\delta_{2k}\leq 0.4348$ instead of the theorem-assumed condition $\delta_{2k}\leq 0.3814$.

There is also a trade-off between $\delta_{2k}$ and $\alpha$. Under the worst case $\|{\mathbf{x}}^{0}_{\mathcal{S}}\|_{\infty}=\|{\mathbf{x}}^{0}_{\mathcal{Z}}\|_{\infty}$, imposing to $\alpha\geq 25\|{\mathbf{x}}^{0}\|_{\infty}$ leads to the relaxed condition $\delta_{2k}\leq 0.4489$.

4 Spherical section property

Next, we derive exact and stable recovery conditions based on the spherical section property (SSP) of ${\mathbf{A}}$, which has the advantage of invariance to left-multiplying nonsingular matrices to the sensing matrix ${\mathbf{A}}$, as pointed out in . On the other hand, more matrices are known to satisfy the RIP than the SSP.

holds for all nonzero ${\mathbf{h}}\in{\mathcal{V}}$.

with probability at least $1-\exp(C_{1}(n-m))$, where $C_{0}$ and $C_{1}$ are universal constants. Hence, $m>4k\Delta$ guarantees (17) to hold, and furthermore, if $\Null({\mathbf{A}})$ is uniformly random, $m=O(k\log(n/m))$ is sufficient for (17) to hold with overwhelming probability . These results can be extended to the augmented model (5).

Suppose $\Null({\mathbf{A}})$ satisfies the $\Delta$-SSP. Let us fix $\|{\mathbf{x}}^{0}\|_{\infty}$ and $\alpha>0$. If

then the null-space condition (18) holds for all ${\mathbf{h}}\in\Null({\mathbf{A}})$ and coordinate sets ${\mathcal{S}}$ of cardinality $|{\mathcal{S}}|\leq k$. By Theorem 1, (36) guarantees that problem (5) recovers any $k$-sparse ${\mathbf{x}}^{0}$ from measurements ${\mathbf{b}}={\mathbf{A}}{\mathbf{x}}^{0}$.

Let ${\mathcal{S}}$ be a coordinate set with $|{\mathcal{S}}|\leq k$. Condition (18) is equivalent to

Since $\|{\mathbf{h}}_{\mathcal{S}}\|_{1}\leq\sqrt{k}\|{\mathbf{h}}_{\mathcal{S}}\|_{2}\leq\sqrt{k}\|{\mathbf{h}}\|_{2}$, (37) holds provided that

which itself holds, in light of (35), provided that (36) holds. ∎

Now we consider the case ${\mathbf{A}}{\mathbf{x}}^{0}={\mathbf{b}}$ where ${\mathbf{x}}^{0}$ is an approximately sparse vector.

then the solution ${\mathbf{x}}^{*}$ of (5) satisfies

where $\|{\mathbf{x}}_{{\mathcal{Z}}}^{0}\|_{1}$ is the best $k$-term approximation error of ${\mathbf{x}}^{0}$.

Let ${\mathbf{h}}={\mathbf{x}}^{*}-{\mathbf{x}}^{0}\in\Null({\mathbf{A}})$. Let

Adding $\|{\mathbf{h}}_{\mathcal{S}}\|_{1}$ to (25) and plugging in (41) gives us

or $(1-2C_{4}\bar{C}^{-1})\|{\mathbf{h}}\|_{1}\leq(1+C_{3})\|{\mathbf{h}}_{\mathcal{S}}\|_{1}$. If $\bar{C}\leq 2C_{4}$, (42) naturally holds. Otherwise, we have $\bar{C}>2C_{4}$ and

Now, combining $\Delta$-SSP and (39), we obtain

5 “RIPless” analysis

The “RIPless” analysis gives non-uniform recovery guarantees for a wide class of compressive sensing matrices such as those with iid subgaussian entries, orthogonal transform ensembles satisfying an incoherence condition, random Toeplitz/circulant ensembles, as well as certain tight and continuous frame ensembles, at $O(k\log(n))$ measurements. This analysis is especially useful in situations where the RIP, as well as NSP and SSP, is difficult to check or does not hold. In this subsection, we describe how to adapt the “RIPless” analysis to model (5).

where $C_{0}$ is a universal constant and $\mu({\mathbf{A}})$ is the incoherence parameter of ${\mathbf{A}}$ (see for its definition and values for various kinds of compressive sensing matrices).

The proof is mostly the same as that of Theorem 1.1 of except we shall adapt Lemma 3.2 of to Lemma 2 below for our model (5). We describe the proof of the theorem very briefly here. For any matrix ${\mathbf{A}}$ satisfying property (46) in Lemma 2, the golfing scheme can be used to construct a dual vector ${\mathbf{y}}$ such that ${\mathbf{A}}^{*}{\mathbf{y}}$ satisfies property (47) in Lemma 2. The properties (46) and (47) and the construction are exactly the same as in . Then Lemma 2 below lets this ${\mathbf{A}}^{*}{\mathbf{y}}$ guarantee the optimality of ${\mathbf{x}}^{0}$ to (12). ∎

and there exists ${\mathbf{y}}$ such that ${\mathbf{v}}={\mathbf{A}}^{*}{\mathbf{y}}$ satisfies

then ${\mathbf{x}}^{0}$ is the unique solution to (5) with ${\mathbf{b}}={\mathbf{A}}{\mathbf{x}}^{0}$ and $\alpha\geq 8\|{\mathbf{x}}^{0}\|_{2}$.

Since the last term of (48) is strictly positive, ${\mathbf{x}}^{0}$ is the unique solution to (5) provided that

Following the proof of Lemma 3.2 in and from (46) and (47) we obtain

which together with $\alpha\geq 8\|{\mathbf{x}}^{0}\|_{2}$ give

Hence, ${\mathbf{x}}^{0}+{\mathbf{h}}$ gives a strictly worse objective (5) than ${\mathbf{x}}^{0}$, so ${\mathbf{x}}^{0}$ is the unique solution to (5). ∎

6 Matrix Recovery Guarantees

It is fairly easy to extend the results above, except the “RIPless” analysis, to the recovery of low-rank matrices. Throughout this subsection, we let $\sigma_{i}({\mathbf{X}}),~i=1,\cdots,m$ denote the $i$th largest singular value of matrix ${\mathbf{X}}$ of rank $m$ or less, and let $\|{\mathbf{X}}\|_{*}:=\sum_{i=1}^{m}\sigma_{i}({\mathbf{X}})$, $\|{\mathbf{X}}\|_{F}:=\left(\sum_{i=1}^{m}\sigma^{2}_{i}({\mathbf{X}})\right)^{1/2}$, and $\|{\mathbf{X}}\|_{2}=\sigma_{1}({\mathbf{X}})$ denote the nuclear, Frobenius, and spectral norms of ${\mathbf{X}}$, respectively.

The extension is based on the following property of unitarily invariant matrix norms.

Let ${\mathbf{X}}$ and ${\mathbf{Y}}$ be two matrices of the same size. Any unitarily invariant norm $\|\cdot\|_{\phi}$ satisfies

In particular, matrices ${\mathbf{X}}$ and ${\mathbf{Y}}$ obey

By applying (51), shows that any sufficient conditions based on RIP and SSP of ${\mathbf{A}}$ for recovering sparse vectors by model (1) can be translated to sufficient conditions based on similar properties of ${\mathcal{A}}$ for recovering low-rank matrices by model (3). We can establish similar translations from model (12) to model (9) using both inequalities (51) and (52). Hence, we present the low-rank matrix recovery results only with the parts that are different from their vector counterparts.

Paper presents the NSP condition for problem (3): all matrices ${\mathbf{X}}^{0}$ of rank $r$ or less can be exactly recovered by problem (3) from measurements ${\mathbf{b}}={\mathcal{A}}({\mathbf{X}}^{0})$ if and only if all ${\mathbf{H}}\in\Null(\mathcal{A})\backslash\{\mathbf{0}\}$ satisfy

We can extend this result to problem (8) by applying inequalities (51) and (52).

Assume that $\|{\mathbf{X}}^{0}\|_{2}$ is fixed. Problem (8) uniquely recovers all matrices ${\mathbf{X}}^{0}$ (with the specified $\|{\mathbf{X}}^{0}\|_{2}$) of rank $r$ or less from measurements ${\mathbf{b}}={\mathcal{A}}({\mathbf{X}}^{0})$ if and only if

holds for all matrices ${\mathbf{H}}\in\Null({\mathcal{A}})$.

where the second inequality follows from (19) by letting ${\mathbf{h}}=-s({\mathbf{H}})$ and ${\mathcal{S}}=\{1,\ldots,r\}$ and noticing ${\mathbf{h}}_{\mathcal{S}}=\sum_{i=1}^{r}\sigma_{i}({\mathbf{H}})$ and ${\mathbf{h}}_{\mathcal{Z}}=\sum_{i=r+1}^{m}\sigma_{i}({\mathbf{H}})$.

For any nonzero ${\mathbf{H}}\in\Null({\mathcal{A}})$, $\|{\mathbf{H}}\|_{F}>0$. Hence, from (55) and (54), it follows that ${\mathbf{X}}^{0}+{\mathbf{H}}$ leads to a strictly worse objective than ${\mathbf{X}}^{0}$. That is, ${\mathbf{X}}^{0}$ is the unique solution to problem (8).

Necessity: For any nonzero ${\mathbf{H}}\in\Null({\mathcal{A}})$ obeying (54), let ${\mathbf{H}}={\mathbf{U}}\Sigma{\mathbf{V}}^{\top}$ be the SVD of ${\mathbf{H}}$. Construct ${\mathbf{X}}^{0}=-{\mathbf{U}}\Sigma_{r}{\mathbf{V}}^{\top}$, where $\Sigma_{r}$ keeps only the largest $r$ diagonal entries of $\Sigma$ and sets the rest to 0. Scale ${\mathbf{X}}^{0}$ so that it has the specified $\|{\mathbf{X}}^{0}\|_{2}$. We have

for any $t>0$. For ${\mathbf{X}}^{0}$ to be the unique solution to (8) given ${\mathbf{b}}={\mathcal{A}}({\mathbf{X}}^{0})$, we must have

for all $t>0$. Hence, (54) is necessary. ∎

Paper introduces the following RIP for matrix recovery.

holds for all ${{\mathbf{X}}}\in{\mathcal{M}}_{r}$.

To uniformly recover all matrices of rank $r$ or less by solving (3), it is sufficient for ${\mathcal{A}}$ to satisfy $\delta_{5r}<0.1$ , which has been improved to the RIP with $\delta_{4r}<\sqrt{2}-1$ in and to $\delta_{2r}<0.307$, as well as ones involving $\delta_{3r}$, $\delta_{4r}$, and $\delta_{5r}$, in . The algorithm SVP provably achieves exact recovery if $\delta_{2r}<1/3$.

Next, we present a stronger RIP-based condition for the unsmoothed problem (3), and then extend it to the smoothed problem (8) without a proof.

Let ${\mathbf{X}}^{0}$ be a matrix with rank $r$ or less. Problem (3) exactly recovers ${\mathbf{X}}^{0}$ from measurements ${\mathbf{b}}={\mathcal{A}}({\mathbf{X}}^{0})$ if ${\mathcal{A}}$ satisfies the RIP with $\delta_{2r}<0.4931$.

The proof is a straightforward extension to the arguments in using arguments in ; the interested reader can find it in Appendix. Next we present the result for the augmented model (8).

Let ${\mathbf{X}}^{0}$ be a matrix with rank $r$ or less. The augmented model (8) exactly recovers ${\mathbf{X}}^{0}$ from measurements ${\mathbf{b}}={\mathcal{A}}({\mathbf{X}}^{0})$ if ${\mathcal{A}}$ satisfies the RIP with $\delta_{2r}<0.4404$ and in (8) $\alpha\geq 10\|{\mathbf{X}}^{0}\|_{2}$.

The proof of Theorem 8 in Appendix establishes that any ${\mathbf{H}}\in\Null({\mathcal{A}})$ satisfies $\|{\mathbf{H}}_{0}\|_{*}\leq\theta_{2r}\|\sum_{i\geq 1}{\mathbf{H}}_{i}\|_{*}.$ Hence, (54) holds if $\left(1+\frac{\|{\mathbf{X}}^{0}\|_{2}}{\alpha}\right)^{-1}\geq\theta_{2r}.$ The rest of the proof is similar to that of Theorem 2. ∎

Skipping a proof similar to that of Theorem 3, we present the stable recovery result as follows.

where $\hat{\sigma}({\mathbf{X}}^{0}):=\sum_{i=r+1}^{\min\{n_{1},n_{2}\}}\sigma_{i}({\mathbf{X}}^{0})$ is the best rank-$r$ approximation error of ${\mathbf{X}}^{0}$, $C_{1}$, $C_{2}$, $\bar{C}_{1}$, and $\bar{C}_{2}$ are given by formulas (33a)–(34b) in which $\theta_{2k}$ shall be replaced by $\theta_{2r}$ (given in (99)), and

Although there are few discussions on SSP for low-rank matrix recovery in the literature (cf. ), we present two SSP-based results without proofs.

Assume that $\|{\mathbf{X}}^{0}\|_{2}$ and $\alpha>0$ are fixed. If

then the null-space condition (54) holds for all ${\mathbf{H}}\in\Null({\mathcal{A}})$. Hence, (60) is sufficient for problem (8) to recover any matrices ${\mathbf{X}}^{0}$ of rank $r$ or less from measurements ${\mathbf{b}}={\mathcal{A}}({\mathbf{X}}^{0})$.

then the solution ${\mathbf{X}}^{*}$ of (8) satisfies

where $\hat{\sigma}({\mathbf{X}}^{0}):=\sum_{i=r+1}^{\min\{n_{1},n_{2}\}}\sigma_{i}({\mathbf{X}}^{0})$ is the best rank-$r$ approximation error of ${\mathbf{X}}^{0}$.

Global Linear Convergence

Now we turn to study the numerical properties of the linearized Bregman algorithm (LBreg) for the augmented model (5). In this section, we show that LBreg, as well as its two fast variants, achieves global linear convergence with no assumptions on the solution sparsity or aforementioned properties of matrix ${\mathbf{A}}$. First, we review its four equivalent forms of LBreg that have appeared in different papers. We start off with the dual gradient descent iteration : give a step size $h>0$, ${\mathbf{y}}^{(0)}=\mathbf{0}$, and $k$ starting from 0,

where ${\mathbf{x}}^{(0)}={\mathbf{p}}^{(0)}=\mathbf{0}$ and the Bregman “distance” $D_{f}^{{\mathbf{p}}}(\cdot,\cdot)$ is defined as

The last two terms of (63f) replace the term $\frac{h}{2}\|{\mathbf{A}}{\mathbf{x}}-{\mathbf{b}}\|_{2}^{2}$ in the original Bregman iteration. Following , one can obtain (63d)-(63e) from (63f)-(63g) by setting ${\mathbf{v}}^{(k)}={\mathbf{p}}^{(k)}+h{\mathbf{A}}^{\top}({\mathbf{b}}-{\mathbf{A}}{\mathbf{x}}^{(k)})+\frac{{\mathbf{x}}^{(k)}}{\alpha}$.

It is most convenient to work with (63a) due to its simplicity and gradient-descent interpretation. In the rest of this section, we let $f({\mathbf{y}})$ be the objective function of (10) and have ${\nabla}f({\mathbf{y}})=-{\mathbf{b}}+\alpha{\mathbf{A}}\shrink({\mathbf{A}}^{\top}{\mathbf{y}})$.

In this subsection, we prove a few key results that will be used to prove the restricted strongly convex property in the next subsection.

Let $\lambda_{\min}^{++}({\mathbf{S}})$ denote the minimum strictly positive eigenvalue of a nonzero symmetric matrix ${\mathbf{S}}$, assuming its existence. Namely,

where $\{\lambda_{i}({\mathbf{S}})\}$ is the set of eigenvalues of ${\mathbf{S}}$.

Let ${\mathbf{A}}$ be a nonzero $m$-by-$n$ matrix. Let ${\mathbf{D}}\succ\mathbf{0}$ be an $n$-by-$n$ diagonal matrix with strictly positive diagonal entries. We have

Next, we show that a constrained eigenvalue problem, which will appear in our proof of restricted strong convexity, has a strictly positive minimum objective.

so $({\mathbf{A}}{\mathbf{D}}{\mathbf{A}}^{\top}+{\mathbf{b}}_{1}{\mathbf{b}}_{1}^{\top})=({\mathbf{A}}_{1}{\mathbf{D}}_{1}{\mathbf{A}}_{1}^{\top})$. Furthermore, drop the constraints ${\mathbf{b}}_{1}^{\top}({\mathbf{A}}{\mathbf{c}}^{*}+{\mathbf{B}}{\mathbf{d}}^{*})\leq 0$ and $d_{1}\geq 0$, and consider the resulting problem

(67) would have the same objective value as (65) if the active constraint ${\mathbf{b}}_{1}^{\top}({\mathbf{A}}{\mathbf{c}}^{*}+{\mathbf{B}}{\mathbf{d}}^{*})=\mathbf{0}$ was present. As (67) does not have this constraint, we conclude

We apply the same argument to (67) and then inductively to the subsequent problems: let

i.e., no constraint is violated. Then, from $\bar{{\mathbf{d}}}\geq\mathbf{0}$, (72), and (71), it follows

Let $\shrink$ be the shrinkage operator $\shrink(s)=\hbox{sign}(s)\max\{|s|-1,0\}$. Then the following inequality

The first inequality in (75) can be proved by elementary case-by-case analysis. The second one is trivial. ∎

2 Globally Linear Convergence

In this subsection, we show that the LBreg iteration (63a), as a fixed-step size gradient descent iteration for (10), generates a globally linearly convergent sequences $\{{\mathbf{y}}^{k}\}$ and $\{{\mathbf{x}}^{k}\}$.

To do this, we need the following theorem from with our modifications for better clarity. Below, we use the notion

Let $f$ denote the objective function of problem (10), and ${\mathbf{x}}^{*}$ denote the solution of (5), which is unique since it has a strictly convex objective. Define coordinate sets ${\mathcal{S}}_{+},{\mathcal{S}}_{-},{\mathcal{S}}_{0}$ as the sets of positive, negative, and zero components of ${\mathbf{x}}^{*}$, respectively. Corresponding to ${\mathcal{S}}_{+},{\mathcal{S}}_{-},{\mathcal{S}}_{0}$, decompose

Then, the set of solutions of (10) is given by

which is a convex set. Furthermore, ${\nabla}f({\mathbf{y}}^{\prime})=\mathbf{0},~\forall{\mathbf{y}}^{\prime}\in{\mathcal{Y}}^{*}$.

Any ${\mathbf{y}}^{\prime}\in{\mathcal{Y}}^{*}$ must satisfy the strong duality condition, namely, the primal objective equal to the dual objective: $-f({\mathbf{y}}^{\prime})=\|{\mathbf{x}}^{*}\|_{1}+\frac{1}{2\alpha}\|{\mathbf{x}}^{*}\|_{2}^{2}$. From this and ${\mathbf{A}}{\mathbf{x}}^{*}={\mathbf{b}}$, it is easy to derive $\alpha\shrink({\mathbf{A}}^{\top}{\mathbf{y}}^{\prime})={\mathbf{x}}^{*}$ using a case-by-case analysis on the sign of $x^{*}_{i}$. Conversely, since $\nabla f({\mathbf{y}})=-{\mathbf{b}}+{\mathbf{A}}(\alpha\shrink({\mathbf{A}}^{\top}{\mathbf{y}}))$ and ${\mathbf{A}}{\mathbf{x}}^{*}={\mathbf{b}}$, any ${\mathbf{y}}^{\prime}$ obeying $\alpha\shrink({\mathbf{A}}^{\top}{\mathbf{y}}^{\prime})={\mathbf{x}}^{*}$ satisfies $\nabla f({\mathbf{y}}^{\prime})=\mathbf{0}$. Then, ${\mathbf{y}}^{\prime}\in{\mathcal{Y}}^{*}$.

By the definition (76b), ${\mathcal{Y}}^{*}$ is a polyhedron, so it is convex. ∎

In general, the two sets of equality equations in (76b) do not define a unique ${\mathbf{y}}^{*}$, so ${\mathcal{Y}}^{*}$ can include multiple solutions.

A typical tool for obtaining global convergence at a linear rate (or, global geometric convergence) is the strong convexity of the objective function. A function $g$ is strongly convex with a constant $c$ if it satisfies

Strong convexity, however, does not hold for our $f({\mathbf{y}})$ since ${\nabla}f({\mathbf{y}}^{*})=\mathbf{0},~\forall{\mathbf{y}}^{*}\in{\mathcal{Y}}^{*}$, while ${\mathcal{Y}}^{*}$ is not necessarily a singleton. Nevertheless, we establish the “restricted” strong convexity (78) below.

and \lambda_{{\mathbf{A}}}=\min\left\{\lambda^{++}_{\min}({{\mathbf{C}}}{{\mathbf{C}}}^{\top}):{{\mathbf{C}}}~\text{is a nonzero submatrix of}~{\mathbf{A}}~\text{ofmrows}\right\}.

Hence, ${\mathbf{y}}^{\prime}$ satisfies the KKT conditions of (81). Using the expression of ${\mathcal{Y}}^{*}$ in (76b), these conditions are

Part 2. Let ${\mathbf{A}}_{\pm}=[{\mathbf{A}}_{+},{\mathbf{A}}_{-}].$ We first argue that ${\mathbf{A}}_{\pm}$ is a nonzero submatrix of ${\mathbf{A}}$. Since ${\mathbf{A}}$ and ${\mathbf{b}}$ are both nonzero, the solution ${\mathbf{x}}^{*}$ to problem (5) is nonzero. If some column ${\mathbf{a}}_{i}$ of ${\mathbf{A}}$ is a zero vector, then $x_{i}$ is free from the constraints ${\mathbf{A}}{\mathbf{x}}={\mathbf{b}}$ and thus $x^{*}_{i}=0$. Hence, all the columns of ${\mathbf{A}}_{\pm}$ are nonzero vectors.

From ${\nabla}f({\mathbf{y}})=-{\mathbf{b}}+\alpha{\mathbf{A}}\shrink({\mathbf{A}}^{\top}{\mathbf{y}})$ and $\mathbf{0}=\nabla f({\mathbf{y}}^{\prime})=-{\mathbf{b}}+\alpha{\mathbf{A}}\shrink({\mathbf{A}}^{\top}{\mathbf{y}}^{\prime})$, we obtain

By definition, every component of $\shrink({\mathbf{A}}_{\pm}^{\top}{\mathbf{y}}^{\prime})=\alpha^{-1}{\mathbf{x}}^{*}_{\pm}$ is nonzero, and all components of $\shrink({\mathbf{A}}_{0}^{\top}{\mathbf{y}}^{\prime})=\alpha^{-1}{\mathbf{x}}_{0}^{*}$ are zero. For this reason, we deal with (84b) and (84c) separately.

Applying inequality (75) to (84b), we can “remove” the “$\shrink$” operators for it as

Equations (86) mean the followings: (i) the projected point ${\mathbf{y}}^{\prime}$ is actively confined by the boundaries of ${\mathcal{Y}}^{*}$ involving $[{\mathbf{A}}_{1}~{\mathbf{A}}_{2}~{\mathbf{A}}_{3}~{\mathbf{A}}_{4}]$ (c.f., the last term of (76b)); (ii) ${\mathbf{A}}_{5}$ does not contribute to ${\mathbf{y}}-{\mathbf{y}}^{\prime}$; (iii) by applying (83) and (86b)–(86e), we get ${\mathbf{A}}_{1}^{\top}{\mathbf{y}}^{\prime}={\mathbf{1}}$, ${\mathbf{A}}_{3}^{\top}{\mathbf{y}}^{\prime}=-{\mathbf{1}}$ and can thus simplify the components of (84c) involving ${\mathbf{A}}_{1}$ and ${\mathbf{A}}_{3}$ as follows:

Now we “drop” the components of (84c) involving ${\mathbf{A}}_{2}$, ${\mathbf{A}}_{4}$, and ${\mathbf{A}}_{5}$ as follows: from (75), it follows that $\langle{\mathbf{A}}_{i}^{\top}{\mathbf{y}}-{\mathbf{A}}_{i}^{\top}{\mathbf{y}}^{\prime},\shrink({\mathbf{A}}_{i}^{\top}{\mathbf{y}})-\shrink({\mathbf{A}}_{i}^{\top}{\mathbf{y}}^{\prime})\rangle\geq 0$ for $i=2,4,5$. Hence,

However, (89) is still not enough to bound (80) from zero since $\bar{{\mathbf{A}}}\bar{{\mathbf{D}}}\bar{{\mathbf{A}}}^{\top}$ may still be rank deficient.

Part 3. To bound (80), we now include the “dropped” parts of ${\mathbf{A}}$ and apply Lemma 5. From (83), we have ${\mathbf{A}}_{2}^{\top}{\mathbf{y}}^{\prime}={\mathbf{1}}$ and ${\mathbf{A}}_{4}^{\top}{\mathbf{y}}^{\prime}=-{\mathbf{1}}$, and further from (86c) and (86e),

Now for the objective of (80), we apply (89) and then Lemma 5 to obtain

Note that under our convention, an $m$-by-$0$ matrix vanishes. Since matrix $\bar{{\mathbf{A}}}$ contains the nonzero matrix ${\mathbf{A}}_{\pm}$ as a submatrix, $\bar{{\mathbf{A}}}\bar{{\mathbf{D}}}\bar{{\mathbf{A}}}^{\top}+\bar{{\mathbf{C}}}\bar{{\mathbf{C}}}^{\top}$ is nonzero. Therefore, we have

If the entries of ${\mathbf{A}}$ are in general positions, i.e., any $m$ distinct columns of ${\mathbf{A}}$ are linearly independent, or in other words, ${\mathbf{A}}$ has completely full rank , then all $m$-by-$m$ submatrices of ${\mathbf{A}}$ have full rank and thus \lambda_{{\mathbf{A}}}=\{\lambda_{\min}({\mathbf{C}}{\mathbf{C}}^{\top}):{\mathbf{C}}~\text{is anm$-by-$msubmatrix of}~{\mathbf{A}}\}. This is often the case when those entries are samples from i.i.d. subgaussian distributions, or the columns of ${\mathbf{A}}$ are data vector independent of one another. In general, the submatrix ${\mathbf{C}}^{*}$ achieving the minimum $\lambda_{{\mathbf{A}}}$ has the maximum number of independent columns, i.e., it contains $r$ columns from ${\mathbf{A}}$ where $r=\rank({\mathbf{A}})$.

With the restricted strong convexity property, we next show the main convergence result with the help of the standard notion of point–to–set distance

where ${\mathbf{z}}$ is a vector and ${\mathcal{Z}}$ is a set of vectors. By convention, the convergence $\dist({\mathbf{z}}^{k},{\mathcal{Z}})\to 0$ is called globally Q-linear if there exists $\mu\in(0,1)$ such that $\dist({\mathbf{z}}^{k+1},{\mathcal{Z}})/\dist({\mathbf{z}}^{k},{\mathcal{Z}})\leq\mu$ for all $k$, and the convergence $s^{k}\to 0$ is called globally R-linear if there exists a globally Q-linear converging sequence $t^{k}\to 0$ such that $|s^{k}|\leq|t^{k}|$. Unlike Q-linear convergence, R-linear convergence does not require $|s^{k}|$ to be monotonic in $k$.

where the strong convexity constant $\nu$ is given in (79), generates a globally Q-linearly converging sequence $\{{\mathbf{y}}^{(k)},k\geq 1\}$

where ${\mathcal{Y}}^{*}$ is given in (76). The objective value sequence converges R-linearly as

Furthermore, $\{{\mathbf{x}}^{(k)}\}$ is a globally R-linear converging sequence since

where we have used the nonexpansive property of the shrinkage operator (cf. ). Hence, we obtain (91).

To get (92), we recall for any convex $f$ with $L$-Lipschitz $\nabla f$, $f({\mathbf{y}})-f({\mathbf{x}})\leq\langle\nabla f({\mathbf{x}}),{\mathbf{y}}-{\mathbf{x}}\rangle+\frac{L}{2}\|{\mathbf{x}}-{\mathbf{y}}\|_{2}^{2}$ (see Theorem 2.1.5 in ). Let ${\mathbf{y}}={\mathbf{y}}^{(k)}$ and ${\mathbf{x}}={\mathbf{y}}^{\prime(k)}$. We have $f({\mathbf{y}}^{\prime(k)})=f^{*}$, $\nabla f({\mathbf{y}}^{\prime(k)})=\mathbf{0}$ and from (91),

which shows (92). When $0<h<2\nu/(\alpha^{2}\|{\mathbf{A}}\|^{4})$, we have $\left(1-2h\nu+h^{2}\alpha^{2}\|{\mathbf{A}}\|^{4}\right)<1$. Due to (63b), (76a), and the non-expansiveness of $\shrink(\cdot)$, we get

If we set $h=\nu/(\alpha^{2}\|{\mathbf{A}}\|_{2}^{4})$, then the geometric decay factor $\left(1-2h\nu+h^{2}\alpha^{2}\|{\mathbf{A}}\|_{2}^{4}\right)=\left(1-\nu^{2}/(\alpha^{2}\|{\mathbf{A}}\|_{2}^{4})\right)$. Hence, we find the convergence rate affected by $\nu$, $\alpha$, and $\|{\mathbf{A}}\|_{2}$. From the definition of $\nu$ in (79), we get

For recovering a sparse vector, recall that both the simulations in Section 3.1 and the analysis in Section 3 show that if ${\mathbf{x}}^{*}$ has faster decaying nonzero entries, $C$ can be set smaller. So, when $r^{*}$ is large, one can choose a small $C$ to counteract.

The proved rate of convergence is quite conservative. The dependence on the solution dynamic range is due to (85), which considers the worst case of (75), yet when this worst case happens, the inequality between (94d) and (94e) can be improved due to properties of the shrinkage operator. In addition, our analysis on the global rate does not exploit the possibility that the algorithm may reach the optimal active set in a finite number of iterations and then exhibit faster linear convergence, typically at a rate depending only on the active set of columns of ${\mathbf{A}}$ and independent of the solution’s dynamic range.

The step size $h\leq 2\nu/(\alpha^{2}\|{\mathbf{A}}\|_{2}^{4})$ is also very conservative. As one will see in the simulation results in the next section, classical techniques for gradient descents such as line search can significantly accelerate the convergence.

3 Extensions to two faster variants of LBreg

We extend the linear convergence results to two variants of LBreg (iteration (63)) that can run significantly faster than LBreg: BB-line-search and kicking . The former dynamically sets the step size $h$ in (63) by the Barzilai-Borwein method with nonmontone line search using techniques from . The latter is a simple add-on to iteration (63) to consolidate a sequence of consecutive iterations in which ${\mathbf{x}}^{k}$ is unchanged. If ${\mathbf{x}}^{k}=\cdots={\mathbf{x}}^{k+j}$, shows that ${\mathbf{y}}^{k},\ldots,{\mathbf{y}}^{k+j}$ stay on the same line, so it is easy to skip all the intermediate iterations and go directly to the end of the line.

Obviously, since kicking only skips certain LBreg iterations, it remains have global linear convergence. On the other hand, given strong convexity, Theorems 3.1 and 3.2 of shows that BB-line-search also enjoys global linear convergence (though the results are weakened to the R-linear convergence of ${\mathbf{A}}{\mathbf{x}}^{(k)}-{\mathbf{b}}$ in our case); it is not difficult to verify that the proof of the theorem remains to hold given only restricted strong convexity In , Theorem 3.1 relies on its inequality (3.4), which in turn require inequalities (3.3) and (3.2) to hold between a current point and its projection to the solution set. The latter is precise our (77). Theorem 3.2 needs (3.12) and in turn (3.11). (3.11) is obtained from (3.1) restricted to between a current point and its projection to the solution set, which can be proved by assuming (3.2) or our (77)..

4 Numerical Demonstration

We present the results of simple tests to demonstrate the convergence of three algorithms: the original LBreg iteration (63), kicking , and BB-line-search . Their numerical efficiency and properties have been previously studied in papers and are not the focus of this paper, so we merely use two examples to illustrate global linear convergence. We generated two compressive sensing tests where both tests had signals ${\mathbf{x}}^{0}$ with 512 entries, out of which $50$ were nonzero and sampled from the standard Gaussian distribution (for Figure 2) or the Bernoulli distribution (for Figure 3). Both tests had the same sensing matrix ${\mathbf{A}}$ with 256 rows and entries sampled from the standard Gaussian distribution. We set $\alpha=10\|{\mathbf{x}}^{0}\|_{\infty}$ in each test and stopped all the three algorithms upon $\|{\nabla}f({\mathbf{y}})\|_{2}<10^{-6}$. The iterative errors $\|{\mathbf{x}}^{k}-{\mathbf{x}}^{*}\|_{2}$ and $\|{\mathbf{y}}^{k}-{\mathbf{y}}^{*}\|_{2}$ of the three algorithms are depicted in Figures 2 and 3.

In both tests, the original version was the slowest. Besides the obvious speed differences, we observe that $\{{\mathbf{x}}^{(k)}\}$ were not monotonic, there were sets of consecutive iterations in which ${\mathbf{x}}^{(k)}$ did not change or fluctuated. Indeed, it is impossible to improve its R-linear convergence to Q-linear convergence. In addition, unlike the other two algorithms, BB-line-search has non-monotonic $\{{\mathbf{y}}^{(k)}\}$, which converges R-linearly instead of Q-linearly.

The convergence appears to have different stages. The early-middle stage has much slower convergence than the final stage.

Comparing the results of two tests, the convergence was faster on the Bernoulli sparse signal than the Gaussian sparse signal. Since the two tests used the same sensing matrix ${\mathbf{A}}$ and the same sparsity, the main reason should be the dynamic range of the signals. A smaller dynamic range leads to faster convergence, which matches our theoretical result on the convergence rate.

Appendix

We establish the theorem by showing that (53) holds for any ${\mathbf{H}}\in\Null({\mathcal{A}})\setminus\{\mathbf{0}\}$.

Based on the SVD ${\mathbf{H}}=\sum_{i=1}^{m}\sigma_{i}({\mathbf{H}}){\mathbf{u}}_{i}{\mathbf{v}}_{i}^{\top}$, where $\sigma_{i}({\mathbf{H}})$ is the $i$-th largest singular value of ${\mathbf{H}}$, we decompose ${\mathbf{H}}={\mathbf{H}}_{0}+{\mathbf{H}}_{1}+{\mathbf{H}}_{2}+\cdots$ where ${\mathbf{H}}_{0}=\sum_{i=1}^{r}\sigma_{i}({\mathbf{H}}){\mathbf{u}}_{i}{\mathbf{v}}_{i}$, ${\mathbf{H}}_{1}=\sum_{i=r+1}^{2r}\sigma_{i}({\mathbf{H}}){\mathbf{u}}_{i}{\mathbf{v}}_{i}$, ${\mathbf{H}}_{2}=\sum_{i=2r+1}^{3r}\sigma_{i}({\mathbf{H}}){\mathbf{u}}_{i}{\mathbf{v}}_{i}$, …. Following these definitions, condition (53) can be equivalently written as

From ${\mathbf{H}}\not=\mathbf{0}$ and the definition of ${\mathbf{H}}_{0}$, we know that ${\mathbf{H}}_{0}\not=\mathbf{0}$ and thus ${\mathcal{A}}({\mathbf{H}}_{0})\not=\mathbf{0}$ due to the RIP of ${\mathcal{A}}$. From ${\mathcal{A}}({\mathbf{H}})=\mathbf{0}$ and ${\mathcal{A}}({\mathbf{H}}_{0})\not=\mathbf{0}$, it follows that ${\mathcal{A}}(\sum_{i\geq 1}{\mathbf{H}}_{i})\not=\mathbf{0}$ and thus $\sum_{i\geq 1}{\mathbf{H}}_{i}\not=\mathbf{0}$. Therefore, $\sum_{i\geq 1}\|{\mathbf{H}}_{i}\|_{*}>0$, and we can define $t:=\|{\mathbf{H}}_{1}\|_{*}/(\sum_{i\geq 1}\|{\mathbf{H}}_{i}\|_{*})>0$ and $\rho:=\|{\mathbf{H}}_{0}\|_{*}/(\sum_{i\geq 1}\|{\mathbf{H}}_{i}\|_{*})>0$.

Next, we present two inequalities without proofs (the interested reader can verify them following the proofs of Lemmas 2.3 and 2.4 in ):

Since ${\mathcal{A}}({\mathbf{H}}_{0}+{\mathbf{H}}_{1})+{\mathcal{A}}\left(\sum_{i\geq 2}{\mathbf{H}}_{i}\right)={\mathcal{A}}({\mathbf{H}})=\mathbf{0}$, the two right-hand sides of (98) equal each other. Hence,

or after a simple calculation of the maximum of $t\in$,

If $\delta_{2r}<(77-\sqrt{1337})/82\approx 0.4931$, then $\theta_{2r}<1$ and thus $\rho<1$. By definition, we get (97) and (53). ∎

Acknowledgements

We thank Hui Zhang, who was visiting Rice from National U of Defense Technology, for his suggestions on the RIPless analysis, as well as Profs. Shiqian Ma and Qing Ling for valuable discussions. We also thank the anonymous referees for numerous suggestions and corrections that have helped improve this manuscript.

References