ALISTA: ANALYTIC WEIGHTS ARE AS GOOD AS LEARNED WEIGHTS IN LISTA
Year: 2019

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.

heading: 1 INTRODUCTION
text: 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 ∈ RM . We observe its noisy linear measurements:
b = M∑ m=1 dmx ∗ m + ε = Dx ∗ + ε, (1)
where b ∈ RN , D = [d1, · · · ,dM ] ∈ RN×M is the dictionary, and ε ∈ RN is additive Gaussian white noise. For simplicity, each column of D, named as a dictionary kernel, is normalized, that is, ‖dm‖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‖22 + λ‖x‖1 (2)
using iterative algorithms such as the iterative shrinkage thresholding algorithm (ISTA):
x(k+1) = ηλ/L ( x(k) + 1
L DT (b−Dx(k))
) , k = 0, 1, 2, . . . (3)
∗These authors contributed equally and are listed alphabetically.
where ηθ is the soft-thresholding function1 and L is usually taken as the largest eigenvalue of DTD.
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:
x(k+1) = ηθ(k)(W (k) 1 b + W (k) 2 x (k)), k = 0, 1, · · · ,K − 1. (4)
If we set W(k)1 ≡ 1LD T , W(k)2 ≡ I − 1LD TD, θ(k) ≡ 1Lλ, 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 Θ
Ex∗,b∼P ∥∥∥x(K)(Θ,b,x(0))− x∗∥∥∥2
2 . (5)
This problem is approximately solved over a training dataset {(x∗i ,bi)}Ni=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).

heading: 1.1 RELATED WORK
text: 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 `0- minimization. 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
1Soft- thresholding function is defined in a component-wise way: ηθ(x) = sign(x)max(0, |x| − θ)
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.

heading: 1.2 MOTIVATION AND CONTRIBUTIONS
text: 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.

heading: 2 ANALYTIC LISTA: CALCULATING WEIGHTS WITHOUT TRAINING
text: We theoretically analyze the LISTA-CPSS model defined in (Chen et al., 2018):
x(k+1) = ηθ(k) ( x(k) − (W(k))T (Dx(k) − b) ) , (6)
where W(k) = [w(k)1 , · · · ,w (k) M ] ∈ RN×M is a linear operator with the same dimensionality with D, x(k) = [ x
(k) 1 , · · · , x (k) M ] is the kth layer node. In (6), Θ = {W(k), θ(k)}k are parameters to train.
Model (6) can be derived from (4) with W(k)1 = (W (k))T ,W (k) 2 = I−W (k) 1 D. (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.
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.

heading: 2.1 RECOVERY ERROR UPPER BOUND
text: 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:
x∗ ∈ X (B, s) , { x∗ ∣∣∣|x∗i | ≤ B, ∀i, ‖x∗‖0 ≤ s}. (7)
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.
2A signal is s-sparse if it has no more than s non-zero entries.
Definition 1. Given D ∈ RN×M with each of its column normalized, we define the generalized mutual coherence:
µ̃(D) = inf W∈RN×M
(W:,i) TD:,i=1,1≤i≤M
{ max i 6=j
1≤i,j≤M
(W:,i) TD:,j } . (8)
Additionally, We define W(D) = { W ∈ RN×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) 6= ∅, 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, 22µ̃s−µ̃+1 ). Using them, define the parameters {W (k), θ(k)}:
W(k) = γ(k)W, θ(k) = γ(k)µ̃(D) sup x∗∈X (B,s)
{ ‖x(k)(x∗)− x∗‖1 } , (9)
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
support(x(k)(x∗)) ⊂ S, ‖x(k)(x∗)− x∗‖2 ≤ sB exp ( − k−1∑ τ=0 c(τ) ) , k = 1, 2, . . . (10)
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 kth 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.

heading: 2.2 RECOVERY ERROR LOWER BOUND
text: 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 PX . Let S = support(x∗). PX 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 −WTD ≈ 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 WTD is bounded away from the identity I. In Appendix D, we discuss the feasibility of (11).
Definition 2. Given D ∈ RN×M , s ≥ 2, σ̄min > 0, we define a set that W(k) are chosen from: W̄(D, s, σ̄min) = { W ∈ RN×M ∣∣∣σmin(I−(W:,S)TD:,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):
T = { {W(k) ∈ W̄(D, s, σ̄min), θ(k)}∞k=0 ∣∣∣support(x(k)(x∗)) ⊂ S, ∀x∗ ∈ X (B, s), ∀k} (12) 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
‖x(k)(x∗)− x∗‖2 ≥ ‖x∗‖2 exp(−c̄k), (13)
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.

heading: 2.3 ANALYTIC LISTA: LESS PARAMETERS TO LEARN
text: Following Theorems 1 and 2, we set W(k) = γ(k)W, where γ(k) is a scalar, and propose Tied LISTA (TiLISTA): x(k+1) = ηθ(k) ( x(k) − γ(k)WT (Dx(k) − b) ) , (14)
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 fromW(D) depends on D only. Hence we propose the analytic LISTA (ALISTA) that decomposes tied-LISTA into two stages:
x(k+1) = ηθ(k) ( x(k) − γ(k)W̃T (Dx(k) − b) ) , (15)
where W̃ is pre-computed by solving the following problem (Stage 1)3:
W̃ ∈ arg min W∈RN×M
∥∥WTD∥∥2 F , s.t. (W:,m)TD:,m = 1, ∀m = 1, 2, · · · ,M, (16)
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 WTD (a quadratic objective), over linear constraints. This is a standard convex quadratic program, which is easier to solve than to solve (8) directly.

heading: 3 CONVOLUTIONAL ANALYTIC LISTA
text: 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; GarciaCardona & 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:
b = M∑ m=1 dm ∗ x∗m + ε, (17)
where each dm is a dictionary kernel (or filter). {dm}Mm=1 is the dictionary of filters, M denotes the number of filters. {x∗m}Mm=1 is the set of coefficient maps that are assumed to have sparse structure,
3Some details and a complexity analysis of Stage 1 are discussed in Appendix E.1
and ∗ is the convolution operator. Now we consider 2D convolution and take4 b ∈ RN2 ,dm ∈ RD2 ,xm ∈ R(N+D−1) 2 . Equation (17) is pointwisely defined as5:
b(i, j) = D−1∑ k=0 D−1∑ l=0 M∑ m=1 dm(k, l)xm(i+ k, j + l) + ε(i, j), 0 ≤ i, j ≤ N − 1. (18)
We concatenate dms and xms: d = [d1, · · · ,dM ]T , x = [x1, · · · ,xM ]T , and rewrite (18) as:
b = M∑ m=1 DNconv,m(dm)xm + ε = D N conv(d)x + ε, (19)
where the matrix DNconv(d) = [D N conv,1(d1), · · · ,DNconv,M (dM )] ∈ RN 2×(N+D−1)2M , 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):
x(k+1)m = ηθ(k) ( x(k)m − ( w(k)m )′ ∗ ( M∑ m̄=1 dm̄ ∗ x(k)m̄ − b )) , m = 1, 2, · · · ,M, (20)
where {w(k)m }Mm=1 share the same sizes with {dm}Mm=1 and (·)′ means a 180 rotation of the filter (Chalasani et al., 2013). We concatenate the filters together: w(k) = [w(k)1 , · · · ,w (k) M ]
T ∈ RD2M . Parameters to train are Θ = {w(k), θ(k)}k.
Let WNconv(w (k)) be the matrix induced by dictionary w(k) with the same dimensionality as DNconv(d). Since convolution can be written as a matrix form (19), (20) is equivalent to
x(k+1) = ηθ(k) ( x(k) − (WNconv(w(k)))T (DNconv(d)x(k) − b) ) . (21)
Then by just substituting D,W(k) with DNconv(d),W N conv(w (k)) respectively, Theorems 1 and 2 can be applied to the convolutional LISTA.
Proposition 1. Let D = DNconv(d) and W(k) = WNconv(w(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.
Similar to the fully connected case (15), based on the results in Proposition 1, we should set w(k)m = γ (k) m w̃m, m = 1, 2, · · · ,M , where w̃ = [w̃1, · · · , w̃M ]T is chosen from
w̃ ∈ WNconv = arg min w∈RD 2M
wm·dm=1, 1≤m≤M
∥∥∥(WNconv(w))TDNconv(d)∥∥∥2 F . (22)
However, (22) is not as efficient to solve as (16). To see that, matrices DNconv(d) and W N conv(w) are both of size N2 × (N + D − 1)2M , the coherence matrix ( WNconv(w) )T DNconv(d) is thus of size (N +D−1)2M × (N +D−1)2M . 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)2M × (N +D − 1)2M ≈ 5× 1011.

heading: 3.1 CALCULATING CONVOLUTIONAL WEIGHTS ANALYTICALLY AND EFFICIENTLY
text: 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 dm(k, l)xm ( (i+k)modN , (j+l)modN ) +ε(i, j), 0 ≤ i, j ≤ N−1, (23)
4Here, 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 .
5Strictly 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.
where b ∈ RN2 ,dm ∈ RD 2 ,xm ∈ RN 2 . Similar to (18), we rewrite (23) in a compact way:
b = M∑ m=1 DNcir,m(dm)xm + ε = D N cir(d)x + ε,
where DNcir(d) : RN 2M → RN2 is a matrix depending on the signal size N and the dictionary d. Then the coherence minimization with the circular convolution is given by
WNcir = arg min w∈RD 2M
wm·dm=1, 1≤m≤M
∥∥∥(WNcir(w))TDNcir(d)∥∥∥2 F . (24)
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. WNcir =W 2D−1 cir ,∀N ≥ 2D − 1.
2. If at least one of the matrices {D2D−1cir,1 , · · · ,D 2D−1 cir,M } is non-singular, W 2D−1 cir involves
only a unique element. Furthermore,
lim N→∞
WNconv =W 2D−1 cir . (25)
The solution set WNcir 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 useW2D−1cir to approximateWNconv. In Appendix E.2, we introduce the algorithm details of solving (24). Based on Proposition 1 and Theorem 3, we obtain the convolutional ALISTA:
x(k+1)m = ηθ(k) ( x(k)m − γ(k)m ( w̃m )′ ∗ ( M∑ m̄=1 dm̄ ∗ x(k)m̄ − b )) , m = 1, 2, · · · ,M, (26)
where w̃ = [w̃1, · · · , w̃M ]T ∈ W2D−1cir and Θ = {{γ (k) m }m,k, {θ(k)}k} are the parameters to train. (26) is a simplified form, compared to the empirically unfolded CSC model recently proposed in (Sreter & Giryes, 2018)

heading: 4 ROBUST ALISTA TO MODEL PERTURBATION
text: 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. 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 mini-
mization 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.

heading: 5 NUMERICAL RESULTS
text: 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.

heading: 5.1 VALIDATION OF THEOREMS 1 AND 2 (ANALYTIC LISTA)
text: 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, Dij ∼ 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 pb = 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:
NMSEdB(x̂,x ∗) = 10 log10
( E‖x̂− x∗‖2/E‖x∗‖2 ) ,
where x∗ is the ground truth and x̂ is the estimated one.
6Our codes are uploaded to https://github.com/xchen-tamu/alista.
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 supx∗ ‖x(k)(x∗) − x∗‖1 (“supx∗” is taken over the test set). Figure 3 reports the average magnitude7 of the false positives and the true positives in xk(x∗) of ALISTA: the “true positives” curve draws the values of E{‖xkS(x∗)‖22/‖xk(x∗)‖22} w.r.t. k (the expectation is taken over the test set), while “false positives” for
E{‖xkSc(x∗)‖22/‖xk(x∗)‖22}. False positives take up small proportion over the positives, which supports the Theorem 1 conclusion that support(xk(x∗)) ⊂ S.

heading: 5.2 VALIDATION OF THEOREM 3 (CONVOLUTIONAL ANALYTIC LISTA)
text: 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 Ns.
In Table 2, we take wNcir ∈ W N cir, w ∗ ∈ W50cir (consider 50 as large enough) For this example,W N cir has only one element. Table 2 shows that wNcir = 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 wNconv ∈ W N conv and w ∗ ∈ w13cir, where W N conv also has only one element. Table 3 shows wNconv → 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∗ ∈ w13cir is displayed in Appendix F.
7The 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.
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.
Table 2: Validation of Conclusion 1 in Theorem 3. D = 7. wNcir ∈ W N cir and w ∗ ∈ W50cir.
‖wNcir −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
Table 3: Validation of Conclusion 2 in Theorem 3. D = 7. wNconv ∈ W N conv and w ∗ ∈ w13cir.

heading: 5.3 VALIDATION OF ROBUST ALISTA
text: 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, σ2max). 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 D̃s 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.

heading: 6 CONCLUSIONS AND FUTURE WORK
text: 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.

heading: A PROOF OF THEOREM 1
text: 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
support(x(k)) ⊂ S, ∀k. (27)
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) TD:,j(x (k) j − x ∗ j ) ) , ∀i /∈ S .
By (9), the thresholds are taken as θ(k) = µ̃γ(k) supx∗{‖x(k) − x∗‖1}. Also, since W ∈ W(D), we have |(W:,i)TD:,j | ≤ µ̃ for all j 6= 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)
TD:,j(x (k) j − x ∗ j ) ∣∣∣,
which implies x(k+1)i = 0,∀i /∈ S by the definition of ηθ(k) , i.e.,
support(x(k+1)) ⊂ S
By induction, (27) is proved.
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) TD:,S(x (k) S − x ∗ S) )
∈ x(k)i − γ (k)(W:,i) TD:,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 6= 0, [−1, 1] if x = 0.
The choice of W ∈ W(D) gives (W:,i)TD:,i = 1. Thus,
x (k) i − γ (k)(W:,i) TD:,S(x (k) S − x ∗ S)
=x (k) i − γ
(k) ∑
j∈S,j 6=i (W:,i)
TD:,j(x (k) j − x ∗ j )− γ(k)(x (k) i − x ∗ i )
=x∗i − γ(k) ∑
j∈S,j 6=i (W:,i)
TD:,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 6=i (W:,i) TD:,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 6=i γ(k) ∣∣∣(W:,i)TD:,j∣∣∣|x(k)j − x∗j |+ θ(k) + |1− γ(k)|∣∣x(k)i − x∗i ∣∣
≤µ̃γ(k) ∑
j∈S,j 6=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 6=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 ∈ Rn, 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.

heading: B PROOF OF THEOREM 2
text: Proof of Theorem 2. We fix D and sample a x∗ ∼ PX . If we can prove
P ( (13) does not hold ∣∣∣support(x∗) = S) ≤ |S|+ |S|, (28)
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− s3/2 − 2) ∑
2≤|S|≤s
P ( support(x∗) = S )
=1− s3/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
3s
)k} .
Let S = support(x∗). For x∗ ∈ X (k)( ), we consider two cases:
1. |x∗i | > ‖x∗‖2(σ̄min/3s)k, ∀i ∈ S.
2. |x∗i | ≤ ‖x∗‖2(σ̄min/3s)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 ) TD:,S(x (k−1) S − x ∗ S)− θ(k−1)sign(x (k) S ).
(29)
Subtracting x∗ from the two sides of (29), we obtain∥∥∥(I− (W(k−1):,S )TD:,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 )TD:,S)(x(k−1)S − x∗S)− θ(k−1)sign(x(k)S )∥∥∥ 2 ≤ ‖x∗‖2(σ̄min/3s)k. (30)
Multiplying both sides of (30) by (I− (W(k−1):,S )TD:,S)−1, we have
‖x(k−1)S − x ∗ S − θ(k−1)(I− (W (k−1) :,S ) TD:,S) −1sign(x(k)S )‖2
≤‖(I− (W(k−1):,S ) TD:,S) −1‖2 · ‖x∗‖2(σ̄min/3s)k ≤ ‖x∗‖2(σ̄min)k−13−ks,
where the last inequality is due to (11). Let x̃(k−1) denote the bias of x(k−1):
x̃(k−1) , θ(k−1)(I− (W(k−1):,S ) TD:,S) −1sign(x(k)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−13−ks}.
If case 2 holds, x∗ belongs to the following set: X̃ (k)( ) = { x∗ ∣∣∣|x∗i | ≤ ‖x∗‖2(σ̄min/3s)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)( ) andX (k−3)( ). Repeating the process, until dividingX (1)( ) into X̃ (1)( ) and X (0)( ).
By induction, we have
X (k)( ) ⊂ X̃ (k)( ) ∪ X̃ (k−1)( ) ∪ X̃ (k−2)( ) ∪ · · · ∪ X̃ (1)( ) ∪ X (0)( ), (31)
where the sets are defined as follows for all j = 0, 1, 2, · · · , k:
X̃ (k−j)( ) = { x∗ ∣∣∣|x∗i + x̃(k−j)i (x∗)| < ‖x∗‖2(σ̄min)k−j3−ks, for some i ∈ S.}, (32)
X (k−j)( ) = { x∗ ∣∣∣‖x(k−j)S (x∗)− x∗S − x̃(k−j)(x∗)‖2 ≤ ‖x∗‖2(σ̄min)k−j3−ks} (33)
and the bias is defined as following for all j = 0, 1, 2, · · · , k:
x̃(k−j)(x∗) = j∑ t=1 ( I− ( W (k−j+t−1) :,S )T D:,S )−t θ(k−j+t−1)sign ( x (k−j+t) S (x ∗) ) . (34)
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) = { (s1, s2, · · · , sn) ∣∣∣si ∈ {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−j3−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−j3−ks ∣∣∣support(x∗) = S) ≤ ∑ i∈S ∑ s(1)∈Si(|S|) ∑ s(2)∈Si(|S|) · · · ∑ s(j)∈Si(|S|) √ |S|B(σ̄min)k−j3−ks B
=|S|3j|S|( √ |S| ( (σ̄min) k−j3−ks ) ≤ |S|3/2(σ̄min)k−j3(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∗‖23−ks ∣∣∣support(x∗) = S) =
∑ s(1)∈Si(|S|) ∑ s(2)∈Si(|S|) · · · ∑ s(k)∈Si(|S|)
P ( ‖x∗ + x̃(0)(x∗)‖2 ≤ ‖x∗‖23−ks, sign(x(1)S ) = s (1), · · · , sign(x(k)S ) = s (k) ∣∣∣support(x∗) = S)
≤3k|S| ( ( 3−ks)|S| ) ≤ |S|.
Then we obtain P (x∗ ∈ X (k)( )|support(x∗) = S)
≤ k−1∑ j=0 |S|3/2(σ̄min)k−j3(j−k)|S| + |S| = k∑ j=1 |S|3/2(σ̄min)j3−j|S| + |S|
= |S|3/2 σ̄min3 −|S| 1− σ̄min3−|S| ( 1− (σ̄min3−|S|)k ) + |S| ≤ |S|3/2 + |S|.
Then (28) is proved.

heading: C PROOF OF THEOREM 3
text: 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 DNcir in detail.
The circular convolution (23) is equivalent with:
b(i, j) = N−1∑ k=0 N−1∑ l=0 M∑ m=1 DNcir(i, j; k, l,m)xm(k, l), 0 ≤ i, j ≤ N − 1,
where the circulant matrix is element-wise defined as:
DNcir(i, j; k, l,m) =
{ dm ( (k − i)modN , (l − j)modN ) , 0 ≤ (k − i)modN , (l − j)modN ≤ D − 1
0, others (35)
Similarly, the corresponding circulant matrix WNcir(i, j; k, l,m) of dictionary w is:
WNcir(i, j; k, l,m) =
{ wm ( (k − i)modN , (l − j)modN ) , 0 ≤ (k − i)modN , (l − j)modN ≤ D − 1
0, others (36)
As we defined in Section 3, b is a vector. With x = [x1, · · · ,xM ]T , x is a vector. Then the operator DNcir 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 DNcir(i, j; k, l,m) and W N cir(i, j; k, l,m): Bcoh = (D N cir)
TWNcir is elementwise defined by:
Bcoh(k1, l1,m1; k2, l2,m2) = N−1∑ i=0 N−1∑ j=0 DNcir(i, j; k1, l1,m1)W N cir(i, j; k2, l2,m2)
= ∑
i∈I(k1,k2) ∑ j∈J (l1,l2) dm1 ( I(i, k1), I(j, l1) ) wm2 ( I(i, k2), I(j, l2) ) .
where
I(k1, k2) = {i|0 ≤ i ≤ N − 1, 0 ≤ I(i, k1) ≤ D − 1, 0 ≤ I(i, k2) ≤ D − 1}, J (l1, l2) = {j|0 ≤ j ≤ N − 1, 0 ≤ I(j, l1) ≤ D − 1, 0 ≤ I(j, l2) ≤ D − 1}.
Lemma 2. Given N ≥ 2D − 1, it holds that:
(a) I(k1, k2) 6= ∅ if and only if “ 0 ≤ (k1−k2)modN ≤ D−1” or “ 0 < (k2−k1)modN ≤ D−1” holds.
(b) J (l1, l2) 6= ∅ if and only if “ 0 ≤ (l1 − l2)modN ≤ D− 1” or “ 0 < (l2 − l1)modN ≤ D− 1” holds.
Proof. Now we prove Conclusion (a). Firstly, we prove “if.” If 0 ≤ (k1 − k2)modN ≤ D − 1 and N ≥ 2D − 1, we have
I(k1, k2) = { (k1 − δ)modN ∣∣δ ∈ Z, (k1 − k2)modN ≤ δ ≤ D − 1} 6= ∅. (37)
If 0 < (k2 − k1)modN ≤ D − 1 and N ≥ 2D − 1, we have I(k1, k2) = { (k2 − δ)modN ∣∣δ ∈ Z, (k2 − k1)modN ≤ δ ≤ D − 1} 6= ∅. (38)
Secondly, we prove “only if.” If I(k1, k2) 6= ∅, we can select an i ∈ I(k1, k2). Let r1 = (k1 − i)modN and r2 = (k2 − i)modN . By the definition of I(k1, k2), we have 0 ≤ r1, r2 ≤ D − 1. Two cases should be considered here. Case 1: r1 ≥ r2. Since 0 ≤ r1 − r2 ≤ D − 1 ≤ N − 1, it holds that r1 − r2 = (r1 − r2)modN . Thus,
r1 − r2 = (r1 − r2)modN = ( (k1 − i)modN − (k2 − i)modN ) modN
= ( (k1 − i)− (k2 − i) ) modN =(k1 − k2)modN .
The equality “0 ≤ r1 − r2 ≤ D − 1” leads to the conclusion “0 ≤ (k1 − k2)modN ≤ D − 1”. In case 2 where r1 < r2, we can obtain 0 < (k2 − k1)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 k1, l1 and consider what values of k2, l2 give I(k1, k2) 6= ∅ and J (l1, l2) 6= ∅. Define four index sets given 0 ≤ k1, l1 ≤ N − 1:
K(k1) ={k|0 ≤ (k1 − k)modN ≤ D − 1} K̄(k1) ={k|0 < (k − k1)modN ≤ D − 1}
L(l1) ={l|0 ≤ (l1 − l)modN ≤ D − 1} L̄(l1) ={l|0 < (l − l1)modN ≤ D − 1}
Lemma 3. If N ≥ 2D − 1, we have:
(a) The cardinality of K(k1), K̄(k1): | K(k1)| = D, | K̄(k1)| = D − 1.
(b) K(k1) ∩ K̄(k1) = ∅.
(c) The cardinality of L(l1), L̄(l1): | L(l1)| = D, | L̄(l1)| = D − 1.
(d) L(l1) ∩ L̄(l1) = ∅.
Proof. Now we prove Conclusion (a). The set K(k1) can be equivalently written as
K(k1) = {(k1 − rk)modN |rk = 0, 1, · · · , D − 1} (39)
Let k(rk) = (k1 − rk)modN . We want to show that k(r1k) 6= k(r2k) as long as r1k 6= r2k. Without loss of generality, we assume 0 ≤ r1k < r2k ≤ D − 1. By the definition of modulo operation, There exist two integers q, q′ such that
k(r1k) = qN + k1 − r1k, k(r2k) = q′N + k1 − r2k.
Suppose k(r1k) = k(r 2 k). Taking the difference between the above two equations, we obtain r 2 k − r1k = (q ′ − q)N , i.e, N divides r2k − r1k. However, 0 ≤ r1k < r2k ≤ D − 1 implies 1 ≤ r2k − r1k ≤ D − 1 ≤ N − 1, which contradicts with “N dividing r2k − r1k.” Thus, it holds that k(r1k) 6= k(r2k). Then we have | K(k1)| = D. In the same way, we have
K̄(k1) = {(k1 + rk)modN |rk = 1, 2, · · · , D − 1} (40) and | K̄(k1)| = D − 1. Conclusion (a) is proved. Now we prove Conclusion (b). Suppose K(k1) ∩ K̄(k1) 6= ∅. Pick a k2 ∈ K(k1) ∩ K̄(k1). Let r3 = (k1−k2)modN and r4 = (k2−k1)modN . Then we have 0 ≤ r3 ≤ D−1 and 0 < r4 ≤ D−1. By the definition of modulo operation, There exist two integers q, q′ such that
k1 − k2 = qN + r3, k2 − k1 = q′N + r4 which imply
r3 + r4 + (q + q ′)N = 0.
However, 0 < r3 +r4 ≤ 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
L(l1) ={(l1 − rl)modN |rl = 0, 1, · · · , D − 1} (41) L̄(l1) ={(l1 + rl)modN |rl = 1, 2, · · · , D − 1} (42)
and | L(l1)| = D, | L̄(l1)| = 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 k1 ∈ {0, 1, · · · , N−1} and consider k2 ∈ K(k1). Let rk = (k1 − k2)modN . Then equation (37) implies that, for any i ∈ I(k1, k2), there exists a δ (rk ≤ δ ≤ D − 1) such that
I(i, k1) = ( k1 − (k1 − δ)modN ) modN = (δ)modN = δ,
I(i, k2) = ( k2 − (k1 − δ)modN ) modN = (δ − rk)modN = δ − rk. (43)
Now we consider another case for k2: k2 ∈ K̄(k1), rk = (k2 − k1)modN . Equation (38) implies that, for any i ∈ I(k1, k2), there exists a δ (rk ≤ δ ≤ D − 1) such that
I(i, k1) = ( k1 − (k2 − δ)modN ) modN
= (δ − rk)modN = δ − rk, I(i, k2) = ( k2 − (k2 − δ)modN ) modN = (δ)modN = δ. (44)
Similarly, for any l1 ∈ {0, 1, · · · , N − 1} and l2 ∈ L(l1), we denote rl = (l1 − l2)modN . For any j ∈ J (l1, l2), there exists a δ (rl ≤ δ ≤ D − 1) such that
I(j, l1) = ( l1 − (l1 − δ)modN ) modN = (δ)modN = δ,
I(j, l2) = ( l2 − (l1 − δ)modN ) modN = (δ − rl)modN = δ − rl. (45)
Another case for l2: l2 ∈ L̄(l1), rl = (l2 − l1)modN . For any j ∈ J (l1, l2), there exists a δ (rl ≤ δ ≤ D − 1) such that
I(j, l1) = ( l1 − (l2 − δ)modN ) modN
= (δ − rl)modN = δ − rl, I(j, l2) = ( l2 − (l2 − δ)modN ) modN = (δ)modN = δ. (46)
Now let us consider the following function. By results in Lemmas 2 and 3, we have
f(k1, l1,m1,m2) = N−1∑ k2=0 N−1∑ l2=0 ( Bcoh(k1, l1,m1; k2, l2,m2) )2 =f1 + f2 + f3 + f4,
where
f1 = ∑
k2∈K(k1) ∑ l2∈L(l1) ( Bcoh(k1, l1,m1; k2, l2,m2) )2 f2 =
∑ k2∈K̄(k1) ∑ l2∈L(l1) ( Bcoh(k1, l1,m1; k2, l2,m2) )2 f3 =
∑ k2∈K(k1) ∑ l2∈L̄(l1) ( Bcoh(k1, l1,m1; k2, l2,m2) )2 f4 =
∑ k2∈K̄(k1) ∑ l2∈L̄(l1) ( Bcoh(k1, l1,m1; k2, l2,m2) )2 .
Combining equations (39), (41), (43) and (45), we obtain
f1 = D−1∑ rk=0 D−1∑ rl=0 D−1∑ δk=rk D−1∑ δl=rl ( dm1(δk, δl)wm2(δk − rk, δl − rl) )2 .
Combining (40), (41), (44) and (45), we obtain
f2 = D−1∑ rk=1 D−1∑ rl=0 D−1∑ δk=rk D−1∑ δl=rl ( dm1(δk − rk, δl)wm2(δk, δl − rl) )2 .
Combining (39), (42), (43) and (46), we obtain
f3 = D−1∑ rk=0 D−1∑ rl=1 D−1∑ δk=rk D−1∑ δl=rl ( dm1(δk, δl − rl)wm2(δk − rk, δl) )2 .
Combining (40), (42), (44) and (46), we obtain
f4 = D−1∑ rk=1 D−1∑ rl=1 D−1∑ δk=rk D−1∑ δl=rl ( dm1(δk − rk, δl − rl)wm2(δk, δl) )2 .
By the above explicit formulas of fi, 1 ≤ i ≤ 4, we have f1, f2, f3, f4 are all independent of k1, l1 and N . They are only related with m1,m2 for fixed d and m. Thus, we are able to denote f(k1, l1,m1,m2) as f(m1,m2) for simplicity. Consequently,
1
N2 ‖(DNcir)TWNcir‖2F =
1
N2 N−1∑ k1=0 N−1∑ l1=0 N−1∑ k2=0 N−1∑ l2=0 M∑ m1=1 M∑ m2=1 ( Bcoh(k1, l1,m1; k2, l2,m2) )2 = 1
N2 N−1∑ k1=0 N−1∑ l1=0 M∑ m1=1 M∑ m2=1 f(k1, l1,m1,m2)
= 1
N2 N−1∑ k1=0 N−1∑ l1=0 M∑ m1=1 M∑ m2=1 f(m1,m2)
= 1
N2 ·N2 · M∑ m1=1 M∑ m2=1 f(m1,m2) = M∑ m1=1 M∑ m2=1 f(m1,m2)
Thus, 1N2 ‖(D N cir) TWNcir‖2F is dependent of N :
1
N2 ‖(DNcir)TWNcir‖2F =
1
(2D − 1)2 ‖(D2D−1cir ) TW2D−1cir ‖ 2 F , ∀N ≥ 2D − 1, (47)
which impliesWNcir =W 2D−1 cir ,∀N ≥ 2D − 1.
C.2 PROOF OF CONCLUSION 2.
Before proving Conclusion 2, let us analyze the relationship between DNconv and D N+D−1 cir .
Similar to Dcir, we use (i, j) as the row index and (k, l,m) as the column index of Dconv. For 0 ≤ i, j ≤ N − 1, 1 ≤ m ≤M ,
DN+D−1cir (i, j; k, l,m) = D N conv(i, j; k, l,m) = { dm(k − i, l − j), 0 ≤ k − i, l − j ≤ D − 1 0, k, l taken as others
(48) Matrix DN+D−1cir is of dimension (N +D− 1)2× (N +D− 1)2M , where 0 ≤ i, j ≤ N +D− 2; matrix DNconv is of dimension (N)
2 × (N +D − 1)2M , where 0 ≤ i, j ≤ N − 1. Thus, DNconv is a block in DN+D−1cir , i.e.,
DN+D−1cir = [ DNconv ∆ND ] .
The matrix ∆ND is of dimension ((N +D − 1)2 −N2)× (N +D − 1)2M : ∆ND = [ DN+D−1cir (i, j; :, :, :) ] , (i, j) ∈ I∆
where I∆ = I1 ∪I2 ∪I3 I1 ={(i, j)|N ≤ i ≤ N +D − 2, 0 ≤ j ≤ N − 1} I2 ={(i, j)|0 ≤ i ≤ N − 1, N ≤ j ≤ N +D − 2} I3 ={(i, j)|N ≤ i ≤ N +D − 2, N ≤ j ≤ N +D − 2}.
Similarly,
WN+D−1cir = [ WNconv ∆NW ] , ∆NW = [ WN+D−1cir (i, j; :, :, :) ] , (i, j) ∈ I∆ .
Then, (DN+D−1cir ) TWN+D−1cir = (D N conv) TWNconv + (∆ N D) T∆NW. (49)
Lemma 4. For any (i, j) ∈ I∆, one has
‖DN+D−1cir (i, j; :, :, :)‖ 2 2 =‖d‖22, (50) ‖WN+D−1cir (i, j; :, :, :)‖ 2 2 =‖w‖22. (51)
Proof. Equation (35) implies that, for (i, j) ∈ I1, 1 ≤ m ≤M ,
DN+D−1cir (i, j; k, l,m) =  dm(k − i, l − j), i ≤ k ≤ N +D − 2, j ≤ l ≤ j +D − 1 dm(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) ∈ I1,
‖DN+D−1cir (i, j; :, :, :)‖ 2 2 = N+D−2∑ k=0 N+D−2∑ l=0 M∑ m=1 ∣∣DN+D−1cir (i, j; k, l,m)∣∣2 = ‖d‖22 Similarly,
‖DN+D−1cir (i, j; :, :, :)‖ 2 2 = ‖d‖22, (i, j) ∈ I2 ∪I3 .
Equation (50) is proved. With the same argument, equation (51) is also proved.
Lemma 5. If N ≥ 2D − 1, we have ‖(∆ND)T∆NW‖2F ≤ ( 2N(D − 1) + (D − 1)2 ) (2D − 1)2‖d‖22‖w‖22. (52)
Proof. For simplicity, we denote two row vectors:
di,j , D N+D−1 cir (i, j; :, :, :) ∈ R 1×(N+D−1)2M
wi,j , W N+D−1 cir (i, j; :, :, :) ∈ R 1×(N+D−1)2M
Then, ‖(∆ND)T∆NW‖2F = ∥∥∥∥ ∑
(i,j)∈I∆
dTi,jwi,j ∥∥∥∥2 F = ∑
(i1,j1)∈I∆ ∑ (i2,j2)∈I∆ 〈 dTi1,j1wi1,j1 ,d T i2,j2wi2,j2 〉 F ,
where〈 dTi1,j1wi1,j1 ,d T i2,j2wi2,j2 〉 F = trace ( wTi1,j1di1,j1d T i2,j2wi2,j2 ) = (di1,j1d T i2,j2) · (wi1,j1w T i2,j2).
Since
di1,j1d T i2,j2 = N−1∑ k=0 N−1∑ l=0 M∑ m=1 = DN+D−1cir (i1, j1; k, l,m)D N+D−1 cir (i2, j2; k, l,m),
with the same argument in Lemma 2, we have: di1,j1d T i2,j2 6= 0 implies
i2 ∈ I ′∆ , {i|0 ≤ (i1 − i)mod(N+D−1) ≤ D − 1 or 0 ≤ (i− i1)mod(N+D−1) ≤ D − 1}
j2 ∈ J ′∆ , {j|0 ≤ (j1 − j)mod(N+D−1) ≤ D − 1 or 0 ≤ (j − j1)mod(N+D−1) ≤ D − 1}
Then ‖(∆ND)T∆NW‖2F = ∑ (i1,j1)∈I∆ ∑ i2∈I′∆ ∑ j2∈J ′∆ (di1,j1d T i2,j2) · (wi1,j1w T i2,j2)
≤ ∑
(i1,j1)∈I∆ ∑ i2∈I′∆ ∑ j2∈J ′∆ ‖d‖22‖w‖22
=| I∆ | · | I ′∆ | · | J ′ ∆ | · ‖d‖22‖w‖22 = ( 2N(D − 1) + (D − 1)2 ) (2D − 1)2‖d‖22‖w‖22,
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 Wnormal = { w ∈ RD 2M ∣∣∣wm · dm = 1, ∀m = 1, · · · ,M}. (53)
Since d ∈ Wnormal, the set is nonempty:
Wnormal 6= ∅. (54)
Define functions FNconv : RD 2M → R, FNcir : RD 2M → R.
FNconv(w) = 1
N +D − 1 ∥∥∥(DNconv(d))TWNconv(w)∥∥∥ F + ιWnormal(w)
FNcir(w) = 1
N ∥∥∥(DNcir(d))TWNcir(w)∥∥∥ F + ιWnormal(w)
By the definitions ofWNconv,W N cir, we have
WNconv = arg min w FNconv(w), W N cir = arg min w FNcir(w)
Step 1: Proving FNconv(w) uniformly converges to F 2D−1 cir (w) on X ∩Wnormal for any compact set X ⊂ RD2M . We arbitrarily choose such a compact set X . Based on (47), (49) and (52), one has, for all w ∈ X ∩Wnormal,
|FNconv(w)− F 2D−1cir (w)| =|F N conv(w)− FN+D−1cir (w)|
= 1
N +D − 1 ∣∣∣∣∥∥∥(DN+D−1cir )TWN+D−1cir ∥∥∥ F − ∥∥∥(DNconv)TWNconv∥∥∥ F ∣∣∣∣ ≤ 1 N +D − 1
∥∥(∆ND)T∆NW∥∥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
|FNconv(w)− F 2D−1cir (w)| ≤ B√ N sup w∈X∩Wnormal ‖w‖, ∀ w ∈ X ∩Wnormal . (55)
Step 2: Proving FNconv(w) epigraphically converges8 to F 2D−1 cir (w).
We want to show, at each point w it holds that
lim inf N→∞
FNconv(w N ) ≥ F 2D−1cir (w) for every sequence w N → w (56)
lim sup N→∞
FNconv(w N ) ≤ F 2D−1cir (w) for some sequence w N → w (57)
Firstly, we prove (56). We arbitrarily pick a sequence {wN}∞N=0 such that wN → w.
If w /∈ Wnormal, F 2D−1cir (w) = +∞. Since Wnormal is a closed set, there exists a N+ such that wN /∈ Wnormal for all N ≥ N+. Thus, one has FNconv(wN ) = +∞ for all N ≥ N+, i.e.,
lim inf N→∞
FNconv(w N ) = F 2D−1cir (w) = +∞.
If w ∈ Wnormal, two cases should be considered. The first case is that any subsequences of {wN}∞N=0 are not kept within Wnormal, i.e., there exists a N+ such that wN /∈ Wnormal for all N ≥ N+. Then we have
lim inf N→∞
FNconv(w N ) = +∞ > F 2D−1cir (w).
The second case is that there exists a subsequence {wNk}∞k=0 ⊂ {wN}∞N=0 such that
wNk ∈ Wnormal, ∀k = 0, 1, 2, · · · . Since wN converges to w, any subsequences should be Cauchy. Given any Cauchy sequence {wNk}∞k=0 in finite dimensional Euclidean space, there exists a compact set X such that
wNk ∈ X, ∀k = 0, 1, 2, · · · Let B′ = supw∈X∩Wnormal ‖w‖. By (55), we obtain
|FNkconv(wNk)− F 2D−1cir (w)| ≤|F Nk conv(w Nk)− F 2D−1cir (w Nk)|+ |F 2D−1cir (w Nk)− F 2D−1cir (w)|
≤ BB ′
√ Nk + |F 2D−1cir (w Nk)− F 2D−1cir (w)|.
8Epigraphic 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).
For any > 0, by the continuity of F 2D−1cir , we are able to find a K > 0 such that |F 2D−1 cir (w Nk)− F 2D−1cir (w)| < for all k ≥ K. Pick a K ′ such that NK′ ≥ (BB′/ )2. Then, for all k ≥ max(K,K ′), we have |FNkconv(wNk)− F 2D−1cir (w)| < 2 , i.e.,
lim k→∞
FNkconv(w Nk) = F 2D−1cir (w).
The above conclusion holds for all subsequences {wNk}∞k=0 ⊂ Wnormal. F 2D−1 cir (w) is an accumulation point of {FNconv(wN )}∞N=0. All the other accumulation points of {FNconv(wN )}∞N=0 must be +∞ because FNconv(w) = F 2D−1cir (w) = +∞ for all w /∈ Wnormal. Thus,
lim inf N→∞
FNconv(w N ) = F 2D−1cir (w) < +∞.
Secondly, we prove (57). We set wN = w for all N = 0, 1, 2, · · · . Then (57) is a direct result of (55).
Step 3: proving (25). Define G(w) = ∥∥(D2D−1cir )TW2D−1cir ∥∥2F .
We want to show that G(w) is strongly convex.
Let w̃i ∈ R(2D−1) 2 be the ith column of W2D−1cir , i.e., W2D−1cir = [ w̃1, w̃2, · · · , w̃(2D−1)2M ] Then
G(w) = (2D−1)2M∑ i=1 (w̃i) T ( D2D−1cir (D 2D−1 cir ) T ) w̃i.
Let w̃ ∈ R(2D−1)4M vectorize W2D−1cir , i.e.,
w̃ = [ (w̃1) T , (w̃2) T , · · · , (w̃(2D−1)2M )T ]T .
Then G(w) can be written as a quadratic form of w̃:
G(w) = w̃TQw̃,
where
Q =  ( D2D−1cir (D 2D−1 cir ) T ) · · · (
D2D−1cir (D 2D−1 cir )
T ) 
︸ ︷︷ ︸ totally (2D−1)2M diagonal blocks
.
As long as at least one of the matrices {D2D−1cir,0 , · · · ,D 2D−1 cir,M−1} is non-singular, D 2D−1 cir is full row rank, which implies that D2D−1cir (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̃ = Tw.
It’s trivial that ‖w̃‖22 = 0 implies ‖W2D−1cir ‖2F = 0. By the definition of W 2D−1 cir , ‖W 2D−1 cir ‖2F = 0 implies ‖w‖22 = 0. Thus, linear operator T is full column rank. Thus, TTQT is positive definite, and
G(w) = wT (TTQT )w is strongly convex. ThenF 2D−1cir (w) = √ G(w)+ιWnormal(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.
1. FNconv e−→ F 2D−1cir . This is proved in Step 2.
2. F 2D−1cir is level bounded. Since G(w) is strongly convex, F 2D−1 cir (w) =
√ G(w) +
ιWnormal(w) must be level bounded.
3. F 2D−1cir 6≡ +∞. SinceWnormal is nonempty (54), dom F 2D−1 cir 6= ∅, F 2D−1 cir is not con-
stantly +∞.
4. All the level set of FNconv are connected. This can be derived from the convexity of F N conv.
5. F 2D−1cir and F N conv are all lower semi-continuous and proper. This condition follows from
the fact that the functions F 2D−1cir and F N conv are all continuous functions defined on a nonempty closed convex domainWnormal.
Applying Proposition 7.32(c), we have {FNconv} is eventually level bounded. If we arbitrarily pick a wN ∈ WNconv and let wcir be the unique point in W 2D−1 cir . Applying Proposition 7.33, we have wN → wcir. By Definition 4.1 in (Rockafellar & Wets, 2009), we obtain the convergence of the sequence of sets {WNconv}: limN→∞W N conv =W 2D−1 cir .

heading: D DISCUSSION OF DEFINITION 2 (11)
text: 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)TD:,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. There are totally(
500
50
) + ( 500
49
) + · · ·+ ( 500
2 ) possible Ss satisfying 2 ≤ |S | ≤ s. It’s impossible to check (11) on all possible Ss. 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)TD:,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.
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
{ σmin ( I− (W:,S)TD:,S )} ≥ σ̄min > 0,
that is, W ∈ W̄(D, s, σ̄min).

heading: E EFFICIENT ALGORITHM TO CALCULATE ANALYTIC WEIGHTS
text: E.1 AN EFFICIENT ALGORITHM TO SOLVE (16)
In this section, we introduce an algorithm to solve (16) (we copy (16) below to facilitate reading):
min W∈RN×M
∥∥WTD∥∥2 F , s.t. (W:,m)TD:,m = 1, ∀m = 1, 2, · · · ,M,
By the definition of the Frobenius norm, it holds that
‖WTD‖2F = ‖(WTD)T ‖2F = ‖DTW‖2F . (58)
Thus, the above problem is equivalent with
min W∈RN×M
∥∥DTW∥∥2 F , s.t. (D:,m)TW:,m = 1, ∀m = 1, 2, · · · ,M.
We apply projected gradient descent (PGD) to solve the above problem. The gradient of ‖DTW‖2F is ∇‖DTW‖2F = DDTW. Denote the set by
W = {W ∈ RN×M |(D:,m)TW:,m = 1, ∀m = 1, 2, · · · ,M.}
Then the projection ontoW can be calculated by ProjW(W) = W + ∆W, ∆W = [ (1− (D:,1)TW:,1)W:,1, · · · , (1− (D:,M )TW:,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 ∈ RN×M . Initialize: Let W0 = D.
1 for j = 0, 1, 2, . . . until convergence do 2 Update W by Wj+1 = ProjW ( Wj − ηD(D)TWj ) .
3 end Output: WJ , where J is the last iterate.
In each step, calculating the gradient has the complexity of O(N2M) because DDT can be precomputed. Calculating the projection takes O(NM) 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/ )N2M). We should note that the bounds given in Table 1 are 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(IBKNM),where I is the number of iterations for training,B is the batch size , andK 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.
E.2 AN EFFICIENT ALGORITHM TO SOLVE (24)
In this section, we introduce an algorithm to solve (24) (we copy (24) below to facilitate reading):
min w∈RD 2M
wm·dm=1, 1≤m≤M
∥∥∥(WNcir(w))TDNcir(d)∥∥∥2 F .
Similarly, by (58), the above problem is equivalent with
min w∈RD 2M
dm·wm=1, 1≤m≤M
∥∥∥(DNcir(d))TWNcir(w)∥∥∥2 F . (59)
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 DNcir(d),W N cir(w) in the frequency domain. To simplify the notation, we denote the operators as DNcir and W N cir respectively. Let F be the FFT operator. Thus, b = DNcirx is equivalent with F b = F DNcir F
H F x. Let b̂ = F b, x̂ = F x be the frequency domain signals, let D̂Ncir = F DNcir F
H be the frequency domain operator. The above equation is:
b̂ = D̂Ncirx̂.
The frequency domain operator D̂Ncir 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 toN×N and do FFT: d̂m = FFT ( zeropad(dm, N −D) ) , then the above operator can be explicitly written as:
b̂ = M∑ m=1 d̂m x̂m,
where ·̄ means complex conjugate. This is due to DNcir 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.
Further, since∥∥∥(D̂Ncir)HŴNcir∥∥∥2 F = ∥∥∥(F DNcir FH )H FWNcir FH ∥∥∥2 F = ∥∥∥F (DNcir)TWNcir FH ∥∥∥2 F
= ∥∥∥(DNcir)TWNcir∥∥∥2
F ,
problem (59) is equivalent with
min w∈RD 2M
dm·wm=1, 1≤m≤M
∥∥∥(D̂Ncir)HŴNcir∥∥∥2 F ,
which can be efficiently solved by the frequency domain ISTA in (Liu et al., 2017). The details are outlined in Algorithm 2.
F VISUALIZATION OF THE ANALYTIC CONVOLUTIONAL WEIGHTS
Fig. 6 visualizes the dictionary d (7× 7× 64) and the weights w̃ ∈ W13cir, used in the convolutional A-LISTA simulation of Section 5.2. It is obtained by Algorithm 2 in Appendix E.2.

heading: G ALGORITHM DETAILS OF TRAINING ROBUST ALISTA
text: G.1 MODEL ARCHITECTURE
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:
arg min W∈RN×M
∥∥Q (DTW − IM )∥∥2F , (60)
Algorithm 2: Frequency-domain ISTA for solving (24)
Input: Dictionary d = [d1, · · · ,dM ]T , dm ∈ RD 2
, m = 1, 2, · · · ,M . Initialize: Let w0 = d.
1 for j = 0, 1, 2, . . . until convergence do 2 Zeropad and FFT:
ŵjm = FFT ( zeropad ( wjm, N −D )) , m = 1, · · · ,M.
3 Compute frequency domain gradient:
(∇f)m = ( M∑ m=1 d̂m ¯̂dm ) ŵjm, m = 1, · · · ,M,
where ·̄ represents the conjugate of a complex number. 4 Compute the next iterate:
wj+1m = ProjWnormal ( IFFT ( ŵjm − η(∇f)m )) , m = 1, · · · ,M,
where the setWnormal is defined in (53). 5 end
Output: wJ , where J is the last iterate.
where is the Hadamard product and Q is a weight matrix that put more penalty on errors on diagonals, because entries on the diagonal will be far smaller than off-diagonal. The above relaxed coherence minimization can be solved using the gradient descent algorithm:
W(k+1) = W(k) − γ(k)D(Q2 (DTW(k) − IM )). (61)
By unfolding (61) and truncate to K steps, and considering the γ(k)DT outside the residual as learnable parameters B, we will have a deep neural network W = E(D) as a coherence minimizer. We call it a Stage 1 encoder as it encodes a dictionary D into a weight matrix, that can be used in the Stage 2 of ALISTA, refered as a decoder. One layer of this model is shown in Fig. 7(a).
The illustration of the whole feed-forward robust model is shown in Fig. 7(b). The two parts, the encoder and the decoder, can be jointly trained to gain the most from data-driven learning. We further adopt pre-training and curriculum learning to stabilize training, as to be discussed below.
G.2 MODEL TRAINING
To stabilize the training process, we train the model in two stages: the pre-training stage and curriculum (joint) training stage.
Pre-Training Stage. We first pre-train the encoder and the decoder individually. The pre-training of the decoder, e.g., ALISTA, follows the standard training procedure in Section 5.1, without bothering the perturbations of D. On the other hand, the encoder will always see perturbed dictionaries D̃ = D + εD, where εD’s entries are sampled from i.i.d. normal distribution with zero mean and σ2pre variance, and update its weight to minimize loss function defined by (60). The σpre is a hyperparameter that we manually select for the pre-training stage, with a default value of 0.01. We use an exponentially decaying learning rate for encoder pre-training with an initial value αpre = 10−4.
Curriculum (Joint) Training Stage. After the pre-training stage, we concatenate these two parts and do joint training. However, a direct end-to-end tuning was observed to cause much instability, due to the randomness in weights. Inspired by the curriculum learning technique, we first perturb the dictionaries with smaller standard deviations and gradually increase the perturbation level during training. Specifically, starting from a small standard deviation σt = σ0, the curriculum joint training procedure repeats the routine below:
• First uniformly sample a batch of standard deviations {σi}BDi=1 from [0, σt], where BD is the batch size for perturbations of the original dictionary D. Use the sampled standard deviations to sample BD perturbations and apply to D and then normalized them to get {D̃}BDi=1.
• Sample a batch of sparse codes {xj}Bxj=1 from a pre-defined Gaussian-Bernoulli distribution; the supports of the sparse codes are decided i.i.d. by a Bernoulli distribution to have around 10% non-zero entries; and the magnitudes are sampled from i.i.d. standard Gaussian. Bx is the batch size for sparse codes in training.
• Measure yi,j = D̃ixj . Then (yi,j ,xj , D̃i) forms a tripelet of training sample. Note that only Robust ALISTA needs D̃i as part of the training samples.
• Feed in the data and update the encoder and decoder with learning rates αe and αd, using the Adam Optimizer, respectively.
• Increase σt to the next larger value, after C training batches.
• Repeat the above steps, until the value of σt exceeds the pre-defined σmax, that represents the maximal standard deviation to sample the dictionary perturbation.
In the experiment, we have BD = 4, Bx = 16, C = 50000, αe = 10−6, αd = 10−4 and σmax ∈ {0.02, 0.03} and σt is obtained by linearly interpolating between [0, σmax] for L − 1 times. We choose L = 5; hence σt takes i5σmax, i = 1, 2, . . . , 5 in order.

heading: H RESULTS OF NATURAL IMAGE DENOISING USING CONV ALISTA
text: The natural image denoising experiment is conducted on the same BSD 500 dataset using the 400- image training set, 50-image validation set and 50-image test set. We convert them all to grayscale, and then add σ = 20 Gaussian i.i.d. noise. We train both Conv LISTA (i.e., model (20)) and Conv ALISTA (i.e., model (26)). Both networks have 5 layers, with the same dictionary D obtained from the training set by solving (24). We reconstruct the denoised images using by convolving the learned feature maps with the original dictionary D. The mean-square-error (MSE) between denoised and clean images are adopted as the network training loss, as inspired by (Zhou et al., 2018).
Six popular benchmark images (adding σ = 20 noise) are tested and reported in Table 4. The APSNR denotes the average PSNR over all images and the A-Times represents the average inference time (in seconds) for denoising one image. We compare Conv LISTA and Conv ALISTA, as well as the classical KSVD denoising algorithm (Elad & Aharon, 2006) and the recent CSC denoising algorithm with gradient regularization (CSC-GR) (Wohlberg, 2018). The results show that Conv LISTA and Conv ALISTA (without heavy tuning done for their optimal performance) can perform comparably with KSVD and outperforms CSC-GR, but with tremendously faster inference speeds than KSVD/CSC-GR. More importantly, Conv LISTA and Conv ALISTA only have marginal performance differences, validating again the analytic weights in convolutional cases.

heading: I RESULTS OF ABLATION STUDIES IN ROBUSTNESS EXPERIMENTS
text: As one anonymous reviewer kindly pointed out, Robust ALISTA has larger parameter space over TiLISTA and ALISTA trained. Therefore, they suggested that we increased the number of layers in TiLISTA, ALISTA and the baseline model LISTA-CPSS in (Chen et al., 2018) to see if their performance in this above evaluation setting can be improved in that way. Note that LISTA-CPSS has tens of layers, hence actually containing more parameters than robust ALISTA. In addition, the reviewers also suggested a set of ablation studies to investigate whether more layers in the encoder can endorse the model better adaptivity to higher level perturbations. We conduct the suggested experiments and present the results in this section.
I.1 NUMBER OF LAYERS IN TILISTA, ALISTA AND LISTA-CPSS
As the reviewers pointed out, the comparison we present in Section. 5.3 might be unfair because robust ALISTA contains much more parameters (because it contains a 4-layer encoder) comparing to ALISTA, which only learns two series of scalars, and TiLISTA which has just one more matrix weight than ALISTA. Therefore, we add the following experiments to consolidate our claim on the effectiveness of the robust ALISTA model:
• we increase the number of layers of TiLISTA and ALISTA which are then trained in the same data augmentation setting as we do in Section. 5.3, to see if they could yield comparitive robustness against dictionary perturbations;
• we also compare the robust ALISTA with the baseline LISTA-CPSS model in Chen et al. (2018), which contains tens of layers of independent weight matrices, thus having even more parameters than robust ALISTA. This comparison could consolidate our claim that the outstanding adaptiveness to dictionary perturbations of robust ALISTA is brought by its encoder-decoder structure rather than its learning capacity alone.
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.
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.
Results are shown in 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.

Review summaries:
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 paper raises many important questions about unrolled iterative optimization algorithms, and answers many questions for the case of iterative soft thresholding algorithm (ISTA, and learned variant LISTA). The authors demonstrate that a major simplification is available for the learned network: instead of learning a matrix for each layer, or even a single (potentially large) matrix, one may obtain the matrix analytically and learn only a series of scalars. These simplifications are not only practically useful but allow for theoretical analysis in the context of optimization theory. On top of this seminal contribution, the results are extended to the convolutional-LISTA setting. Finally, yet another fascinating result is presented, namely that the analytic weights can  be determined from a Gaussian-perturbed version of the dictionary. Experimental validation of all results is presented.