Sparse Recovery via Differential Inclusions

Introduction

Such a problem has been widely studied in applied mathematics , engineering, and statistics , see for example surveys in . In these works, convex regularization or relaxation approach has been exploited to overcome the combinatorial explosion of searching the best sparse signals using subset least squares. However, it has been known since that all convex regularization approaches lead to biased estimators whose expectation does not meet the true signal, which motivates the exploration of using nonconvex regularization yet it may suffer from a computational hurdle of locating the global optima .

To address this dilemma between statistical accuracy and computational hurdle, in this paper we introduce some dynamics from the Inverse Scale Space (ISS) method, which first appeared in the image restoration literature in and analyzed and implemented carefully in . The name refers to the observation there that large-scale (image) features are recovered before small-scale ones. Our goal here is to show that such dynamics provides a surprisingly simple way to statistically accurate (unbiased and sign-consistent) estimator if equipped with a new type regularization – early stopping. Our results also extend those early error analysis on ISS to statistical consistency, establishing model selection consistency as well as minimax optimal $l_{2}$ error bounds under comparable conditions to LASSO, etc.

The first one, called Bregman ISS here, is given by the nonlinear differential inclusions:

A damping version of the first one, called Linearized Bregman ISS, has its solution path $\{\rho_{t},\beta_{t}\}_{t\geq 0}$ governed by the nonlinear differential inclusions:

where $\kappa>0$ is a constant. Compared to (1.2a), equation (1.3a) has the additional term $\frac{1}{\kappa}\dot{\beta}$. As $\kappa\to\infty$, (1.3) is reduced to (1.2), and the solution path of (1.3) may converge to that of (1.2) exponentially fast as $\kappa$ increases. We will see that (1.3) has a unique solution path $\rho_{t}$ and $\beta_{t}$, $t\geq 0$, which are both continuous for all $\kappa>0$. Alternatively, (1.3) can be obtained as a differential inclusion replacing the $l_{1}$-norm in (1.2b) by the Elastic Net penalty $\|\beta_{t}\|_{1}+\frac{1}{\kappa}\|\beta_{t}\|_{2}^{2}$ which will be discussed later.

The discretizations of (1.2) and (1.3) are known as Bregman Iteration (equation (3.7) of ) and Linearized Bregman Iteration (equations (5.19-20) of ), respectively. They were introduced in the literature of variational imaging and compressive sensing before (1.2) and (1.3). Through a change of variable, Bregman Iteration becomes the iteration of the Augmented Lagrangian Method . On the other hand, Linearized Bregman Iteration is a simple two-line iteration:

which is evidently a forward Euler discretization to (1.3), where $\alpha_{k}>0$ is a step size. Define $z_{k}=\rho_{k}+\frac{1}{\kappa}\beta_{k}$. Then (1.4) can be simplified to:

which is called LASSO in statistics literature .

To see this, consider the general LASSO problem ,

where for the convenience of comparison we replace the regularization parameter $\lambda$ by $t=1/\lambda$ in the following equivalent form

Aside from the obvious relation $t=1/\lambda$, solution $\beta$ is piece-wise linear in $\lambda$ though not so in $t$. Despite this, $t$ will be convenient to our analysis by reflecting a nature of time evolution of the solution.

Since (1.7) is a convex program, $\hat{\beta}_{t}$ is a solution to (1.7) if and only if it obeys the first-order optimality conditions

which are obtained by taking the subdifferential of the objective in (1.7).

It is well-known that LASSO solution $\hat{\beta}_{t}$ is biased . For example, considering the simple case that $n=p=1$, $X$ is the identity and $y\geq 0$, then (1.8) yields

Moreover, the Linearized Bregman ISS (1.3) has the solution,

which converges to the unbiased Bregman ISS estimator exponentially fast.

In reality we are not given the support set $S$, so the following two properties are used to evaluate the performance of an estimator $\hat{\beta}$.

Asymptotic normality: $\sqrt{n}(\hat{\beta}-\beta^{\ast})\to\mathcal{N}(0,\Sigma^{*})$, where

Since these properties hold for the oracle estimator, they are often referred to as the oracle properties.

The bias can be removed by a simple differentiation of LASSO solution. To see this, by multiplying $t$ on both sides of (1.8a) and differentiating it with respect to $t$, any point on the LASSO path satisfies

With path consistency assumed at time $t=\tau$, we have $\beta_{\tau,i}=0,~\forall~i\not\in S$, and from (1.14) we have

Generically, sign consistency occurs in a neighborhood and thus $\dot{\rho}_{\tau,S}=0$. Therefore,

which is the oracle estimator without bias! This motivates us to replace $(\hat{\beta}_{t}+t\dot{\hat{\beta}}_{t})$ in (1.14) by just $\beta_{t}$, which gives the differential inclusions (1.2a) of Bregman ISS. Later we will show that the resulting $\beta_{t}$ in (1.2) indeed reaches sign-consistency under nearly the same condition as LASSO and hence gives the unbiased oracle estimator.

Compared to LASSO, our dynamic approach has great advantages in algorithmic simplicity and estimate quality. In practice, while LASSO is solved for a sequence of regularization parameters (i.e. regularization path), with or without extra debiasing steps such as subset least squares, a single run of our algorithms gives the entire path or, in the case of (discrete) linearized Bregman, a (discrete) approximate path. In addition, the linearized Bregman algorithms such as the simple iterative scheme in (1.5) are readily parallelizable. Such regularization paths can be unbiased, or, in the case of (discrete) linearized Bregman, have less bias than LASSO paths. Generally speaking, our algorithms return regularization paths of improved quality at just a fraction of cost by LASSO.

In addition, points out that it is impossible to achieve unbiased estimator with convex regularization. To avoid bias in regularized least square problem, it is thus necessary to introduce non-convex penalties (e.g. SCAD etc.) which however suffers from the computational difficulty (typically NP-hardness ) on locating the global optima. In a contrast, the dynamic approach studied in this paper, without optimizing any objective function, will be seen to play the same role as non-convex regularization but using a new regularization – early stopping. Such a debiasing ability naturally inherits from the dynamic solution paths, without suffering the cost of finding global optima in nonconvex optimization.

Therefore, in addition to giving the basic solution properties such as existence, uniqueness, and (dis)continuity, we also attempt to explain the good behaviors of the new solution paths and sequence by establishing their statistical path consistency property. Basically we argue that

Under nearly the same conditions for LASSO that the covariates $x_{i}$ are sufficiently uncorrelated and the signal $\beta_{S}^{*}$ is strong enough, Bregman ISS (1.2) with a proper early stopping rule will return the oracle estimator;

Sign consistency and $l_{2}$-error bounds of minimax rates can be generalized to the Linearized Bregman iteration (1.4) and its limit dynamics (1.3), under similar conditions.

2 Notation and assumptions

Throughout the paper, given two numbers $a$ and $b$, let $a\vee b:=\max(a,b)$.

3 Outline

In the rest of this paper, we establish basic solution properties of Bregman and Linearized Bregman ISS in Section 2. Section 3 and Section 4 describe statistical consistency properties of Bregman ISS and their generalizations to Linearized Bregman ISS/discretization, respectively. Section 5 is dedicated to the ideas of proofs. Section 6 presents some preliminary data-dependent stopping rules and Section 7 collects some comments on related works. Section 8 provides some preliminary numerical results. Conclusions are summarized in Section 9.

Bregman and Linearized Bregman solution paths

It has been pointed out in that the solution to Bregman ISS (1.2) is a piece-wise regularization path given iteratively by the following nonnegative least squares, starting with $k=0$, $t_{0}=0$, and $\rho_{0}=\beta_{0}=0$:

set $t_{k+1}:=\sup\{t>t_{k}:\rho_{t_{k}}+\frac{t-t_{k}}{n}X^{T}(y-X\beta_{t_{k}})\in\partial\|\beta_{t_{k}}\|_{1}\}$; if $t_{k+1}=\infty$, then exit;

set $\rho_{t_{k+1}}:=\rho_{t_{k}}+\frac{t_{k+1}-t_{k}}{n}X^{T}(y-X\beta_{t_{k}})$;

set $S_{k+1}:=\{i:|(\rho_{t_{k+1}})_{i}|=1\}$ and $T_{k+1}=\{1,\ldots,p\}\setminus S_{k+1}$;

Hence the solution path to (1.2) is established by

where $\rho_{t}$ is piece-wise linear and $\beta_{t}$ is piece-wise constant. As the nonnegative least square in (2.1) is a convex quadratic programming with linear constraints, the solution always exists but may not be unique, especially in high dimensional setting $n<p$ which fails the strong convexity in least square. The following theorem presents some general conditions to ensure both the existence and uniqueness of solution path.

Let $\rho_{t}$ be right continuously differentiable and $\beta_{t}$ be right continuous. Then (1.3) has a unique solution.

The proofs of these theorems are collected in Appendix.

Consistency of Bregman ISS Dynamics

In this section necessary and sufficient conditions are established for noisy sparse signal recovery with Bregman ISS (1.2).

Restricted Strong Convexity: there is a $\gamma\in(0,1]$,

Irrepresentable Condition: there is a $\eta\in(0,1)$,

where $X_{S}^{\dagger}:=X_{S}\left(\frac{1}{n}X^{T}_{S}X_{S}\right)^{-1}$.

Condition A1 says that the Hessian matrix of the empirical risk $\frac{1}{2n}\|y-X\beta\|^{2}_{2}$ restricted on the index set $S\times S$ is strictly positive definitive, so the empirical risk is strongly convex when restricted on the support set $S$. Such a condition is necessary in the sense that once it fails, $X_{S}$ will be linearly dependent and no unique representation is possible under the basis $X_{S}$.

so in this sense one cannot represent the irrelevant covariates $X_{T}$ by the relevant ones $X_{S}$ effectively.

Neither A1 nor A2 can be checked when the support set $S$ of signal is not known. Alternatively we can use a more strict but checkable condition proposed in .

It can be shown that once A3 holds, then A1 and A2 simultaneously hold with

since $(1-\mu(s-1))I_{S}\leq X_{S}^{\ast}X_{S}\leq(1+\mu(s-1))I_{S}$, and

We note that condition A3 is shown to be sharp in the noisy case in . With these one can translate all the theoretical results with condition A1 and A2 into condition A3.

With these assumptions, the following stopping time is crucial throughout this section

In applications, we often normalize the measurement matrix $X$ such that $\|X_{j}\|_{n}=1$. So the crucial dependence is $\overline{\tau}\sim\eta\sqrt{n/\log p}/\sigma$, which is equivalent to the optimal choice of LASSO parameters .

We will examine two scenarios for establishing these results: the first is an interesting mean Bregman ISS path, which is another biased path, distinct to LASSO, yet qualitatively equivalent; the second is the unbiased Bregman ISS path itself, which meets the consistency results above under nearly the same conditions as LASSO.

2 Mean Bregman ISS Path versus LASSO Path

As we have seen in Section 1.1 near equation (1.14), Bregman ISS (1.2) can be derived by differentiating LASSO’s KKT conditions. Such a relation can be seen precisely by considering the consistency conditions of LASSO on the following temporal mean path of Bregman ISS:

A connection between Bregman ISS and LASSO lies in the same condition under which their paths from start to time $t$ are supported within the true support $S$. In addition, the Bregman ISS mean path $\bar{\beta}(t)$ is identical to the LASSO path if the Bregman ISS path is incremental with only adding variables, but without dropping. In general, the two paths are distinct.

Let $(\beta_{t},\rho_{t})$ be either the Bregman ISS path (1.2) or the LASSO path (1.8) with $\rho(t)\in\partial\|\beta_{t}\|_{1}$. Assume that for all $t\leq\tau$,

the Bregman ISS path, its mean path, and the LASSO path all have supports in $S$;

the mean Bregman ISS path $\bar{\beta}(1/\lambda)$ is piecewise linear with $\lambda=1/t$;

In particular in noiseless setting, $\epsilon=0$, (3.3) becomes

which is sufficient and necessary to guarantee that both Bregman ISS, LASSO, and OMP recovers the sparse signal in noiseless setting; once it is violated there is some $S$-sparse signal for which these methods fail.

From (3.4a) one gets the Bregman ISS solution

which leads to the following equation by plugging into (3.4b)

Integration on both sides of this equation and setting

which ensures that $\beta_{T}(t)=0$. So is the mean path.

On the other hand, LASSO starts from the KKT condition (1.8) which splits into

Following the same trick above one can see the same condition (3.7) is met for LASSO to ensure $\hat{\beta}_{T}(t)=0$. This finishes the proof of part A.

As to part B, for $t\leq\tau$, the mean path is obtained by integration on (3.4a)

Equation (2.2) implies that $\frac{1}{t}\rho_{t}=\frac{1}{t}\rho_{t_{k}}+\frac{1-t_{k}/t}{t_{k+1}-t_{k}}\rho_{t_{k+1}}$, which is piecewise linear with respect to $\lambda=1/t$.

Despite of the difference to the LASSO path, the mean Bregman ISS path may reach statistical model-selection consistency under the same conditions as LASSO.

Assume that both (A.1) and (A.2) hold. Then the following holds.

(Sign-Consistency) moreover if the signal is strong enough such that $\beta^{*}_{\min}>c_{1}/\bar{\tau}$,

Under the same conditions as LASSO with $\lambda^{*}=1/\bar{\tau}$ , the mean path $\bar{\beta}$ of Bregman ISS reaches sign-consistency. These conditions are sufficient and necessary in the sense that once violated, there exists an instance such that the probability of failure will be larger than $1/2$ due to noise. In this sense, the mean path estimator $\bar{\beta}(\bar{\tau})$ is “statistically equivalent” to the LASSO estimator.

The mean Bregman ISS path geometrically sheds light on why LASSO incurs bias while Bregman ISS can avoid it. The LASSO path, likes the mean Bregman ISS path, involves some kind of averaging that ensures the path continuity but causes bias. The Bregman ISS path is piecewise constant, allows it to be bias-free.

Now we need to answer the following question: what are conditions to ensure the sign consistency of the Bregman ISS path?

3 Consistency of Bregman ISS

The following theorem tells us that under the irrepresentable (incoherence) condition, the Bregman ISS dynamics always evolves in the support of true signals in the early stage; furthermore if the signal is strong enough then the dynamics will pick up all the true variables before selecting any incorrect ones. When such a sign consistency is reached, Bregman ISS returns the oracle estimator which is unbiased.

Assume that both (A.1) and (A.2) hold. Then Bregman ISS (1.2) has paths satisfying:

(Sign-consistency) moreover if the signal is strong enough such that

To have sign consistency, Theorem 3.3 makes a strong signal condition with a lower bound on $\beta^{*}_{\min}$. However even without such a strong signal assumption, the minimax optimal $l_{2}$-error rates can be achieved disregarding sign consistency.

Assume that both (A1) and (A2) hold. There is a $\tau\in[0,\overline{\tau}]$ such that with probability at least $1-\frac{2}{p\sqrt{\pi\log p}}$,

The existence of such $\tau$ does not ensure us to find it easily. However one can use $\bar{\tau}$ at a cost of enlarging the constants by a square root of condition number of $\Sigma_{S}=X^{*}_{S}X_{S}$.

Under the same condition of Theorem 3.4 and assuming an upper eigenvalue bound $X^{*}_{S}X_{S}\leq\gamma_{\max}I_{S}$, then the following holds for all $t\in[\tau,\bar{\tau}]$ with probability at least $1-\frac{2}{p\sqrt{\pi\log p}}$

where $\mathcal{K}(X_{S}^{*}X_{S})=\gamma_{\max}/\gamma$ is the condition number of $X_{S}^{*}X_{S}$.

All the results in this subsection follow from the more general results on Linearized Bregman ISS (1.3) in the next section by taking $\kappa\to\infty$, whose proofs will be summarized in Section 5.

Generalizations to Linearized Bregman ISS and Its Discretization

In this section, we state a general consistency result for Linearized Bregman ISS (1.3) and Linearized Bregman Iterations (1.4) whose proofs will be given in the next section.

The following new stopping time replaces (3.1) throughout this section.

Clearly when $\kappa\to\infty$, it reduces to (3.1).

The following theorem establishes general conditions for statistical consistency of Linearized Bregman ISS (LBISS) (1.3).

Assume (A1), (A2), and $\kappa$ is big enough such that

(No-false-negative for Mean Path) moreover if the signal is strong enough such that $\beta^{*}_{\min}>c_{1}/\bar{\tau}$,

(Sign-consistency for LBISS) moreover if the smallest magnitude $\beta^{*}_{\min}$ is strong enough and $\kappa$ big enough such that

($l_{2}$-bound) for some constant $C$ and $\kappa$ large enough to satisfy

there is a $\tau\in[0,\overline{\tau}]$ such that $\|\beta_{\tau}-\beta^{*}\|_{2}\leq(C+\frac{2\sigma}{\gamma^{1/2}})\sqrt{\frac{s\log p}{n}}$ with probability at least $1-\frac{2}{p\sqrt{\pi\log p}}-\frac{1}{n\sqrt{\pi\log n}}$.

is enough to guarantee the existence of $\kappa$.

is enough to guarantee the existence of $\kappa$.

Taking $\kappa=\infty$, we get the Theorem 3.3 for Bregman ISS.

where $\mathcal{K}(X_{S}^{*}X_{S})$ is the condition number of $X_{S}^{*}X_{S}$.

2 Consistency of Linearized Bregman iterations

The following theorem establishes statistical consistency conditions for Linearized Bregman Iteration (1.4).

Let $t_{n}=\sum_{k=0}^{n-1}\alpha_{k}$. Assume (A1), (A2), and $\kappa$ is big enough such that

and step size $\alpha$ is small such that $\kappa\alpha\|X_{S}X_{S}^{*}\|<2$. Then any solution path of (1.3) satisfies

(Sign-consistency) moreover if the smallest magnitude $\beta^{*}_{\min}$ is strong enough and $\kappa$ is big enough to ensure

($l_{2}$-bound) for some large enough constants $\kappa$ and $C$ such that

with probability at least $1-\frac{2}{p\sqrt{\pi\log p}}-\frac{1}{n\sqrt{\pi\log n}}$, there is a $k^{*}$, $t_{k^{*}}\leq\bar{\tau}$, such that $\|\beta_{k^{*}}-\beta^{*}\|_{2}\leq(C+\frac{2\sigma}{\gamma^{1/2}})\sqrt{\frac{s\log p}{n}}$.

Analysis of ISS Dynamics

The general idea to analyze differential inclusions in (1.2) and (1.3) is to associate these dynamics with some potential or Lyapunov functions, which control a fast convergence of solutions to the oracle estimator. When the solution path $\beta(t)$ evolves in the support set $S$, a suitable choice of potential functions should be expected with exponentially fast decay, which enables us to estimate the stopping time of reaching sign consistency and small $l_{2}$-error.

The difficulty lies in that ISS dynamics are differential inclusions, hence we exploit differential inequalities of such a potential function to derive the bounds.

One would like to study the dynamics of the following differential inclusion

2 Differential inequality with restricted exponential decay of potential

Define the following Oracle Dynamics as if an oracle discloses the true variable set $S$ such that we restrict our attention on a subspace defined by $S$,

Here $X^{\ast}_{S}X_{S}$ is a $s\times s$ symmetric matrix satisfying the strong convexity $X^{\ast}_{S}X_{S}\geq\gamma I_{s}$, which will lead to exponentially fast decay of potential function.

To reach this goal, our key treatment here is a differential inequality associated differential inclusion in Oracle Dynamics which is tight enough to ensure the exponential decay of potential function. This is a Bihari’s type nonlinear differential inequality, which generalizes the linear cases of Grönwall-Bellman inequalities . In our treatment, a piecewise continuous bound is given which leads to the tight rates in this paper.

The potential $\Psi$ of the Oracle Dynamics above satisfies the following differential inequality

where $F^{-1}$ is the right-continuous inverse of the following strictly increasing function

and hence an exponential decay of Bregman distance. As we shall see in the proof, such a fast rate is crucial to ensuring all the strong signals selected before wrong components. Therefore one can achieve the tight stopping rules below for sign-consistency under nearly the same conditions as LASSO.

Such an inequality ensures a decrease of the potential function at a fast enough speed which leads to the following tight estimates on stopping time.

We are concerned with the following stopping time reaching sign-consistency and $l_{2}$-consistency of Oracle Dynamics, respectively. Define

Equipped with the generalized Bihari’s inequality, one can build up the following bounds for stopping time on sign-consistency and $l_{2}$-consistency, respectively.

The following bounds hold for the Oracle Dynamics (5.4)

3 Sign-consistency and l2l_{2}-error bound

Data-dependent Stopping Rules for Bregman ISS

All the previous results enable us to select $\bar{\tau}$ as a stopping time which however depends on unknown parameters $\gamma$, $\eta$, and noise level $\sigma$, hence is not a data-dependent stopping rule. In this section we present two preliminary results with early stopping rules comparable to , which only depend on the noise level $\sigma$ and thus can be estimated from data. We leave it our future work to explore fully adaptive stopping rules.

In the following, define the residue $r(t):=y-X\beta(t)$. The first theorem adopts the stopping rule based on $\|r(t)\|_{2}$ and the second theorem is based on $\|Xr(t)\|_{\infty}$.

Then Bregman ISS with the stopping rule $\|r(t)\|_{2}\leq\sigma\sqrt{n+2\sqrt{n\log{n}}}$ selects the true subset $S$ with probability at least $1-O(1/n)$.

This result is comparable to Theorem 7 in .

The first condition on the minimum of magnitude of signals ensures the model selection consistency of the Bregman ISS path and thus indicates that one can find some $t$ along the path so that the residual term satisfies $\|r(t)\|_{2}\leq\sigma\sqrt{n+2\sqrt{n\log{n}}}$. Once the path achieves sign consistency, the Bregman ISS must stop.

The second condition $\beta^{*}_{\min}\geq\frac{2\sigma}{\sqrt{\gamma}}\left(\sqrt{1+2\sqrt{\frac{\log n}{n}}}+\sqrt{\frac{\log s}{n}}\right)$ guarantees that one can not stop earlier before Bregman ISS achieves a full recovery. Note that as $n\to\infty$, one needs $\beta^{*}_{\min}\geq 2\sigma/\sqrt{\gamma}$ which is a constant.

Then Bregman ISS with the stopping rule $\|X^{T}r(t)\|_{\infty}\leq 2\sigma\sqrt{\max_{i}\|X_{i}\|\log p}$ ($\delta>0$) selects the true subset $S$ with probability at least $1-O(1/p+1/n)$.

This result is comparable to Theorem 8 in , though the lower bound $\beta^{*}_{\min}\geq O(\sigma\sqrt{s\log p/n})$ loses a factor $\sqrt{s}$ here. As $n\to\infty$, the lower bound can be arbitrarily small.

The remaining of this section presents the proofs of the theorems above.

Lemma 3 in or Lemma 5.2 in shows that with probability at least $1-1/n$, $\epsilon$ is essentially $l_{2}$ upper bounded

We have now shown that the Bregman ISS stops once the path acheives sign consistency.

Next we are going to show that the algorithm will not stop whenever there is some $i\in S$ such that $\beta_{i}(t)=0$. By Lemma A.5,

so it suffices to have $\beta^{*}_{\min}\geq\frac{2\sigma(\sqrt{n+2\sqrt{n\log n}}+\sqrt{\log s})}{\sqrt{n\gamma}}$. ∎

Hence, according to Theorem 4.2, the Bregman ISS achieves the sign consistency with high probability. Assume that at time $\tau^{*}$, $\beta(\tau^{*})$ has the same sign as the underlying sparse signal $\beta$. For each $t$,

where $s_{t}=(I-P_{S(t)})X_{S}\beta_{S}$ is the signal part of the residual and $n_{t}=(I-P_{S(t)})\epsilon$ is the noise part of the residual. Then $r_{\tau^{*}}=n_{\tau^{*}}$. Let $b_{\infty}=\sigma\sqrt{2(1+c)\max_{i}\|X_{i}\|\log p}$.

which means the algorithm stops at $\tau^{*}$.

Next we are going to show that the algorithm will not stop whenever there is some $i\in A_{t}\subseteq S$ such that $\beta_{i}(t)=0$. By Lemma A.5,

with probability at least $1-O(p^{-1}+n^{-1})$. ∎

Related work

For general penalized least square problems, has shown that no convex penalty functions can fully achieve the oracle properties and thus one has to resort to non-convex regularization, whose global minimizer is, however, algorithmically difficult to locate. Alternatively, one can apply LASSO for variable selection and then remove the bias in LASSO by solving a subset least squares in the second stage. On the other hand, noticed that Bregman iteration may reduce bias, also known as contrast loss, in the context of Total Variation image denoising. In this paper, we shall see that dynamics (1.2) can automatically remove bias without any non-convexity or second-stage subset least squares. It is a different kind of regularization via early stopping.

Early stopping regularization has been studied widely in linear inverse problems, e.g. , and recently in Boosting, e.g. . In fact, Linearized Bregman iterations can be viewed as an extension of Landweber iteration (also called $L_{2}$-Boost in statistics),

which follows the primal path $\beta_{t}$ as a gradient descent method solving least square problem. To have solution sparsity, Linearized Bregman iterations (1.4) adds the dual path $\rho_{t}$ in favor of sparse solutions. For ISS, notices that early stopping regularization is needed as the Bregman distance between the signal $\beta^{*}$ and the path $\beta_{t}$ will first decrease and then increase after the prediction error $\|X\beta-y\|$ drops below the noise level. A further quantification of such early stopping regularization is given in under a source condition.

Linearized Bregman iteration (1.5) is shown in equivalent to the gradient ascent iteration applied to the Lagrange dual of the problem

Such a combined $l_{1}$ and $l_{2}$ penalty is called Elastic Net in statistics . In particular, $\beta_{k}$ converges to the unique solution of (7.1) at a linear rate (as long as $X\not=0$ and $X\beta=y$ has a solution); see . In addition, for sufficiently large $\kappa$, the solution to (7.1) is a solution to the basis pursuit model , which is (7.1) without $\frac{1}{2\kappa}\|\beta\|_{2}^{2}$. In noisy settings, early stopping regularization is necessary for signal recovery. The introduction of Elastic Net in statistics is due to a limitation of LASSO that can select at most $s=n$ variables from $p\gg n$ candidates, where the additional $l_{2}$-penalty ($\|\beta\|_{2}^{2}$) enables one to select $s>n$ variables which might be highly correlated, at the cost of a biased estimator. This scenario is beyond the scope of this paper with the assumption $s\leq n<p$ and is left to be explored in the future. However we note that although the Inverse Scale Space (1.3) can be equivalently viewed as differential inclusions (with a discretization (1.5)) associated with the Elastic Net penalty, its dynamics does not follow the regularization paths of Elastic Net. The results in this paper basically say that under nearly the same condition as LASSO, Bregman ISS (1.2) with early stopping regularization may recover the signal without bias, while the bias in (1.3) and (1.5) can be controlled to be arbitrarily small by increasing $\kappa$ with the same sign-consistency. Finally, we note that such iterative algorithms can be easily extended to general settings with differentiable convex loss and non-differentiable convex penalty, e.g. Linearized Bregman iteration in matrix completion .

One should not confuse Linearized Bregman iteration (1.5) with iterative soft-thresholding algorithm (ISTA), which has appeared under different names in the literature (for example, see ),

Both have an iterative thresholded dynamics with similar computational costs. However by moving the shrinkage operator to a different place in (1.5), Linearized Bregman iteration generates a sparse solution path, while ISTA treats $\lambda_{k}$ as the regularization parameter and its iterates converge to a LASSO solution with a regularization parameter $\lambda=\lim_{k\to\infty}\lambda_{k}$. Most ISTA-based LASSO solvers simply use a fixed $\lambda_{k}$ through out the iteration. Although the others update $\lambda_{k}$ over the iterations, they do so not aiming to provide a full regularization path but to accelerate convergence; this technique is known as “continuation” or homotopy method .

2 Parallel and distributed computing

It is very easy to implement iteration (1.5) in parallel and distributed manners and apply it to very large-scale datasets. Suppose

where the all-reduce step collects inputs from and then returns the sum to all the $L$ workstations. It is the sum of $L$ $n$-dimensional vectors, so no matter how the all-reduce step is implemented, the communication cost is independent of $p$. It is important to note that the algorithm is not changed at all. In particular, distributing the data into more computing units, i.e., increasing $L$, does not increase the number of iterations. Therefore, the parallel implementation is nearly embarrassingly parallel and truly scalable. In addition, it is also possible to develop implementations for data divided into blocks of rows of $X$ or even smaller subblocks that split both rows and columns. Recently, (1.5) has also been extended in to a decentralized setting where not only data and computation are distributed but communication is restricted to computing units with direct communication links so there is no data fusion center or long distance communication. The scheme fits sensor network or multi-party regression over the internet, where long-distance communication incurs long delays and high costs.

Experiments

In this section we provide some experimental results to illustrate the relations among LASSO, Bregman ISS (ISS) and Linearized Bregman iteration (LB). The LASSO paths in comparison are computed by R-package ‘lars’, while LB paths are computed with our R-package ‘libra’.

Figure 1 is an example of regularization path of three methods. As $\kappa$ goes bigger, the LB path becomes closer to that of ISS. For LB we choose $\kappa\alpha=1/10$ such that the step size of gradient decent is $1/10$, to satisfy the convergence condition. Note that if $\alpha$ is too big, the solution is oscillating.

To compare the performance of three methods quantitatively, we choose the AUC of ROC curve, to measure the goodness of three regularization paths. ROC (receiver-operating-characteristic) curve is plotted by thresholding the regularization parameter $\lambda$ in LASSO, $t$ in ISS, or $k$ in LB at different levels which create different true positive rates (TPR) and false positive rates (FPR):

ROC is a curve from $(0,0)$ to $(1,1)$. AUC (Area Under the Curve) means the area under the ROC curve. Large AUC values indicate that the signals are picked out earlier than noise on regularization paths. Repeating the experiments for 100 times, in Table 1 we report the mean AUC with standard deviations for the three methods at different noise levels. It shows that all the three methods work reasonably well in this example, while Bregman ISS performs slightly better than LASSO. As $\kappa$ becomes bigger, the performance of LB gets closer to that of Bregman ISS. Notice that as noise level $\sigma$ gets larger, all the methods have their performance decay since signal and noise get confused.

Conclusion and Future Directions

In this paper, noisy sparse signal recovery is approached via dynamics, called Bregman ISS, which can be viewed as a dual gradient descent derived from LASSO KKT conditions. A damped version of this dynamics, Linearized Bregman ISS, can be viewed as a dual gradient descent associated with Elastic Nets. A discretization of Linearized Bregman ISS leads to the widely used Linearized Bregman Iteration algorithm. Equipped with an early stopping regularization, Bregman ISS can simultaneously achieve model selection consistency and unbiased estimation, under nearly the same conditions as LASSO whose estimators are biased though. As a discretization of Linearized Bregman ISS paths, model selection consistency and minimax optimal $l_{2}$-error bounds for Linearized Bregman Iteration are also established. Some data-dependent stopping rules are given for Bregman ISS solution paths.

Future directions of our study include fully data-dependent stopping rules and generalization of our results in nonlinear settings.

A Proofs

The existence part follows from , noticing the nonnegative least squares always have solutions.

We show that the uniqueness part. Define $f(\beta):=\frac{1}{2n}\|y-X\beta\|^{2}$. Then, the differential inclusion (1.2) is equivalent to

Let $S^{+}_{t}:=\{i:(\rho_{t})_{i}=1\}$, $S^{-}_{t}:=\{i:(\rho_{t})_{i}=-1\}$, and $S_{t}=S_{t}^{+}\cup S_{t}^{-}$. By (1.2b), in the case of $S_{t}=\emptyset$, we have $\beta_{t}=0$, so $-\nabla f(\beta_{t})=-\nabla f(0)$ is unique. In the case of $S_{t}\not=\emptyset$, we show below that $X\beta_{t}$ and $-\nabla f(\beta_{t})$ are both unique. The uniqueness of $\rho_{t}$ follows from these results and (A.1a).

In fact, (1.2a) and (1.2b) impose the following constraints on $\beta_{t}$:

To see how $\nabla f(\beta_{t})$ is involved, notice that $\left(\nabla f(\beta_{t})\right)_{i}\geq 0$ must hold for $\forall i\in S_{t}^{+}$ since $(\rho_{t})_{i}\in$ is already at its maximal value 1 and $\nabla f(\beta_{t})<0$ is forbidden as it would further increase $(\rho_{t})_{i}$ to an impossible value. The same argument holds for $\left(\nabla f(\beta_{t})\right)_{i}\leq 0$ for $\forall i\in S_{t}^{-}$.

Furthermore, we will have $(\beta_{t})_{i}\cdot\left(\nabla f(\beta_{t})\right)_{i}=0$ for all $i$. To see this, assume $(\beta_{t})_{i}\not=0$. Then by the right continuity assumption, there exists an interval $[t,t+\epsilon)$ in which $\beta_{i}$ remains nonzero with the same sign. By (A.1b), $(\rho_{t})_{i}$ will remain either $+1$ or $-1$ in the same interval, so $\left(\nabla f(\beta_{t})\right)_{i}=0$. On the other hand, assume $\left(\nabla f(\beta_{t})\right)_{i}\not=0$. Then by (A.1a), $\rho_{i}$ will change and thus it cannot stay either $+1$ or $-1$. By the right continuity of $\beta$, it must hold that $(\beta_{t})_{i}=0$. Therefore, we have the addition constraints

Conditions (A.2) and (A.3) are precisely the KKT optimality conditions for

which is identify to (2.1) except (2.1) specifies the time $t_{k+1}$. Let $\beta_{t}$ be the solution to problem (A.4).

In general, if $f$ is strictly convex, then the solution $\beta_{t}$ is unique. In our case, $f$ is not necessarily strictly convex, but $f=g(X\beta)$ for a strictly convex function $g$. Therefore, $X\beta_{t}$ is unique, and thus so is $\nabla f(\beta_{t})=X^{T}\nabla g(X\beta_{t})$. Lastly, $\beta_{t}$ is unique if the columns of $X$ corresponding to nonzero entries of $\beta_{t}$ are linearly independent since $X\beta_{t}$ is unique. ∎

Obviously, $g(x)$ is Lipschitz continuous. Therefore, the Picard-Lindelöf Theorem implies that there exists a unique solution to this ODE, which leads to the solution of (1.3). ∎

We note that the solution of (1.3), though not piece-wise linear or constant, can still be computed in a piece-wise closed form where on each piece, the signs of $\beta_{t}$ remain unchanged. This is left to the reader.

A.2 Proof of Consistency of LBISS

Assume that $X_{S}$ has full column rank.

For all $t\leq\tau$, solution of (1.3) $\beta(t)$ contains no false positive if

where $P_{T}=I-X_{S}^{\dagger}X_{S}^{\ast}$ is the projection operator onto the column space of $X_{T}$.

Mean path $\bar{\beta}(\tau)$ is sign-consistent if

where $\Phi_{S}=X^{\ast}_{S}X_{S}=\frac{1}{n}X^{T}_{S}X_{S}$.

No-false-positivity and the sign-consistency for mean path in Theorem 4.1, directly follow this lemma.

Consider the differential inclusion (1.3)

From (A.5) one gets $-(\beta_{S}-\beta^{*}_{S})=(X^{*}_{S}X_{S})^{-1}(\dot{\rho}_{S}+\dot{\beta}_{S}/\kappa)-(X^{*}_{S}X_{S})^{-1}X_{S}^{\ast}\epsilon$, which leads to the following equation by plugging into (A.6)

The second part is obtained by integration on (A.5)

Suppose $\epsilon\sim\mathcal{N}(0,\sigma^{2}I_{n})$, and $X\in R^{n\times p}$

From the Gaussian tail probability bound,

The first inequality is directly the union bound of index $j$. The second inequality is obtained by the fact

then according to the definition of $\Psi$ and $F$, we have

Combining the following result from right continuous differentiability

and the strong convexity conditions of $X^{*}_{s}X_{s}$, we have

Note $tr(X_{S}(X_{S}^{\ast}X_{S})^{-1}X_{S}^{\ast})=s,(X_{S}^{\ast}X_{S})^{-1}X_{S}^{\ast}\cdot X_{S}(X_{S}^{\ast}X_{S})^{-1}=(X_{S}^{\ast}X_{S})^{-1}\preceq 1/\gamma,$ and $X_{T}^{\ast}P_{T}\cdot P_{T}X_{T}\preceq X_{T}^{\ast}X_{T}$, using Lemma A.2, we have

(no-false-positivity for $\beta(t)$ up to $\tau$) First consider the LB-ISS

using $\left<d\rho_{S}(t)/dt,d\beta_{S}(t)/dt\right>=0$ from the assumption of Bregman ISS paths. On the set $\mathcal{A}^{c}\bigcup\mathcal{B}^{c}$,

Denote this upper bound as $B$. Returning to the original problem, by Lemma A.1, it suffices to have for all $t\leq\tau$,

which leads to that on the set $\mathcal{C}^{c}$, $t\|X^{\ast}_{T}P_{T}\epsilon\|_{\infty}<(1-B/\kappa\eta)\eta$.

(no-false-negativity for the mean path) it suffices to ensure

where $\Phi_{S}=X_{S}^{\ast}X_{S}$. The second part on the right hand side is $\|\frac{1}{\tau}\Phi_{S}^{-1}(\rho_{S}+\beta_{S}/\kappa)\|_{\infty}\leq\frac{1}{\tau}\|\Phi_{S}^{-1}\|_{\infty}(1+B/\kappa)$. The first part is bounded on the set $\mathcal{B}^{c}$.

($l_{2}$-error bound) Lemma 5.2 implies if $C>\frac{8\sigma\left(\max_{j\in T}\|X_{j}\|_{n}\right)}{\eta\gamma}$, when $\kappa$ is big enough, we have

(Sign Consistency for $\beta_{t}$) The condition

which is ensured by $\kappa$ big enough and

A.3 Proof of Consistency of Linearized Bregman Iterations

First of all, we give a discrete version of generalized Bihari’s inequality which is useful for Linearized Bregman iterations (1.4).

where $X^{\ast}_{S}X_{S}\geq\gamma I$. Let the potential (or Lyapunov) function be

Then the following difference inequality holds

Note that for $i\in S$, $(\rho^{(i)}_{k+1}-\rho^{(i)}_{k})\beta^{(i)}_{k+1}=|\beta^{(i)}_{k+1}|-\rho^{(i)}_{k}\beta^{(i)}_{k+1}\geq 0$

Next we present a discrete stopping time bound from the inequality above.

where $X^{\ast}_{S}X_{S}\geq\gamma I$ and $\alpha_{k}\leq\alpha$, for all $k>0$.

Taking $\alpha\to 0$, it recovers the stopping time bounds in continuous case, Lemma 5.2.

For a uniform upper bound on step sizes $\alpha_{t}\leq\alpha$, by the discrete Bihari’s inequality in Lemma A.3

Note that this implies that $\|r_{t}\|:=\|y-X\beta_{t}\|$ is monotonically nonincreasing for all $t\in(0,\bar{\tau})$. The following lemma makes it precise.

For $t\in[0,\bar{\tau}]$, the residue admits an orthogonal decomposition

Acknowledgements

We thank Dr. Ming Yan for helps on fast Matlab codes for computing Bregman ISS paths. The research of Stanley Osher was supported in part by NSF grant 1118971 and ONR grant N000141210838. The research of Yuan Yao was supported in part by National Basic Research Program of China under grant 2012CB825501 and 2015CB856000, as well as NSFC grant 61071157 and 11421110001. The research of Wotao Yin was supported in part by NSF grants DMS-1349855 and DMS-1317602 and ARO MURI grant W911NF-09-1-0383.

Supplementary Material

Supplement A: Matlab Linearized Bregman codes (http://www.math.ucla.edu/˜wotaoyin/software.html).

Supplement B: R Package of Linearized Bregman algorithms (https://cran.r-project.org/web/packages/Libra/index.html).

References