Title: Understanding Trainable Sparse Coding via matrix factorization

Abstract: Sparse coding is a core building block in many data analysis and machine learning pipelines. Typically it is solved by relying on generic optimization techniques, such as the Iterative Soft Thresholding Algorithm and its accelerated version (ISTA, FISTA). These methods are optimal in the class of first-order methods for non-smooth, convex functions. However, they do not exploit the particular structure of the problem at hand nor the input data distribution. An acceleration using neural networks, coined LISTA, was proposed in Gregor & Le Cun (2010), which showed empirically that one could achieve high quality estimates with few iterations by modifying the parameters of the proximal splitting appropriately.
In this paper we study the reasons for such acceleration. Our mathematical analysis reveals that it is related to a specific matrix factorization of the Gram kernel of the dictionary, which attempts to nearly diagonalise the kernel with a basis that produces a small perturbation of the 1 ball. When this factorization succeeds, we prove that the resulting splitting algorithm enjoys an improved convergence bound with respect to the non-adaptive version. Moreover, our analysis also shows that conditions for acceleration occur mostly at the beginning of the iterative process, consistent with numerical experiments. We further validate our analysis by showing that on dictionaries where this factorization does not exist, adaptive acceleration fails.

1. Introduction:

Feature selection is a crucial point in high dimensional data analysis. Different techniques have been developed to tackle this problem efficiently, and amongst them sparsity has emerged as a leading paradigm. In statistics, the LASSO estimator (Tibshirani, 1996) provides a reliable way to select features and has been extensively studied in the last two decades (Hastie et al. (2015) and references therein). In machine learning and signal processing, sparse coding has made its way into several modern architectures, including large scale computer vision (Coates & Ng, 2011) and biologically inspired models (Cadieu & Olshausen, 2012) . Also, Dictionary learning is a generic unsupervised learning method to perform nonlinear dimensionality reduction with efficient computational complexity (Mairal et al., 2009) . All these techniques heavily rely on the resolution of 1 -regularized least squares.

The 1 -sparse coding problem is defined as solving, for a given input x ∈ R n and dictionary D ∈ R n×m , the following problem:

EQUATION

)

This problem is convex and can therefore be solved using convex optimization machinery. Proximal splitting methods (Beck & Teboulle, 2009 ) alternate between the minimization of the smooth and differentiable part using the gradient information and the minimization of the non-differentiable part using a proximal operator (Combettes & Bauschke, 2011) . These methods can also be accelerated by considering a momentum term, as it is done in FISTA (Beck & Teboulle, 2009; Nesterov, 2005) . Coordinate descent (Friedman et al., 2007; Osher & Li, 2009) leverages the closed formula that can be derived for optimizing the problem (1) for one coordinate z i given that all the other are fixed. At each step of the algorithm, one coordinate is updated to its optimal value, which yields an inexpensive scheme to perform each step. The choice of the coordinate to update at each step is critical for the performance of the optimization procedure. Least Angle Regression (LARS) (Hesterberg et al., 2008) is another method that computes the whole LASSO regularization path. These algorithms all provide an optimization procedure that leverages the local properties of the cost function iteratively. They can be shown to be optimal among the class of first-order methods for generic convex, non-smooth functions (Bubeck, 2014) .

But all these results are given in the worst case and do not use the distribution of the considered problem. One can thus wonder whether a more efficient algorithm to solve (1) exists for a fixed dictionary D and generic input x drawn from a certain input data distribution. In Gregor & Le Cun (2010) , the authors introduced LISTA, a trained version of ISTA that adapts the parameters of the proximal splitting algorithm to approximate the solution of the LASSO using a finite number of steps. This method exploits the common structure of the problem to learn a better transform than the generic ISTA step. As ISTA is composed of a succession of linear operations and piecewise non linearities, the authors use the neural network framework and the backpropagation to derive an efficient procedure solving the LASSO problem. In Sprechmann et al. (2012) , the authors extended LISTA to more generic sparse coding scenarios and showed that adaptive acceleration is possible under general input distributions and sparsity conditions.

In this paper, we are interested in the following question: Given a finite computational budget, what is the optimum estimator of the sparse coding? This question belongs to the general topic of computational tradeoffs in statistical inference. Randomized sketches (Alaoui & Mahoney, 2015; Yang et al., 2015) reduce the size of convex problems by projecting expensive kernel operators into random subspaces, and reveal a tradeoff between computational efficiency and statistical accuracy. Agarwal (2012) provides several theoretical results on perfoming inference under various computational constraints, and Chandrasekaran & Jordan (2013) considers a hierarchy of convex relaxations that provide practical tradeoffs between accuracy and computational cost. More recently, Oymak et al. (2015) provides sharp time-data tradeoffs in the context of linear inverse problems, showing the existence of a phase transition between the number of measurements and the convergence rate of the resulting recovery optimization algorithm. Giryes et al. (2016) builds on this result to produce an analysis of LISTA that describes acceleration in conditions where the iterative procedure has linear convergence rate. Finally, Xin et al. (2016) also studies the capabilities of Deep Neural networks at approximating sparse inference. The authors show that unrolled iterations lead to better approximation if one allows the weights to vary at each layer, contrary to standard splitting algorithms. Whereas their focus is on relaxing the convergence hypothesis of iterative thresholding algorithms, we study a complementary question, namely when is speedup possible, without assuming strongly convex optimization. Their results are consistent with ours, since our analysis also shows that learning shared layer weights is less effective.

Inspired by the LISTA architecture, our mathematical analysis reveals that adaptive acceleration is related to a specific matrix factorization of the Gram matrix of the dictionary B = D T D as B = A T SA -R ,where A is unitary, S is diagonal and the residual is positive semidefinite: R 0. Our factorization balances between near diagonalization by asking that R is small and small perturbation of the 1 norm, i.e. Az 1 -z 1 is small. When this factorization succeeds, we prove that the resulting splitting algorithm enjoys a convergence rate with improved constants with respect to the non-adaptive version. Moreover, our analysis also shows that acceleration is mostly possible at the beginning of the iterative process, when the current estimate is far from the optimal solution, which is consistent with numerical experiments. We also show that the existence of this factorization is not only sufficient for acceleration, but also necessary. This is shown by constructing dictionaries whose Gram matrix diagonalizes in a basis that is incoherent with the canonical basis, and verifying that LISTA fails in that case to accelerate with respect to ISTA.

In our numerical experiments, we design a specialized version of LISTA called FacNet, with more constrained parameters, which is then used as a tool to show that our theoretical analysis captures the acceleration mechanism of LISTA. Our theoretical results can be applied to FacNet and as LISTA is a generalization of this model, it always performs at least as well, showing that the existence of the factorization is a sufficient certificate for acceleration by LISTA. Reciprocally, we show that for cases where no acceleration is possible with FacNet, the LISTA model also fail to provide acceleration, linking the two speedup mechanisms. This numerical evidence suggest that the existence of our proposed factorization is sufficient and somewhat necessary for LISTA to show good results.

The rest of the paper is structured as follows. Section 2 presents our mathematical analysis and proves the convergence of the adaptive algorithm as a function of the quality of the matrix factorization. Finally, Section 3 presents the generic architectures that will enable the usage of such schemes and the numerical experiments, which validate our analysis over a range of different scenarios.

2 Accelerating Sparse Coding with Sparse Matrix Factorizations

2.1. Unitary Proximal Splitting:

In this section we describe our setup for accelerating sparse coding based on the Proximal Splitting method. Let Ω ⊂ R n be the set describing our input data, and D ∈ R n×m be a dictionary, with m > n. We wish to find fast and accurate approximations of the sparse coding z * (x) of any x ∈ Ω, defined in (1) For simplicity, we denote B = D T D and y = D † x to rewrite (1) as

z * (x) = arg min z F x (z) = 1 2 (y -z) T B(y -z) E(z) + λ z 1 G(z)

.

(2)

For clarity, we will refer to F x as F and to z * (x) as z * . The classic proximal splitting technique finds z * as the limit of sequence (z k ) k , obtained by successively constructing a surrogate loss F k (z) of the form

EQUATION

satisfying F k (z) ≥ F (z) for all z ∈ R m . Since F k is separable in each coordinate of z, z k+1 = arg min z F k (z) can be computed efficiently. This scheme is based on a majoration of the quadratic form (y -z) T B(y -z) with an isotropic quadratic form L k z k -z 2 2 . The convergence rate of the splitting algorithm is optimized by choosing L k as the smallest constant satisfying F k (z) ≥ F (z), which corresponds to the largest singular value of B.

The computation of z k+1 remains separable by replacing the quadratic form L k I by any diagonal form. However, the Gram matrix B = D T D might be poorly approximated via diagonal forms for general dictionaries. Our objective is to accelerate the convergence of this algorithm by finding appropriate factorizations of the matrix B such that B ≈ A T SA , and Az 1 ≈ z 1 , where A is unitary and S is diagonal positive definite. Given a point z k at iteration k, we can rewrite F (z) as

EQUATION

with

Q B (v, w) := 1 2 (v -w) T B(v -w) + λ v 1 .

For any diagonal positive definite matrix S and unitary matrix A, the surrogate loss

F (z, z k ) := E(z k ) + (z k -y) T B(z -z k ) + Q S (Az, Az k ) can be explicitly minimized, since arg min z F (z, z k ) = A T arg min u (z k -y) T BA T (u -Az k ) + Q S (u, Az k ) = A T arg min u Q S u, Az k -S -1 AB(z k -y) (5)

where we use the variable change u = Az. As S is diagonal positive definite, ( 5) is separable and can be computed easily, using a linear operation followed by a point-wise non linear soft-thresholding. Thus, any couple (A, S) ensures an computationally cheap scheme. The question is then how to factorize B using S and A in an optimal manner, that is, such that the resulting proximal splitting sequence converges as fast as possible to the sparse coding solution.

2.2. Non-asymptotic Analysis:

We will now establish convergence results based on the previous factorization. These bounds will inform us on how to best choose the factors A k and S k in each iteration.

For that purpose, let us define

EQUATION

The quantity δ A (z) thus measures how invariant the 1 norm is to the unitary operator A, whereas R corresponds to the residual of approximating the original Gram matrix B by our factorization A T SA . Given a current estimate z k , we can rewrite

EQUATION

By imposing that R is a positive semidefinite residual one immediately obtains the following bound.

Proposition 2.1. Suppose that R = A T SA -B is positive definite, and define

EQUATION

Proof. By definition of z k+1 and using the fact that R 0 we have

F (z k+1 ) -F (z * ) ≤ F (z k+1 ) -F (z k+1 , z k ) + F (z * , z k ) -F (z * ) = - 1 2 (z k+1 -z k ) T R(z k+1 -z k ) -δ A (z k+1 ) + 1 2 (z * -z k ) T R(z * -z k ) + δ A (z * ) ≤ 1 2 (z * -z k ) T R(z * -z k ) + δ A (z * ) -δ A (z k+1 ) .

where the first line results from the definition of z k+1 and the third line makes use of R positiveness.

This simple bound reveals that to obtain fast approximations to the sparse coding it is sufficient to find S and A such that R is small and that the 1 commutation term δ A is small. These two conditions will be often in tension: one can always obtain R ≡ 0 by using the Singular Value Decomposition of B = A T 0 S 0 A 0 and setting A = A 0 and S = S 0 . However, the resulting A 0 might introduce large commutation error δ A0 . Similarly, as the absolute value is non-expansive, i.e. |a| -|b| ≤ a -b , we have that

|δ A (z)| = λ Az 1 -z 1 ≤ λ (A -I)z 1 (10) ≤ λ 2 max( Az 0 , z 0 ) • A -I • z 2 ,

where we have used the Cauchy-Schwartz inequality x 1 ≤ x 0 x 2 in the last equation. In particular, (10) shows that unitary matrices in the neighborhood of I with A -I small have small 1 commutation error δ A but can be inappropriate to approximate general B matrix.

The commutation error also depends upon the sparsity of z and Az . If both z and Az are sparse then the commutation error is reduced, which can be achieved if A is itself a sparse unitary matrix. Moreover, since

|δ A (z) -δ A (z )| ≤ λ| z 1 -z 1 | + λ| Az 1 -Az 1 | and | z 1 -z 1 | ≤ z -z 1 ≤ z -z 0 z -z 2 it

results that δ A is Lipschitz with respect to the Euclidean norm; let us denote by L A (z) its local Lipschitz constant in z, which can be computed using the norm of the subgradient in z1 . An uniform upper bound for this constant is (1 + A 1 )λ √ m, but it is typically much smaller when z and Az are both sparse. Equation ( 8) defines an iterative procedure determined by the pairs {(A k , S k )} k . The following theorem uses the previous results to compute an upper bound of the resulting sparse coding estimator. Theorem 2.2. Let A k , S k be the pair of unitary and diagonal matrices corresponding to iteration k, chosen such that

EQUATION

)

with α = k-1 i=1 2LA i (zi+1) z * -zi+1 2 + (z * -zi) T (Ri-1 -Ri)(z * -zi) , β = k-1 i=0 (i + 1) (zi+1 -zi) T Ri(zi+1 -zi) + 2δA i (zi+1) -2δA i (zi) ,

where L A (z) denote the local lipschitz constant of δ A at z.

Remarks: If one sets A k = I and S k = B I for all k ≥ 0, (11) corresponds to the bound of the ISTA algorithm (Beck & Teboulle, 2009) .

We can specialize the theorem in the case when A 0 , S 0 are chosen to minimize the bound (9) and

EQUATION

This corollary shows that by simply replacing the first step of ISTA by the modified proximal step detailed in (5), one can obtain an improved bound at fixed k as soon as

EQUATION

)

More generally, given a current estimate z k , searching for a factorization (A k , S k ) will improve the upper bound when

EQUATION

)

We emphasize that this is not a guarantee of acceleration, since it is based on improving an upper bound. However, it provides a simple picture on the mechanism that makes non-asymptotic acceleration possible.

2.3. Interpretation:

In this section we analyze the consequences of Theorem 2.2 in the design of fast sparse coding approximations, and provide a possible explanation for the behavior observed numerically.

2.3.1. 'Phase Transition" and Law of Diminishing Returns:

(14) reveals that the optimum matrix factorization in terms of minimizing the upper bound depends upon the current scale of the problem, that is, of the distance z * -z k . At the beginning of the optimization, when z * -z k is large, the bound ( 14) makes it easier to explore the space of factorizations (A, S) with A further away from the identity. Indeed, the bound tolerates larger increases in L A (z k+1 ), which is dominated by i.e. the sparsity of both z 1 and A 0 (z 1 ). On the other hand, when we reach intermediate solutions z k such that z * -z k is small with respect to L A (z k+1 ), the upper bound is minimized by choosing factorizations where A is closer and closer to the identity, leading to the non-adaptive regime of standard ISTA (A = Id).

L A (z k+1 ) ≤ λ( z k+1 0 + Az k+1 0 ) , W e + X Z W g (a) ISTA -Recurrent Neural Network X W (0) e W (1) g + W (1) e W (2) g + W (2) e Z (b) LISTA -Unfolded network

This is consistent with the numerical experiments, which show that the gains provided by learned sparse coding methods are mostly concentrated in the first iterations. Once the estimates reach a certain energy level, section 3 shows that LISTA enters a steady state in which the convergence rate matches that of standard ISTA.

The natural follow-up question is to determine how many layers of adaptive splitting are sufficient before entering the steady regime of convergence. A conservative estimate of this quantity would require an upper bound of z * -z k from the energy bound

F (z k ) -F (z * ).

Since in general F is convex but not strongly convex, such bound does not exist unless one can assume that F is locally strongly convex (for instance for sufficiently small values of F ).

2.3.2. Improving the factorization to particular input distributions:

Given an input dataset D = (x i , z

EQUATION

and sparse coding solutions z * i , the factorization adapted to D is defined as min

EQUATION

)

Therefore, adapting the factorization to a particular dataset, as opposed to enforcing it uniformly over a given ball B(z * ; R) (where the radius R ensures that the initial value z 0 ∈ B(z * ; R)), will always improve the upper bound (9). Studying the gains resulting from the adaptation to the input distribution will be let for future work.

3. Numerical Experiments:

This section provides numerical arguments to analyse adaptive optimization algorithms and their performances, and relates them to the theoretical properties developed in the previous section. All the experiments were run using Python and Tensorflow. For all the experiments, the training is performed using Adagrad (Duchi et al., 2011) . The code to reproduce the figures is available online2 .

3.1. Adaptive Optimization Networks Architectures:

LISTA/LFISTA In Gregor & Le Cun (2010), the authors introduced LISTA, a neural network constructed by considering ISTA as a recurrent neural net. At each step, ISTA performs the following 2-step procedure :

1. u k+1 = z k - 1 L D T (Dz k -x) = (I - 1 L D T D) Wg z k + 1 L D T We x , 2. z k+1 = h λ L (u k+1 ) where h θ (u) = sign(u)(|u| -θ) + ,            step k of ISTA (16)

This procedure combines a linear operation to compute u k+1 with an element-wise non linearity. It can be summarized as a recurrent neural network, presented in Figure 1a ., with tied weights. The autors in Gregor & Le Cun (2010) considered the architecture Φ K Θ with parameters Θ = (W

EQUATION

e , θ (k) ) k=1,...K obtained by unfolding K times the recurrent network, as presented in Figure 1b . The layers φ k Θ are defined as

EQUATION

If

EQUATION

= D T L and θ (k) = λ L are fixed for all the K layers, the output of this neural net is exactly the vector z K resulting from K steps of ISTA. With LISTA, the parameters Θ are learned using back propagation to minimize the cost function:

f (Θ) = E x F x (Φ K Θ (x)) .

A similar algorithm can be derived from FISTA, the accelerated version of ISTA to obtain LFISTA (see Figure 5 in Appendix A ). The architecture is very similar to LISTA, now with two memory tapes:

z k+1 = h θ (W g z k + W m z k-1 + W e x) .

Factorization network Our analysis in Section 2 suggests a refactorization of LISTA in more a structured class of parameters. Following the same basic architecture, and using ( 5), the network FacNet, Ψ K Θ is formed using layers such that:

EQUATION

with S diagonal and A unitary, the parameters of the k-th layer. The parameters obtained after training such a network with back-propagation can be used with the theory developed in Section 2. Up to the last linear operation A T of the network, this network is a re-parametrization of LISTA in a more constrained parameter space. Thus, LISTA is a generalization of this proposed network and should have performances at least as good as FacNet, for a fixed number of layers.

The optimization can also be performed using backpropagation. To enforce the unitary constraints on A (k) , the cost function is modified with a penalty:

EQUATION

with Θ = (A (k) , S (k) ) k=1...K the parameters of the K layers and µ a scaling factor for the regularization. The resulting matrix A (k) is then projected on the Stiefel Manifold using a SVD to obtain final parameters, coherent with the network structure.

Linear model Finally, it is important to distinguish the performance gain resulting from choosing a suitable starting point and the acceleration from our model. To highlights the gain obtain by changing the starting point, we considered a linear model with one layer such that z out = A (0) x. This model is learned using SGD with the convex cost function

f (A (0) ) = (I -DA (0) )x 2 2 + λ A (0) x 1 .

It computes a tradeoff between starting from the sparsest point 0 and a point with minimal reconstruction error y . Then, we observe the performance of the classical iteration of ISTA using z out as a stating point instead of 0 .

3.2. Synthetic problems with known distributions:

Gaussian dictionary In order to disentangle the role of dictionary structure from the role of data distribution structure, the minimization problem is tested using a synthetic generative model with no structure in the weights distribution. First, m atoms d i ∈ R n are drawn iid from a multivariate Gaussian with mean 0 and covariance I n and the dictionary D is defined as di / di 2 i=1...m . The data points are generated from its sparse codes following a Bernoulli-Gaussian model. The coefficients z = (z 1 , . . . , z m ) are constructed with z i = b i a i , where b i ∼ B(ρ) and a i ∼ N (0, σI m ) , where ρ controls the sparsity of the data. The values are set to m=100, n=64 for the dictionary dimension, ρ = 5 /m for the sparsity level and σ=10 for the activation coefficient generation parameters. The sparsity Figure 2 displays the cost performance for methods ISTA/FISTA/Linear relatively to their iterations and for methods LISTA/LFISTA/FacNet relatively to the number of layers used to solve our generated problem. Linear has performances comparable to learned methods with the first iteration but a gap appears as the number of layers increases, until a point where it achieves the same performances as non adaptive methods. This highlights that the adaptation is possible in the subsequent layers of the networks, going farther than choosing a suitable starting point for iterative methods. The first layers permit to achieve a large gain over the classical optimization strategy, by leveraging the structure of the problem. This appears even with no structure in the sparsity patterns of input data, in accordance with the results in the previous section. We also observe diminishing returns as the number of layers increases. This results from the phase transition described in Subsubsection 2.3.1, as the last layers behave as ISTA steps and do not speed up the convergence. The 3 learned algorithms are always performing at least as well as their classical counterpart, as it was stated in Theorem 2.2. We also explored the effect of the sparsity level in the training and learning of adaptive networks. In the denser setting, the arbitrage between the 1 -norm and the squared error is easier as the solution has a lot of non zero coefficients. Thus in this setting, the approximate method is more precise than in the very sparse setting where the approximation must perform a fine selection of the coefficients. But it also yield lower gain at the beggining as the sparser solution can move faster.

There is a small gap between LISTA and FacNet in this setup. This can be explained from the extra constraints on the weights that we impose in the FacNet, which effectively reduce the parameter space by half. Also, we implement the unitary constraints on the matrix A by a soft regularization (see ( 19)), involving an extra hyper-parameter µ that also contributes to the small performance gap. In any case, these experiments show that our analysis accounts for most of the acceleration provided by LISTA, as the performance of both methods are similar, up to optimization errors.

Adversarial dictionary:

The results from Section 2 show that problems with a gram matrix composed of large eigenvalues associated to non sparse eigenvectors are harder to accelerate. Indeed, it is not possible in this case to find a quasi diagonalization of the matrix B that does not distort the 1 norm. It is possible to generate such a dictionary using Harmonic Analysis. The Discrete Fourier Transform (DFT) distorts a lot the 1 ball, since a very sparse vector in the temporal space is transformed in widely spread spectrum in the Fourier domain. We can thus design a dictionary for which LISTA and FacNet performances should be degraded.

D = di / di 2 i=1...m is constructed such that d j,k = e -2πijζ k , with ζ k k≤n

randomly selected from 1 /m, . . . , m /2 /m without replacement.

The resulting performances are reported in Figure 3 . The first layer provides a big gain by changing the starting point of the iterative methods. It realizes an arbitrage of the tradeoff between starting from 0 and starting from y . But the next layers do not yield any extra gain compared to the original ISTA algorithm. After 4 layers, the cost performance of both adaptive methods and ISTA are equivalent. It is clear that in this case, FacNet does not accelerate efficiently the sparse coding, in accordance with our result from Section 2. LISTA also displays poor performances in this setting. This provides further evidence that FacNet and LISTA share the same acceleration mechanism as adversarial dictionaries for FacNet are also adversarial for LISTA.

3.3. Sparse coding with over complete dictionary on images:

Wavelet encoding for natural images A highly structured dictionary composed of translation invariant Haar wavelets is used to encode 8x8 patches of images from the PASCAL VOC 2008 dataset. The network is used to learn an efficient sparse coder for natural images over this family. 500 images are sampled from dataset to train the encoder. Training batches are obtained by uniformly sampling patches from the training image set to feed the stochastic optimization of the network. The encoder is then tested with 10000 patches sampled from 100 new images from the same dataset.

Learned dictionary for MNIST To evaluate the performance of LISTA for dictionary learning, LISTA was used to encode MNIST images over an unconstrained dictionary, learned a priori using classical dictionary learning techniques. The dictionary of 100 atoms was learned from 10000 MNIST images in grayscale rescaled to 17x17 using the implementation of Mairal et al. (2009) proposed in scikit-learn, with λ = 0.05. Then, the networks were trained through backpropagation using all the 60000 images from the training set of MNIST. Finally, the perfornance of these encoders were evaluated with the 10000 images of the training set of MNIST.

The Figure 4 displays the cost performance of the adaptive procedures compared to nonadaptive algorithms. In both scenario, FacNet has performances comparable to the one of LISTA and their behavior are in accordance with the theory developed in Section 2. The gains become smaller for each added layer and the initial gain is achieved for dictionary either structured or unstructured. The MNIST case presents a much larger gain compare to the experiment with natural images. This results from the difference of structure of the input distribution, as the MNIST digits are much more constrained than patches from natural images and the network is able to leverage it to find a better encoder. In the MNIST case, a network composed of 12 layers is sufficient to achieve performance comparable to ISTA with more than 1000 iterations.

4. Conclusions:

In this paper we studied the problem of finite computational budget approximation of sparse coding. Inspired by the ability of neural networks to accelerate over splitting methods on the first few iterations, we have studied which properties of the dictionary matrix and the data distribution lead to such acceleration. Our analysis reveals that one can obtain acceleration by finding approximate matrix factorizations of the dictionary which nearly diagonalize its Gram matrix, but whose orthogonal transformations leave approximately invariant the 1 ball. By appropriately balancing these two conditions, we show that the resulting rotated proximal splitting scheme has an upper bound which improves over the ISTA upper bound under appropriate sparsity.

In order to relate this specific factorization property to the actual LISTA algorithm, we have introduced a reparametrization of the neural network that specifically computes the factorization, and incidentally provides reduced learning complexity (less parameters) from the original LISTA. Numerical experiments of Section 3 show that such reparametrization recovers the same gains as the original neural network, providing evidence that our theoretical analysis is partially explaining the behavior of the LISTA neural network. Our acceleration scheme is inherently transient, in the sense that once the iterates are sufficiently close to the optimum, the factorization is not effective anymore. This transient effect is also consistent with the performance observed numerically, although the possibility remains open to find alternative models that further exploit the particular structure of the sparse coding. Finally, we provide evidence that successful matrix factorization is not only sufficient but also necessary for acceleration, by showing that Fourier dictionaries are not accelerated.

Despite these initial results, a lot remains to be understood on the general question of optimal tradeoffs between computational budget and statistical accuracy. Our analysis so far did not take into account any probabilistic consideration (e.g. obtain approximations that hold with high probability or in expectation). Another area of further study is the extension of our analysis to the FISTA case, and more generally to other inference tasks that are currently solved via iterative procedures compatible with neural network parametrizations, such as inference in Graphical Models using Belief Propagation or other ill-posed inverse problems.

A Learned Fista:

A similar algorithm can be derived from FISTA, the accelerated version of ISTA to obtain LFISTA (see Figure 5 ). The architecture is very similar to LISTA, now with two memory taps: It introduces a momentum term to improve the convergence rate of ISTA as follows:

1.

y k = z k + t k-1 -1 t k (z k -z k-1 ) , 2. z k+1 = h λ L y k - 1 L ∇E(y k ) = h λ L (I - 1 L B)y k + 1 L D T x , 3. t k+1 = 1 + 1 + 4t 2 k 2 .

By substituting the expression for y k into the first equation, we obtain a generic recurrent architecture very similar to LISTA, now with two memory taps, that we denote by LFISTA: This model is equivalent to running K-steps of FISTA when its parameters are initialized with

EQUATION

W (k) g = 1 + t k-1 -1 t k I - 1 L B , W (k) m = 1 -t k-1 t k I - 1 L B , W (k) e = 1 L D T .

The parameters of this new architecture, presented in Figure 5 , are trained analogously as in the LISTA case.

B Proofs:

Lemma B.1. Suppose that R = A T SA -B is positive definite, and define

EQUATION

δ A (z) = Az 1 -z 1 . Then we have

EQUATION

Proof. We define f (t) = F tz k+1 + (1 -t)z * , t ∈ [0, 1] .

Since F is convex, f is also convex in [0, 1]. Since f (0) = F (z * ) is the global minimum, it results that f (t) is increasing in (0, 1], and hence

F (z k+1 ) -F (z * ) = f (1) -f (0) = f (t)dt ≤ f (1) ,

where f (1) is any element of ∂f (1). Since δ A (z) is a difference of convex functions, its subgradient can be defined as a limit of infimal convolutions Hiriart-Urruty (1991). We have ∂f Adding these inequalities for n = 0 . . . k -1 we obtain

  k-1 n=0 F (z n )   -kF (z * ) ≤ k-1 n=0

∇δ An (z n+1 ), (z * -z n+1 ) + ( 25)

+ 1 2 (z * -z 0 ) T R 0 (z * -z 0 ) -(z * -z k ) T R k-1 (z * -z k ) + + 1 2 k-1 n=1 (z * -z n ) T (R n-1 -R n )(z * -z n ) .

On the other hand, we also have

F (z n ) -F (z n+1 ) ≥ F (z n ) -F (z n , z n ) + F (z n+1 , z n ) -F (z n+1 ) = -δ An (z n ) + δ An (z n+1 ) + 1 2 (z n+1 -z n ) T R n (z n+1 -z n ) ,

which results in (n + 1) 1 2 (z n+1 -z n ) T R n (z n+1 -z n ) + δ An (z n+1 ) -δ An (z n ) .

Combining ( 25) and ( 26) we obtain (28)

EQUATION

Proof: We verify that in that case, R n-1 -R n ≡ 0 and for n > 1 and δ An ≡ 0 for n > 0 .

:

This quantity exists as δA is a difference of convex. See proof of ?? in appendices for precisions.

The code can be found at https://github.com/tomMoral/AdaptiveOptim

Figures and Tables:

Figure 1: Network architecture for ISTA/LISTA. The unfolded version (b) is trainable through backpropagation and permits to approximate the sparse coding solution efficiently.

Figure 2: Evolution of the cost function F (z k ) -F (z * ) with the number of layers or the number of iteration k for different sparsity level. (left) ρ = 1 /20 and (right)ρ = 1 /4 .

Figure 3: Evolution of the cost function F (z k ) -F (z * ) with the number of layers or the number of iteration k for a problem generated with an adversarial dictionary.

Figure 4: Evolution of the cost function F (z k ) -F (z * ) with the number of layers or the number of iteration k for two image datasets.

Figure5: Network architecture for LFISTA. This network is trainable through backpropagation and permits to approximate the sparse coding solution efficiently.

1) = ∂F (z k+1 ), z k+1 -z * , and since∂F (z) = ∂ F (z, z k ) -R(z -z k ) -∂δ A (z) and 0 ∈ ∂ F (z k+1 , z k ) it results that ∂F (z k+1 ) = -R(z k+1 -z k ) -∂δ A (z k+1) , and thusF (z k+1 ) -F (z * ) ≤ (z * -z k+1 ) T R(z k+1 -z k ) + ∂δ A (z k+1 ), (z * -z k+1 ) . (22) (21) is obtained by observing that (z * -z k+1 ) T R(z k+1 -z k ) ≤ 1 2 (z * -z k ) T R(z * -z k ) -(z * -z k+1 ) T R(z * -z k+1 ) ,(23)thanks to the fact that R 0.Theorem B.2. Let A k , S k be the pair of unitary and diagonal matrices corresponding to iteration k, chosen such thatR k = A T k S k A k -B 0. It results that F (z k ) -F (z * ) ≤ (z * -z 0 ) T R 0 (z * -z 0 ) + 2 ∇δ A0 (z 1 ), (z * -∇δ An (z n+1 ), (z * -z n+1 ) + (z * -z n ) T (R n-1 -R n )(z * -z n ) (z n+1 -z n ) T R n (z n+1 -z n ) + 2δ An (z n+1 ) -2δ An (z n ) .Proof: The proof is adapted from(Beck & Teboulle, 2009), Theorem 3.1. From Lemma B.1, we start by using (21) to bound terms of the form F (z n ) -F (z * ):F (zn)-F (z * ) ≤ ∇δA n (zn+1), (z * -zn+1) + 1 2 (z * -zn) T Rn(z * -zn) -(z * -zn+1) T Rn(z * -zn+1) .

(z n+1 -z n ) T R n (z n+1 -z n ) 1) δ An (z n+1 ) -δ An (z n )

∇δ An (z n+1 ), (z * -z n+1 ) + (z * -z n ) T (R n-1 -R n )(z * -z n ) z n+1 -z n ) T R n (z n+1 -z n ) + 2δ An (z n+1 ) -2δ An (z n ) . Corollary B.3. If A k = I, S k = B I for k > 0 then F (z k )-F (z * ) ≤ (z * -z 0 ) T R 0 (z * -z 0 ) + 2L A0 (z 1 )( z * -z 1 + z 1 -z 0 ) + (z * -z 1 ) T R 0 (z * -z 1 ) T 2k .

Peer review (iclr2017):
This paper performs theoretical analysis to understand how sparse coding could be accelerated by neural networks. The neural networks are generated by unfolding the ISTA/FISTA iterations. Based on the results, the authors proposed a reparametrization approach for the neural network architecture to enforce the factorization property and recovered the original gain of LISTA, which justified the theoretical analysis. My comments are listed below. It is not clear about the purpose of Section 2.3.2. Adapting the factorization to the input distribution based on (15) would be time consuming because the overhead of solving (15) may not save the total time. In fact, the approach does not use (15) but back propagation to learn the factorization parameters. Minor comments: - E(z_k) in (3) and (4) are not defined. - E_x in (19) is not defined. - Forward referencing (“Equation (20) defines…”) in the paragraph above Theorem 2.2. needs to be corrected.