UNDERSTANDING APPROXIMATE AND UNROLLED DICTIONARY LEARNING FOR PATTERN RECOVERY
Year: 2022

Abstract:
Dictionary learning consists of finding a sparse representation from noisy data and is a common way to encode data-driven prior knowledge on signals. Alternating minimization (AM) is standard for the underlying optimization, where gradient descent steps alternate with sparse coding procedures. The major drawback of this method is its prohibitive computational cost, making it unpractical on large real-world data sets. This work studies an approximate formulation of dictionary learning based on unrolling and compares it to alternating minimization to find the best trade-off between speed and precision. We analyze the asymptotic behavior and convergence rate of gradients estimates in both methods. We show that unrolling performs better on the support of the inner problem solution and during the first iterations. Finally, we apply unrolling on pattern learning in magnetoencephalography (MEG) with the help of a stochastic algorithm and compare the performance to a state-of-the-art method.

text: Dictionary learning consists of finding a sparse representation from noisy data and is a common way to encode data-driven prior knowledge on signals. Alternating minimization (AM) is standard for the underlying optimization, where gradient descent steps alternate with sparse coding procedures. The major drawback of this method is its prohibitive computational cost, making it unpractical on large real-world data sets. This work studies an approximate formulation of dictionary learning based on unrolling and compares it to alternating minimization to find the best trade-off between speed and precision. We analyze the asymptotic behavior and convergence rate of gradients estimates in both methods. We show that unrolling performs better on the support of the inner problem solution and during the first iterations. Finally, we apply unrolling on pattern learning in magnetoencephalography (MEG) with the help of a stochastic algorithm and compare the performance to a state-of-the-art method.

heading: 1 INTRODUCTION
text: Pattern learning provides insightful information on the data in various biomedical applications. Typical examples include the study of magnetoencephalography (MEG) recordings, where one aims to analyze the electrical activity in the brain from measurements of the magnetic field around the scalp of the patient (Dupré la Tour et al., 2018). One may also mention neural oscillations study in the local field potential (Cole & Voytek, 2017) or QRS complex detection in electrocardiograms (Xiang et al., 2018) among others.
Dictionary learning (Olshausen & Field, 1997; Aharon et al., 2006; Mairal et al., 2009) is particularly efficient on pattern learning tasks, such as blood cells detection (Yellin et al., 2017) and MEG signals analysis (Dupré la Tour et al., 2018). This framework assumes that the signal can be decomposed into a sparse representation in a redundant basis of patterns – also called atoms. In other words, the goal is to recover a sparse code Z ∈ Rn×T and a dictionary D ∈ Rm×n from noisy measurements Y ∈ Rm×T which are obtained as the linear transformation DZ, corrupted with noise B ∈ Rm×T : Y = DZ + B. Theoretical elements on identifiability and local convergence have been proven in several studies (Gribonval et al., 2015; Haeffele & Vidal, 2015; Agarwal et al., 2016; Sun et al., 2016). Sparsity-based optimization problems related to dictionary learning generally rely on the usage of the `0 or `1 regularizations. In this paper, we study Lasso-based (Tibshirani, 1996) dictionary learning where the dictionary D is learned in a set of constraints C by solving
min Z∈Rn×T ,D∈C
F (Z,D) , 1
2 ‖DZ − Y ‖22 + λ ‖Z‖1 . (1)
Dictionary learning can be written as a bi-level optimization problem to minimize the cost function with respect to the dictionary only, as mentioned in Mairal et al. (2009),
min D∈C G(D) , F (Z∗(D),D) with Z∗(D) = arg min Z∈Rn×T F (Z,D) . (2)
Computing the data representation Z∗(D) is often referred to as the inner problem, while the global minimization is the outer problem. Classical constraint sets include the unit norm, where each atom is normalized to avoid scale-invariant issues, and normalized convolutional kernels to perform Convolutional Dictionary Learning (Grosse et al., 2007).
Classical dictionary learning methods solve this bi-convex optimization problem through Alternating Minimization (AM) (Mairal et al., 2009). It consists in minimizing the cost function F over Z with a fixed dictionary D and then performing projected gradient descent to optimize the dictionary with a fixed Z. While AM provides a simple strategy to perform dictionary learning, it can be inefficient on large-scale data sets due to the need to resolve the inner problems precisely for all samples. In recent years, many studies have focused on algorithm unrolling (Tolooshams et al., 2020; Scetbon et al., 2021) to overcome this issue. The core idea consists of unrolling the algorithm, which solves the inner problem, and then computing the gradient with respect to the dictionary with the help of back-propagation through the iterates of this algorithm. Gregor & LeCun (2010) popularized this method and first proposed to unroll ISTA (Daubechies et al., 2004) – a proximal gradient descent algorithm designed for the Lasso – to speed up the computation of Z∗(D). The N + 1-th layer of this network – called LISTA – is obtained as ZN+1 = ST λ
L (W 1Y + W 2ZN ), with ST be-
ing the soft-thresholding operator. This work has led to many contributions aiming at improving this method and providing theoretical justifications in a supervised (Chen et al., 2018; Liu & Chen, 2019) or unsupervised (Moreau & Bruna, 2017; Ablin et al., 2019) setting. For such unrolled algorithms, the weights W 1 and W 2 can be re-parameterized as functions of D – as illustrated in Figure A in appendix – such that the output ZN (D) matches the result of N iterations of ISTA, i.e.
W 1D = 1
L D> and W 2D =
( I − 1
L D>D
) , where L = ‖D‖2 . (3)
Then, the dictionary can be learned by minimizing the loss F (ZN (D),D) over D with backpropagation. This approach is generally referred to as Deep Dictionary Learning (DDL). DDL and variants with different kinds of regularization (Tolooshams et al., 2020; Lecouat et al., 2020; Scetbon et al., 2021), image processing based on metric learning (Tang et al., 2020), and classification tasks with scattering (Zarka et al., 2019) have been proposed in the literature, among others. While these techniques have achieved good performance levels on several signal processing tasks, the reasons they speed up the learning process are still unclear.
In this work, we study unrolling in Lasso-based dictionary learning as an approximate bi-level optimization problem. What makes this work different from Bertrand et al. (2020), Ablin et al. (2020) and Tolooshams & Ba (2021) is that we study the instability of non-smooth bi-level optimization and unrolled sparse coding out of the support, which is of major interest in practice with a small number of layers. In Section 2, we analyze the convergence of the Jacobian computed with automatic differentiation and find out that its stability is guaranteed on the support of the sparse codes only. De facto, numerical instabilities in its estimation make unrolling inefficient after a few dozen iterations. In Section 3, we empirically show that unrolling leads to better results than AM only with a small number of iterations of sparse coding, making it possible to learn a good dictionary in this setting. Then we adapt a stochastic approach to make this method usable on large data sets, and we apply it to pattern learning in magnetoencephalography (MEG) in Section 4. We do so by adapting unrolling to rank one convolutional dictionary learning on multivariate time series (Dupré la Tour et al., 2018). We show that there is no need to unroll more than a few dozen iterations to obtain satisfying results, leading to a significant gain of time compared to a state-of-the-art algorithm.

heading: 2 BI-LEVEL OPTIMIZATION FOR APPROXIMATE DICTIONARY LEARNING
text: As Z∗(D) does not have a closed-form expression, G cannot be computed directly. A solution is to replace the inner problem Z∗(D) by an approximation ZN (D) obtained through N iterations of a numerical optimization algorithm or its unrolled version. This reduces the problem to minimizing GN (D) , F (ZN (D),D). The first question is how sub-optimal global solutions of GN are
compared to the ones of G. Proposition 2.1 shows that the global minima of GN converge as fast as the numerical approximation ZN in function value.
Proposition 2.1 Let D∗ = arg minD∈C G(D) and D∗N = arg minD∈C GN (D), where N is the number of unrolled iterations. We denote by K(D∗) a constant depending on D∗, and by C(N) the convergence speed of the algorithm, which approximates the inner problem solution. We have
GN (D ∗ N )−G(D∗) ≤ K(D∗)C(N) .
The proofs of all theoretical results are deferred to Appendix C. Proposition 2.1 implies that when ZN is computed with FISTA (Beck & Teboulle, 2009), the function value for global minima of GN converges with speed C(N) = 1N2 towards the value of the global minima of F . Therefore, solving the inner problem approximately leads to suitable solutions for equation 2, given that the optimization procedure is efficient enough to find a proper minimum of GN . As the computational cost of zN increases with N , the choice of N results in a trade-off between the precision of the solution and the computational efficiency, which is critical for processing large data sets.
Moreover, learning the dictionary and computing the sparse codes are two different tasks. The loss GN takes into account the dictionary and the corresponding approximation ZN (D) to evaluate the quality of the solution. However, the dictionary evaluation should reflect its ability to generate the same signals as the ground truth data and not consider an approximate sparse code that can be recomputed afterward. Therefore, we should distinguish the ability of the algorithm to recover a good dictionary from its ability to learn the dictionary and the sparse codes at the same time. In this work, we use the metric proposed in Moreau & Gramfort (2020) for convolutions to evaluate the quality of the dictionary. We compare the atoms using their correlation and denote as C the cost matrix whose entry i, j compare the atom i of the first dictionary and j of the second. We define a sign and permutation invariant metric S(C) = maxσ∈Sn 1 n ∑n i=1 |Cσ(i),i|, where Sn is the group of permutations of [1, n]. This metric corresponds to the best linear sum assignment on the cost matrix C, and it can be computed with the Hungarian algorithm. Note that doing so has several limitations and that evaluating the dictionary is still an open problem. Without loss of generality, let T = 1 and thus z ∈ Rn in the rest of this section.
Gradient estimation in dictionary learning. Approximate dictionary learning is a non-convex problem, meaning that good or poor local minima of GN may be reached depending on the initialization, the optimization path, and the structure of the problem. Therefore, a gradient descent onGN has no guarantee to find an adequate minimizer of G. While complete theoretical analysis of these problems is arduous, we propose to study the correlation between the gradient obtained with GN and the actual gradient of G, as a way to ensure that the optimization dynamics are similar. Once z∗(D) is known, Danskin (1967, Thm 1) states that g∗(D) = ∇G(D) is equal to∇2F (z∗(D),D), where∇2 indicates that the gradient is computed relatively to the second variable in F . Even though the inner problem is non-smooth, this result holds as long as the solution z∗(D) is unique. In the following, we will assume that D>D is invertible on the support of z∗(D), which implies the uniqueness of z∗(D). This occurs with probability one if D is sampled from a continuous distribution (Tibshirani, 2013). AM and DDL differ in how they estimate the gradient of G. AM relies on the analytical formula of g∗ and uses an approximation zN of z∗, leading to the approximate gradient g1N (D) = ∇2F (zN (D),D). We evaluate how well g1N approximates g∗ in Proposition 2.2.
Proposition 2.2 Let D ∈ Rm×n. Then, there exists a constant L1 > 0 such that for every number of iterations N ∥∥g1N − g∗∥∥ ≤ L1 ‖zN (D)− z∗(D)‖ . Proposition 2.2 shows that g1N converges as fast as the iterates of ISTA converge. DDL computes the gradient automatically through zN (D). As opposed to AM, this directly minimizes the loss GN (D). Automatic differentiation yields a sub-gradient g2N (D) such that
g2N (D) ∈ ∇2F (zN (D),D) + J + N ( ∂1F (zN (D),D) ) , (4)
where JN : Rm×n → Rn is the weak Jacobian of zN (D) with respect to D and J+N denotes its adjoint. The product between J+N and ∂1F (zN (D),D) is computed via automatic differentiation.
Proposition 2.3 Let D ∈ Rm×n. Let S∗ be the support of z∗(D), SN be the support of zN and S̃N = SN ∪ S∗. Let f(z,D) = 12 ‖Dz − y‖ 2 2 be the data-fitting term in F . Let R(J, S̃) =
J+ ( ∇21,1f(z∗,D) 1S̃ ) +∇22,1f(z∗,D) 1S̃ . Then there exists a constant L2 > 0 and a subsequence of (F)ISTA iterates zφ(N) such that for all N ∈ N:
∃ g2φ(N) ∈ ∇2f(zφ(N),D) + J + φ(N) ( ∇1f(zφ(N),D) + λ∂‖·‖1(zφ(N)) ) s.t. :∥∥∥g2φ(N) − g∗∥∥∥ ≤ ∥∥∥R(Jφ(N), S̃φ(N))∥∥∥∥∥zφ(N) − z∗∥∥+ L22 ∥∥zφ(N) − z∗∥∥2 .
This sub-sequence zφ(N) corresponds to iterates on the support of z∗.
Proposition 2.3 shows that g2N may converge faster than g 1 N once the support is reached.
Ablin et al. (2020) and Tolooshams & Ba (2021) have studied the behavior of strongly convex functions, as it is the case on the support, and found similar results. This allowed Tolooshams & Ba (2021) to focus on support identification and show that automatic differentiation leads to a better gradient estimation in dictionary learning on the support under minor assumptions.
However, we are also interested in characterizing the behavior outside of the support, where the gradient estimation is difficult because of the sub-differential. In practice, automatic differentiation uses the sign operator as a sub-gradient of ‖·‖1. The convergence behavior of g2N is also driven by R(JN , S̃N ) and thus by the weak Jacobian computed via back-propagation. We first compute a closed-form expression of the weak Jacobian of z∗(D) and zN (D). We then show that R(JN , S̃N ) ≤ L ‖JN − J∗‖ and we analyze the convergence of JN towards J∗.
Study of the Jacobian. The computation of the Jacobian can be done by differentiating through ISTA. In Theorem 2.4, we show that JN+1 depends on JN and the past iterate zN , and converges towards a fixed point. This formula can be used to compute the Jacobian during the forward pass, avoiding the computational cost of back-propagation and saving memory.
Theorem 2.4 At iteration N + 1 of ISTA, the weak Jacobian of zN+1 relatively to Dl, where Dl is the l-th row of D, is given by induction:
∂(zN+1)
∂Dl = 1|zN+1|>0
( ∂(zN )
∂Dl − 1 L
( Dlz > N + (D > l zN − yl)In + D>D ∂(zN )
∂Dl
)) .
∂(zN ) ∂Dl will be denoted by JNl . It converges towards the weak Jacobian J ∗ l of z ∗ relatively to Dl, whose values are
J∗l S∗ = −(D > :,S∗D:,S∗) −1(Dlz ∗> + (D>l z ∗ − yl)In)S∗ , on the support S∗ of z∗, and 0 elsewhere. Moreover, R(J∗, S∗) = 0.
This result is similar to Bertrand et al. (2020) where the Jacobian of z is computed over λ to perform hyper-parameter optimization in Lasso-type models. Using R(J∗, S∗) = 0, we can write∥∥∥R(JN , S̃N )∥∥∥ ≤ ∥∥∥R(JN , S̃N )−R(J∗, S∗)∥∥∥ ≤ L ‖JN − J∗‖ , (5) as ∥∥∇21,1f(z∗, D)∥∥2 = L. If the back-propagation were to output an accurate estimate JN of the
weak Jacobian J∗, ∥∥∥R(JN , S̃N )∥∥∥ would be 0, and the convergence rate of g2N could be twice as fast as the one of g1N . To quantify this, we now analyze the convergence of JN towards J ∗. In
Proposition 2.5, we compute an upper bound of ∥∥JNl − J∗l ∥∥ with possible usage of truncated backpropagation (Shaban et al., 2019). Truncated back-propagation of depth K corresponds to an initial estimate of the Jacobian JN−K = 0 and iterating the induction in Theorem 2.4.
Proposition 2.5 Let N be the number of iterations and K be the back-propagation depth. We assume that ∀n ≥ N−K, S∗ ⊂ Sn. Let ĒN = Sn\S∗, let L be the largest eigenvalue ofD>:,S∗D:,S∗ , and let µn be the smallest eigenvalue of D>:,SnD:,Sn−1 . Let Bn = ∥∥∥PEn −D>:,ĒnD†>:,S∗PS∗∥∥∥, where PS is the projection on RS and D† is the pseudo-inverse of D. We have∥∥JNl − J∗l ∥∥ ≤ K∏
k=1
( 1− µN−k
L
) ‖J∗l ‖+ 2
L ‖Dl‖ K−1∑ k=0 k∏ i=1 (1−µN−i L ) (∥∥zN−kl − z∗l ∥∥+BN−k ‖z∗l ‖) .
Proposition 2.5 reveals multiple stages in the Jacobian estimation. First, one can see that if all iterates used for the back-propagation lie on the support S∗, the Jacobian estimate has a quasi-linear convergence, as shown in the following corollary.
Corollary 2.6 Let µ > 0 be the smallest eigenvalue of D>:,S∗D:,S∗ . Let K ≤ N be the backpropagation depth and let ∆N = F (zN , D) − F (z∗, D) + L2 ‖zN − z
∗‖. Suppose that ∀n ∈ [N −K,N ]; Sn ⊂ S∗. Then, we have∥∥J∗l − JNl ∥∥ ≤ (1− µL)K ‖J∗l ‖+K (1− µL)K−1 ‖Dl‖ 4∆N−KL2 . Once the support is reached, ISTA also converges with the same linear rate (1 − µL ). Thus the gradient estimate g2N converges almost twice as fast as g 1 N in the best case – with optimal subgradient – as O(K(1− µL ) 2K). This is similar to Ablin et al. (2020, Proposition.5) and Tolooshams
& Ba (2021). Second, Proposition 2.5 shows that ∥∥J∗l − JNl ∥∥ may increase when the support is not well-estimated, leading to a deterioration of the gradient estimate. This is due to an accumulation of errors materialized by the sum in the right-hand side of the inequality, as the term BN ‖z∗‖ may not vanish to 0 as long as SN 6⊂ S∗. Interestingly, once the support is reached at iteration S < N , the errors converge linearly towards 0, and we recover the fast estimation of g∗ with g2. Therefore, Lasso-based DDL should either be used with a low number of steps or truncated back-propagation to ensure stability. These results apply for all linear dictionaries, including convolutions.
Numerical illustrations. We now illustrate these theoretical results depending on the number N of unrolled iterations. The data are generated from a random Gaussian dictionary D of size 30×50, with Bernoulli-Gaussian sparse codes z (sparsity 0.3, σ2z = 1), and Gaussian noise (σ2noise = 0.1) – more details in Appendix A.
Figure 1 confirms the linear convergence of JNl once the support is reached. However, the convergence might be unstable when the number of iteration grows, leading to exploding gradient, as illustrated in the second case. When this happens, using a small number of iterations or truncated back-propagation becomes necessary to prevent accumulating errors. It is also of interest to look at the proportion of unstable Jacobians (see Figure 2). We recover behaviors observed in the first and second case in Figure 1. 40% samples suffer from numerical instabilities in this example. This has a negative impact on the gradient estimation outside of the support.
We display the convergence behavior of the gradients estimated by AM and by DDL with different back-propagation depths (20, 50, full) for simulated data and images in Figure 3. We unroll FISTA instead of ISTA to make the convergence faster. We observed similar

heading: BP depth
text: behaviors for both algorithms in early iterations but using ISTA required too much memory to reach full convergence. As we optimize using a line search algorithm, we are mainly interested in the ability of the estimate to provide an adequate descent direction. Therefore, we display the convergence in angle defined as the cosine similarity 〈g, g∗〉 = Tr(g
T g∗) ‖g‖‖g∗‖ . The angle provides a good metric to
assert that the two gradients are correlated and thus will lead to similar optimization paths. We also provide the convergence in norm in appendix. We compare g1N and g 2 N with the relative difference of their angles with g∗, defined as 〈g 2 N ,g ∗〉−〈g1N ,g
∗〉 1−〈g1N ,g∗〉 . When its value is positive, DDL provides the
best descent direction. Generally, when the back-propagation goes too deep, the performance of g2N decreases compared to g1N , and we observe large numerical instabilities. This behavior is coherent with the Jacobian convergence patterns studied in Proposition 2.5. Once on the support, g2N reaches back the performance of g1N as anticipated. In the case of a real image, unrolling beats AM by up to 20% in terms of gradient direction estimation when the number of iterations does not exceed 50, especially with small back-propagation depth. This highlights that the principal interest of unrolled algorithms is to use them with a small number of layers – i.e., a small number of iterations.

heading: 3 APPROXIMATE DICTIONARY LEARNING IN PRACTICE
text: This section introduces practical guidelines on Lasso-based approximate dictionary learning with unit norm constraint, and we provide empirical justifications for its ability to recover the dictionary. We also propose a strategy to scale DDL with a stochastic optimization method. We provide a full description of all our experiments in Appendix A. We optimize with projected gradient descent combined to a line search to compute high-quality steps sizes. The computations have been performed on a GPU NVIDIA Tesla V100-DGXS 32GB using PyTorch (Paszke et al., 2019).1
Improvement of precision. As stated before, a low number of iterations allows for efficient and stable computations, but this makes the sparse code less precise. One can learn the steps sizes of (F)ISTA to speed up convergence and compensate for imprecise representations, as proposed by Ablin et al. (2019) for LISTA. To avoid poor results due to large degrees of freedom in unsupervised learning, we propose a method in two steps to refine the initialization of the dictionary before relaxing the constraints on the steps sizes:
1. We learn the dictionary with fixed steps sizes equal to 1L where L = ‖D‖ 2, given by convergence
conditions. Lipschitz constants or upper bounds are computed at each gradient step with norms, or the FFT for convolutions, outside the scope of the network graph.
2. Then, once convergence is reached, we jointly learn the step sizes and the dictionary. Both are still updated using gradient descent with line search to ensure stable optimization.
1Code is available at https://github.com/bmalezieux/unrolled_dl.
The use of LISTA-like algorithms with no ground truth generally aims at improving the speed of sparse coding when high precision is not required. When it is the case, the final sparse codes can be computed separately with FISTA (Beck & Teboulle, 2009) or coordinate descent (Wu et al., 2008) to improve the quality of the representation.

heading: 3.1 OPTIMIZATION DYNAMICS IN APPROXIMATE DICTIONARY LEARNING
text: In this part, we study empirical properties of approximate dictionary learning related to global optimization dynamics to put our results on gradient estimation in a broader context.
Unrolling v. AM. In Figure 4, we show the number of gradient steps before reaching convergence, the behavior of the loss FN , and the recovery score defined at the beginning of the section for synthetic data generated by a Gaussian dictionary. As a reminder, S(C) = maxσ∈Sn 1 n ∑n i=1 |Cσ(i),i| where C is the correlation matrix between the columns of the true dictionary and the estimate. The number of iterations corresponds to N in the estimate zN (D). First, DDL leads to fewer gradient steps than AM in the first iterations. This suggests that automatic differentiation better estimates the directions of the gradients for small depths. However, computing the gradient requires backpropagating through the algorithm, and DDL takes 1.5 times longer to perform one gradient step than AM on average for the same number of iterationsN . When looking at the loss and the recovery score, we notice that the advantage of DDL for the minimization of FN is minor without learning the steps sizes, but there is an increase of performance concerning the recovery score. DDL better estimates the dictionary for small depths, inferior to 50. When unrolling more iterations, AM performs as well as DDL on the approximate problem and is faster.
Approximate DL. Figure 4 shows that high-quality dictionaries are obtained before the convergence of FN , either with AM or DDL. 40 iterations are sufficient to reach a reasonable solution
concerning the recovery score, even though the loss is still very far from the optimum. This suggests that computing optimal sparse codes at each gradient step is unnecessary to recover the dictionary. Figure 5 illustrates that by showing the PSNR of a noisy image reconstruction depending on the number of iterations, compared to full AM dictionary learning with 103 iterations. As for synthetic data, optimal performance is reached very fast. In this particular case, the model converges after 80 seconds with approximate DL unrolled for 20 iterations of FISTA compared to 600 seconds in the case of standard DL. Note that the speed rate highly depends on the value of λ. Higher values of λ tend to make FISTA converge faster, and unrolling becomes unnecessary in this case. On the contrary, unrolling is more efficient than AM for lower values of λ.
Loss landscape. The ability of gradient descent to find adequate local minima strongly depends on the structure of the problem. To quantify this, we evaluate the variation of PSNR depending on the Signal to Noise Ratio (SNR) (10 log10 (σ
2/σ2b) where σ2b is the variance of the noise) for 50 random initializations in the context of convolutional dictionary learning on a task of image denoising, with 20 unrolled iterations. Figure 5 shows that approximate CDL is robust to random initialization when the level of noise is not too high. In this case, all local minima are similar in terms of reconstruction quality. We provide a visualization of the loss landscape with the help of ideas presented in Li et al. (2018). The algorithm computes a minimum, and we chose two properly rescaled vectors to create a plan from this minimum. The 3D landscape is displayed on this plan in Figure B using the Python library K3D-Jupyter2. We also compare in Figure 5 (right) the shapes of local minima in 1D by computing the values of the loss along a line between two local minima. These visualizations confirm that dictionary learning locally behaves like a convex function with similar local minima.
3.2 STOCHASTIC DDL
In order to apply DDL in realistic settings, it is tempting to adapt Stochastic Gradient Descent (SGD), commonly used for neural networks. The major advantage is that the sparse coding is not performed on all data at each forward pass, leading to significant time and memory savings. The issue is that the choice of gradient steps is critical to the optimization process in dictionary learning, and SGD methods based on simple heuristics like rate decay are difficult to tune in this context. We propose to leverage a new optimization scheme introduced in Vaswani et al. (2019), which consists of performing a stochastic line search. The algorithm computes a good step size at each epoch, after which a heuristic decreases the maximal step. Figure 6 displays the recovery score function of the time for various mini-
batch sizes on a problem with 105 samples. The data were generated as in Figure 1 but with a larger dictionary (50×100). The algorithm achieves good performance with small mini-batches and thus limited memory usage. We also compare this method with Online dictionary learning (Mairal et al., 2009) in Figure E. It shows that our method speeds up the dictionary recovery, especially for lower values of λ. This strategy can be adapted very easily for convolutional models by taking sub-windows of the full signal and performing a stochastic line search, as demonstrated in Section 4. See Tolooshams et al. (2020) for another unrolled stochastic CDL algorithm applied to medical data.

heading: 4 APPLICATION TO PATTERN LEARNING IN MEG SIGNALS
text: In magnetoencephalography (MEG), the measurements over the scalp consist of hundreds of simultaneous recordings, which provide information on the neural activity during a large period. Convolutional dictionary learning makes it possible to learn cognitive patterns corresponding to physiological activities (Dupré la Tour et al., 2018). As the electromagnetic waves propagate through the brain at the speed of light, every sensor measures the same waveform simultaneously but not at the same
2Package available at https://github.com/K3D-tools/K3D-jupyter.
Spatial pattern 0 Spatial pattern 1 Spatial pattern 2
0.0 0.5 1.0 Time (s)
0.5
0.0
Temporal pattern 0
0.0 0.5 1.0 Time (s)
0.2
0.0
Temporal pattern 1
0.0 0.5 1.0 Time (s)
0.0
0.2
0.4 Temporal pattern 2
Figure 7: Stochastic Deep CDL on 6 minutes of MEG data (204 channels, sampling rate of 150Hz). The algorithm uses 40 atoms, 30 unrolled iterations and 100 iterations with batch size 20. We recover heartbeat (0), blinking (1) artifacts, and an auditory evoked response (2) among others.
Minibatch Time window Steps learning Corr. u Corr. v Mean corr. Time 5 20 s True 0.85 ± 0.02 0.84 ± 0.06 0.845 110 s 5 20 s False 0.88 ± 0.02 0.78 ± 0.06 0.83 57 s 5 10 s True 0.83 ± 0.01 0.82 ± 0.09 0.825 56 s
20 10 s True 0.85 ± 0.01 0.75 ± 0.09 0.80 163 s
Table 1: Stochastic Deep CDL on MEG data (as in Figure 7). We compare u and v to 12 important atoms output by alphacsc (correlation averaged on 5 runs), depending on several hyperparameters, with 30 layers, 10 epochs and 10 iterations per epochs. λrescaled = 0.3λmax, λmax = ∥∥DT y∥∥∞. The best setups achieve 80% – 90% average correlation with alphacsc in around 100 sec. compared to around 1400 sec. Our method is also faster than convolutional K-SVD (Yellin et al., 2017).
intensity. The authors propose to rely on multivariate convolutional sparse coding (CSC) with rank1 constraint to leverage this physical property and learn prototypical patterns. In this case, space and time patterns are disjoint in each atom: Dk = ukvTk where u gathers the spatial activations on each channel and v corresponds to the temporal pattern. This leads to the model
min zk∈RT ,uk∈RS ,vk∈Rt
1
2 ∥∥∥∥∥ n∑ k=1 (ukv > k ) ∗ zk − y ∥∥∥∥∥ 2
2
+ λ n∑ k=1 ‖zk‖1 , (6)
where n is the number of atoms, T is the total recording time, t is the kernel size, and S is the number of sensors. We propose to learn u and v with Stochastic Deep CDL unrolled for a few iterations to speed up the computations of the atoms. Figure 7 reproduces the multivariate CSC experiments of alphacsc3 (Dupré la Tour et al., 2018) on the dataset sample of MNE (Gramfort et al., 2013) – 6 minutes of recordings with 204 channels sampled at 150Hz with visual and audio stimuli.
The algorithm recovers the main waveforms and spatial patterns with approximate sparse codes and without performing the sparse coding on the whole data set at each gradient iteration, which leads to a significant gain of time. We are able to distinguish several meaningful patterns as heartbeat and blinking artifacts or auditive evoked response. As this problem is unsupervised, it is difficult to provide robust quantitative quality measurements. Therefore, we compare our patterns to 12 important patterns recovered by alpahcsc in terms of correlation in Table 1. Good setups achieve between 80% and 90% average correlation ten times faster.

heading: 5 CONCLUSION
text: Dictionary learning is an efficient technique to learn patterns in a signal but is challenging to apply to large real-world problems. This work showed that approximate dictionary learning, which consists in replacing the optimal solution of the Lasso with a time-efficient approximation, offers a valuable trade-off between computational cost and quality of the solution compared to complete Alternating Minimization. This method, combined with a well-suited stochastic gradient descent algorithm, scales up to large data sets, as demonstrated on a MEG pattern learning problem. This work provided a theoretical study of the asymptotic behavior of unrolling in approximate dictionary learning. In particular, we showed that numerical instabilities make DDL usage inefficient when too many iterations are unrolled. However, the super-efficiency of DDL in the first iterations remains unexplained, and our first findings would benefit from theoretical support.
3Package and experiments available at https://alphacsc.github.io

heading: ETHICS STATEMENT
text: The MEG data conform to ethic guidelines (no individual names, collected under individual’s consent, . . . ).

heading: REPRODUCIBILITY STATEMENT
text: Code is available at https://github.com/bmalezieux/unrolled_dl. We provide a full description of all our experiments in Appendix A, and the proofs of our theoretical results in Appendix C.

heading: ACKNOWLEDGMENTS
text: This work was supported by grants from Digiteo France.

heading: A FULL DESCRIPTION OF THE EXPERIMENTS
text: This section provides complementary information on the experiments presented in the paper.

heading: A.1 CONVERGENCE OF THE JACOBIANS - FIGURE 1 AND FIGURE 2
text: We generate a normalized random Gaussian dictionary D of dimension 30 × 50, and sparse codes z from a Bernoulli Gaussian distribution of sparsity 0.3 and σ2 = 1. The signal to process is y = Dz + b where b is an additive Gaussian noise with σ2noise = 0.1. The Jacobians are computed for a random perturbation D+bD of D where bD is a Gaussian noise of scale 0.5σ2D. J N l corresponds to the approximate Jacobian with N iterations of ISTA with λ = 0.1. J∗l corresponds the true Jacobian computed with sparse codes obtained after 104 iterations of ISTA with λ = 0.1.
In Figure 2, the norm ∥∥JNl − J∗l ∥∥ is computed for 50 samples.

heading: A.2 CONVERGENCE OF THE GRADIENT ESTIMATES - FIGURE 3
text: Synthetic data. We generate a normalized random Gaussian dictionary D of dimension 30× 50, and 1000 sparse codes z from a Bernoulli Gaussian distribution of sparsity 0.3 and σ2 = 1. The signal to process is y = Dz + b where b is an additive Gaussian noise with σ2noise = 0.1. The gradients are computed for a random perturbation D + bD of D where bD is a Gaussian noise of scale 0.5σ2D.
Noisy image. A 128× 128 black-and-white image is degraded by a Gaussian noise with σ2noise = 0.1 and normalized. We processed 1000 patches of dimension 10 × 10 from the image, and we computed the gradients for a dictionary composed of 128 random patches.
gN corresponds to the gradient for N iterations of FISTA with λ = 0.1. g∗ corresponds to the true gradient computed with a sparse code obtained after 104 iterations of FISTA.

heading: A.3 OPTIMIZATION DYNAMICS ON SYNTHETIC DATA - FIGURE 4
text: We generate a normalized random Gaussian dictionary D of dimension 30 × 50, and sparse codes z from a Bernoulli Gaussian distribution of sparsity 0.3 and σ2 = 1. The signal to process is y = Dz + b where b is an additive Gaussian noise with σ2noise = 0.1. The initial dictionary
is taken as a random perturbation D + bD of D where bD is a Gaussian noise of scale 0.5σ2D. N corresponds to the number of unrolled iterations of FISTA. F ∗ is the value of the loss for 103 iterations minus 10−3. S∗ is the score obtained after 103 iterations plus 10−3. The optimization is done with λ = 0.1. We compare the number of gradient steps (left), the loss values (center), and the recovery scores (right) for 50 different dictionaries. DDL with steps sizes learning is evaluated on 100 iterations only due to memory and optimization time issues.

heading: A.4 OPTIMIZATION DYNAMICS AND LOSS LANDSCAPES ON IMAGES - FIGURE 5
text: A 128× 128 black-and-white image is degraded by a Gaussian noise and normalized.
Left. In this experiment, σ2noise = 0.1. We learn a dictionary composed of 128 atoms on 10 × 10 patches with FISTA and λ = 0.1 in all cases. The PSNR is obtained with sparse codes output by the network. The results are compared to the truth with the Peak Signal to Noise Ratio. Dictionary learning denoising with 1000 iterations of FISTA is taken as a baseline.
Center. We learn 50 dictionaries from 50 random initializations in convolutional dictionary learning with 50 kernels of size 8 × 8 with 20 unrolled iterations of FISTA and λ = 0.1. The PSNR is obtained with sparse codes output by the network. We compare the average, minimal and maximal PSNR, and recovery scores with all other dictionaries to study the robustness to random initialization depending on the level of noise (SNR).
Right. In this experiment, σ2noise = 0.1. We learn 2 dictionaries from 2 random initializations in convolutional dictionary learning with 50 kernels of size 8 × 8 with 20 unrolled iterations of FISTA and λ = 0.1. We display the loss values on the line between these two dictionaries. The experiment is repeated on 10 different random initializations.

heading: A.5 STOCHASTIC DDL ON SYNTHETIC DATA - FIGURE 6
text: We generate a normalized random Gaussian dictionary D of dimension 50 × 100, and 105 sparse codes z from a Bernoulli Gaussian distribution of sparsity 0.3 and σ2 = 1. The signal to process is y = Dz + b where b is an additive Gaussian noise with σ2noise = 0.1. The initial dictionary is taken as a random gaussian dictionary. We compare stochastic and full-batch line search projected gradient descent with 30 unrolled iterations of FISTA and λ = 0.1, without steps sizes learning. Stochastic DDL is run for 10 epochs with a maximum of 100 iterations for each epoch.

heading: A.6 PATTERN LEARNING IN MEG - FIGURE 7
text: Stochastic Deep CDL on 6 minutes of recordings of MEG data with 204 channels and a sampling rate of 150Hz. We remove the powerline artifacts and high-pass filter the signal to remove the drift which can impact the CSC technique. The signal is also resampled to 150 Hz to reduce the computational burden. This preprocessing procedure is presented in alphacsc, and available in the code in the supplementary materials. The algorithm learns 40 atoms of 1 second on mini batches of 10 seconds, with 30 unrolled iterations of FISTA, λscaled = 0.3, and 10 epochs with 10 iterations per epoch. The number of mini-batches per iteration is 20, with possible overlap.

heading: B EXTRA FIGURES AND EXPERIMENTAL RESULTS
text: LISTA - Figure A. Illustration of LISTA for Dictionary Learning with initialization Z0 = 0 for N = 3. W 1D = 1 L (D) >, W 2D = (I − 1L (D) >D), where L = ‖D‖2. The result ZN (D) output by the network is an approximation of the solution of the LASSO.
Loss landscape in 2D - Figure B. We provide a visualization of the loss landscape with the help of ideas presented in Li et al. (2018). The algorithm computes a minimum, and we chose two properly rescaled vectors to create a plan from this minimum. The 3D landscape is displayed on this plan in the appendix using the Python library K3D-Jupyter. This visualization and the visualization in 1D confirm that (approximate) dictionary learning locally behaves like a convex function with smooth local minima.
y
W 1D
W 2D
W 1D
W 2D
W 1D
zN (D)
Figure A: LISTA
Figure B: Loss landscape in approximate CDL
Gradient convergence in norm - Figure C. Gradient estimates convergence in norm for synthetic data (left) and patches from a noisy image (right). The setup is similar to Figure 3. Both gradient estimates converge smoothly in early iterations. When the back-propagation goes too deep, the performance of g2N decreases compared to g 1 N , and we observe large numerical instabilities. This behavior is coherent with the Jacobian convergence patterns studied in Proposition 2.5. Once on the support, g2N reaches back the performance of g 1 N .
100 101 102 103
Iterations N
10-1
101
103
||g ∗ − g| |
Gaussian dictionary
BP depth 20 50 full AM
100 101 102 103
Iterations N
101
103
||g ∗ − g| |
Noisy image
Figure C: Gradient estimates convergence in norm for synthetic data (left) and patches from a noisy image (right). Both gradient estimates converge smoothly in early iterations, after what DDL gradient becomes unstable. The behavior returns to normal once the algorithm reaches the support.
Computation time to reach 0.95 recovery score - Figure D. The setup is similar to Figure 6. A random Gaussian dictionary of size 50× 100 generates the data from 105 sparse codes with sparsity 0.3. The approximate sparse coding is solved with λ = 0.1 and 30 unrolled iterations of FISTA. The algorithm achieves good performances with small mini-batches and thus limited memory usage. Stochastic DDL can process large amounts of data and recovers good quality dictionaries faster than full batch DDL.
Sto DDL vs. Online DL - Figure E. We compare the time Online DL from spams4 (Mairal et al., 2009) and Stochastic DDL need to reach a recovery score of 0.95 with a batch size of 2000. Online DL is run with 10 threads. We repeat the experiment 10 times for different values of λ from
4package available at http://thoth.inrialpes.fr/people/mairal/spams/
100 500 2000 10000 DDL Oracle DL Minibatch size
0
10
20
30
40
50
Ti m
e (s
)
Ti m
eo ut
Time to reach a score of 0.95
Figure D: Time to reach a recovery score of 0.95. Intermediate batch sizes offer a good trade-off between speed and memory usage compared to full-batch DDL.
0.1 to 1.0. The setup is similar to Figure D, and we initialize both methods randomly. Stochastic DDL is more efficient for smaller values of λ, due to the fact that sparse coding is slower in this case. For higher values of λ, both methods are equivalent. Another advantage of Stochastic DDL is its modularity. It works on various kinds of dictionary parameterization thanks to automatic differentiation, as illustrated on 1-rank multivariate convolutional dictionary learning in Figure 7.
0.2 0.4 0.6 0.8 1.0 λ
101
102
Ti m
e (s
)
Time to reach a score of 0.95
Sto. DDL Online DL
Figure E: Comparison between Online DL and Stochastic DDL. Stochastic DDL is more efficient for smaller values of λ, due to the fact that sparse coding is slower in this case.

heading: C PROOFS OF THEORETICAL RESULTS
text: This section gives the proofs for the various theoretical results in the paper.

heading: C.1 PROOF OF PROPOSITION 2.1.
text: Proposition 2.1 Let D∗ = arg minD∈C G(D) and D∗N = arg minD∈C GN (D), where N is the number of unrolled iterations. We denote by K(D∗) a constant depending on D∗, and by C(N) the convergence speed of the algorithm, which approximates the inner problem solution. We have
GN (D ∗ N )−G(D∗) ≤ K(D∗)C(N) .
Let G(D) , F (Z∗(D),D) and GN (D) , F (ZN (D),D) where Z∗(D) = arg minZ∈Rn×T F (Z,D) and ZN (D) = FISTA(D, N). Let D
∗ = arg minD∈C G(D) and D∗N = arg minD∈C GN (D). We have
GN (D ∗ N )−G(D∗) = GN (D∗N )−GN (D∗) +GN (D∗)−G(D∗) (7)
= F (ZN (DN ),DN )− F (ZN (D∗),D∗) (8) + F (ZN (D ∗),D∗)− F (Z(D∗),D∗) (9)
By definition of D∗N F (ZN (D ∗ N ),D ∗ N )− F (ZN (D∗),D∗) ≤ 0 (10)
The convergence rate of FISTA in function value for a fixed dictionary D is
F (ZN (D),D)− F (ZN (D),D) ≤ K(D)
N2 (11)
Therefore
F (ZN (D ∗),D∗)− F (Z(D∗),D∗) ≤ K(D
∗)
N2 (12)
Hence
GN (D ∗ N )−G(D∗) ≤
K(D∗)
N2 (13)

heading: C.2 PROOF OF PROPOSITION 2.2
text: Proposition 2.2 Let D ∈ Rm×n. Then, there exists a constant L1 > 0 such that for every number of iterations N ∥∥g1N − g∗∥∥ ≤ L1 ‖zN (D)− z∗(D)‖ . We have
F (z,D) = 1
2 ‖Dz − y‖22 + λ ‖z‖1 (14)
∇2F (z,D) = (Dz − y)z> (15)
z0(D) = 0 and the iterates (zN (D))N∈N converge towards z∗(D). Hence, they are contained in a closed ball around z∗(D). As ∇2F (·,D) is continuously differentiable, it is locally Lipschitz on this closed ball, and there exists a constant L1(D) depending on D such that∥∥g1N − g∗∥∥ = ‖∇2F (zN (D),D)−∇2F (z∗(D),D)‖ (16)
≤ L1(D) ‖zN (D)− z∗(D)‖ (17)

heading: C.3 PROOF OF PROPOSITION 2.3.
text: Proposition 2.3 Let D ∈ Rm×n. Let S∗ be the support of z∗(D), SN be the support of zN and S̃N = SN ∪ S∗. Let f(z,D) = 12 ‖Dz − y‖ 2 2 be the data-fitting term in F . Let R(J, S̃) =
J+ ( ∇21,1f(z∗,D) 1S̃ ) +∇22,1f(z∗,D) 1S̃ . Then there exists a constant L2 > 0 and a subsequence of (F)ISTA iterates zφ(N) such that for all N ∈ N:
∃ g2φ(N) ∈ ∇2f(zφ(N),D) + J + φ(N) ( ∇1f(zφ(N),D) + λ∂‖·‖1(zφ(N)) ) s.t. :∥∥∥g2φ(N) − g∗∥∥∥ ≤ ∥∥∥R(Jφ(N), S̃φ(N))∥∥∥∥∥zφ(N) − z∗∥∥+ L22 ∥∥zφ(N) − z∗∥∥2 .
This sub-sequence zφ(N) corresponds to iterates on the support of z∗.
We have g2N (D) ∈ ∇2f(zN (D),D) + J + N ( ∇1f(zN (D),D) + λ∂‖·‖1(zN ) ) (18)
We adapt equation (6) in Ablin et al. (2020)
g2N = g ∗ +R(JN , S̃N )(zN − z∗) +RD,zN + J + NR z,z N (19)
where
R(J, S̃) = J+ ( ∇21,1f(z∗,D) 1S̃ ) +∇22,1f(z∗,D) 1S̃ (20)
RD,zN = ∇2f(zN ,D)−∇2f(z ∗,D)−∇22,1f(z∗,D)(zN − z∗) (21)
Rz,zN ∈ ∇1f(zN ,D) + λ∂‖·‖1(zN )−∇ 2 1,1f(z ∗,D)(zN − z∗) (22)
As zN and z∗ are on S̃N ∇22,1f(z∗,D)(zN − z∗) = ( ∇22,1f(z∗,D) 1S̃N ) (zN − z∗) (23)
J+ ( ∇21,1f(z∗,D)(zN − z∗) ) = J+ ( ∇21,1f(z∗,D) 1S̃N (zN − z ∗) )
(24)
As stated in Proposition 2.2, ∇2f(·,D) is locally Lipschitz, and RD,zN is the Taylor rest of ∇2f(·,D). Therefore, there exists a constant LD,z such that
∀N ∈ N, ∥∥∥RD,zN ∥∥∥ ≤ LD,z2 ‖zN (D)− z∗(D)‖2 (25)
We know that 0 ∈ ∇1f(z∗,D)+λ∂‖·‖1(z ∗). In other words, ∃u∗ ∈ λ∂‖·‖1(z ∗) s.t.∇1f(z∗,D)+ u∗ = 0. Therefore we have:
Rz,zN ∈ ∇1f(zN ,D)−∇1f(z ∗,D)−∇21,1f(z∗, x)(zN − z∗) + λ∂ ‖zN‖1 − u ∗ (26)
Let Lz,z be the Lipschitz constant of ∇1f(·,D). (F)ISTA outputs a sequence such that there exists a sub-sequence (zφ(N))N∈N which has the same support as z∗. For this sub-sequence, u∗ ∈ λ∂‖·‖1(zφ(N)). Therefore, there exists R z,z φ(N) such that
1. Rz,zφ(N) ∈ ∇1f(zφ(N),D) + λ∂‖·‖1(zφ(N))−∇ 2 1,1f(z ∗, x)(zφ(N) − z∗) 2. ∥∥∥Rz,zφ(N)∥∥∥ ≤ Lz,z2 ∥∥zφ(N) − z∗∥∥2
For this sub-sequence, we can adapt Proposition 2 from Ablin et al. (2020). Let L2 = LD,z +Lz,z , we have
∃ g2φ(N) ∈ ∇2f(zφ(N),D) + Jφ(N) ( ∇1f(zφ(N),D) + λ∂ ∥∥zφ(N)∥∥1 ), s.t. : (27)∥∥∥g2φ(N) − g∗∥∥∥ ≤ ∥∥∥R(Jφ(N), S̃φ(N))∥∥∥∥∥zφ(N) − z∗∥∥+ L22 ∥∥zφ(N) − z∗∥∥2 (28)

heading: C.4 PROOF OF THEOREM 2.4.
text: Theorem 2.4 At iteration N + 1 of ISTA, the weak Jacobian of zN+1 relatively to Dl, where Dl is the l-th row of D, is given by induction:
∂(zN+1)
∂Dl = 1|zN+1|>0
( ∂(zN )
∂Dl − 1 L
( Dlz > N + (D > l zN − yl)In + D>D ∂(zN )
∂Dl
)) .
∂(zN ) ∂Dl will be denoted by JNl . It converges towards the weak Jacobian J ∗ l of z ∗ relatively to Dl, whose values are
J∗l S∗ = −(D > :,S∗D:,S∗) −1(Dlz ∗> + (D>l z ∗ − yl)In)S∗ ,
on the support S∗ of z∗, and 0 elsewhere. Moreover, R(J∗, S∗) = 0.
We start by recalling a Lemma from Deledalle et al. (2014).
Lemma C.1 The soft-thresholding STµ defined by STµ(z) = sgn(z) (|z| − µ)+ is weakly differentiable with weak derivative dSTµ(z)dz = 1|z|>µ.
Coordinate-wise, ISTA corresponds to the following equality:
zN+1 = STµ((I − 1
L D>D)zN +
1 L D>y) (29)
(zN+1)i = STµ((zN )i − 1
L m∑ p=1 ( n∑ j=1 DjiDjp)(zN )p + 1 L n∑ j=1 Djiyj) (30)
The Jacobian is computed coordinate wise with the chain rule:
∂(zN+1)i ∂Dlk = 1|(zN+1)i|>0 · ( ∂(zN )i ∂Dlk − 1 L ∂ ∂Dlk ( m∑ p=1 ( n∑ j=1 DjiDjp)(zN )p) + 1 L ∂ ∂Dlk n∑ j=1 Djiyj))
(31) Last term:
∂
∂Dlk n∑ j=1 Djiyj = δikyl (32)
Second term:
∂
∂Dlk m∑ p=1 n∑ j=1 DjiDjp(zN )p = m∑ p=1 n∑ j=1 DjiDjp ∂(zN )p ∂Dlk + m∑ p=1 n∑ j=1 ∂DjiDjp ∂Dlk (zN )p (33)
∂DjiDjp ∂Dlk =  2Dlk if j = l and i = p = k Dlp if j = l and i = k and p 6= k Dli if j = l and i 6= k and p = k 0 else
(34)
Therefore: m∑ p=1 n∑ j=1 ∂DjiDjp ∂Dlk (zN )p = m∑ p=1 (2Dlkδipδik +Dliδpk1i6=k +Dlpδik1k 6=p)(zN )p (35)
= 2Dlk(zN )kδik +Dli(zN )k1i6=k + m∑ p=1 p 6=k Dlp(zN )pδik (36)
= Dli(zN )k + δik m∑ p=1 Dlp(zN )p (37)
Hence:
∂(zN+1)i ∂Dlk = 1|(zN+1)i|>0 · (∂(zN )i ∂Dlk − 1 L (Dli(zN )k+ (38)
δik( m∑ p=1 Dlp(zN )p) + m∑ p=1 n∑ j=1 ∂(zN )p ∂Dlk DjiDjp − δikyl) )
This leads to the following vector formulation:
∂(zN+1)
∂Dl = 1|zN+1|>0
( ∂(zN )
∂Dl − 1 L
( Dlz > N + (D > l zN − yl)Im + D>D ∂(zN )
∂Dl
)) (39)
On the support of z∗, denoted by S∗, this quantity converges towards the fixed point:
J∗l = −(D>:,S∗D:,S∗)−1(Dlz∗ > + (D>l z ∗ − yl)Im)S∗ (40)
Elsewhere, J∗l is equal to 0. To prove that R(J ∗, S∗) = 0, we use the expression given by equation 39
J∗ = 1S∗ (
J∗ − 1 L
( ∇22,1f(z∗,Dl)> +∇21,1f(z∗,D)>J ∗)) (41) J∗ − 1S∗ J∗ = 1
L 1S∗ ∇22,1f(z∗,Dl)> + 1S∗ ∇21,1f(z∗,D)>J ∗ (42)
0 = J∗+ ( ∇21,1f(z∗,D) 1S∗ ) +∇22,1f(z∗,D) 1S∗ (43) 0 = R(J∗, S∗) (44)

heading: C.5 PROOF OF PROPOSITION 2.5 AND COROLLARY 2.6
text: Proposition 2.5 Let N be the number of iterations and K be the back-propagation depth. We assume that ∀n ≥ N−K, S∗ ⊂ Sn. Let ĒN = Sn\S∗, let L be the largest eigenvalue ofD>:,S∗D:,S∗ , and let µn be the smallest eigenvalue of D>:,SnD:,Sn−1 . Let Bn = ∥∥∥PEn −D>:,ĒnD†>:,S∗PS∗∥∥∥, where PS is the projection on RS and D† is the pseudo-inverse of D. We have ∥∥JNl − J∗l ∥∥ ≤ K∏ k=1 ( 1− µN−k L ) ‖J∗l ‖+ 2 L ‖Dl‖ K−1∑ k=0 k∏ i=1 (1−µN−i L ) (∥∥zN−kl − z∗l ∥∥+BN−k ‖z∗l ‖) .
We denote by G the matrix (I − 1LD >D). For zN with support SN and z∗ with support S∗, we have with the induction in Theorem 2.4
JNl,SN = ( GJN−1l + u N−1 l ) SN
(45) J∗l,S∗ = ( GJ∗l + u ∗ l ) S∗ (46)
where uNl = − 1L ( Dlz > N + (D > l zN − yl)I ) and the other terms on S̄N and S̄∗ are 0. We can thus decompose their difference as the sum of two terms, one on the support S∗ and one on this complement EN = SN \ S∗
J∗l − JNl = (J∗l − JNl )S∗ + (J∗l − JNl )EN .
Recall that we assume S∗ ⊂ SN . Let’s study the terms separately on S∗ and EN = SN \S∗. These two terms can be decompose again to constitute a double recursion system,
(JNl − J∗l )S∗ = GS∗(JN−1l − J ∗ l ) + (u N−1 l − u ∗ l )S∗ (47)
= GS∗,S∗(JN−1l − J ∗ l )S∗ +GS∗,EN−1(J N−1 l − J ∗)EN−1 + (u N−1 l − u ∗ l )S∗ ,
(48)
(JNl − J∗l )EN = (J N l )EN = GEN (J N−1 l − J ∗ l ) +GEN ,S∗J ∗ l + (u N−1 l )EN (49)
= GEN ,S∗(J N−1 l − J ∗ l )S∗ +GEN ,EN−1(J N−1 l − J ∗ l )EN−1 (50)
+ (uN−1l − u ∗ l )EN + ( (u∗l )EN −D > :,EN D:,S∗(D > :,S∗D:,S∗) −1(u∗l )S∗ ) .
We define as PSN ,EN the operator which projects a vector from EN on (SN , EN ) with zeros on SN . As S∗ ∪ EN = SN , we get by combining these two expressions,
(JNl − J∗l )SN =GSN ,SN−1(JN−1l − J ∗ l )SN−1 + (u N−1 l − u ∗ l )SN (51) + PSN ,EN ( (u∗l )EN −D > :,EN D:,S∗(D > :,S∗D:,S∗) −1(u∗l )S∗ )
Taking the norm yields to the following inequality,∥∥JNl − J∗l ∥∥ ≤∥∥GSN ,SN−1∥∥∥∥JN−1l − J∗l ∥∥+ ∥∥uN−1l − u∗l ∥∥ (52) + ∥∥∥(u∗l )EN −D>:,END:,S∗(D>:,S∗D:,S∗)−1(u∗l )S∗∥∥∥ .
Denoting by µN the smallest eigenvalue of D>:,SND:,SN−1 , then ∥∥GSN ,SN−1∥∥ = (1 − µNL ) and we get that ∥∥JNl − J∗l ∥∥ ≤ K∏ k=1 (1− µN−k L ) ∥∥JN−Kl − J∗l ∥∥ (53)
+ K−1∑ k=0 k∏ i=1 (1− µN−i L ) (∥∥uN−kl − u∗l ∥∥+ ∥∥∥(u∗l )EN−k −D>:,EN−kD†>:,S∗(u∗l )S∗∥∥∥) .
The back-propagation is initialized as JN−Kl = 0. Therefore ∥∥JN−Kl − J∗l ∥∥ = ‖J∗l ‖. More-
over ∥∥uN−kl − u∗l ∥∥ ≤ 2L ‖Dl‖ ∥∥zN−kl − z∗l ∥∥. Finally, ∥∥∥(u∗l )EN−k −D>:,EN−kD†>:,S∗(u∗l )S∗∥∥∥ can be rewritten with projection matrices PEN−k and PS̄∗ to obtain∥∥∥(u∗l )EN−k −D>:,EN−kD†>:,S∗(u∗l )S∗∥∥∥ ≤∥∥∥PEN−ku∗l −D>:,EN−kD†>:,S∗PS∗u∗l ∥∥∥ (54) ≤ ∥∥∥PEN−k −D>:,EN−kD†>:,S∗PS∗∥∥∥ ‖u∗l ‖ (55) ≤ ∥∥∥PEN−k −D>:,EN−kD†>:,S∗PS∗∥∥∥ 2L ‖Dl‖ ‖z∗l ‖ . (56)

heading: Let BN−k =
text: ∥∥∥PEN−k −D>:,EN−kD†>:,S∗PS∗∥∥∥. We have∥∥JNl − J∗l ∥∥ ≤ K∏ k=1 (1−µN−k L ) ‖J∗l ‖+ 2 L ‖Dl‖ K−1∑ k=0 k∏ i=1 (1−µN−i L ) (∥∥zN−kl − z∗l ∥∥+BN−k ‖z∗l ‖) . (57) We now suppose that the support is reached at iteration N − s, with s ≥ K. Therefore, ∀n ∈ [N − s,N ] Sn = S∗. Let ∆n = F (zn, D) − F (z∗, D) + L2 ‖zn − z
∗‖. On the support, F is a µ-strongly convex function and the convergence rate of (zN ) is
‖z∗ − zN‖ ≤ ( 1− µ
L )s 2∆N−s L
(58)
Thus, we obtain∥∥JNl − J∗l ∥∥ ≤ K∏ k=1 (1− µN−k L ) ‖J∗l ‖ (59)
+ 2
L ‖Dl‖ K−1∑ k=0 k∏ i=1 (1− µN−i L ) (∥∥zN−kl − z∗l ∥∥+BN−k ‖u∗l ‖)
≤ K∏ k=1 (1− µN−k L ) ‖J∗l ‖ (60)
+ 2
L ‖Dl‖ s−1∑ k=0 (1− µ L )k (∥∥zN−kl − z∗l ∥∥)
+ 2
L ‖Dl‖ (1−
µ L )s K−1∑ k=s−1 k∏ i=s−1 (1− µN−i L ) (∥∥zN−kl − z∗l ∥∥+BN−k ‖(u∗l )‖)
≤ K∏ k=1 (1− µN−k L ) ‖J∗l ‖ (61)
+ 2
L ‖Dl‖ s−1∑ k=0 (1− µ L )k ( 1− µ L )s−1−k 2∆N−s L
+ 2
L ‖Dl‖ (1−
µ L )s K−1∑ k=s−1 k∏ i=s−1 (1− µN−i L ) (∥∥zN−kl − z∗l ∥∥+BN−k ‖(u∗l )‖)
≤ K∏ k=1 (1− µN−k L ) ‖J∗l ‖ (62)
+ ‖Dl‖ (1− µ
L )s−1s 4∆N−s L2
+ 2
L ‖Dl‖ (1−
µ L )s K−1∑ k=s−1 k∏ i=s−1 (1− µN−i L ) (∥∥zN−kl − z∗l ∥∥+BN−k ‖(u∗l )‖)
(63)
Corollary 2.6 Let µ > 0 be the smallest eigenvalue of D>:,S∗D:,S∗ . Let K ≤ N be the backpropagation depth and let ∆N = F (zN , D) − F (z∗, D) + L2 ‖zN − z
∗‖. Suppose that ∀n ∈ [N −K,N ]; Sn ⊂ S∗. Then, we have∥∥J∗l − JNl ∥∥ ≤ (1− µL)K ‖J∗l ‖+K (1− µL)K−1 ‖Dl‖ 4∆N−KL2 . The term 2L ‖Dl‖ (1− µ L ) s ∑K−1 k=s−1 ∏k i=s−1(1− µN−i L )
(∥∥zN−kl − z∗l ∥∥+BN−k ‖(u∗l )‖) vanishes when the algorithm is initialized on the support. Otherwise, it goes to 0 as s,K → N and N →∞ because ∀n > N − s, µn = µ < 1.
D ITERATIVE ALGORITHMS FOR SPARSE CODING RESOLUTION.
ISTA. Algorithm to solve minz 12 ‖y −Dz‖ 2 2 + λ ‖z‖1
Algorithm 1 ISTA y, D, λ, N z0 = 0, n = 0 Compute the Lipschitz constant L of D>D while n < N do un+1 ← zN − 1LD
>(Dzn − y) zn+1 ← ST λ
L (un+1)
n← n+ 1 end while
FISTA. Algorithm to solve minz 12 ‖y −Dz‖ 2 2 + λ ‖z‖1
Algorithm 2 FISTA y, D, λ, N z0 = x0 = 0, n = 0, t0 = 1 Compute the Lipschitz constant L of D>D while n < N do un+1 ← zn − 1LD
>(Dzn − y) xn+1 ← ST λ
L (un+1)
tn+1 ← 1+ √ 1+4t2n 2 zn+1 ← xn+1 + tn−1tn+1 (xn+1 − xn) n← n+ 1
end while

Review summaries:
The authors theoretically study the performance of dictionary learning using "unrolling" based methods. As opposed to alternating minimization (AM) which switches back and forth between dictionary estimation and sparse recovery, the paper writes down the target dictionary as the solution to a bi-level optimization, where the "inner" optimization is approximated by unrolling with N steps. The main contribution is an approach (along with careful analysis) of computing the subgradient for the outer optimization, and experiments to show that this method works on synthetic and real datasets.

The paper studies dictionary learning where assumes that data can be represented as a linear combination of a few atoms of a matrix called dictionary. Traditionally, one way to approach the problem is to set up a min-min (bi-convex) optimization problem known as lasso or basis pursuit and solve it through alternating minimization (alternate between a sparse coding step and a dictionary update). The paper compares gradient based alternating minimization which uses an analytic gradient (given the code estimate, compute the gradient) to unrolled-based dictionary learning which uses backpropagation (automatic differentiation) through an iterative algorithm estimating the code (inner problem) to compute the gradient for update of the dictionary. This paper borrows results from (Ablin et al., 2020) that had studied min-min optimization problems when the objective functions are smooth, differentiable, and strongly convex. Specifically, solving lasso iteratively through ISTA or FISTA, after support selection, with some assumption on the dictionary, the problem of this paper is reduced to the strongly convex case of (Ablin et al., 2020). Hence, the results follow. Their contribution that makes their paper different from (Ablin et al., 2020) is the study of the Jacobian and instability of the convergence prior to support selection. This very same model (dictionary learning through unrolled algorithms) have been already studied theoretically by [1] (which is missing in the citations) and empirically in the context of convolutional dictionary learning by this paper (Tolooshams et al., 2020) that they cite.

The authors investigate the asymptotic behavior of unrolling applied to dictionary learning. The applicability of unrolling to dictionary learning results from the circumstance that dictionary learning can be reformulated in terms of bilevel optimization, where the lower-level (or inner) problem is a sparse coding problem (in the LASSO case considered here). Unrolling means to replace the argmin in the lower-level problem with the $N$-th iterate of a suitable optimization algorithm. Gradients of the upper-level loss can then be computed by means of backpropagation through algorithmic iterates. The authors study the behavior of resulting gradients depending on $N$ and draw a comparison with alternating minimization. They find that unrolling constitutes a scalable alternative to alternating minimization, where unrolling a relatively small number of iterations or using truncated backpropagation is favorable to ensure stable approximate gradients.