Tradeoffs between Convergence Speed and Reconstruction Accuracy in Inverse Problems

I Introduction

Often the recovery of ${\bf x}$ from ${\bf y}$ is an ill-posed problem. For example, when ${\bf M}$ has fewer rows than columns ($m<n$), rendering (1) an underdetermined linear system of equations. In this case, it is impossible to recover ${\bf x}$ without introducing additional assumptions on its structure. A popular strategy is to assume that ${\bf x}$ resides in a low dimensional set $\mathcal{K}$, e.g., sparse vectors or a Gaussian Mixture Model (GMM) . The natural by-product minimization problem then becomes

This can be reformulated in an unconstrained form as

A popular technique for solving $\eqref{eq:min_const}$ and $\eqref{eq:min_unconst}$ is using iterative programs such as proximal methods that include the iterative shrinkage-thresholding algorithm (ISTA) and the alternating direction method of multipliers (ADMM) . This strategy is particularly useful for large dimensions $d$.

Many applications impose time constraints, which limit the number of computations that can be performed to recover ${\bf x}$ from the measurements. One way to minimize time and computations is to reduce the number of iterations without increasing the computational cost of each iteration. A different approach is to use momentum methods or random projections to accelerate convergence. Another alternative is to keep the number of iterations fixed while reducing the cost of each iteration. For example, since the complexity of iterative methods rely, among other things, on $m$, a common technique to save computations is to sub-sample the measurements ${\bf y}$, removing “redundant information,” to an amount that still allows reconstruction of ${\bf x}$. A series of recent works suggest that by obtaining more measurements one can benefit from simple efficient methods that cannot be applied with a smaller number of measurements.

In the generalization properties of large-scale learning systems have been studied showing a tradeoff between the number of measurements and the target approximation. The work in showed how it is possible to make the run-time of SVM optimization decrease as the size of the training data increases. In , it is shown that the problem of supervised learning of halfspaces over 3-sparse vectors with trinary values $\{-1,1,0\}$ may be solved with efficient algorithms only if the number of training examples exceeds a certain limit. Similar phenomena are encountered in the context of sparse recovery, where efficient algorithms are guaranteed to reconstruct the sparsest vector only if the number of samples is larger than a certain quantity . In it was shown that by having a larger number of training examples it is possible to design more efficient optimization problems by projecting onto simpler sets. This idea is further studied in by changing the amount of smoothing applied in convex optimization. In the authors show that more measurements may allow increasing the step-size in the projected gradient algorithm (PGD) and thus accelerating its convergence.

While these works studied a tradeoff between convergence speed and the number of available measurements, this paper takes a different route. Consider the case in which due to time constraints we need to stop the iterations before we achieve the desired reconstruction accuracy. For the original algorithm, this can result in the recovery being very far from the optimum. An important question is whether we can modify the original iterations (e.g., those dictated by the shrinkage or ADMM techniques), such that the method convergences to an improved solution with fewer iterations without adding complexity to them. This introduces a tradeoff between the recovery error we are willing to tolerate and the computational cost. As we demonstrate, this goes beyond the trivial relationship between the approximation error and the number of iterations that exists for various iterative methods .

Such a tradeoff is experimentally demonstrated by the success of learned ISTA (LISTA) for sparse recovery with $f(\cdot)=\left\|\cdot\right\|_{1}$. This technique learns a neural network with only several layers, where each layer is a modified version of the ISTA iteration.ISTA and its variants is one of the most powerful optimization techniques for sparse coding. It achieves virtually the same accuracy as the original ISTA using one to two orders of magnitude less iterations. The acceleration of iterative algorithms with neural networks is not unique only to the sparse recovery problem and $f(\cdot)=\left\|\cdot\right\|_{1}$. This behavior was demonstrated for other models such as the analysis cosparse and low-rank matrix models , Poisson noise , acceleration of Eulerian fluid simulation , and feature learning . However, a proper theoretical justification to this phenomena is still lacking.

Contribution. In this work, we provide theoretical foundations elucidating the tradeoff between the allowed minimization error and the number of simple iterations used for solving inverse problems. We formally show that if we allow a certain reconstruction error in the solution, then it is possible to change iterative methods by modifying the linear operations applied in them such that each iteration has the same complexity as before but the number of steps required to attain a certain error is reduced.

Such a tradeoff seems natural when working with real data, where both the data and the assumed models are noisy or approximate; searching for the exact solution of an optimization problem, where all the variables are affected by measurement or model noise may be an unnecessary use of valuable computational resources. We formally prove this relation for iterative projection algorithms. Interestingly, a related tradeoff exists also in the context of sampling theory, where by allowing some error in the reconstruction we may use fewer samples and/or quantization levels . We argue that the tradeoff we analyze may explain the smaller number of iterations required in LISTA compared to ISTA.

Parallel efforts to our work also provide justification for the success of LISTA. In , the fast convergence of LISTA is justified by connecting between the convergence speed and the factorization of the Gram matrix of ${\bf M}$. In , the convergence speed of ISTA and LISTA is analyzed using the restricted isometry property (RIP) , showing that LISTA may reduce the RIP, which leads to faster convergence. A relation between LISTA and approximate message passing (AMP) strategies is drawn in .

Our paper differs from previous contributions in three main points: (i) it goes beyond the case of standard LISTA with sparse signals and considers variants that apply to general low-dimensional models; (ii) our theory relies on the concept of inexact projections and their relation to the tradeoff between convergence-speed and recovery accuracy, which differs significantly from other attempts to explain the success of LISTA; and (iii) besides exploring LISTA, we provide acceleration strategies to other programs such as model-based compressed sensing and sparse recovery with side-information.

Organization. This paper is organized as follows. In Section II we present preliminary notation and definitions, and describe the ISTA, LISTA and PGD techniques. Section III introduces a new theory for PGD for non-convex cones. Section IV shows how it is possible to tradeoff between convergence speed and reconstruction accuracy by introducing the inexact projected gradient descent (IPGD) method using spectral compressed sensing as a motivating example. The reconstruction error of IPGD is analyzed as a function of the iterations in Section V. Section VI discusses the relation between our theory and model-based compressed sensing and sparse recovery with side information . Section VII proposes a LISTA version of IPGD, the learned IPGD (LIPGD), and demonstrates its usage in the task of image super-resolution. Section VIII relates the approximation of minimization problems studied here with neural networks and deep learning, providing a theoretical foundation for the success of LISTA and suggesting a “mixture-model” extension of this technique. Section IX concludes the paper.

II Preliminaries and Background

A popular iterative technique for minimizing (3) is ISTA. Each of its iterations is composed of a gradient step with step size $\mu$, obeying $\frac{1}{\mu}\geq\left\|{\bf M}\right\|$ to ensure convergence , followed by a proximal mapping $\mathcal{S}_{f,\lambda}(\cdot)$ of the function $f$, defined as

where $\lambda$ is a parameter of the mapping. The resulting ISTA iteration can be written as

where ${\bf z}_{t}$ is an estimate of ${\bf x}$ at iteration $t$. Note that the step size $\mu$ multiplies the parameter of the proximal mapping.

The proximal mapping has a simple form for many functions $f$. For example, when $f(\cdot)=\left\|\cdot\right\|_{1}$, it is an element-wise shrinkage function,

Therefore, the advantage of ISTA is that its iterations require only the application of matrix multiplications and then a simple non-linear function. Nonetheless, the main drawback of ISTA is the large number of iterations that is typically required for convergence.

Many acceleration techniques have been proposed to speed up convergence of ISTA (see as a partial list of such works). A prominent strategy is LISTA, which has the same structure as ISTA but with different linear operations in (5). Empirically, it is observed that it is able to attain a solution very close to that of ISTA with a significantly smaller fixed number of iterations $T$. The LISTA iterations are given byWe present the more general version that can be used for any signal model and not only for sparsity.

While other acceleration techniques for ISTA have been proposed together with a thorough theoretical analysis, the powerful LISTA method has been introduced without mathematical justification for its success. In this work, we focus on the PGD algorithm, whose iterations are almost identical to the ones of ISTA but with an orthogonal projection instead of a proximal mapping. We propose an acceleration technique for it, which is very similar to the one of LISTA, accompanied by a theoretical analysis.

II-B Projected gradient descent (PGD)

PGD is a generalization of the iterative hard thresholding (IHT) algorithm, which was developed for $\mathcal{K}$ being the set of sparse vectors . This important method has been analyzed in various works. For example, for standard sparsity in , for sparsity patterns that belong to a certain model in , for a general union of subspaces in , for nonlinear measurements in , and more recently in for a set of the form

The formulation (8) generalizes the special cases above. For example, if $f(\cdot)=\left\|\cdot\right\|_{0}$ and $R$ is the sparsity level then we have the IHT method from ; when $f$ counts the number of non-zeros of only certain sparsity patterns, which are bounded by $R$, we have the model-based IHT of . PGD may also be applied to non-linear inverse problems .

Theorem II.5 below provides convergence guarantees on PGD (it is the noiseless version of Theorem 1.2 in ). Before presenting the result, we introduce several properties of the set $\mathcal{K}$ and some basic lemmas.

The descent set of the function $f$ at a point ${\bf x}$ is defined as

The tangent cone $C_{f}({\bf x})$ at a point ${\bf x}$ is the conic hull of $D_{f}({\bf x})$, i.e., the smallest closed cone $C_{f}({\bf x})$ satisfying $D_{f}({\bf x})\subseteq C_{f}({\bf x})$.

For concise writing, below we denote $D_{f}({\bf x})$ and $C_{f}({\bf x})$ as $D$ and $C$, respectively.

where $\kappa_{f}=1$ if $D$ is a convex set and $\kappa_{f}=2$ otherwise.

We now introduce the convergence rate provided in for PGD. For brevity, we present only its noiseless version.

where $\kappa_{f}$ is defined in Lemma II.4, and

II-C Gaussian mean width

When ${\bf M}$ is a random matrix with i.i.d. Gaussian distributed entries $\mathcal{N}(0,1)$, it has been shown in that the convergence rate $\rho(C)$ is tightly related to the dimensionality of the set (model) ${\bf x}$ resides in. A very useful expression for measuring the “intrinsic dimensionality” of sets is the (Gaussian) mean width.

The Gaussian mean width of a set $\Upsilon$ is defined as

Two variants of this measure are generally used. The cone Gaussian mean width, $\omega_{C}=\omega(C)$, which measures the dimensionality of the tangent cone $C=C_{f}({\bf x})$; and the set Gaussian mean width, $\omega_{\mathcal{K}}=\omega(\mathcal{K}-\mathcal{K})$, which is related directly to the set $\mathcal{K}$ through its Minkowski difference $\mathcal{K}-\mathcal{K}=\left\{{\bf z}-{\bf v}:{\bf z},{\bf v}\in\mathcal{K}\right\}$. The cone Gaussian mean width relies on both the set $\mathcal{K}$ (through $f$) and a specific target point ${\bf x}$, while the set Gaussian mean width considers only $\mathcal{K}$. On the other hand, the dependence of $\omega_{C}$ on $\mathcal{K}$ is indirect via the descent set at the point ${\bf x}$. There is a series of works, which developed convergence and reconstruction guarantees for various methods based on $\omega_{C}$ , and others that rely on $\omega_{\mathcal{K}}$ . The first ($\omega_{C}$) is mainly employed in the case of convex functions $f$, which are used to relax the non-convex set in which ${\bf x}$ resides. In this setting, often $D$ is convex and ${\bf x}\in\mathcal{K}$.

II-D PGD convergence rate and the cone Gaussian mean width

In , it has been shown that the smaller $\omega_{C}$, the faster the convergence. More specifically, if $m$ is very close to $\omega_{C}$, then we may apply PGD with a step-size $\mu=\frac{1}{(\sqrt{d}+\sqrt{m})^{2}}\simeq\frac{1}{d}$ and have a convergence rate of (Theorem 2.4 in )

If $\omega_{C}$ is smaller than $\sqrt{m}$ by a certain constant factor, then we may apply PGD with a larger step size $\mu\simeq\frac{1}{m}$, which leads to improved convergence (Theorem 2.2 in )

These relationships rely on the fact that with larger $m$ the eigenvalues of ${\bf I}-\mu{\bf M}^{T}{\bf M}$ (after projection onto $C$, see (14)) are better positioned such that it is possible to improve convergence by increasing $\mu$.

The connections in (17) and (18) between $\rho(C)$ and $\omega_{C}$ are not unique only to the case that ${\bf M}$ is a random Gaussian matrix. Similar relationships hold for many other types of matrices .

III PGD Theory based on the Projection Set

A similar phenomenon also occurs with the set of sparse vectors with a tree structure $\hat{\mathcal{K}}$ (see (16)), where again $\omega_{C}=\sqrt{d}$ implying $m=O(d)$. Yet, from the work in , we know that in this setting it is sufficient to choose $m=O(k)$. Note that for $\hat{\mathcal{K}}$, the set Gaussian mean width is $\omega_{\hat{\mathcal{K}}}=\sqrt{k}$. If we would have relied on it instead of on $\omega_{C}$ in the bound for the required size of $m$, it would have coincided with .

In order to address these deficiencies in the convergence rate, we provide a variant of Theorem II.5 that relies on the set $\mathcal{K}$ directly through the Minkowski difference $\mathcal{K}-\mathcal{K}$ in lieu of $C$. For simplicity we present only the noiseless case but the extension to the noisy setting can be performed using the strategy in .

where $\kappa_{\mathcal{K}}=1$ if $\mathcal{K}$ is convex and $\kappa_{\mathcal{K}}=2$ otherwise, and

Proof: We repeat similar steps to the ones in the proof of Theorem 1.2 in .

We start by noting that the PGD error at iteration $t+1$ is,

where the last inequality is due to Lemma II.3 and the fact that ${\bf y}={\bf M}{\bf x}$. Since $\mathcal{K}$ is a closed cone, also the Minkowski difference $\mathcal{K}-\mathcal{K}$ is a closed cone. Moreover, $D\subset\mathcal{K}-\mathcal{K}$ as ${\bf x}\in\mathcal{K}$. Thus, following Lemma II.4 we have

where the second inequality is due to the fact that ${\bf z}_{t}$ is of the form $\mathcal{P}_{\mathcal{K}}({\bf v})$ for some vector ${\bf v}$, and the last inequality follows from Lemma II.3.

where the last inequality follows from Lemma II.2. Using the definition of $\rho(\mathcal{K})$ and applying the inequality in (23) recursively leads to the desired result. $\Box$

When ${\bf M}$ is a random Gaussian matrix, the relationships in (17) and (18) hold with $\rho(\mathcal{K})$ and $\omega_{\mathcal{K}}$ replacing $\rho(C)$ and $\omega_{C}$ respectively. This implies that we need $m=O(\omega_{\mathcal{K}}^{2})$ for convergence. This result is in line with the conditions on $m$ that appear in previous works for $k$-sparse vectors , for which $\omega_{\mathcal{K}}=O(\sqrt{k\log(d/k)})$, and for sparse vectors with tree structure , where $\omega_{\mathcal{K}}=O(\sqrt{k})$.

As discussed in Section II-C, the measure $\omega_{\mathcal{K}}$ is related directly to the set $\mathcal{K}$ (may be non-convex) in which ${\bf x}$ resides. Thus, it provides a better measure for the complexity of $\mathcal{K}$ when it is unbounded or has some specific structure as is the case for sparsity with tree structure . In such settings, Theorem III.1 should be favored over Theorem II.5.

Notice that if ${\bf x}\in\mathcal{K}$, then we have $D\subset\mathcal{K}-\mathcal{K}$. Thus, in the settings that $D$ is convex and ${\bf x}\in\mathcal{K}$, we have $C=D\subset\mathcal{K}-\mathcal{K}$ implying that $\omega_{C}\leq\omega_{\mathcal{K}}$; when ${\bf M}$ is random Gaussian, this also implies $\rho(C)\leq\rho(\mathcal{K})$. Therefore, in this scenario Theorem II.5 has an advantage over Theorem III.1.

IV Inexact projected gradient descent (IPGD)

It may happen that the function $f$ or the set $\mathcal{K}$ are too loose for describing ${\bf x}$. Instead, we may select a set $\hat{\mathcal{K}}$ that better characterizes ${\bf x}$ and therefore leads to a smaller $\omega$, resulting in faster convergence. This improvement can be very significant; smaller $\omega$ both improves the convergence rate and allows using a larger step-size (see Section II-D).

A related study showed that it is enough to use a small number of Gaussians to represent all the patches in natural images instead of using a dictionary that spans a much larger union of subspaces. This work relied on Gaussian Mixture Models (GMM), whose mean width scales proportionally to the number of Gaussians used, which is significantly smaller than the mean width of the sparse model.

A difficulty often encountered is that the projection onto $\hat{\mathcal{K}}$, which may even be unknown, is more complex to implement than the projection onto $\mathcal{K}$. The latter can be easier to project onto but provides a lower convergence rate.

Thus, in this work we introduce a technique that compromises between the reconstruction error and convergence speed by using PGD with an inexact “projection” that projects onto a set that is approximately as small as $\hat{\mathcal{K}}$ but yet is as computationally efficient as the projection onto $\mathcal{K}$. In this way, the computational complexity of each projected gradient descent iteration remains the same while the convergence rate becomes closer to that of the more complex PGD with a projection onto $\hat{\mathcal{K}}$.

The “projection” we propose is composed of a simple operator $p$ (e.g., a linear or an element-wise function) and the projection onto $\mathcal{K}$, $\mathcal{P}_{\mathcal{K}}$, such that it introduces only a slight distortion into ${\bf x}$. In particular, we require the following:

If $\mathcal{K}$ is convex, then we require

From the fact that $\left\|\mathcal{P}_{D}\left({\bf x}-p({\bf x})\right)\right\|\leq\left\|{\bf x}-p({\bf x})\right\|$, it is sufficient that

to ensure (24). Examples for projections that satisfy condition (24) are given hereafter in sections IV-B and VI-B.

IV-A2 The projection condition for non-convex 𝒦𝒦\mathcal{K}

In the case that $\mathcal{K}$ is non-convex, we require

Due to Lemma II.3 and a simple change of variables, (27) is equivalent to

which by another simple change of variables is the same as

An example for a projection that satisfies condition (27) is provided in Section VI-A.

IV-B Inexact PGD

Plugging the inexact projection into the PGD step results in the proposed inexact PGD (IPGD) iteration (compare to (8))

To motivate this algorithm consider the problem of spectral compressed sensing , in which one wants to recover a sparse representation in a dictionary that has high local coherence. It has been shown that if the non-zeros in the representation are far from each other then it is easier to obtain good recovery .

V IPGD Convergence Analysis

We turn to analyze the performance of IPGD. For simplicity of the discussion, we analyze the convergence of this technique only for a linear operator $p$ and the noiseless setting, i.e., ${\bf e}=0$. The extension to other types of operators and the noisy case is straightforward by arguments similar to those used in for treating the noise term and other classes of matrices.

We present two theorems on the convergence of IPGD. The first result provides a bound in terms of $\rho(C)$ (i.e., depends on $\omega_{C}$ if ${\bf M}$ is a random Gaussian matrix) for the case that $D$ is convex corresponding to $\kappa_{f}=1$ in Theorem II.5; the second provides a bound in terms of $\rho(K)$ (i.e., depends on $\omega_{\mathcal{K}}$ if ${\bf M}$ is a random Gaussian matrix) when $\mathcal{K}$ is a closed cone but not necessarily convex. The proofs of both theorems are deferred to appendices A and B.

is the “effective convergence rate” of IPGD for small $\epsilon$.

where $\kappa_{\mathcal{K}}$ and $\rho(K)$ are defined in Theorem III.1,

is the “effective convergence rate” of IPGD for small $\epsilon$.

Theorems V.1 and V.2 imply that if $\epsilon$ is small enough (compared to $\rho_{p}^{t}$, where $t$ is the iteration number and $\rho_{p}$ is defined in (V.2)) then IPGD has an effective convergence rate of $\rho_{p}=\rho_{p}(C)$ when $D$ is convex, and $\rho_{p}=\kappa_{\mathcal{K}}\rho_{p}(\mathcal{K})$ in the case that $\mathcal{K}$ is a closed cone but not necessarily convex. Note that if $p={\bf I}$ then $\epsilon=0$ and our results coincide with theorems II.5 and III.1.

As we shall see hereafter, for some operators $p$ the rate $\rho_{p}$ may be significantly smaller than $\rho(C)$ and $\rho(\mathcal{K})$. The smaller the set that $p$ maps to, the smaller $\rho_{p}$ becomes. At the same time, when $p$ maps to smaller sets it usually provides a “coarser estimate” and thus the approximation error $\epsilon$ in (25) and (27) increases. Thus, IPGD allows us to tradeoff approximation error $\epsilon$ and improved convergence $\rho_{p}$.

The error term in theorems V.1 and V.2 at iteration $t$ is comprised of two components. The first goes to zero as $t$ increases while the second increases with iterations and is on the order of $\epsilon$. The fewer iterations we perform the larger $\epsilon$ we may allow. An alternative perspective is that the larger the reconstruction error we can tolerate, the larger $\epsilon$ may be and thus we require fewer iterations. Therefore, the projection $p$ introduces a tradeoff. On the one hand, it leads to an increase in the reconstruction error. On the other hand, it simplifies the projected set, which leads to faster convergence (to a solution with larger error).

The works in use a similar concept of near-optimal projection (compared to that assumes only exact projections). The main difference between these contributions and ours is that these papers focus on specific models, while we present a general framework that is not specific to a certain low-dimensional prior. In addition, in these papers the projection is performed to make it possible to recover a vector from a certain low-dimensional set, while in this work the main purpose of our inexact projections is to accelerate the convergence within a limited number of iterations. For a larger number of iterations these projections may not lead to a good reconstruction error.

VI Examples

This section presents examples of IPGD with an operator $p$ that accelerates the convergence of PGD for a given set $\mathcal{K}$.

The best way to recover ${\bf x}$ is by using a projection onto the set $\hat{\mathcal{K}}$ in (16), which is the strategy proposed in the context of model-based compressed sensing . Yet, this projection requires some additional computations at each iteration . Our technique suggests to approximate it by a linear projection onto the first levels of the tree (a simple operation) followed by a projection onto $\mathcal{K}=\{{\bf z}:\left\|{\bf z}\right\|_{0}\leq k\}$.

The more levels we add in the projection $p$, the smaller the approximation error $\epsilon$ turns out to be. More specifically, it is easy to show that $\epsilon$ in (27) is bounded by two times the energy of the entries eliminated from ${\bf x}$ divided by the total energy of ${\bf x}$, i.e., by $2\frac{\left\|p({\bf x})-{\bf x}\right\|}{\left\|{\bf x}\right\|}$. Clearly, the more layers we add the smaller $\epsilon$ becomes. Yet, assuming that all nodes in each layer are selected with equal probability, the probability of selecting a node at layer $l$ is equal to $\prod_{i=1}^{l}0.5^{i-1}$, where we take into account the fact that a node can be selected only if all its forefathers have been chosen. Thus, the upper layers have more significant impact on the values of $\epsilon$.

On the other hand, the convergence rate $\rho_{p}(\mathcal{K})$ for a projection with $l$ layers is equivalent to the convergence rate for the set of vectors of size $2^{l}$ (denoted by $\mathcal{K}_{l}$). Thus, we get that $\rho_{p}(\mathcal{K})=\rho(\mathcal{K}_{l})$, which is dependent on the Gaussian mean width $\omega_{\mathcal{K}_{l}}$ that scales as $\max(kl,k\log(2^{l}/k))$. Clearly, when we take all the layers $l=\log(d)$ and we have $\omega_{\mathcal{K}_{l}}=\omega_{\mathcal{K}}=O(k\log(d/k))$.

Figure 3(a) presents the signal reconstruction error ($\left\|{\bf x}-{\bf z}_{t}\right\|_{2}$) as a function of the number of iterations for PGD with the sets $\mathcal{K}$ (IHT ) and $\hat{\mathcal{K}}$ (model-based IHT )For demonstration purposes we plot only the cases where model-based IHT converges to zero. and for the proposed IPGD with $p$ that projects onto a different number of levels (1-5) of the tree. All algorithms use step size $\mu=\frac{1}{(\sqrt{d}+\sqrt{m})^{2}}$. It is interesting to note that if $p$ projects only onto the first layer, then the algorithm does not converge as the resulting approximation error $\epsilon$ is too large. However, starting from the second layer, we get a faster convergence at the first iterations with $p$ that projects onto a smaller set, which yields a smaller $\rho$. As the number of iterations increases, the more accurate projections achieve a lower reconstruction error, where the plateau attained is proportional to the approximation error of $p$ as predicted by our theory.

This tradeoff can be used to further accelerate the convergence by changing the projection in IPGD over the iterations. Thus, in the first iterations we enjoy the fast convergence of the coarser projections and in the later ones we use more accurate projections that allow achieving a lower plateau. The last line in Fig. 3 demonstrates this strategy, where at the first iteration $p$ is set to be a projection onto the first two levels, and then every four iterations another tree level is added to the projection until it becomes a projection onto all the tree levels (in this case IPGD coincides with PGD). Note that IPGD converges faster than PGD also when the projection in it becomes onto all the tree levels. This can be explained by the fact that typically convergence of non-linear optimization techniques depends on the initialization point .

While here we arbitrarily chose to add another level every fixed number of iterations, in general, a control set can be used for setting the number of iterations to be performed in each training level. We demonstrate this strategy in Section VII.

Since PGD with $\hat{\mathcal{K}}$ does not introduce an error in its projection and projects onto a precise set, it achieves the smallest recovery error throughout all iterations. Yet, as its projection is computationally demanding, it converges slower than IPGD if we take into account the run time of each iteration, as can been seen in Fig. 3(b). This clearly demonstrates the advantage of using simple projections with IPGD compared to accurate but more complex projections with PGD.

VI-B Sparse recovery with side information

Another possible strategy to improve reconstruction that relates to our framework is using side information about the recovered signal, e.g., from estimates of similar signals. This approach was applied to improve the quality of MRI and CT scans , and also in the general context of sparse recovery .

Assume that someone gives us oracle side information on the set of Haar columns corresponding to the largest coefficients that contain $95\%$ of the energy in a patch ${\bf D}{\bf x}$. While there are many ways to incorporate the side information in the recovery, we show here how IPGD can be used for this purpose. Denoting by ${\bf P}_{{\bf x},95\%}^{oracle}$ the linear projection onto this set of columns, one may apply IPGD with $p={\bf D}{\bf P}_{{\bf x},95\%}^{oracle}{\bf D}^{*}$ and $\mathcal{K}$. As $\left\|{\bf D}{\bf x}\right\|=\left\|{\bf x}\right\|$ (since ${\bf D}$ is unitary), we have that $\epsilon=0.05$ in (24). Figure 5 compares between PGD with $\mathcal{K}$ and IPGD with $p$ and $\mathcal{K}$. We average over $100$ different randomly selected sensing matrices and patches.

Since projections onto smaller sets lead to faster convergence we suggest as in the previous example to apply PGD with an oracle projection that uses less columns from the Haar basis at the first iteration (i.e., has larger $\epsilon$) and then adds columns gradually throughout the iterations. The third (red) line in Fig. 5 demonstrates this option, where the first iterations use a projection onto the columns that contain $50\%$ of the energy of ${\bf x}$ and then every $5$ iterations the next $50$ columns correspoding to the coefficients with the largest energy are added. We continue until the columns span $95\%$ of the energy of the signal. Thus, IPGD with changing projections converges faster than IPGD with a constant $p={\bf D}{\bf P}_{{\bf x},95\%}^{oracle}{\bf D}^{*}$ but reaches the same plateau.

Typically, oracle information on the coefficients of ${\bf x}$ in the Haar basis is not accessible. Even though, it is still possible to use common statistics of the data to accelerate convergence. For example, in our case it is known that most of the energy of the signal is concentrated in the low-resolution Haar filters. Therefore, we propose to use IPGD with a projection $p$ that projects onto the first $512$ columns of the Haar basis. As before, it is possible to accelerate convergence by projecting first on a smaller number of columns and then increasing the number as the iterations proceed (in this case we add columns till IPGD coincides with PGD). These two options are presented in the fourth and fifth line of Fig. 5, respectively. Both of these options provide faster convergence, where IPGD with a fixed projection $p$ incurs a higher error as it uses less accurate projections in the last iterations compared to PGD and IPGD with changing projections. The plateau of the latter is the same one of the regular PGD (which is not attained in the graph due to its early stop) but is achieved with a much smaller number of iterations.

VII Learning the Projection – Learned IPGD (LIPGD)

In many scenarios, we may not know what type of simple operator $p$ causes $\mathcal{P}_{\mathcal{K}}(p(\cdot))$ to approximate $\hat{\mathcal{K}}$ in the best possible way. Therefore, a useful strategy is to learn $p$ for a given dataset. Assuming a linear $p$, we may rewrite (30) as

Instead of learning $p$ directly, we may learn two matrices ${\bf A}$ and ${\bf U}$, where the first replaces $p\mu{\bf M}^{*}$ and the second $p\left({\bf I}-\mu{\bf M}^{*}{\bf M}\right)$. This results in the iterations

which is very similar to those of LISTA in (7). The only difference between (35) and LISTA is the non-linear part, which is an orthogonal projection in the first and a proximal mapping in the second.

We apply this method to replace the sparse coding step in the super-resolution algorithm proposed in , where a pair of low and high resolution dictionaries is used to reconstruct the patches of the high-resolution image from the low-resolution one. In the code provided by the authors of , orthogonal matching pursuit (OMP) with sparsity $3$ is used. The complexity of this strategy corresponds to IHT with $3$ iterations. The target sparsity we use with IHT is higher ($k=40$) as it was observed to provide better reconstruction results. Note that in IHT, unlike OMP, the number of iterations may be different than the sparsity level. For optimal hyperparameter selection (such as choosing the target sparsity level), we use the training set used for the training of the dictionary in , which contains $91$ images.

Since IHT does not converge with only $3$ iterations, we apply LIPGD to accelerate convergence. We use the same dictionary dimension as in ($30\times 1000$ and $81\times 1000$ for the low and high resolution dictionaries, respectively), and train an LIPGD network to infer the sparse code of the image patches in the low-resolution dictionary. Training of the weights is performed by stochastic gradient descent with batch-size $1000$ and Nesterov momentum for adaptively setting the learning rate. We train the network using only the first $85$ images in the training set, keeping the last $6$ as a validation set. We reduce the training rate by a factor of $2$ if the validation error stops decreasing. The initial learning rate is set to $0.001$ and the Nesterov parameter to $0.9$. We use the sparse representations of the training data calculated by IHT or LIPGD to generate the high-resolution dictionary as in .

Table I summarizes the reconstruction results of regular bicubic interpolation, the OMP-based super-resolution technique of (with $3$ iterations) and its version with IHT and LIPGD (replacing OMP). It can be seen clearly that IHT leads to inferior results compared to OMP since it does not converge in $3$ iterations. LIPGD improves over both IHT and OMP as the training of the network allows it to provide good sparse approximation with only $3$ iterations. This demonstrates the efficiency of the proposed LIPGD technique, which has the same computational complexity of both OMP and IHT.

VIII Learning the Projection – LISTA Mixture Model

Though the theory in this paper applies directly only to (35) (with some constraints on ${\bf A}$ and ${\bf U}$ that stem from the constraints on $p$), the fast convergence of LISTA may be explained by the resemblance of the two methods. The success of LISTA may be interpreted as learning to approximate the set $\hat{\mathcal{K}}$ in an indirect way by learning the linear operators ${\bf A}$ and ${\bf U}$. In other words, it can be viewed as a method for learning a linear operator that together with the proximal mapping $\mathcal{S}_{f,\lambda}$ approximates a more accurate proximal mapping of a true unknown function $\hat{f}$ that leads to much faster convergence.

With this understanding, we argue that using multiple inexact projections may lead to faster convergence as each can approximate in a more accurate way different parts of the set $\hat{\mathcal{K}}$. In order to show this, we propose a LISTA mixture model (MM), similar to the Gaussian mixture model proposed in , in which we train several LISTA networks, one for each part of the dataset. Then, once we get a new vector, we apply all the networks on it in parallel (and therefore with negligible impact on the latency, which is very important in many applications) and chose the one that attains the smallest value in the objective of the minimization problem (3).

We test this strategy on the house image by extracting from it patches of size $5\times 5$, adding random Gaussian noise to each of them with variance $\sigma^{2}=25$ and then removing the DC and normalizing each. We take $7/9$ of the patches for training and $1/9$ for validation and testing. We train LISTA to minimize directly the objective (3) as in and stop the optimization after the error of the validation set increases. For the LISTA-MM we use $6$ LISTA networks such that we train the first one on the whole data. We then remove $1/6$ of the data whose objective value in (3) is the closest to the one ISTA attains after $1000$ iterations. We use this LISTA network as the initialization of the next one that is trained on the rest of the data. We repeat this process by removing in the same way the part of the data with the smallest relative error and then train the next network. After training $6$ networks we cluster the data points by selecting for each patch the network that leads to the smallest objective error for it in (3) and fine tune each network for its corresponding group of patches. We repeat this process $5$ times. The objective error of (3) as a function of the number of iterations/depth of the networks is presented in Fig. 6. Indeed, it can be seen that partitioning the data, which leads to a better approximation, accelerates convergence.

Our proposed LISTA-MM strategy bears some resemblence to the recently proposed rapid and accurate image super resolution (RAISR) algorithm . In this method, different filters are trained for different types of patches in natural images. This leads to improved quality in the attained up-scaled images with only minor overhead in the computational cost, leading to a very efficient super-resolution technique.

IX Conclusion

In this work we suggested an approach to trade-off between approximation error and convergence speed. This is accomplished by approximating complicated projections by inexact ones that are computationally efficient. We provided theory for the convergence of an iterative algorithm that uses such an approximate projection and showed that at the cost of an error in the projection one may achieve faster convergence in the first iterations. The larger the error the smaller the number of iterations that enjoy fast convergence. This suggests that if we have a budget for only a small number of iterations (with a given complexity), then it may be worthwhile to use inexact projections which can result in a worse solution in the long term but make better use of the given computational constraints. Moreover, we showed that even when we can afford a larger number of iterations, it may be worthwhile to use inexact projections in the first iterations and then change to more accurate ones at latter stages.

Appendix A Proof of Theorem V.1

The proof of Theorem V.1 relies on the following lemma.

Proof: Since ${\bf z}_{t}=\mathcal{P}_{\mathcal{K}}(p{\bf v})$ for a certain vector ${\bf v}$, we have

where the last equality follows from Lemma II.3. Using the triangle inequality with (37) leads to

We turn now to bound the first and second terms in the right-hand-side (rhs) of (38). For the second term, note that since $\mathcal{P}_{D}({\bf x}-p{\bf x})\in D$, we have

For the first term in the rhs of (38) we use the inverse triangle inequality $\left\|\mathcal{P}_{D}(p{\bf v}-p{\bf x})\right\|-\left\|\mathcal{P}_{D}(p{\bf x}-{\bf x})\right\|\leq\left\|\mathcal{P}_{D}(p{\bf v}-{\bf x})\right\|$ and (25). Combining the results leads to

Proof: The IPGD error at iteration $t+1$ is,

Using Lemma II.3 and the fact that ${\bf y}={\bf M}{\bf x}$ we have

where $(a)$ follows from the convexity of $D$ and the triangle inequality; and $(b)$ from (25) and Lemma II.4. Using Lemma A.1 with (42) leads to

Applying the inequality in (43) recursively provides the desired result. $\Box$

Appendix B Proof of Theorem V.2

The proof Theorem V.2 relies on the following lemma.

Under the same conditions of Theorem V.2,

Using (29) and the same steps of the proof of Theorem III.1, we may bound the second term in the rhs of (45) by $2\epsilon\rho(\mathcal{K})\left\|{\bf x}\right\|$. This leads to

From the inverse triangle inequality together with (29), we have that $\left\|\mathcal{P}_{D}(p{\bf v}-{\bf x})\right\|\geq\left\|\mathcal{P}_{D}(p{\bf v}-p{\bf x})\right\|-\epsilon\left\|{\bf x}\right\|$. Thus,

where the last inequality follows from the same line of argument used for deriving (40) in Lemma A.1 (with $\mathcal{K}-\mathcal{K}$ instead of $C$). $\Box$

where $(c)$ follows from the triangle inequality; and $(d)$ from (28) and Lemma II.4. Using Lemma B.1 with (48), we get

Applying (49) recursively leads to the desired result. $\Box$

Acknowledgments

RG is partially supported by GIF grant no. I-2432-406.10/2016 and ERC-StG grant no. 757497 (SPADE). YE is partially supported by the ERC grant no. 646804-ERC-COG-BNYQ. AB is partially supported by ERC-StG RAPID. GS is partially supported by ONR, NSF, NGA, and ARO. We thank Dr. Pablo Sprechmann for early work and insights into this line of research, and Prof. Ron Kimmel and Prof. Gilles Blanchard for insightful comments.

References