Title: ALISTA: ANALYTIC WEIGHTS ARE AS GOOD AS LEARNED WEIGHTS IN LISTA

Abstract: Deep neural networks based on unfolding an iterative algorithm, for example, LISTA (learned iterative shrinkage thresholding algorithm), have been an empirical success for sparse signal recovery. The weights of these neural networks are currently determined by data-driven "black-box" training. In this work, we propose Analytic LISTA (ALISTA), where the weight matrix in LISTA is computed as the solution to a data-free optimization problem, leaving only the stepsize and threshold parameters to data-driven learning. This significantly simplifies the training. Specifically, the data-free optimization problem is based on coherence minimization. We show our ALISTA retains the optimal linear convergence proved in (Chen et al., 2018) and has a performance comparable to LISTA. Furthermore, we extend ALISTA to convolutional linear operators, again determined in a data-free manner. We also propose a feed-forward framework that combines the data-free optimization and ALISTA networks from end to end, one that can be jointly trained to gain robustness to small perturbations in the encoding model.

1. INTRODUCTION:

Sparse vector recovery, or sparse coding, is a classical problem in source coding, signal reconstruction, pattern recognition and feature selection. There is an unknown sparse vector

x * = [x * 1 , • • • , x * M ] T ∈ R M .

We observe its noisy linear measurements:

EQUATION

where

b ∈ R N , D = [d 1 , • • • , d M ] ∈ R N ×M

is the dictionary, and ε ∈ R N is additive Gaussian white noise. For simplicity, each column of D, named as a dictionary kernel, is normalized, that is,

d m 2 = D :,m 2 = 1, m = 1, 2, • • • , M

. Typically, we have N M , so Equation ( 1) is an under-determined system. However, when x * is sufficiently sparse, it can be recovered faithfully. A popular approach is to solve the LASSO problem below (where λ is a scalar):

minimize x 1 2 b -Dx 2 2 + λ x 1 (2)

using iterative algorithms such as the iterative shrinkage thresholding algorithm (ISTA):

EQUATION

where η θ is the soft-thresholding function1 and L is usually taken as the largest eigenvalue of D T D.

Inspired by ISTA, the authors of (Gregor & LeCun, 2010) proposed to learn the weights in the matrices in ISTA rather than fixing them. Their methods is called Learned ISTA (LISTA) and resembles a recurrent neural network (RNN). If the iteration is truncated to K iterations, LISTA becomes a K-layer feed-forward neural network with side connections. Specifically, LISTA is:

EQUATION

If we set

EQUATION

≡ I -1 L D T D, θ (k) ≡ 1 L λ, then LISTA recovers ISTA. Given each pair of sparse vector and its noisy measurements (x * , b), applying (4) from some initial point x (0) and using b as the input yields x (k) . Our goal is to choose the parameters

Θ = {W (k) 1 , W (k)

w , θ (k) } k=0,1,...,K-1 such that x (k) is close to x * for all sparse x * following some distribution P. Therefore, given the distribution P, all parameters in Θ are subject to learning:

minimize Θ E x * ,b∼P x (K) Θ, b, x (0) -x * 2 2 .

(

EQUATION

This problem is approximately solved over a training dataset {(x * i , b i )} N i=1 sampled from P. Many empirical results, e.g., (Gregor & LeCun, 2010; Sprechmann et al., 2015; Wang et al., 2016b) , show that a trained K-layer LISTA (with K usually set to 10 ∼ 20) or its variants can generalize more than well to unseen samples (x , b ) from the same distribution and recover x from b to the same accuracy within one or two order-of-magnitude fewer iterations than the original ISTA. Additionally, the accuracies of the outputs {x (k) } of the layers k = 1, .., K gradually improve. However, such networks will generalize worse when the input deviates from the training distribution (e.g., when D varies), in contrast to the classical iterative algorithms such as ISTA that are trainingfree and thus agnostic to the input distribution. The Analysis-Synthesis model (Rubinstein & Elad, 2014; Yang et al., 2016) could also be viewed as a special LISTA model with only one layer (K = 1).

More recently, the convolutional sparse coding (CSC), an extension of the sparse coding (1), gains increasingly attention in the machine learning area. (Sreter & Giryes, 2018) showed that the CSC could be similarly approximated and accelerated by a LISTA-type feed-forward network. (Tolooshams et al., 2018) designed a structure of sparse auto-encoder inspired by multi-layer CSC. (Papyan et al., 2016; Sulam et al., 2017) also revealed CSC as a potentially useful tool for understanding general convolutional neural networks (CNNs).

1.1. RELATED WORK:

Despite the empirical success (Sprechmann et al., 2015; Wang et al., 2016a; b; c; d; Zhang & Ghanem, 2018; Zhou et al., 2018; Ito et al., 2018) in constructing fast trainable regressors for approximating iterative sparse solvers, the theoretical understanding of such approximations remains limited.

A handful of recent works have been investigating the theory of LISTA. (Moreau & Bruna, 2017) re-factorized the Gram matrix of dictionary, by trying to nearly diagonalize the Gram matrix with a basis, subject to a small 1 perturbation. They thus re-parameterized LISTA a new factorized architecture that achieved similar acceleration gain to LISTA, hence ending up with an "indirect" proof. They concluded that LISTA can converge faster than ISTA, but still sublinearly. (Giryes et al., 2018) interpreted LISTA as a projected gradient descent descent (PGD) where the projection step was inaccurate, which enables a trade-off between approximation error and convergence speed. The latest work (Chen et al., 2018) presented the more related results to ours: they introduced necessary conditions for the LISTA weight structure in order to achieve asymptotic linear convergence of LISTA, which also proved to be a theoretical convergence rate upper bound. They also introduced a thresholding scheme for practically improving the convergence speed. Note that, none of the above works extended their discussions to CSC and its similar LISTA-type architectures.

Several other works examined the theoretical properties of some sibling architectures to LISTA. (Xin et al., 2016) studied the model proposed by (Wang et al., 2016b) , which unfolded/truncated the iterative hard thresholding (IHT) algorithm instead of ISTA, for approximating the solution to 0minimization. They showed that the learnable fast regressor can be obtained by using a transformed dictionary with improved restricted isometry property (RIP). However, their discussions are not applicable to LISTA directly, although IHT is linearly convergent (Blumensath & Davies, 2009) under rather strong assumptions. Their discussions were also limited to linear sparse coding and resulting fully-connected networks only. (Borgerding et al., 2017; Metzler et al., 2017) studied a similar learning-based model inspired from another LASSO solver, called approximated message passing (AMP). (Borgerding et al., 2017) showed the MMSE-optimality of an AMP-inspired model, but not accompanied with any convergence rate result. Also, the popular assumption in analyzing AMP algorithms (called "state evolution") does not hold when analyzing ISTA.

1.2. MOTIVATION AND CONTRIBUTIONS:

This paper presents multi-fold contributions in advancing the theoretical understanding of LISTA, beyond state-of-the-art results. Firstly, we show that the layer-wise weights in LISTA need not being learned from data. That is based on decoupling LISTA training into a data-free analytic optimization stage followed by a lighter-weight data-driven learning stage without compromising the optimal linear convergence rate proved in (Chen et al., 2018) . We establish a minimum-coherence criterion between the desired LISTA weights and the dictionary D, which leads to an efficient algorithm that can analytically solve the former from the latter, independent of the distribution of x. The data-driven training is then reduced to learning layer-wise step sizes and thresholds only, which will fit the distribution of x. The new scheme, called Analytic LISTA (ALISTA), provides important insights into the working mechanism of LISTA. Experiments shows ALISTA to perform comparably with previous LISTA models (Gregor & LeCun, 2010; Chen et al., 2018) with much lighter-weight training. Then, we extend the above discussions and conclusions to CSC, and introduce an efficient algorithm to solve the convolutional version of coherence minimization. Further, we introduce a new robust LISTA learning scheme benefiting from the decoupled structure, by adding perturbations to D during training. The resulting model is shown to possess much stronger robustness when the input distribution varies, even when D changes to some extent, compared to classical LISTA models that learn to (over-)fit one specific D.

2. ANALYTIC LISTA: CALCULATING WEIGHTS WITHOUT TRAINING:

We theoretically analyze the LISTA-CPSS model defined in (Chen et al., 2018) :

EQUATION

where

EQUATION

M ] ∈ R N ×M is a linear operator with the same dimensionality with

D, x (k) = x (k) 1 , • • • , x (k)

M is the k th layer node. In (6), Θ = {W (k) , θ (k) } k are parameters to train. Model (6) can be derived from (4) with (Chen et al., 2018) showed that (6) has the same representation capability with (4) on the sparse recovery problem, with a specifically light weight structure.

W (k) 1 = (W (k) ) T , W (k) 2 = I -W (k) 1 D.

Our theoretical analysis will further define and establish properties of "good" parameters Θ in (6), and then discuss how to analytically compute those good parameters rather than relying solely on black-box training. In this way, the LISTA model could be further significantly simplified, with little performance loss. The proofs of all the theorems in this paper are provided in the appendix.

2.1. RECOVERY ERROR UPPER BOUND:

We start with an assumption on the "ground truth" signal x * and the noise ε. Assumption 1 (Basic assumptions). Signal x * is sampled from the following set:

EQUATION

In other words, x * is bounded and s-sparse2 (s ≥ 2). Furthermore, we assume ε = 0.

The zero-noise assumption is for simplicity of the proofs. Our experiments will show that our models are robust to noisy cases.

The mutual coherence of the dictionary D is a significant concept in compressive sensing (Donoho & Elad, 2003; Elad, 2007; Lu et al., 2018) . A dictionary with small coherence possesses better sparse recovery performance. Motivated by this point, we introduce the following definition.

Definition 1. Given D ∈ R N ×M with each of its column normalized, we define the generalized mutual coherence:

EQUATION

Additionally, We define W(D) = W ∈ R N ×M : W attains the infimum given (8) . A weight matrix W is "good" if W ∈ W(D).

In the above definition, problem ( 8) is feasible and attainable, i.e., W(D) = ∅, which was proven in Lemma 1 of (Chen et al., 2018) .

Theorem 1 (Recovery error upper bound). Take any x * ∈ X (B, s), any W ∈ W(D), and any sequence γ (k) ∈ (0, 2 2μs-μ+1 ). Using them, define the parameters {W (k) , θ (k) }:

EQUATION

while the sequence {x (k) (x * )} ∞ k=1 is generated by (6) using the above parameters and x (0) = 0 (Note that each x (k) (x * ) depends only on θ (k-1) , θ (k-2) , . . . and defines θ (k) ). Let Assumption 1 hold with any B > 0 and s < (1 + 1/μ)/2. Then, we have

EQUATION

)

where S is the support of x * and c

(k) = -log (2μs -μ)γ (k) + |1 -γ (k) | is a positive constant.

In Theorem 1, Eqn. ( 9) defines the properties of "good" parameters:

• The weights W (k) can be separated as the product of a scalar γ (k) and a matrix W independent of layer index k, where W has small coherence with D.

• γ (k) is bounded in an interval.

• θ (k) /γ (k) is proportional to the 1 error of the output of the k th layer.

The factor c (k) takes the maximum at γ (k) = 1. If γ (k) ≡ 1, the recovery error converges to zero in a linear rate (Chen et al., 2018) :

x (k) (x * ) -x * 2 ≤ sB exp -ck ,

where c = -log(2μs -μ) ≥ c (k) . Although γ (k) ≡ 1 gives the optimal theoretical upper bound if there are infinitely many layers k = 0, 1, 2, • • • , it is not the optimal choice for finite k. Practically, there are finitely many layers and γ (k) obtained by learning is bounded in an interval.

2.2. RECOVERY ERROR LOWER BOUND:

In this subsection, we introduce a lower bound of the recovery error of LISTA, which illustrates that the parameters analytically given by ( 9) are optimal in the convergence order (linear). Assumption 2. The signal x * is a random variable following the distribution P X . Let S = support(x * ). P X satisfies: 2 ≤ | S | ≤ s; S uniformly distributes on the whole index set; nonzero part x * S satisfies the uniform distribution with bound B: |x * i | ≤ B, ∀i ∈ S. Moreover, the observation noise ε = 0.

Theorem 1 tells that an ideal weight W ∈ W(D) satisfies I -W T D ≈ 0. But this cannot be met exactly in the overcomplete D case, i.e., N < M . Definition 2 defines the set of matrices W such that W T D is bounded away from the identity I. In Appendix D, we discuss the feasibility of (11).

Definition 2. Given D ∈ R N ×M , s ≥ 2, σmin > 0, we define a set that W (k) are chosen from:

W(D, s, σmin ) = W ∈ R N ×M σ min I-(W :,S ) T D :,S ≥ σmin , ∀ S with 2 ≤ | S | ≤ s . (11)

Based on Definition 2, we define a set that Θ = {W (k) , θ (k) } ∞ k=0 are chosen from:

Definition 3. Let {x (k) (x * )} ∞ k=1 be generated by (6) with {W (k) , θ (k) } ∞ k=0 and x (0) = 0. Then we define T as the set of parameters that guarantee there is no false positive in x (k) :

EQUATION

The conclusion (10) demonstrates that T is nonempty because "support(x (k) (x * )) ⊂ S" is satisfied as long as θ (k-1) large enough. Actually, T contains almost all "good" parameters because considerable false positives lead to large recovery errors. With T defined, we have:

Theorem 2 (Recovery error lower bound). Let the sequence {x (k) (x * )} ∞ k=1 be generated by ( 6) with {W (k) , θ (k) } ∞ k=0 and x (0) = 0. Under Assumption 2, for all parameters {W (k) , θ (k) } ∞ k=0 ∈ T and any sufficient small > 0, we have

EQUATION

with probability at least (1 -s3/2 -2 ), where c = s log(3) -log(σ min ).

This theorem illustrates that, with high probability, the convergence rate of LISTA cannot be faster than a linear rate. Thus, the parameters given in (9), that leads to the linear convergence if γ k is bounded within an interval near 1, are optimal with respect to the order of convergence of LISTA.

2.3. ANALYTIC LISTA: LESS PARAMETERS TO LEARN:

Following Theorems 1 and 2, we set W (k) = γ (k) W, where γ (k) is a scalar, and propose Tied LISTA (TiLISTA):

EQUATION

where Θ = {γ (k) } k , {θ (k) } k , W are parameters to train. The matrix W is tied over all the layers. Further, we notice that the selection of W from W(D) depends on D only. Hence we propose the analytic LISTA (ALISTA) that decomposes tied-LISTA into two stages:

EQUATION

)

where W is pre-computed by solving the following problem (Stage 1) 3 :

EQUATION

Then with W fixed, {γ (k) , θ (k) } k in (15) are learned from end to end (Stage 2). ( 16) reformulates (8) to minimizing the Frobenius norm of W T D (a quadratic objective), over linear constraints. This is a standard convex quadratic program, which is easier to solve than to solve (8) directly. Table 1 : Summary: variants of LISTA and the number of parameters to learn.

Vanilla LISTA (4) LISTA-CPSS (6) TiLISTA ( 14) ALISTA ( 15)

O(KM 2 + K + M N ) O(KN M + K) O(N M + K) O(K)

3. CONVOLUTIONAL ANALYTIC LISTA:

We extend the analytic LISTA to the convolutional case in this section, starting from discussing the convolutional sparse coding (CSC). Many works studied CSC and proposed efficient algorithms for that (Bristow et al., 2013; Heide et al., 2015; Wohlberg, 2014; 2016; Papyan et al., 2017; Garcia-Cardona & Wohlberg, 2018; Wang et al., 2018; Liu et al., 2017; 2018) . In CSC, the general linear transform is replaced by convolutions in order to learn spatially invariant features:

EQUATION

where each d m is a dictionary kernel (or filter). {d m } M m=1 is the dictionary of filters, M denotes the number of filters. {x * m } M m=1 is the set of coefficient maps that are assumed to have sparse structure, and * is the convolution operator. Now we consider 2D convolution and take 17) is pointwisely defined as5 :

4 b ∈ R N 2 , d m ∈ R D 2 , x m ∈ R (N +D-1) 2 . Equation (

EQUATION

We concatenate d m s and x m s: and rewrite (18) as:

d = [d 1 , • • • , d M ] T , x = [x 1 , • • • , x M ] T ,

EQUATION

where the matrix

D N conv (d) = [D N conv,1 (d 1 ), • • • , D N conv,M (d M )] ∈ R N 2 ×(N +D-1) 2 M

, depending on the signal size N and the dictionary d, is defined in detail in (48) in Appendix C.2.

From (17), the convolutional LISTA becomes a natural extension of the fully-connected LISTA (6):

EQUATION

where {w

(k)

m } M m=1 share the same sizes with {d m } M m=1 and (•) means a 180 rotation of the filter (Chalasani et al., 2013) . We concatenate the filters together: k) ) be the matrix induced by dictionary w (k) with the same dimensionality as D N conv (d). Since convolution can be written as a matrix form ( 19), ( 20) is equivalent to

w (k) = [w (k) 1 , • • • , w (k) M ] T ∈ R D 2 M . Parameters to train are Θ = {w (k) , θ (k) } k . Let W N conv (w (

EQUATION

Then by just substituting D, W (k) with D N conv (d), W N conv (w (k) ) respectively, Theorems 1 and 2 can be applied to the convolutional LISTA. k) ). With Assumption 1 and other settings the same with those in Theorem 1, (10) holds. With Assumption 2 and other settings the same with those in Theorem 2, (13) holds.

Proposition 1. Let D = D N conv (d) and W (k) = W N conv (w (

Similar to the fully connected case (15), based on the results in Proposition 1, we should set w

EQUATION

However, ( 22) is not as efficient to solve as ( 16). To see that, matrices

D N conv (d) and W N conv (w) are both of size N 2 × (N + D -1) 2 M , the coherence matrix W N conv (w) T D N conv (d) is thus of size (N + D -1) 2 M × (N + D -1) 2 M .

In the typical application setting of CSC, b is usually an image rather than a small patch. For example, if the image size is 100 × 100, dictionary size is 7

× 7 × 64, N = 100, D = 7, M = 64, then (N + D -1) 2 M × (N + D -1) 2 M ≈ 5 × 10 11 .

3.1. CALCULATING CONVOLUTIONAL WEIGHTS ANALYTICALLY AND EFFICIENTLY:

To overcome the computational challenge of solving (22), we exploit the following circular convolution as an efficient approximation:

b(i, j) = D-1 k=0 D-1 l=0 M m=1 d m (k, l)x m (i+k) modN , (j+l) modN +ε(i, j), 0 ≤ i, j ≤ N -1, (23) where b ∈ R N 2 , d m ∈ R D 2 , x m ∈ R N 2

. Similar to (18), we rewrite (23) in a compact way:

b = M m=1 D N cir,m (d m )x m + ε = D N cir (d)x + ε, where D N cir (d) : R N 2 M → R N 2

is a matrix depending on the signal size N and the dictionary d. Then the coherence minimization with the circular convolution is given by

EQUATION

)

The following theorem motivates us to use the solution to (24) to approximate that of ( 22). Theorem 3. The solution sets of ( 22) and ( 24) satisfy the following properties:

1. W N cir = W 2D-1 cir , ∀N ≥ 2D -1.

2.. If at least one of the matrices {D:

EQUATION

)

The solution set W N cir is not related with the image size N as long as N ≥ 2D -1, thus one can deal with a much smaller-size problem (let N = 2D -1). Further, (25) indicates that as N gets (much) larger than D, the boundary condition becomes less important. Thus, one can use W 2D-1 cir to approximate W N conv . In Appendix E.2, we introduce the algorithm details of solving (24). Based on Proposition 1 and Theorem 3, we obtain the convolutional ALISTA:

EQUATION

where 26) is a simplified form, compared to the empirically unfolded CSC model recently proposed in (Sreter & Giryes, 2018) 4 ROBUST ALISTA TO MODEL PERTURBATION Many applications, such as often found in surveillance video scenarios (Zhao et al., 2011; Han et al., 2013) , can be formulated as sparse coding models whose dictionaries are subject to small dynamic perturbations (e.g, slowly varied over time). Specifically, the linear system model (1) may have uncertain D: D = D + ε D , where ε D is some small stochastic perturbation. Classical LISTA entangles the learning of all its parameters, and the trained model is tied to one static D. The important contribution of ALISTA is to decompose fitting W w.r.t. D, from adapting other parameters {γ (k) , θ (k) } k to training data.

w = [ w1 , • • • , wM ] T ∈ W 2D-1 cir and Θ = {{γ (k) m } m,k , {θ (k) } k } are the parameters to train. (

In this section, we develop a robust variant of ALISTA that is a fast regressor not only for a given D, but all its randomly perturbations D to some extent. Up to our best knowledge, this approach is new. Robust ALISTA can be sketched as the following empirical routine (at each iteration):

• Sample a perturbed dictionary D. Sample x and ε to generate b w.r.t. D.

• Apply Stage 1 of ALISTA w.r.t. D and obtain W; however, instead of an iterative minimization algorithm, we use a neural network that unfolds that algorithm to produce W. • Apply Stage 2 of ALISTA w.r.t. W, D, x, and b to obtain {γ (k) , θ (k) } k .

In Robust ALISTA above, D becomes a part of the data for training the neural network that generates W. This neural network is faster to apply than the minimization algorithm. One might attempt to use D in the last step, rather than D, but D makes training less stable, potentially because of larger weight variations between training iterations due to the random perturbations in D. We observe that using D stabilizes training better and empirically achieves a good prediction. More details of training Robust ALISTA are given in Appendix G.

5. NUMERICAL RESULTS:

In this section, we conduct extensive experiments on both synthesized and real data to demonstrate:6 

• We experimentally validate Theorems 1 and 2, and show that ALISTA is as effective as classical LISTA (Gregor & LeCun, 2010; Chen et al., 2018) but is much easier to train. • Similar conclusions can be drawn for convolutional analytic LISTA.

• The robust analytic LISTA further shows remarkable robustness in sparse code prediction, given that D is randomly perturbed within some extent.

Notation For brevity, we let LISTA denote the vanilla LISTA model ( 4) in (Gregor & LeCun, 2010) ; LISTA-CPSS refers to the lately-proposed fast LISTA variant (Chen et al., 2018) with weight coupling and support selection; TiLISTA is the tied LISTA ( 14); and ALISTA is our proposed Analytic LISTA ( 15). If the model is for convolutional case, then we add "Conv" as the prefix for model name, such as "Conv ALISTA" that represents the convolutional analytic LISTA.

5.1. VALIDATION OF THEOREMS 1 AND 2 (ANALYTIC LISTA):

We follow the same N = 250, M = 500 setting as (Chen et al., 2018) by default. We sample the entries of D i.i.d. from the standard Gaussian distribution, D ij ∼ N (0, 1/N ) and then normalize its columns to have the unit 2 norm. We fix a dictionary D in this section. To generate sparse vectors x * , we decide each of its entry to be non-zero following the Bernoulli distribution with p b = 0.1.

The values of the non-zero entries are sampled from the standard Gaussian distribution. A test set of 1000 samples generated in the above manner is fixed for all tests in our simulations. The analytic weight W that we use in the ALISTA is obtained by solving ( 16).

All networks used (vanilla LISTA, LISTA-CPSS, TiLISTA and ALISTA) have the same number of 16 layers. We also include two classical iterative solvers: ISTA and FISTA. We train the networks with four different levels of noises: SNR (Signal-to-Noise Ratio) = 20, 30, 40, and ∞. While our theory mainly discussed the noise-free case (SNR = ∞), we hope to empirically study the algorithm performance under noise too. As shown in Figure 1 , the x-axes denotes the indices of layers for the networks, or the number of iterations for the iterative algorithms. The y-axes represent the NMSE (Normalized Mean Squared Error) in the decibel (dB) unit:

NMSE dB (x, x * ) = 10 log 10 E x -x * 2 /E x * 2 ,

where x * is the ground truth and x is the estimated one. In Figure 1 (a) noise-less case, all four learned models apparently converge much faster than two iterative solvers (ISTA/FISTA curves almost overlap in this y-scale, at the small number of iterations). Among the four networks, classical-LISTA is inferior to the other three by an obvious margin. LISTA-CPSS, TiLISTA and ALISTA perform comparably: ALISTA is observed to eventually achieve the lowest NMSE. Figure 1 (a) also supports Theorem 2, that all networks have at most linear convergence, regardless of how freely their parameters can be end-to-end learned.

Figure 1 (b) -(d) further show that even in the presence of noise, ALISTA can empirically perform comparably with LISTA-CPSS and TiLISTA, and stay clearly better than LISTA and ISTA/FISTA. Always note that ALISTA the smallest amount of parameters to learn from the end-to-end training (Stage 2). The above results endorse that: i) the optimal LISTA layer-wise weights could be structured as W (k) = γ (k) W; and ii) W could be analytically solved rather than learned from data, without incurring performance loss. We also observe the significant reduction of training time for ALISTA: while LISTA-CPSS of the same depth took ∼1.5 hours to train, ALISTA was trained within only 6 minutes (0.1 hours) to achieve comparable performance, on the same hardware (one 1080 Ti on server). We further supply Figures 2 and 3 to justify Theorem 1 from different perspectives. Figure 2 plots the learned parameters {γ (k) , θ (k) } in ALISTA (Stage 2), showing that they satisfy the properties proposed in Theorem 1: γ (k) bounded; θ (k) and γ (k) is proportional to sup x * x (k) (x * ) -x * 1 ("sup x * " is taken over the test set). Figure 3 reports the average magnitude7 of the false positives and the true positives in x k (x * ) of ALISTA: the "true positives" curve draws the values of E{ x k S (x * ) 2 2 / x k (x * ) 2 2 } w.r.t. k (the expectation is taken over the test set), while "false positives" for

E{ x k S c (x * ) 2 2 / x k (x * ) 2 2 }.

False positives take up small proportion over the positives, which supports the Theorem 1 conclusion that support(x k (x * )) ⊂ S.

5.2. VALIDATION OF THEOREM 3 (CONVOLUTIONAL ANALYTIC LISTA):

For convolutional cases, we use real image data to verify Theorem 3. We train a convolutional dictionary d with D = 7, M = 64 on the BSD500 training set (400 images), using the Algorithm 1 in (Liu et al., 2018) . We then use it for problems ( 22) and ( 24) and solve them with different N s.

In Table 2 , we take w N cir ∈ W N cir , w * ∈ W 50 cir (consider 50 as large enough) For this example, W N cir has only one element. Table 2 shows that w N cir = w * for N ≥ 13, i.e., the solution of the problem (24) is independent of N if N ≥ 2D -1, justifying the first conclusion in Theorem 3. In Table 3 , we take w N conv ∈ W N conv and w * ∈ w 13 cir , where W N conv also has only one element. 

. w N cir -w * 2 / w * 2 N = 10 N = 11 N = 12 N = 13 N = 15 N = 20 2.0 × 10 -2 9.

3 × 10 -3 3.9 × 10 -3 1.4 × 10 -12 8.8 × 10 -13 5.9 × 10 -13 

5.3. VALIDATION OF ROBUST ALISTA:

We empirically verify the effectiveness of Robust ALISTA, by sampling the dictionary perturbation ε D entry-wise i.i.d. from another Gaussian distribution N (0, σ 2 max ). We choose σ max = 0.02 and 0.03. Other simulation settings are by default the same as in Section 5.1. We then build the Robust ALISTA model, following the strategy in Section 4 and using a 4-layer encoder for approximating its second step (see Appendix G for details). Correspondingly, we compare Robust ALISTA with TiLISTA and ALISTA with specific data augmentation: we straightforwardly augment their training sets, by including all data generated with randomly perturbed Ds when training Robust ALISTA. We also include the data-free FISTA algorithm into the comparison.

Figure 4 plots the results when the trained models are applied on the testing data, generated with the same dictionary and perturbed by N (0, σ t ). We vary σ t from zero to slightly above σ max . Not surprisingly, FISTA is unaffected, while the other three data-driven models all slight degrade as σ t increases. Compared to the augmented TiLISTA and ALISTA whose performance are both inferior to FISTA, the proposed Robust ALISTA appears to be much more favorable in improving robustness to model perturbations. In both σ max cases, it consistently achieves much lower NMSE than FISTA, even when σ t has slightly surpassed σ max . Although the NMSE of ALISTA may decrease faster if σ t continues growing larger, such decrease could be alleviated by improving σ max in training, e.g., by comparing σ max = 0.02 and 0.03. Robust ALISTA demonstrates remarkable robustness and maintains the best NMSE performance, within at least the [0, σ max ] range.

6. CONCLUSIONS AND FUTURE WORK:

Based on the recent theoretical advances of LISTA, we have made further steps to reduce the training complexity and improve the robustness of LISTA. Specifically, we no longer train any matrix for LISTA but directly use the solution to an analytic minimization problem to solve for its layer-wise weights. Therefore, only two scalar sequences (stepsizes and thresholds) still need to be trained. Excluding the matrix from training is backed by our theoretical upper and lower bounds. The resulting method, Analytic LISTA or ALISTA, is not only faster to train but performs as well as the state-of-the-art variant of LISTA by (Chen et al., 2018) . This discovery motivates us to further replace the minimization algorithm by its unfolding neural network, and train this neural network to more quickly produce the weight matrix. The resulting algorithm is used to handle perturbations in the model dictionary -we only train once for a dictionary with all its small perturbations. Our future work will investigate the theoretical sensitivity of ALISTA (and its convolutional version) to noisy measurements.

A PROOF OF THEOREM 1:

In this proof, we use the notion x (k) to replace x (k) (x * ) for simplicity. We fix D in the proof, μ(D) can be simply written as μ.

Before proving Theorem 1, we present and prove a lemma. Lemma 1. With all the settings the same with those in Theorem 1, we have

EQUATION

In another word, there are no false positives in x (k) : x

(k) i = 0, ∀i / ∈ S, ∀k.

Proof. Take arbitrary x * ∈ X (B, s). We prove Lemma 1 by induction. As k = 0, (27) is satisfied since x (0) = 0. Fixing k, and assuming support(x (k) ) ⊂ S, we have

x (k+1) i =η θ (k) x (k) i -γ (k) (W :,i ) T (Dx (k) -b) =η θ (k) -γ (k) j∈S (W :,i ) T D :,j (x (k) j -x * j ) , ∀i / ∈ S .

By ( 9), the thresholds are taken as

θ (k) = μγ (k) sup x * { x (k) -x * 1 }. Also, since W ∈ W(D), we have |(W :,i ) T D :,j | ≤ μ for all j = i. Thus, for all i / ∈ S, θ (k) ≥μγ (k) x (k) -x * 1 = j∈support(x (k) ) μγ (k) x (k) j -x * j = j∈S μγ (k) x (k) j -x * j ≥ -γ (k) j∈S (W :,i ) T D :,j (x (k) j -x * j ) ,

which implies x With Lemma 1, we are able to prove Theorem 1 now.

Proof of Theorem 1. Take arbitrary x * ∈ X (B, s). For all i ∈ S, by ( 27), we obtain

x (k+1) i = η θ (k) x (k) i -γ (k) (W :,i ) T D :,S (x (k) S -x * S ) ∈ x (k) i -γ (k) (W :,i ) T D :,S (x (k) S -x * S ) -θ (k) ∂ 1 (x (k+1) i ),

where ∂ 1 (x) is the sub-gradient of |x|, x ∈ R:

∂ 1 (x) = {sign(x)} if x = 0, [-1, 1] if x = 0.

The choice of W ∈ W(D) gives (W :,i ) T D :,i = 1. Thus,

x (k) i -γ (k) (W :,i ) T D :,S (x (k) S -x * S ) =x (k) i -γ (k) j∈S,j =i (W :,i ) T D :,j (x (k) j -x * j ) -γ (k) (x (k) i -x * i ) =x * i -γ (k) j∈S,j =i (W :,i ) T D :,j (x (k) j -x * j ) + (1 -γ (k) )(x (k) i -x * i ).

Then the following inclusion formula holds for all i ∈ S,

x (k+1) i -x * i ∈ -γ (k) j∈S,j =i (W :,i ) T D :,j (x (k) j -x * j ) -θ (k) ∂ 1 (x (k+1) i ) + (1 -γ (k) )(x (k) i -x * i ).

By the definition of ∂ 1 , every element in ∂ 1 (x), ∀x ∈ R has a magnitude less than or equal to 1.

Thus, for all i ∈ S,

|x (k+1) i -x * i | ≤ j∈S,j =i γ (k) (W :,i ) T D :,j |x (k) j -x * j | + θ (k) + |1 -γ (k) | x (k) i -x * i ≤μγ (k) j∈S,j =i |x (k) j -x * j | + θ (k) + |1 -γ (k) | x (k) i -x * i . Equation (27) implies x (k) -x * 1 = x (k) S -x * S 1 for all k. Then x (k+1) -x * 1 = i∈S |x (k+1) i -x * i | ≤ i∈S μγ (k) j∈S,j =i |x (k) j -x * j | + θ (k) + |1 -γ (k) ||x (k) i -x * i | =μγ (k) (|S| -1) i∈S |x (k) i -x * i | + θ (k) |S| + |1 -γ (k) | x (k) -x * 1 =μγ (k) (|S| -1) x (k) -x * 1 + θ (k) |S| + |1 -γ (k) | x (k) -x * 1 .

Taking supremum of the above inequality over x * ∈ X (B, s), by |S| ≤ s,

sup x * { x (k+1) -x * 1 } ≤ μγ (k) (s -1) + |1 -γ (k) | sup x * { x (k) -x * 1 } + θ (k) s.

By the value of θ (k) given in ( 9), we have

sup x * { x (k+1) -x * 1 } ≤ γ (k) (2μs -μ) + |1 -γ (k) | sup x * { x (k) -x * 1 }. Let c (τ ) = -log (2μs -μ)γ (τ ) + |1 -γ (τ ) | . Then, by induction, sup x * { x (k+1) -x * 1 } ≤ exp - k τ =0 c (τ ) sup x * { x (0) -x * 1 } ≤ exp - k τ =0 c (τ ) sB.

Since x 2 ≤ x 1 for any x ∈ R n , we can get the upper bound for 2 norm:

sup x * { x (k+1) -x * 2 } ≤ sup x * { x (k+1) -x * 1 } ≤ sB exp - k τ =0 c (τ ) .

The assumption s < (1

+ 1/μ)/2 gives 2μs -μ < 1. If 0 < γ (k) ≤ 1, we have c (k) > 0. If 1 < γ (k) < 2/(1 + 2μs -μ), we have (2μs -μ)γ (k) + |1 -γ (k) | = (2μs -μ)γ (k) + γ (k) -1 < 1, which implies c (k) > 0.

Theorem 1 is proved.

B PROOF OF THEOREM 2:

Proof of Theorem 2. We fix D and sample a x * ∼ P X .

If we can prove P (13) does not hold support(

EQUATION

then the lower bound (13) in Theorem 2 is proved by

P (13) holds = S,2≤|S|≤s P (13) holds support(x * ) = S P support(x * ) = S ≥(1 -s 3/2 -2 ) 2≤|S|≤s P support(x * ) = S =1 -s 3/2 -2 .

Now we fix k and prove inequality (28) by three steps:

Step 1: If (13) does not hold, then what condition x * should satisfy?

Fixing k, we define a set X (k) ( ), which involves all the x * that does not satisfy (13):

X (k) ( ) = {(13) does not hold} = x * x (k) (x * ) -x * 2 < x * 2 σmin 3 s k .

Let S = support(x * ). For x * ∈ X (k) ( ), we consider two cases:

1. |x * i | > x * 2 (σ min /3 s ) k , ∀i ∈ S. 2. |x * i | ≤ x * 2 (σ min /3 s ) k , for some i ∈ S.

If case 1 holds, we obtain that the support of x (k) is exactly the same with that of x * :

support(x (k) (x * )) = S .

Then the relationship between x (k) and x (k-1) can be reduced to an affine transform:

x (k) S =η θ (k) x (k-1) S -(W (k-1) :,S ) T (Dx (k-1) -b) =x (k-1) S -(W (k-1) :,S

) T D :,S (x

EQUATION

Subtracting x * from the two sides of (29), we obtain

I -(W (k-1) :,S ) T D :,S (x (k-1) S -x * S ) -θ (k-1) sign(x (k) S ) 2 = x (k) S -x * S 2 = x (k) -x * 2 ,

where the last equality is due to Definition 3. Thus, for all x * ∈ X (k) ( ), if case 1 holds, we have

I -(W (k-1) :,S

) T D :,S (x

EQUATION

)

Multiplying both sides of ( 30) by (I -(W (k-1) :,S

) T D :,S ) -1 , we have

x (k-1) S -x * S -θ (k-1) (I -(W (k-1) :,S ) T D :,S ) -1 sign(x (k) S ) 2 ≤ (I -(W (k-1) :,S ) T D :,S ) -1 2 • x * 2 (σ min /3 s ) k ≤ x * 2 (σ min ) k-1 3 -ks ,

where the last inequality is due to (11). Let x(k-1) denote the bias of x (k-1) :

EQUATION

S ), then we get a condition that x * satisfies if case 1 holds:

X (k-1) ( ) = x * x (k-1) S (x * ) -x * S -x(k-1) (x * ) 2 ≤ x * 2 (σ min ) k-1 3 -ks .

If case 2 holds, x * belongs to the following set:

X (k) ( ) = x * |x * i | ≤ x * 2 σmin /3 s k , for some i ∈ S .

Then for any x * ∈ X (k) ( ), either x * ∈ X (k-1) ( ) or x * ∈ X (k) ( ) holds. In another word,

X (k) ( ) ⊂ X (k) ( ) ∪ X (k-1) ( ).

Step 2: By imitating the construction of X (k) ( ), we construct

X (k-2) ( ), X (k-3) ( ), • • • .

Similar to:

Step 1, we divide X (k-1) ( ) into two sets: X (k-1) ( ) and X (k-2) ( ), then we divide

X (k-2) ( ) into X (k-2) (

) and X (k-3) ( ). Repeating the process, until dividing X (1) ( ) into X (1) ( ) and X (0) ( ).

By induction, we have

EQUATION

where the sets are defined as follows for all j = 0, 1, 2, • • • , k:

EQUATION

X (k-j) ( ) = x * x (k-j) S (x * ) -x * S -x(k-j) (x * ) 2 ≤ x * 2 σmin k-j 3 -ks (33)

and the bias is defined as following for all j = 0, 1, 2, • • • , k:

EQUATION

)

Step 3: Estimating the probabilities of all the sets in (31).

By (31), we have

P x * ∈ X (k) ( ) support(x * ) = S ≤ k-1 j=1 P x * ∈ X (k-j) ( ) support(x * ) = S + P x * ∈ X (0) ( ) support(x * ) = S .

Now we have to prove that each of the above terms is small, then P (x * ∈ X (k) ( )|support(x * ) = S) is small and (28) will be proved.

Define a set of n-dimensional sign numbers

Si(n) = (s 1 , s 2 , • • • , s n ) s i ∈ {0, -1, 1}, ∀i = 1, • • • , n . Since sign x (k-j+t) S ∈ Si(|S|) for all t = 1, 2, • • • , j, {sign(x (k-j+t) S

)} j t=1 has finitely possible values. Let sign(x

(k-j+t) S ) = s (t) for t = 1, 2, • • • , j. Then x(k-j) i (x *

) is independent of x * and can be written as

x(k-j) i (s (1) , s (2) , • • • , s (j)

). Thus, we have

P (x * ∈ X (k-j) ( )|support(x * ) = S) = i∈S s (1) ∈Si(|S|) s (2) ∈Si(|S|) • • • s (j) ∈Si(|S|) P |x * i + x(k-j) i (x * )| < x * 2 (σ min ) k-j 3 -ks , sign(x (k) S ) = s (1) , • • • , sign(x (k-j+1) S ) = s (j) support(x * ) = S ≤ i∈S s (1) ∈Si(|S|) s (2) ∈Si(|S|) • • • s (j) ∈Si(|S|) P |x * i + x(k-j) i (s (1) , s (2) , • • • , s (j) )| < |S|B(σ min ) k-j 3 -ks support(x * ) = S ≤ i∈S s (1) ∈Si(|S|) s (2) ∈Si(|S|) • • • s (j) ∈Si(|S|) |S|B(σ min ) k-j 3 -ks B =|S|3 j|S| ( |S| (σ min ) k-j 3 -ks ≤ |S| 3/2 (σ min ) k-j 3 (j-k)|S|

where the second inequality comes from the uniform distribution of x * S (Assumption 2), the last inequality comes from |S| ≤ s.

The last term, due to the uniform distribution of x * S and x (0) = 0, can be bounded by P

(x * ∈ X (0) ( )|support(x * ) = S) =P x * + x(0) (x * ) 2 ≤ x * 2 3 -ks support(x * ) = S = s (1) ∈Si(|S|) s (2) ∈Si(|S|) • • • s (k) ∈Si(|S|) P x * + x(0) (x * ) 2 ≤ x * 2 3 -ks , sign(x (1) S ) = s (1) , • • • , sign(x (k) S ) = s (k) support(x * ) = S ≤3 k|S| ( 3 -ks ) |S| ≤ |S| .

Then we obtain P 28) is proved.

(x * ∈ X (k) ( )|support(x * ) = S) ≤ k-1 j=0 |S| 3/2 (σ min ) k-j 3 (j-k)|S| + |S| = k j=1 |S| 3/2 (σ min ) j 3 -j|S| + |S| = |S| 3/2 σmin 3 -|S| 1 -σmin 3 -|S| 1 -(σ min 3 -|S| ) k + |S| ≤ |S| 3/2 + |S| . Then (

C PROOF OF THEOREM 3:

There are two conclusions in Theorem 3. We prove the two conclusions in the following two subsections respectively.

C.1 PROOF OF CONCLUSION 1.

Before proving Conclusion 1, we analyze the operator D N:

cir in detail. The circular convolution ( 23) is equivalent with:

b(i, j) = N -1 k=0 N -1 l=0 M m=1 D N cir (i, j; k, l, m)x m (k, l), 0 ≤ i, j ≤ N -1,

where the circulant matrix is element-wise defined as:

EQUATION

Similarly, the corresponding circulant matrix W N cir (i, j; k, l, m) of dictionary w is:

EQUATION

As we defined in Section 3, b is a vector. With

x = [x 1 , • • • , x M ] T ,

x is a vector. Then the operator D N cir is a matrix, where (i, j) is its row index and (k, l, m) is its column index. Define a function measuring the difference between i and k:

I(i, k) (k -i) modN , 0 ≤ i, k ≤ N -1.

The coherence between D N cir (i, j; k, l, m) and W N cir (i, j; k, l, m): B coh = (D N cir ) T W N cir is elementwise defined by:

B coh (k 1 , l 1 , m 1 ; k 2 , l 2 , m 2 ) = N -1 i=0 N -1 j=0 D N cir (i, j; k 1 , l 1 , m 1 )W N cir (i, j; k 2 , l 2 , m 2 ) = i∈I(k1,k2) j∈J (l1,l2)

d m1 I(i, k 1 ), I(j, l 1 ) w m2 I(i, k 2 ), I(j, l 2 ) .

where

I(k 1 , k 2 ) = {i|0 ≤ i ≤ N -1, 0 ≤ I(i, k 1 ) ≤ D -1, 0 ≤ I(i, k 2 ) ≤ D -1}, J (l 1 , l 2 ) = {j|0 ≤ j ≤ N -1, 0 ≤ I(j, l 1 ) ≤ D -1, 0 ≤ I(j, l 2 ) ≤ D -1}. Lemma 2. Given N ≥ 2D -1, it holds that: (a) I(k 1 , k 2 ) = ∅ if and only if " 0 ≤ (k 1 -k 2 ) modN ≤ D -1" or " 0 < (k 2 -k 1 ) modN ≤ D -1" holds. (b) J (l 1 , l 2 ) = ∅ if and only if " 0 ≤ (l 1 -l 2 ) modN ≤ D -1" or " 0 < (l 2 -l 1 ) modN ≤ D -1" holds.

Proof. Now we prove Conclusion (a). Firstly, we prove "if.

EQUATION

EQUATION

)

Secondly, we prove "only if.

" If I(k 1 , k 2 ) = ∅, we can select an i ∈ I(k 1 , k 2 ). Let r 1 = (k 1 - i) modN and r 2 = (k 2 -i) modN . By the definition of I(k 1 , k 2 ), we have 0 ≤ r 1 , r 2 ≤ D -1.

Two cases should be considered here. Case 1:

r 1 ≥ r 2 . Since 0 ≤ r 1 -r 2 ≤ D -1 ≤ N -1, it holds that r 1 -r 2 = (r 1 -r 2 ) modN . Thus, r 1 -r 2 = (r 1 -r 2 ) modN = (k 1 -i) modN -(k 2 -i) modN modN = (k 1 -i) -(k 2 -i) modN =(k 1 -k 2 ) modN .

The equality "0 ≤ r 1 -r 2 ≤ D -1" leads to the conclusion "0 ≤ (k 1 -k 2 ) modN ≤ D -1". In case 2 where r 1 < r 2 , we can obtain 0 < (k 2 -k 1 ) modN ≤ D -1 with the similar arguments.

Conclusion (b) can be proved by the same argument with the proof of (a). Lemma 2 is proved.

Now we fix k 1 , l 1 and consider what values of k 2 , l 2 give I(k 1 , k 2 ) = ∅ and J (l 1 , l 2 ) = ∅. Define four index sets given 0 ≤ k 1 , l 1 ≤ N -1:

K(k 1 ) ={k|0 ≤ (k 1 -k) modN ≤ D -1} K(k 1 ) ={k|0 < (k -k 1 ) modN ≤ D -1} L(l 1 ) ={l|0 ≤ (l 1 -l) modN ≤ D -1} L(l 1 ) ={l|0 < (l -l 1 ) modN ≤ D -1} Lemma 3. If N ≥ 2D -1, we have: (a) The cardinality of K(k 1 ), K(k 1 ): | K(k 1 )| = D, | K(k 1 )| = D -1. (b) K(k 1 ) ∩ K(k 1 ) = ∅. (c) The cardinality of L(l 1 ), L(l 1 ): | L(l 1 )| = D, | L(l 1 )| = D -1. (d) L(l 1 ) ∩ L(l 1 ) = ∅.

Proof. Now we prove Conclusion (a). The set K(k 1 ) can be equivalently written as

K(k 1 ) = {(k 1 -r k ) modN |r k = 0, 1, • • • , D -1} (39) Let k(r k ) = (k 1 -r k ) modN . We want to show that k(r 1 k ) = k(r 2 k ) as long as r 1 k = r 2 k . Without loss of generality, we assume 0 ≤ r 1 k < r 2 k ≤ D -1.

By the definition of modulo operation, There exist two integers q, q such that

k(r 1 k ) = qN + k 1 -r 1 k , k(r 2 k ) = q N + k 1 -r 2 k . Suppose k(r 1 k ) = k(r 2 k ).

Taking the difference between the above two equations, we obtain

r 2 k - r 1 k = (q -q)N , i.e, N divides r 2 k -r 1 k . However, 0 ≤ r 1 k < r 2 k ≤ D -1 implies 1 ≤ r 2 k -r 1 k ≤ D -1 ≤ N -1, which contradicts with "N dividing r 2 k -r 1 k ." Thus, it holds that k(r 1 k ) = k(r 2 k ). Then we have | K(k 1 )| = D.

In the same way, we have

K(k 1 ) = {(k 1 + r k ) modN |r k = 1, 2, • • • , D -1} (40) and | K(k 1 )| = D -1. Conclusion (a) is proved. Now we prove Conclusion (b). Suppose K(k 1 ) ∩ K(k 1 ) = ∅. Pick a k 2 ∈ K(k 1 ) ∩ K(k 1 ). Let r 3 = (k 1 -k 2 ) modN and r 4 = (k 2 -k 1 ) modN . Then we have 0 ≤ r 3 ≤ D -1 and 0 < r 4 ≤ D -1.

By the definition of modulo operation, There exist two integers q, q such that

k 1 -k 2 = qN + r 3 , k 2 -k 1 = q N + r 4 which imply r 3 + r 4 + (q + q )N = 0. However, 0 < r 3 +r 4 ≤ 2D -2 contradicts with "q ∈ Z, q ∈ Z, N ∈ Z, N ≥ 2D -1." Conclusion (b) is proved.

Conclusions (c) and (d) are actually the same with Conclusions (a) and (b) respectively. Thus, it holds that

EQUATION

and

| L(l 1 )| = D, | L(l 1 )| = D -1. Lemma 3 is proved.

With the preparations, we can prove Conclusion 1 of Theorem 3 now.

Proof of Theorem 3, Conclusion 1. Firstly we fix k:

1 ∈ {0, 1, • • • , N -1} and consider k 2 ∈ K(k 1 ). Let r k = (k 1 -k 2 ) modN .

Then equation (37) implies that, for any i ∈ I(k 1 , k 2 ), there exists a δ (r k ≤ δ ≤ D -1) such that

EQUATION

Now we consider another case for k 2 :

EQUATION

Similarly, for any

EQUATION

Another case for l 2 : l 2 ∈ L(l 1 ), r l = (l 2 -l 1 ) modN . For any j ∈ J (l 1 , l 2 ), there exists a δ (r l ≤ δ ≤ D -1) such that

EQUATION

Now let us consider the following function. By results in Lemmas 2 and 3, we have

f (k 1 , l 1 , m 1 , m 2 ) = N -1 k2=0 N -1 l2=0 B coh (k 1 , l 1 , m 1 ; k 2 , l 2 , m 2 ) 2 =f 1 + f 2 + f 3 + f 4 ,

where

f 1 = k2∈K(k1) l2∈L(l1) B coh (k 1 , l 1 , m 1 ; k 2 , l 2 , m 2 ) 2 f 2 = k2∈ K(k1) l2∈L(l1) B coh (k 1 , l 1 , m 1 ; k 2 , l 2 , m 2 ) 2 f 3 = k2∈K(k1) l2∈ L(l1) B coh (k 1 , l 1 , m 1 ; k 2 , l 2 , m 2 ) 2 f 4 = k2∈ K(k1) l2∈ L(l1) B coh (k 1 , l 1 , m 1 ; k 2 , l 2 , m 2 ) 2 .

Combining equations ( 39), ( 41), ( 43) and ( 45), we obtain

f 1 = D-1 r k =0 D-1 r l =0 D-1 δ k =r k D-1 δ l =r l d m1 (δ k , δ l )w m2 (δ k -r k , δ l -r l ) 2 .

Combining ( 40), ( 41), ( 44) and ( 45), we obtain

f 2 = D-1 r k =1 D-1 r l =0 D-1 δ k =r k D-1 δ l =r l d m1 (δ k -r k , δ l )w m2 (δ k , δ l -r l ) 2 .

Combining ( 39), ( 42), ( 43) and ( 46), we obtain

f 3 = D-1 r k =0 D-1 r l =1 D-1 δ k =r k D-1 δ l =r l d m1 (δ k , δ l -r l )w m2 (δ k -r k , δ l ) 2 .

Combining (40), ( 42), ( 44) and ( 46), we obtain

f 4 = D-1 r k =1 D-1 r l =1 D-1 δ k =r k D-1 δ l =r l d m1 (δ k -r k , δ l -r l )w m2 (δ k , δ l ) 2 .

By the above explicit formulas of f i , 1 ≤ i ≤ 4, we have f 1 , f 2 , f 3 , f 4 are all independent of k 1 , l 1 and N . They are only related with m 1 , m 2 for fixed d and m. Thus, we are able to denote

EQUATION

which implies Similar to D cir , we use (i, j) as the row index and (k, l, m) as the column index of D conv . For

W N cir = W 2D-1 cir , ∀N ≥ 2D -1.

0 ≤ i, j ≤ N -1, 1 ≤ m ≤ M , D N +D-1 cir (i, j; k, l, m) = D N conv (i, j; k, l, m) = d m (k -i, l -j), 0 ≤ k -i, l -j ≤ D -1 0, k, l taken as others (48) Matrix D N +D-1 cir is of dimension (N + D -1) 2 × (N + D -1) 2 M , where 0 ≤ i, j ≤ N + D -2; matrix D N conv is of dimension (N ) 2 × (N + D -1) 2 M , where 0 ≤ i, j ≤ N -1. Thus, D N conv is a block in D N +D-1 cir , i.e., D N +D-1 cir = D N conv ∆ N D . The matrix ∆ N D is of dimension ((N + D -1) 2 -N 2 ) × (N + D -1) 2 M : ∆ N D = D N +D-1 cir (i, j; :, :, :) , (i, j) ∈ I ∆ where I ∆ = I 1 ∪ I 2 ∪ I 3 I 1 ={(i, j)|N ≤ i ≤ N + D -2, 0 ≤ j ≤ N -1} I 2 ={(i, j)|0 ≤ i ≤ N -1, N ≤ j ≤ N + D -2} I 3 ={(i, j)|N ≤ i ≤ N + D -2, N ≤ j ≤ N + D -2}.

Similarly,

W N +D-1 cir = W N conv ∆ N W , ∆ N W = W N +D-1

cir (i, j; :, :, :) , (i, j) ∈ I ∆ .

Then,

EQUATION

Lemma 4. For any (i, j) ∈ I ∆ , one has

EQUATION

W N +D-1 cir (i, j; :, :, :

EQUATION

Proof. Equation ( 35) implies that, for (i, j)

∈ I 1 , 1 ≤ m ≤ M , D N +D-1 cir (i, j; k, l, m) =    d m (k -i, l -j), i ≤ k ≤ N + D -2, j ≤ l ≤ j + D -1 d m (k -i + N + D -1, l -j), 0 ≤ k ≤ i -N, j ≤ l ≤ j + D -1 0,

k, l taken as others Thus, for any (i, j) ∈ I 1 ,

D N +D-1 cir (i, j; :, :, :) 2 2 = N +D-2 k=0 N +D-2 l=0 M m=1 D N +D-1 cir (i, j; k, l, m) 2 = d 2 2

Similarly, D N +D-1 cir (i, j; :, :, : 50) is proved. With the same argument, equation ( 51) is also proved.

) 2 2 = d 2 2 , (i, j) ∈ I 2 ∪ I 3 . Equation (

EQUATION

Proof. For simplicity, we denote two row vectors:

d i,j D N +D-1 cir (i, j; :, :, :) ∈ R 1×(N +D-1) 2 M w i,j W N +D-1 cir (i, j; :, :, :) ∈ R 1×(N +D-1) 2 M

Then,

(∆ N D ) T ∆ N W 2 F = (i,j)∈I∆ d T i,j w i,j 2 F = (i1,j1)∈I∆ (i2,j2)∈I∆ d T i1,j1 w i1,j1 , d T i2,j2 w i2,j2 F , where d T i1,j1 w i1,j1 , d T i2,j2 w i2,j2 F = trace w T i1,j1 d i1,j1 d T i2,j2 w i2,j2 = (d i1,j1 d T i2,j2 ) • (w i1,j1 w T i2,j2

).

Since

d i1,j1 d T i2,j2 = N -1 k=0 N -1 l=0 M m=1 = D N +D-1 cir (i 1 , j 1 ; k, l, m)D N +D-1 cir (i 2 , j 2 ; k, l, m),

with the same argument in Lemma 2, we have:

d i1,j1 d T i2,j2 = 0 implies i 2 ∈ I ∆ {i|0 ≤ (i 1 -i) mod(N +D-1) ≤ D -1 or 0 ≤ (i -i 1 ) mod(N +D-1) ≤ D -1} j 2 ∈ J ∆ {j|0 ≤ (j 1 -j) mod(N +D-1) ≤ D -1 or 0 ≤ (j -j 1 ) mod(N +D-1) ≤ D -1} Then (∆ N D ) T ∆ N W 2 F = (i1,j1)∈I∆ i2∈I ∆ j2∈J ∆ (d i1,j1 d T i2,j2 ) • (w i1,j1 w T i2,j2 ) ≤ (i1,j1)∈I∆ i2∈I ∆ j2∈J ∆ d 2 2 w 2 2 =| I ∆ | • | I ∆ | • | J ∆ | • d 2 2 w 2 2 = 2N (D -1) + (D -1) 2 (2D -1) 2 d 2 2 w 2 2

, where the inequality in the second line follows from ( 50) and (51). Inequality (52) is proved.

With these preparations, we can prove Theorem 3, Conclusion 2 now.

Proof of Theorem 3, Conclusion 2. Define set:

EQUATION

Since d ∈ W normal , the set is nonempty:

EQUATION

) Define functions F N conv : R D 2 M → R, F N cir : R D 2 M → R. F N conv (w) = 1 N + D -1 D N conv (d) T W N conv (w) F + ι W normal (w) F N cir (w) = 1 N (D N cir (d)) T W N cir (w) F + ι W normal (w)

By the definitions of W N conv , W N cir , we have

W N conv = arg min w F N conv (w), W N cir = arg min w F N cir (w)

Step 1: Proving F N conv (w) uniformly converges to F 2D-1 cir (w) on X ∩ W normal for any compact set X ⊂ R D 2 M .

We arbitrarily choose such a compact set X. Based on (47), ( 49) and ( 52), one has, for all w ∈ X ∩ W normal ,

|F N conv (w) -F 2D-1 cir (w)| =|F N conv (w) -F N +D-1 cir (w)| = 1 N + D -1 (D N +D-1 cir ) T W N +D-1 cir F -(D N conv ) T W N conv F ≤ 1 N + D -1 (∆ N D ) T ∆ N W F ≤ 2N (D -1) + (D -1) 2 (2D -1) N + D -1 d 2 w 2 ≤ (2D -1) 2(D -1) √ N + D -1 d 2 w 2 .

Thus, there exists a constant B > 0, which is independent of N , such that

EQUATION

)

Step 2: Proving F N conv (w) epigraphically converges8 to F 2D-1 cir (w).

We want to show, at each point w it holds that lim inf

EQUATION

Firstly, we prove (56). We arbitrarily pick a sequence {w N } ∞ N =0 such that w N → w. If w / ∈ W normal , F 2D-1 cir (w) = +∞. Since W normal is a closed set, there exists a N + such that w N / ∈ W normal for all N ≥ N + . Thus, one has F N conv (w N ) = +∞ for all N ≥ N + , i.e., lim inf

N →∞ F N conv (w N ) = F 2D-1 cir (w) = +∞.

If w ∈ W normal , two cases should be considered. The first case is that any subsequences of {w N } ∞ N =0 are not kept within W normal , i.e., there exists a N + such that w N / ∈ W normal for all N ≥ N + . Then we have lim inf

N →∞ F N conv (w N ) = +∞ > F 2D-1 cir (w).

The second case is that there exists a subsequence {w

N k } ∞ k=0 ⊂ {w N } ∞ N =0 such that w N k ∈ W normal , ∀k = 0, 1, 2, • • • .

Since w N converges to w, any subsequences should be Cauchy. Given any Cauchy sequence {w N k } ∞ k=0 in finite dimensional Euclidean space, there exists a compact set X such that w N k ∈ X, ∀k = 0, 1, 2, • • • Let B = sup w∈X∩W normal w . By (55), we obtain

|F N k conv (w N k ) -F 2D-1 cir (w)| ≤|F N k conv (w N k ) -F 2D-1 cir (w N k )| + |F 2D-1 cir (w N k ) -F 2D-1 cir (w)| ≤ BB √ N k + |F 2D-1 cir (w N k ) -F 2D-1 cir (w)|.

For any > 0, by the continuity of F 2D-1 cir , we are able to find a K > 0 such that |F 2D-1 cir

(w N k ) - F 2D-1 cir (w)| < for all k ≥ K. Pick a K such that N K ≥ (BB / ) 2 . Then, for all k ≥ max(K, K ), we have |F N k conv (w N k ) -F 2D-1 cir (w)| < 2 , i.e., lim k→∞ F N k conv (w N k ) = F 2D-1 cir (w).

The above conclusion holds for all subsequences {w

N k } ∞ k=0 ⊂ W normal . F 2D-1 cir (w) is an accu- mulation point of {F N conv (w N )} ∞ N =0 . All the other accumulation points of {F N conv (w N )} ∞ N =0 must be +∞ because F N conv (w) = F 2D-1 cir (w) = +∞ for all w / ∈ W normal . Thus, lim inf N →∞ F N conv (w N ) = F 2D-1 cir (w) < +∞.

Secondly, we prove (57). We set w N = w for all N = 0, 1, 2, • • • . Then ( 57) is a direct result of (55).

Step 3: proving (25). Define

G(w) = (D 2D-1 cir ) T W 2D-1 cir 2

F . We want to show that G(w) is strongly convex.

Let wi ∈ R (2D-1) 2 be the i th column of W 2D-1 cir , i.e.,

W 2D-1 cir = w1 , w2 , • • • , w(2D-1) 2 M Then G(w) = (2D-1) 2 M i=1 ( wi ) T D 2D-1 cir (D 2D-1 cir ) T wi .

Let w ∈ R (2D-1) 4 M vectorize W 2D-1 cir , i.e., w = ( w1 ) T , ( w2 ) T , • • • , ( w(2D-1) 2 M ) T T .

Then G(w) can be written as a quadratic form of w:

G(w) = wT Q w,

where

Q =     D 2D-1 cir (D 2D-1 cir ) T • • • D 2D-1 cir (D 2D-1 cir ) T     totally (2D-1) 2 M diagonal blocks .

As long as at least one of the matrices

{D 2D-1 cir,0 , • • • , D 2D-1 cir,M -1 } is non-singular, D 2D-1 cir is full row rank, which implies that D 2D-1 cir (D 2D-1 cir ) T is non-singular. Then Q is positive definite.

The transform between w and w is linear. We denote the transform as T , i.e., w = T w. . This is proved in Step 2.

2.. F:

2D-1 cir is level bounded. Since G(w) is strongly convex, F 2D-1 cir (w) = G(w) + ι W normal (w) must be level bounded. 3. F 2D-1 cir ≡ +∞. Since W normal is nonempty (54), dom F 2D-1 cir = ∅, F 2D-1 cir is not con- stantly +∞.

4. All the level set of F N conv are connected. This can be derived from the convexity of F N conv .

5. F 2D-1 cir and F N conv are all lower semi-continuous and proper. This condition follows from the fact that the functions F In this section, we want to numerically show that, given typical D and s, there is a σmin > 0 such that a random generated matrix W ∈ W(D, s, σmin ). However, given D and W, it's intractable to completely check (11):

σ min I -(W :,S ) T D :,S ≥ σmin , ∀ S with 2 ≤ | S | ≤ s.

The reason is that there are extremely large amount of possible S s. For example, we take M = 250, N = 500, s = 50. Instead of checking all possible Ss, we sample 5000 Ss from the whole set:

S ⊂ S = {S : S ⊂ {1, 2, • • • , 500}|2 ≤ |S| ≤ s},

where S is the set of all the samples. Then we estimate σmin with the following quantity:

σ (D, W) = min S∈S σ min I -(W :,S ) T D :,S

Furthermore, we use the same D as that in Section 5 and generate 1000 Ws with each entry i.i.d sampled from the normal distribution. Then we normalize each column of the generated Ws. This technique is commonly used in sparse coding. Finally, we report the distribution of σ (D, W) with the fixed D and the 1000 sampled Ws in Figure 5 . In this section, we introduce an algorithm to solve (16) (we copy (16) below to facilitate reading): By the definition of the Frobenius norm, it holds that

min W∈R N ×M W T D 2 F , s.t. (W :,m ) T D :,m = 1, ∀m = 1, 2, • • • , M,

EQUATION

Thus, the above problem is equivalent with

min W∈R N ×M D T W 2 F , s.t. (D :,m ) T W :,m = 1, ∀m = 1, 2, • • • , M.

We apply projected gradient descent (PGD) to solve the above problem. The gradient of

D T W 2 F is ∇ D T W 2 F = DD T W. Denote the set by W = {W ∈ R N ×M |(D :,m ) T W :,m = 1, ∀m = 1, 2, • • • , M.}

Then the projection onto W can be calculated by

Proj W (W) = W + ∆W, ∆W = (1 -(D :,1 ) T W :,1 )W :,1 , • • • , (1 -(D :,M ) T W :,M )W :,M

With these formulas, we are able to write down the PGD, which is listed in Algorithm 1.

Algorithm 1: Projected gradient descent for solving (16)

Input: Dictionary D ∈ R N ×M . Initialize: Let W 0 = D. 1 for j = 0, 1, 2, . . . until convergence do 2 Update W by W j+1 = Proj W W j -ηD(D) T W j .

3. end:

Output: W J , where J is the last iterate.

In each step, calculating the gradient has the complexity of O(N 2 M ) because DD T can be precomputed. Calculating the projection takes O(N M ) time consumptions. Due to the objective function to minimize in ( 16) is restricted strongly convex, Algorithm 1 is linear convergent (Zhang & Cheng, 2015) . To get an -accurate solution, PGD takes O(log(1/ )) steps. Thus, the complexity of Algorithm 1 is O(log(1/ )N 2 M ). We should note that the bounds given in In this section, we introduce an algorithm to solve (24) (we copy (24) below to facilitate reading):

min w∈R D 2 M wm•dm=1, 1≤m≤M W N cir (w) T D N cir (d) 2 F .

Similarly, by (58), the above problem is equivalent with

EQUATION

Since the circular convolution is very efficient to calculate in the frequency domain, we consider solving (59) utilizing the fast Fourier transform (FFT).

Firstly, we introduce the operators D N cir (d), W N cir (w) in the frequency domain. To simplify the notation, we denote the operators as D N cir and W N cir respectively. Let F be the FFT operator. Thus, The frequency domain operator DN cir is much cheaper to calculate than the operator D N cir in the spacial domain because it is block diagonal (Wohlberg, 2016) . Specifically, we zero pad d to N × N and do FFT: dm = FFT zeropad(d m , N -D) , then the above operator can be explicitly written as: b = M m=1 dm xm , where • means complex conjugate. This is due to D N cir is actually cross-correlation, not convolution (see ( 18)). Cross-correlation is equal to the transpose of convolution. Thus, there should be complex conjugate in the frequency domain. which can be efficiently solved by the frequency domain ISTA in (Liu et al., 2017) . The details are outlined in Algorithm 2. Inspired by the a similar unrolling and truncating fashion in LISTA, we can approximately solve the coherence minimization problem (16) using a similar finite-layer neural network that is unfolded from iterative algorithms. Because the linear constraints in ( 16) are hard to enforce in deep neural networks, we first relax it to the following form:

b = D N cir x is equivalent with F b = F D N cir F H F x. Let b = F b, x = F x

F VISUALIZATION OF THE ANALYTIC CONVOLUTIONAL WEIGHTS:

EQUATION

The results are shown in Table . 5, where the performances are measured with NMSE in dB, which is defined in Section 5.1. The "Augmented" prefix means the models are trained in the data augmentation setting. σ is the standard deviation of the Gaussian distribution that is used to generate the dictioanry perturbations. T stands for the number of layers (in the case of robust ALISTA it means the nubmer of layers of the ALISTA decoder, with a 4-layer encoder). We follow the training strategy and settings explained in Appendix G, with σ max = 0.02 during training.

On one hand, the comparison of performances of ALISTA, TiLISTA and LISTA-CPSS shows results that are consistent to the intuition that larger parameter space yields larger learning capacity, and therefore, better adaptiveness (LISTA-CPSS > TiLISTA > ALISTA). On the other hand, we can also find that ALISTA with more layers has worse performance. We think this observation is also reasonable for two reasons: 1) adding more layers in ALISTA does not enlarge the parameter volume significantly because it has only two scalar parameters in each layer; noting that ALISTA uses a fixed, analytically solved weight matrix, if this weight matrix is not compatible with the perturbed dictionary, more layers can even hurt the performance instead of improving. Lastly, it's clearly shown that robust ALISTA outperforms LISTA-CPSS, even if it contains less parameters. This proves that the encoding process that adaptively transforms the perturbed dictiories is necessary to achieve good robustness against perturbations in dictionaries.

σ of perturbations during testing 0.0001 0.001 0.01 0.015 0.02 0.025 Augmented ALISTA Augmented LISTA-CPSS Table 5 : The results (recovery NMSE in dB) of ablation study on the influence of model capacity towards the model robustness against dictionary perturbations.

I.2 ABLATION STUDY ON THE DEPTHS OF ENCODERS IN ROBUST ALISTA:

Another constructive suggestion from the reviewers is to design an ablation study to investigate the influence of the depth of encoders in robust ALISTA on its adaptivity to dictionary perturbations. A natural intuition is that, adding more layers to the encoder can increase its ability to sustain larger perturbation levels. But is this true? Therefore, we train another two robust ALISTA models, with a 5-layer and a 6-layer encoders respectively and 16-layer ALISTA decoders for both, and compare them with the originally reported robust ALISTA model with a 4-layer encoder and a 16-layer ALISTA decoders. All three models use one pretrained decoder, and pretrain their encoders using the same method (see Appendix G). We only use σ max = 0.02 during training. One thing to notice is that we observe unstable training process if we use default initial learning rates α pre = 10 -4 in the pre-training stage and α e = 10 -6 in the joint training stage when encoders have 5 or 6 layers. Therefore, we use α pre = 10 -5 for the 6-layer encoder in the pre-training stage, and in the joint training stage use a decreased and uniform initial learning rate α e = 10 -9 for the three encoders while keeping the default initial learning rate α d = 10 -4 for decoders. The other settings remain the same. -68.57 -68.56 -67.94 -66.86 -64.84 -56.63 5 -69.34 -69.34 -69.02 -68.49 -67.20 -65.55 6 -70.38 -70.33 -69.92 -69.22 -67.72 -65.60 Table 6 : The results of ablation study on the influence of the depths of encoder towards the model robustness against dictionary perturbations.

Results are shown in:



:

Soft-thresholding function is defined in a component-wise way: η θ (x) = sign(x) max(0, |x| -θ)

A signal is s-sparse if it has no more than s non-zero entries.

Some details and a complexity analysis of Stage 1 are discussed in Appendix E.1

Here, b, dm, xm are vectors. The notion b(i, j) means the (iN + j) th entry of b. Additionally, dm, xm are defined in the same way for all m = 1, • • • , M .

Strictly speaking, (18) is the cross-correlation rather than convolution. However in TensorFlow, that operation is named as convolution, and we follow that convention to be consistent with the learning community.

Our codes are uploaded to https://github.com/xchen-tamu/alista.

The number and proportion of false alarms are a more straightforward performance metric. However, they are sensitive to the threshold. We found that, although using a smaller threshold leads to more false alarms, the final recovery quality is better and those false alarms have small magnitudes and are easy to remove by thresholding during post-processing. That's why we chose to show their magnitudes, implying that we get easy-to-remove false alarms.

Epigraphic convergence is a standard tool to prove the convergence of a sequence of minimization problems. The definition of epigraphic convergence refers to Definition 7.1 and Proposition 7.2 in(Rockafellar & Wets, 2009).

Figures and Tables:

Figure 1: Justification of Theorems 1 and 2: comparision among LISTA variants.

θ k /γ k

Figure 2: Justification of Theorem 1 (noiseless case): parameters obtained by training satisfy (9).

Figure 3: Justification of Theorem 1 (noiseless case): Proportion of false positives vs true positives in x k (x * ).

Figure 4: Validation of Robust ALISTA. 5.3 VALIDATION OF ROBUST ALISTA

∀i / ∈ S by the definition of η θ (k) , i.e., support(x (k+1) ) ⊂ S By induction, (27) is proved.

PROOF OF CONCLUSION 2. Before proving Conclusion 2, let us analyze the relationship between D N conv and D N +D-1 cir .

Thus, linear operator T is full column rank. Thus, T T QT is positive definite, andG(w) = w T (T T QT )w is strongly convex. Then F 2D-1 cir (w) = G(w)+ι W normal (w) has only one minimizer, i.e., W 2D-1 cir involves only a unique element. Now we check the conditions of Propositions 7.32(c) and 7.33 in (Rockafellar & Wets, 2009) to apply them.

Figure 5 demonstrates that, with the fixed D, most of the random generated Ws have a σ (D, W) within the interval [0.25, 0.35]. Thus, the numerical results support our claim: with high probability, a random generated W satisfies min S∈S

Figure 5: Discussion of Definition 2: distribution of σ (D, W) on random generated Ws.

be the frequency domain signals, let DN cir = F D N cir F H be the frequency domain operator. The above equation is: b = DN cir x.

Fig.6visualizes the dictionary d (7 × 7 × 64) and the weights w ∈ W 13 cir , used in the convolutional A-LISTA simulation of Section 5.2. It is obtained by Algorithm 2 in Appendix E.2.

Table 3 shows w N conv → w * , i.e., the solution of the problem (22) converges to that of (24) as N increases, validating the second conclusion of Theorem 3. Visualized w * ∈ w 13 cir is displayed in Appendix F. Besides validating Theorem 3, we also present a real image denoising experiment to verify the effectiveness of Conv ALISTA. The detailed settings and results are presented in Appendix H. Validation of Conclusion 1 in Theorem 3. D = 7. w N cir ∈ W N cir and w * ∈ W 50 cir

Table content: <table/>

Applying Proposition 7.32(c), we have {F N conv } is eventually level bounded. If we arbitrarily pick a w N ∈ W N conv and let w cir be the unique point in W 2D-1 cir . Applying Proposition 7.33, we have w N → w cir . By Definition 4.1 in (Rockafellar & Wets, 2009), we obtain the convergence of the sequence of sets {W N conv }: lim N →∞ W N

Table content: <table><tr><td colspan="2">2D-1 cir nonempty closed convex domain W normal . and F N conv are all continuous functions defined on a</td></tr><tr><td>conv = W 2D-1 cir</td><td>.</td></tr><tr><td>D DISCUSSION OF DEFINITION 2 (11)</td><td/></tr></table>

Table1are the number of parameters to train, not the training complexity. The training complexity can be estimated by "Number of iterations × complexity of back-propagation", i.e., O(IBKN M ),where I is the number of iterations for training, B is the batch size , and K is the number of layers. Actually, Algorithm 1 (Stage 1) only takes a few seconds on an example of D : 250 × 500, while the training process (Stage 2) of, for example, ALISTA, takes around 0.1 hours.

Table content: <table/>

Table. 6. The performances are measured with NMSE in dB, defined in Section 5.1. From the table we can see that encoders do show better robustness when they have more layers, i.e. larger learning capacity.

Table content: <table><tr><td># Encoder Layers</td><td>σ of perturbations during testing 0.0001 0.001 0.01 0.015 0.02</td><td>0.025</td></tr><tr><td>4</td><td/><td/></tr></table>

Peer review (iclr2019):
 The papers studies neural network-based sparse signal recovery, and derives many new theoretical insights into the classical LISTA model. The authors proposed Analytic LISTA (ALISTA), where the weight matrix in LISTA is pre-computed with a data-free coherence minimization, followed by a separate data-driven learning step for merely (a very small number of) step-size and threshold parameters. Their theory is extensible to convolutional cases. The two-stage decomposed pipeline was shown to keep the optimal linear convergence proved in (Chen et al., 2018). Experiments observe that ALISTA has almost no performance loss compared to the much heavier parameterized LISTA, in contrast to the common wisdom that (brutal-force) “end-to-end” always outperforms stage-wise training. Their contributions thus manifest in both novel theory results, and the practical impacts of simplifying/accelerating LISTA training. Besides, they also proposed an interesting new strategy called Robust ALISTA to overcome the small perturbations on the encoding basis, which also benefits from this decomposed problems structure. The proofs and conclusions are mathematically correct to my best knowledge. I personally worked on similar sparse unfolding problems before so this work looks particularly novel and interesting to me. My intuition then was that, it should not be really necessary to use heavily parameterized networks to approximate a simple linear sparse coding form (LISTA idea). Similar accelerations could have been achieved with line search for something similar to steepest descent (also computational expensive, but need learn step-sizes only, and agnostic to input distribution). Correspondingly, there should exist a more elegant network solution with very light learnable weights. This work perfectly coincides with the intuition, providing very solid guidance on how a LISTA model could be built right. Given in recent three years, many application works rely on unfold-truncating techniques (compressive sensing, reconstruction, super resolution, image restoration, clustering…), I envision this paper to generate important impacts for practitioners pursuing those ideas. Additionally, I like Theorem 3 in Section 3.1, on the provable efficient approximation of general convolution using circular convolution. It could be useful for many other problems such as filter response matching. I therefore hold a very positive attitude towards this paper and support for its acceptance. Some questions I would like the authors to clarify & improve in revision: 1. Eqn (7) assumes noise-free case. The author stated “The zero-noise assumption is for simplicity of the proofs.” Could the authors elaborate which part of current theory/proof will fail in noisy case? If so, can it be overcome (even by less “simpler” way)? How about convolutional case, the same? Could the authors at least provide some empirical results for ALISTA’s performance under noise? 2. Section 5.3. It is unclear to me why Robust ALISTA has to work better than the data augmented ALISTA. Is it potentially because that in the data augmentation baseline, the training data volume is much amplified, and one ALISTA model might become underfitting? It would be interesting to create a larger-capacity ALISTA model (e.g., by increasing unfolded layer numbers), train it on the augmented data, and see if it can compare more favorably against Robust ALISTA? 3. The writeup is overall very good, mature, and easy to follow. But still, typos occur from time to time, showing a bit rush. For example, Section 5.1, “the x-axes denotes is the indices of layers” should remove “is”. Please make sure more proofreading will be done.