Title: NOODL: PROVABLE ONLINE DICTIONARY LEARNING AND SPARSE CODING

Abstract: We consider the dictionary learning problem, where the aim is to model the given data as a linear combination of a few columns of a matrix known as a dictionary, where the sparse weights forming the linear combination are known as coefficients. Since the dictionary and coefficients, parameterizing the linear model are unknown, the corresponding optimization is inherently non-convex. This was a major challenge until recently, when provable algorithms for dictionary learning were proposed. Yet, these provide guarantees only on the recovery of the dictionary, without explicit recovery guarantees on the coefficients. Moreover, any estimation error in the dictionary adversely impacts the ability to successfully localize and estimate the coefficients. This potentially limits the utility of existing provable dictionary learning methods in applications where coefficient recovery is of interest. To this end, we develop NOODL: a simple Neurally plausible alternating Optimization-based Online Dictionary Learning algorithm, which recovers both the dictionary and coefficients exactly at a geometric rate, when initialized appropriately. Our algorithm, NOODL, is also scalable and amenable for large scale distributed implementations in neural architectures, by which we mean that it only involves simple linear and non-linear operations. Finally, we corroborate these theoretical results via experimental evaluation of the proposed algorithm with the current state-of-the-art techniques. Table 1: Comparison of provable algorithms for dictionary learning.

1. INTRODUCTION:

Sparse models avoid overfitting by favoring simple yet highly expressive representations. Since signals of interest may not be inherently sparse, expressing them as a sparse linear combination of a few columns of a dictionary is used to exploit the sparsity properties. Of specific interest are overcomplete dictionaries, since they provide a flexible way of capturing the richness of a dataset, while yielding sparse representations that are robust to noise; see Mallat and Zhang (1993) ; Chen et al. (1998) ; Donoho et al. (2006) . In practice however, these dictionaries may not be known, warranting a need to learn such representations -known as dictionary learning (DL) or sparse coding (Olshausen and Field, 1997) . Formally, this entails learning an a priori unknown dictionary A ∈ R n×m and sparse coefficients x * (j) ∈ R m from data samples y (j) ∈ R n generated as y (j) = A * x * (j) , x * (j) 0 ≤ k for all j = 1, 2, . . .

This particular model can also be viewed as an extension of the low-rank model (Pearson, 1901) .

Here, instead of sharing a low-dimensional structure, each data vector can now reside in a separate low-dimensional subspace. Therefore, together the data matrix admits a union-of-subspace model. As a result of this additional flexibility, DL finds applications in a wide range of signal processing and machine learning tasks, such as denoising (Elad and Aharon, 2006) , image inpainting (Mairal et al., 2009) , clustering and classification (Ramirez et al., 2010; Rambhatla and Haupt, 2013; Rambhatla et al., 2016; 2017; 2019b; a) , and analysis of deep learning primitives (Ranzato et al., 2008; Gregor and LeCun, 2010) ; see also Elad (2010) , and references therein.

Notwithstanding the non-convexity of the associated optimization problems (since both factors are unknown), alternating minimization-based dictionary learning techniques have enjoyed significant success in practice. Popular heuristics include regularized least squares-based (Olshausen and Field, 1997; Lee et al., 2007; Mairal et al., 2009; Lewicki and Sejnowski, 2000; Kreutz-Delgado et al., 2003) , and greedy approaches such as the method of optimal directions (MOD) (Engan et al., 1999) and k-SVD (Aharon et al., 2006) . However, dictionary learning, and matrix factorization models in general, are difficult to analyze in theory; see also Li et al. (2016a) .

To this end, motivated from a string of recent theoretical works (Gribonval and Schnass, 2010; Jenatton et al., 2012; Geng and Wright, 2014) , provable algorithms for DL have been proposed recently to explain the success of aforementioned alternating minimization-based algorithms (Agarwal et al., 2014; Arora et al., 2014; 2015) . However, these works exclusively focus on guarantees for dictionary recovery. On the other hand, for applications of DL in tasks such as classification and clusteringwhich rely on coefficient recovery -it is crucial to have guarantees on coefficients recovery as well.

Contrary to conventional prescription, a sparse approximation step after recovery of the dictionary does not help; since any error in the dictionary -which leads to an error-in-variables (EIV) (Fuller, 2009) model for the dictionary -degrades our ability to even recover the support of the coefficients (Wainwright, 2009) . Further, when this error is non-negligible, the existing results guarantee recovery of the sparse coefficients only in 2 -norm sense (Donoho et al., 2006) . As a result, there is a need for scalable dictionary learning techniques with guaranteed recovery of both factors.

1.1. SUMMARY OF OUR CONTRIBUTIONS:

In this work, we present a simple online DL algorithm motivated from the following regularized least squares-based problem, where S(•) is a nonlinear function that promotes sparsity. (P1)

Although our algorithm does not optimize this objective, it leverages the fact that the problem (P1) is convex w.r.t A, given the sparse coefficients {x (j) }. Following this, we recover the dictionary by choosing an appropriate gradient descent-based strategy (Arora et al., 2015; Engan et al., 1999) . To recover the coefficients, we develop an iterative hard thresholding (IHT)-based update step (Haupt and Nowak, 2006; Blumensath and Davies, 2009) , and show that -given an appropriate initial estimate of the dictionary and a mini-batch of p data samples at each iteration t of the online algorithmalternating between this IHT-based update for coefficients, and a gradient descent-based step for the dictionary leads to geometric convergence to the true factors, i.e., x (j) →x * (j) and A

i →A * i as t→∞. In addition to achieving exact recovery of both factors, our algorithm -Neurally plausible alternating Optimization-based Online Dictionary Learning (NOODL) -has linear convergence properties. Furthermore, it is scalable, and involves simple operations, making it an attractive choice for practical DL applications. Our major contributions are summarized as follows:

• Provable coefficient recovery: To the best of our knowledge, this is the first result on exact recovery of the sparse coefficients {x * (j) }, including their support recovery, for the DL problem. The proposed IHT-based strategy to update coefficient under the EIV model, is of independent interest for recovery of the sparse coefficients via IHT, which is challenging even when the dictionary is known; see also Yuan et al. (2016) and Li et al. (2016b) .

• Unbiased estimation of factors and linear convergence: The recovery guarantees on the coefficients also helps us to get rid of the bias incurred by the prior-art in dictionary estimation. Furthermore, our technique geometrically converges to the true factors. • Online nature and neural implementation: The online nature of algorithm, makes it suitable for machine learning applications with streaming data. In addition, the separability of the coefficient update allows for distributed implementations in neural architectures (only involves simple linear and non-linear operations) to solve large-scale problems. To showcase this, we also present a prototype neural implementation of NOODL.

In addition, we also verify these theoretical properties of NOODL through experimental evaluations on synthetic data, and compare its performance with state-of-the-art provable DL techniques.

1.2. RELATED WORKS:

With the success of the alternating minimization-based techniques in practice, a push to study the DL problem began when Gribonval and Schnass (2010) showed that for m = n, the solution pair (A * , X * ) lies at a local minima of the following non-convex optimization program, where X = [x (1) , x (2) , . . . , x (p) ] and Y = [y (1) , y (2) , . . . , y (p) ], with high probability over the randomness of the coefficients, min

EQUATION

Algorithm 1: NOODL: Neurally plausible alternating Optimization-based Online Dictionary Learning.

Input: Fresh data samples y (j) ∈ R n for j ∈ [p] at each iteration t generated as per (1), where

EQUATION

x and τ (r) chosen as per A.5 and A.6. No. of iterations T = Ω(log(1/ T )) and R = Ω(log(1/δ R )), for target tolerances T and δ R . Output: The dictionary A (t) and coefficient estimates x (t) (j) for j ∈ [p] at each iterate t. Initialize: Estimate A (0) , which is ( 0 , 2)-near to A * for 0 = O * (1/ log(n)) for t = 0 to T -1 do Predict: (Estimate Coefficients) for j = 1 to p do Initialize:

x (0) (j) = T C/2 (A (t) y (j) ) (3) for r = 0 to R -1 do Update: x (r+1) (j)

= T τ (r) (x (r) (j) -η (r) x A (t) (A (t) x (r) (j) -y (j) )) (4) end end x (t) (j) := x (R) (j) for j ∈ [p] Learn: (Update Dictionary) Form empirical gradient estimate:

EQUATION

Take a gradient descent step:

EQUATION

Normalize:

A (t+1) i = A (t+1) i / A (t+1) i ∀ i ∈ [m] end

Due to the streaming nature of the incoming data, NOODL takes a mini-batch of p data samples at the t-th iteration of the algorithm, as shown in Algorithm 1. It then proceeds by alternating between two update stages: coefficient estimation ("Predict") and dictionary update ("Learn") as follows.

Predict Stage: For a general data sample y = A * x * , the algorithm begins by forming an initial coefficient estimate x (0) based on a hard thresholding (HT) step as shown in (3), where T τ (z) := z • 1 |z|≥τ for a vector z. Given this initial estimate x (0) , the algorithm iterates over R = Ω(log(1/δ R ))

IHT-based steps (4) to achieve a target tolerance of δ R , such that (1 -η x ) R ≤ δ R . Here, η (r)

x is the learning rate, and τ (r) is the threshold at the r-th iterate of the IHT. In practice, these can be fixed to some constants for all iterations; see A.6 for details. Finally at the end of this stage, we have estimate x (t) 

:= x (R) of x * .

Learn Stage: Using this estimate of the coefficients, we update the dictionary at t-th iteration A (t) by an approximate gradient descent step (6), using the empirical gradient estimate (5) and the learning rate η A = Θ(m/k); see also A.5. Finally, we normalize the columns of the dictionary and continue to the next batch. The running time of each step t of NOODL is therefore O(mnp log(1/δ R )). For a target tolerance of T and δ T , such that

A (T ) i -A * i ≤ T , ∀i ∈ [m] and | x (T ) i -x * i | ≤ δ T we choose T = max(Ω(log(1/ T )), Ω(log( √ k/δ T ))).

NOODL uses an initial HT step and an approximate gradient descent-based strategy as in Arora et al. (2015) . Following which, our IHT-based coefficient update step yields an estimate of the coefficients at each iteration of the online algorithm. Coupled with the guaranteed progress made on the dictionary, this also removes the bias in dictionary estimation. Further, the simultaneous recovery of both factors also avoids an often expensive post-processing step for recovery of the coefficients.

3. MAIN RESULT:

We start by introducing a few important definitions. First, as discussed in the previous section we require that the initial estimate A (0) of the dictionary is ( 0 , 2)-close to A * . In fact, we require this closeness property to hold at each subsequent iteration t, which is a key ingredient in our analysis. This initialization achieves two goals. First, the σ(i)A π(i) -A * i ≤ 0 condition ensures that the signed-support of the coefficients are recovered correctly (with high probability) by the hard thresholding-based coefficient initialization step, where signed-support is defined as follows.

Definition 2. The signed-support of a vector x is defined as sign(x) • supp(x).

Next, the A-A * ≤ 2 A * condition keeps the dictionary estimates close to A * and is used in our analysis to ensure that the gradient direction (5) makes progress. Further, in our analysis, we ensure t (defined as A (t) i -A * i ≤ t ) contracts at every iteration, and assume 0 , t = O * (1/ log(n)). Also, we assume that the dictionary A is fixed (deterministic) and µ-incoherent, defined as follows.

Definition 3. A matrix A ∈ R n×m with unit-norm columns is µ-incoherent if for all i = j the inner-product between the columns of the matrix follow

| A i , A j | ≤ µ/ √ n.

The incoherence parameter measures the degree of closeness of the dictionary elements. Smaller values (i.e., close to 0) of µ are preferred, since they indicate that the dictionary elements do not resemble each other. This helps us to effectively tell dictionary elements apart (Donoho and Huo, 2001; Candès and Romberg, 2007) . We assume that µ = O(log(n)) (Donoho and Huo, 2001) . Next, we assume that the coefficients are drawn from a distribution class D defined as follows.

Definition 4 (Distribution class D). The coefficient vector x * belongs to an unknown distribution D, where the support

S = supp(x * ) is at most of size k, Pr[i ∈ S] = Θ(k/m) and Pr[i, j ∈ S] = Θ(k 2 /m 2 ). Moreover, the distribution is normalized such that E[x * i |i ∈ S] = 0 and E[x * 2 i |i ∈ S] = 1, and when i ∈ S, |x * i | ≥ C for some constant C ≤ 1.

In addition, the non-zero entries are sub-Gaussian and pairwise independent conditioned on the support.

The randomness of the coefficient is necessary for our finite sample analysis of the convergence. Here, there are two sources of randomness. The first is the randomness of the support, where the non-zero elements are assumed to pair-wise independent. The second is the value an element in the support takes, which is assumed to be zero mean with variance one, and bounded in magnitude. Similar conditions are also required for support recovery of sparse coefficients, even when the dictionary is known (Wainwright, 2009; Yuan et al., 2016) . Note that, although we only consider the case |x * i | ≥ C for ease of discussion, analogous results may hold more generally for x * i s drawn from a distribution with sufficiently (exponentially) small probability of taking values in [-C, C].

Recall that, given the coefficients, we recover the dictionary by making progress on the least squares objective (P1) (ignoring the term penalizing S(•)). Note that, our algorithm is based on finding an appropriate direction to ensure descent based on the geometry of the objective. To this end, we adopt a gradient descent-based strategy for dictionary update. However, since the coefficients are not exactly known, this results in an approximate gradient descent-based approach, where the empirical gradient estimate is formed as (5). In our analysis, we establish the conditions under which both the empirical gradient vector (corresponding to each dictionary element) and the gradient matrix concentrate around their means. To ensure progress at each iterate t, we show that the expected gradient vector is (Ω(k/m), Ω(m/k), 0)-correlated with the descent direction, defined as follows.

Definition 5. A vector g (t) is (ρ -, ρ + , ζ t )-correlated with a vector z * if g (t) , z (t) -z * ≥ ρ -z (t) -z * 2 + ρ + g (t) 2 -ζ t .

This can be viewed as a local descent condition which leads to the true dictionary columns; see also Candès et al. (2015) , Chen and Wainwright (2015) and Arora et al. (2015) . In convex optimization literature, this condition is implied by the 2ρ --strong convexity, and 1/2ρ + -smoothness of the objective. We show that for NOODL, ζ t = 0, which facilitates linear convergence to A * without incurring any bias. Overall our specific model assumptions for the analysis can be formalized as: 

A.1 A * is µ-incoherent (Def. 3), where µ = O(log(n)), A * = O( m/n) and m = O(n); A.2

(r) x < c 1 ( t , µ, n, k) = Ω(k/ √ n) < 1 and τ (r) = c 2 ( t , µ, k, n) = Ω(k 2 /n) for small constants c 1 and c 2 .

We are now ready to state our main result. A summary of the notation followed by a details of the analysis is provided in Appendix A and Appendix B, respectively.

Theorem 1 (Main Result). Suppose that assumptions A.1-A.6 hold, and Algorithm 1 is provided with p = Ω(mk 2 ) new samples generated according to model (1) at each iteration t. Then, with probability at least (1 -δ

(t)

alg ) for some small constant δ

(t)

alg , given R = Ω(log(n)), the coefficient estimate x (t) i at t-th iteration has the correct signed-support and satisfies ( x

(t) i -x * i ) 2 = O(k(1 -ω) t/2 A (0) i -A * i )

, for all i ∈ supp(x * ). Furthermore, for some 0 < ω < 1/2, the estimate A (t) at (t)-th iteration satisfies

A (t) i -A * i 2 ≤ (1 -ω) t A (0) i -A * i 2

, for all t = 1, 2, . . . ..

Our main result establishes that when the model satisfies A.1∼A.3, the errors corresponding to the dictionary and coefficients geometrically decrease to the true model parameters, given appropriate dictionary initialization and learning parameters (step sizes and threshold); see A.4∼A.6. In other words, to attain a target tolerance of T and δ T , where

A (T ) i -A * i ≤ T , | x (T ) i -x * i | ≤ δ T , we require T = max(Ω(log(1/ T )), Ω(log( √ k/δ T ))

) outer iterations and R = Ω(log(1/δ R )) IHT steps per outer iteration. Here, δ R ≥ (1 -η x ) R is the target decay tolerance for the IHT steps. An appropriate number of IHT steps, R, remove the dependence of final coefficient error (per outer iteration) on the initial x (0) . In Arora et al. (2015) , this dependence in fact results in an irreducible error, which is the source of bias in dictionary estimation. As a result, since (for NOODL) the error in the coefficients only depends on the error in the dictionary, it can be made arbitrarily small, at a geometric rate, by the choice of T , δ T , and δ R . Also, note that, NOODL can tolerate i.i.d. noise, as long as the noise variance is controlled to enable the concentration results to hold; we consider the noiseless case here for ease of discussion, which is already highly involved.

Intuitively, Theorem 1 highlights the symbiotic relationship between the two factors. It shows that, to make progress on one, it is imperative to make progress on the other. The primary condition that allows us to make progress on both factors is the signed-support recovery (Def. 2). However, the introduction of IHT step adds complexity in the analysis of both the dictionary and coefficients. To analyze the coefficients, in addition to deriving conditions on the parameters to preserve the correct signed-support, we analyze the recursive IHT update step, and decompose the noise term into a component that depends on the error in the dictionary, and the other that depends on the initial coefficient estimate. For the dictionary update, we analyze the interactions between elements of the coefficient vector (introduces by the IHT-based update step) and show that the gradient vector for the dictionary update is (Ω(k/m), Ω(m/k), 0)-correlated with the descent direction. In the end, this leads to exact recovery of the coefficients and removal of bias in the dictionary estimation. Note that our analysis pipeline is standard for the convergence analysis for iterative algorithms. However, the introduction of the IHT-based strategy for coefficient update makes the analysis highly involved as compared to existing results, e.g., the simple HT-based coefficient estimate in Arora et al. (2015) .

NOODL has an overall running time of O(mnp log(1/δ R ) max(log(1/ T ), log( √ k/δ T )) to achieve target tolerances T and δ T , with a total sample complexity of p•T = Ω(mk 2 ). Thus to remove bias, the IHT-based coefficient update introduces a factor of log(1/δ R ) in the computational complexity as compared to Arora et al. (2015) (has a total sample complexity of p • T = Ω(mk)), and also does not have the exponential running time and sample complexity as Barak et al. (2015) ; see Table 1 .

4. NEURAL IMPLEMENTATION OF NOODL:

The neural plausibility of our algorithm implies that it can be implemented as a neural network. This is because, NOODL employs simple linear and non-linear operations (such as inner-product and hard-thresholding) and the coefficient updates are separable across data samples, as shown in (4) of Algorithm 1. To this end, we present a neural implementation of our algorithm in Fig. 1 , which showcases the applicability of NOODL in large-scale distributed learning tasks, motivated from the implementations described in (Olshausen and Field, 1997) and (Arora et al., 2015) .

The neural architecture shown in Fig. 1 (a) has three layers -input layer, weighted residual evaluation layer, and the output layer. The input to the network is a data and step-size pair (y (j) , η x ) to each input node. Given an input, the second layer evaluates the weighted residuals as shown in Fig. 1 . Finally, the output layer neurons evaluate the IHT iterates x (r+1) (j) (4). We illustrate the operation of this architecture using the timing diagram in Fig. 1(b ). The main stages of operation are as follows. is given by (4). Output:

x ← 0 0 x (0) (j) = Tτ (A (t) y(j)) x (0) (j)

x

(1) (j)

x

(1) (j)

. . .

x (R) (j)

Residual: 0 y(j) y(j) ηx(y(j) -

A (t) x (0) (j) ) ηx(y(j) -A (t) x (0) (j) ) ηx(y(j) -A (t) x (1) (j) ) . . . ηx(y(j) -A (t) x (R-1) (j)

) Input: (y(j), 1)

.

(y(j), ηx) . . . . . . (y(j), 1)

(b) The timing sequence of the neural implementation.

Initial Hard Thresholding Phase: The coefficients initialized to zero, and an input (y (j) , 1) is provided to the input layer at a time instant = 0, which communicates these to the second layer. Therefore, the residual at the output of the weighted residual evaluation layer evaluates to y (j) at = 1. Next, at = 2, this residual is communicated to the output layer, which results in evaluation of the initialization x (0) (j) as per (3). This iterate is communicated to the second layer for the next residual evaluation. Also, at this time, the input layer is injected with (y (j) , η x ) to set the step size parameter η x for the IHT phase, as shown in Fig. 1(b) .

Iterative Hard Thresholding (IHT) Phase: Beginning = 3, the timing sequence enters the IHT phase. Here, the output layer neurons communicate the iterates x

(r+1) (j)

to the second layer for evaluation of subsequent iterates as shown in Fig. 1(b ). The process then continues till the time instance = 2R + 1, for R = Ω(log(1/δ R )) to generate the final coefficient estimate x

(t) (j) := x (R) (j)

for the current batch of data. At this time, the input layer is again injected with (y (j) , 1) to prepare the network for residual sharing and gradient evaluation for dictionary update.

Dictionary Update Phase: The procedure now enters the dictionary update phase, denoted as "Hebbian Learning" in the timing sequence. In this phase, each output layer neuron communicates the final coefficient estimate x (t) (j) = x (R) (j) to the second layer, which evaluates the residual for one last time (with η x = 1), and shares it across all second layer neurons ("Hebbian learning"). This allows each second layer neuron to evaluate the empirical gradient estimate (5), which is used to update the current dictionary estimate (stored as weights) via an approximate gradient descent step. This completes one outer iteration of Algorithm 1, and the process continues for T iterations to achieve target tolerances T and δ T , with each step receiving a new mini-batch of data.

5. EXPERIMENTS:

We now analyze the convergence properties and sample complexity of NOODL via experimental evaluations2 . The experimental data generation set-up, additional results, including analysis of computational time, are shown in Appendix E.

5.1. CONVERGENCE ANALYSIS:

We compare the performance of our algorithm NOODL with the current state-of-the-art alternating optimization-based online algorithms presented in Arora et al. (2015) , and the popular algorithm presented in Mairal et al. (2009) (denoted as Mairal '09). First of these, Arora15(''biased''), is a simple neurally plausible method which incurs a bias and has a sample complexity of Ω(mk). The other, referred to as Arora15(''unbiased''), incurs no bias as per Arora et al. (2015) , but the sample complexity results were not established.

Discussion: Fig. 2 panels (a-i), (b-i), (c-i), and (d-i) show the performance of the aforementioned methods for k = 10, 20, 50, and 100, respectively. Here, for all experiments we set η x = 0.2 and τ = 0.1. We terminate NOODL when the error in dictionary is less than 10 -10 . Also, for coefficient update, we terminate when change in the iterates is below 10 -12 . For k = 10, 20 and k = 50,  and (d-i ) show the convergence of NOODL, Arora15(''biased''), Arora15(''unbiased'') and Mairal '09, for different sparsity levels for n = 1000, m = 1500 and p = 5000. Since NOODL also recovers the coefficients, we show the corresponding recovery of the dictionary, coefficients, and overall fit in panels (a-ii), (b-ii), (c-ii), and (d-ii), respectively. Further, panels (e-i) and (e-ii) show the phase transition in samples p (per iteration) with the size of the dictionary m averaged across 10 Monte Carlo simulations for the two factors. Here, n = 100, k = 3, ηx = 0.2, τ = 0.1, 0 = 2/ log(n), ηA is chosen as per A.5. A trial is considered successful if the relative Frobenius error incurred by A and X is below 5 × 10 -7 after 50 iterations.

we note that Arora15(''biased'') and Arora15(''unbiased'') incur significant bias, while NOODL converges to A * linearly. NOODL also converges for significantly higher choices of sparsity k, i.e., for k = 100 as shown in panel (d), beyond k = O( √ n), indicating a potential for improving this bound. Further, we observe that Mairal '09 exhibits significantly slow convergence as compared to NOODL. Also, in panels (a-ii), (b-ii), (c-ii) and (d-ii) we show the corresponding performance of NOODL in terms of the error in the overall fit ( Y -AX F / Y F ), and the error in the coefficients and the dictionary, in terms of relative Frobenius error metric discussed above. We observe that the error in dictionary and coefficients drops linearly as indicated by our main result.

5.2. PHASE TRANSITIONS:

Fig. 2 panels (e-i ) and (e-ii), shows the phase transition in number of samples with respect to the size of the dictionary m. We observe a sharp phase transition at p m = 1 for the dictionary, and at p m = 0.75 for the coefficients. This phenomenon is similar to that observed by Agarwal et al. (2014) (however, theoretically they required p = O(m 2 )). Here, we confirm number of samples required by NOODL are linearly dependent on the dictionary elements m.

6. FUTURE WORK:

We consider the online DL setting in this work. We note that, empirically NOODL works for the batch setting also. However, analysis for this case will require more sophisticated concentration results, which can address the resulting dependence between iterations of the algorithm. In addition, our experiments indicate that NOODL works beyond the sparsity ranges prescribed by our theoretical results. Arguably, the bounds on sparsity can potentially be improved by moving away from the incoherence-based analysis. We also note that in our experiments, NOODL converges even when initialized outside the prescribed initialization region, albeit it achieves the linear rate once it satisfies the closeness condition A.4. These potential directions may significantly impact the analysis and development of provable algorithms for other factorization problems as well. We leave these research directions, and a precise analysis under the noisy setting, for future explorations.

7. CONCLUSIONS:

We present NOODL, to the best of our knowledge, the first neurally plausible provable online algorithm for exact recovery of both factors of the dictionary learning (DL) model. NOODL alternates between: (a) an iterative hard thresholding (IHT)-based step for coefficient recovery, and (b) a gradient descent-based update for the dictionary, resulting in a simple and scalable algorithm, suitable for large-scale distributed implementations. We show that once initialized appropriately, the sequence of estimates produced by NOODL converge linearly to the true dictionary and coefficients without incurring any bias in the estimation. Complementary to our theoretical and numerical results, we also design an implementation of NOODL in a neural architecture for use in practical applications.

In essence, the analysis of this inherently non-convex problem impacts other matrix and tensor factorization tasks arising in signal processing, collaborative filtering, and machine learning.

A SUMMARY OF NOTATION:

We summarizes the definitions of some frequently used symbols in our analysis in Table 2 . In addition, we use D (v) as a diagonal matrix with elements of a vector v on the diagonal. Given a matrix M, we use M -i to denote a resulting matrix without i-th column. Also note that, since we show that A (t)

i -A * i ≤ t contracts in every step, therefore we fix t , 0 = O * (1/ log(n)) in our analysis. i-th column of the dictionary estimate at the t-th iterate.

t A (t) i -A * i ≤ t = O * ( 1 log(n) )

Upper-bound on column-wise error at the t-th iterate.

µ t µt √ n = µ √ n + 2 t

Incoherence between the columns of A (t) ; See Claim 1.

λ (t) j λ (t) j := | A (t) j -A * j , A * j | ≤ 2 t 2

Inner-product between the error and the dictionary element.

Λ (t) S (i, j) Λ (t) S (i, j) = λ (t) j , for j = i, i ∈ S 0, otherwise. A diagonal matrix of size |S| × |S| with λ (t)

j on the diagonal for j ∈ S.

Coefficient Related Symbol Definition x (r) i i-th element the coefficient estimate at the r-th IHT iterate. C |x * i | ≥ C for i ∈ supp(x * ) and C ≤ 1 Lower-bound on x * i s. S S := supp(x * ) where |S| ≤ k Support of x * δ R δ R := (1 -η x + η x µt √ n ) R ≥ (1 -η x ) R Decay parameter for coefficients. δ T | x (T ) i -x * i | ≤ δ T ∀i ∈ supp(x * ) Target coefficient element error tol- erance. C ( ) i C ( ) i := |x * i -x ( ) i | for i ∈ supp(x * )

Error in non-zero elements of the coefficient vector.

Probabilities Symbol Definition Symbol Definition:

q i q i = Pr[i ∈ S] = Θ( k m ) q i,j q i,j = Pr[i, j ∈ S] = Θ( k 2 m 2 ) p i p i = E[x * i sign(x * i )|x * i = 0] δ (t) T δ (t) T = 2m exp(-C 2 /O * ( 2 t )) δ (t) β δ (t) β = 2k exp(-1/O( t )) δ (t) HW δ (t) HW = exp(-1/O( t )) δ (t) gi δ (t) gi = exp(-Ω(k)) δ (t) g δ (t) g = (n + m) exp(-Ω(m log(n))

Other terms Symbol Definition ξ

(r+1) j ξ (r+1) j := i =j ( A (t) j -A * j , A * i + A * j , A * i )x * i - i =j A (t) j , A (t) i x (r) i β (t) j β (t) j := i =j ( A * j , A * i -A (t) i + A * j -A (t) j , A (t) i + A (t) j -A * j , A * i )x * i t β t β = O( √ k t ) is an upper-bound on β (t) j with probability at least (1 -δ (t) β ) ξ (r+1) j ξ (r+1) j := β (t) j + i =j | A (t) j , A (t) i | |x * i -x (r) i | ∆ (t) j ∆ (t) j := E[A (t) S ϑ (R) S sign(x * j )] ϑ (R) i ϑ (R) i := R r=1 η x ξ (r) i (1 -η x ) R-r + γ (R) i γ (R) i γ (R) i := (1 -η x ) R (x (0) i -x * i (1 -λ (t) i )) γ γ := E[(A (t) x -y)sign(x * j )1 F x * ]; See † below. x i x i := x (R) i = x * i (1 -λ (t) i ) + ϑ (R) i † 1 F x * is the indicator function corresponding to the event that sign(x * ) = sign( x),

denoted by F x * , and similarly for the complement F x *

B PROOF OF THEOREM 1:

We now prove our main result. The detailed proofs of intermediate lemmas and claims are organized in Appendix C and Appendix D, respectively. Furthermore, the standard concentration results are stated in Appendix F for completeness. Also, see Table 3 for a map of dependence between the results.

OVERVIEW:

Given an ( 0 , 2)-close estimate of the dictionary, the main property that allows us to make progress on the dictionary is the recovery of the correct sign and support of the coefficients. Therefore, we first show that the initial coefficient estimate (3) recovers the correct signed-support in Step I.A. Now, the IHT-based coefficient update step also needs to preserve the correct signed-support. This is to ensure that the approximate gradient descent-based update for the dictionary makes progress. Therefore, in Step I.B, we derive the conditions under which the signed-support recovery condition is preserved by the IHT update.

To get a handle on the coefficients, in Step II.A, we derive an upper-bound on the error incurred by each non-zero element of the estimated coefficient vector, i.e., | x i -x * i | for i ∈ S for a general coefficient vector x * , and show that this error only depends on t (the column-wise error in the dictionary) given enough IHT iterations R as per the chosen decay parameter δ R . In addition, for analysis of the dictionary update, we develop an expression for the estimated coefficient vector in

Step II.B.

We then use the coefficient estimate to show that the gradient vector satisfies the local descent condition (Def. 5). This ensures that the gradient makes progress after taking the gradient descent-based step (6). To begin, we first develop an expression for the expected gradient vector (corresponding to each dictionary element) in Step III.A. Here, we use the closeness property Def 1 of the dictionary estimate. Further, since we use an empirical estimate, we show that the empirical gradient vector concentrates around its mean in Step III.B. Now using Lemma 15, we have that descent along this direction makes progress.

Next in Step IV.A and Step IV.B, we show that the updated dictionary estimate maintains the closeness property Def 1. This sets the stage for the next dictionary update iteration. As a result, our main result establishes the conditions under which any t-th iteration succeeds.

Our main result is as follows.

Theorem 1 (Main Result) Suppose that assumptions A.1-A.6 hold, and Algorithm 1 is provided with p = Ω(mk 2 ) new samples generated according to model (1) at each iteration t. Then, with probability at least (1 -δ

(t) alg ), given R = Ω(log(n)), the coefficient estimate x (t)

i at t-th iteration has the correct signed-support and satisfies

( x (t) i -x * i ) 2 = O(k(1 -ω) t/2 A (0) i -A * i ), for all i ∈ supp(x * ). Furthermore, for some 0 < ω < 1/2, the estimate A (t) at (t)-th iteration satisfies A (t) i -A * i 2 ≤ (1 -ω) t A (0) i -A * i 2 , for all t = 1, 2, . . . ..

Here, δ (t):

alg is some small constant, where δ

(t) alg = δ (t) T + δ (t) β + δ HW + δ (t) gi + δ (t) g , δ (t) T = 2m exp(-C 2 /O * ( 2 t )), δ (t) β = 2k exp(-1/O( t )), δ (t) HW = exp(-1/O( t )), δ (t) gi = exp(-Ω(k)), δ (t) g = (n + m) exp(-Ω(m log(n)), and A (t) i -A * i ≤ t .

STEP I: COEFFICIENTS HAVE THE CORRECT SIGNED-SUPPORT:

As a first step, we ensure that our coefficient estimate has the correct signed-support (Def. 2). To this end, we first show that the initialization has the correct signed-support, and then show that the iterative hard-thresholding (IHT)-based update step preserves the correct signed-support for a suitable choice of parameters.

• Step I.A: Showing that the initial coefficient estimate has the correct signed-support-Given an ( 0 , 2)-close estimate A (0) of A * , we first show that for a general sample y the initialization step (3) recovers the correct signed-support with probability at least (1-δ

(t) T ), where δ (t) T = 2m exp(-C 2 O * ( 2 t )

). This is encapsulated by the following lemma. Lemma 1 (Signed-support recovery by coefficient initialization step). Suppose

A (t) is t -close to A * . Then, if µ = O(log(n)), k = O * ( √ n/µ log(n)), and t = O * (1/ log(m)), with probability at least (1 -δ (t) T ) for each random sample y = A * x * : sign(T C/2 ((A (t) ) y) = sign(x * ),

where δ

(t) T = 2m exp(-C 2 O * ( 2 t ) ).

Note that this result only requires the dictionary to be column-wise close to the true dictionary, and works for less stringent conditions on the initial dictionary estimate, i.e., requires

t = O * (1/ log(m)) instead of t = O * (1/ log(m))

; see also (Arora et al., 2015) .

• Step I.B: The iterative IHT-type updates preserve the correct signed support-Next, we show that the IHT-type coefficient update step (4) preserves the correct signed-support for an appropriate choice of step-size parameter η (r)

x and threshold τ (r) . The choice of these parameters arises from the analysis of the IHT-based update step. Specifically, we show that at each iterate r, the step-size η (r)

x should be chosen to ensure that the component corresponding to the true coefficient value is greater than the "interference" introduced by other non-zero coefficient elements. Then, if the threshold is chosen to reject this "noise", each iteration of the IHT-based update step preserves the correct signed-support.

Lemma 2 (IHT update step preserves the correct signed-support). Suppose

A (t) is t - close to A * , µ = O(log(n)), k = O * ( √ n/µ log(n)), and t = O * (1/ log(m)) Then, with probability at least (1 -δ (t) β -δ (t)

T ), each iterate of the IHT-based coefficient update step shown in (4) has the correct signed-support, if for a constant c (r)

1 ( t , µ, k, n) = Ω(k 2 /n), the step size is chosen as η (r) x ≤ c (r)

1 , and the threshold τ (r) is chosen as

τ (r) = η (r) x (t β + µt √ n x (r-1) -x * 1 ) := c (r) 2 ( t , µ, k, n) = Ω(k 2 /n),

for some constants c 1 and c 2 . Here,

EQUATION

T = 2m exp(-C 2 O * ( 2 t )

) ,and δ

(t) β = 2k exp(-1 O( t) ).

Note that, although we have a dependence on the iterate r in choice of η (r)

x and τ (r) , these can be set to some constants independent of r. In practice, this dependence allows for greater flexibility in the choice of these parameters.

STEP II: ANALYZING THE COEFFICIENT ESTIMATE:

We now derive an upper-bound on the error incurred by each non-zero coefficient element. Further, we derive an expression for the coefficient estimate at the t-th round of the online algorithm x (t) := x (R) ; we use x instead of x (t) for simplicity.

•:

Step II.A: Derive a bound on the error incurred by the coefficient estimate-Since Lemma 2 ensures that x has the correct signed-support, we now focus on the error incurred by each coefficient element on the support by analyzing x. To this end, we carefully analyze the effect of the recursive update (4), to decompose the error incurred by each element on the support into two components -one that depends on the initial coefficient estimate x (0) and other that depends on the error in the dictionary. We show that the effect of the component that depends on the initial coefficient estimate diminishes by a factor of (1 -η x + η x Lemma 3 (Upper-bound on the error in coefficient estimation). With probability at least (1 -δ

(t) β -δ (t)

T ) the error incurred by each element i 1 ∈ supp(x * ) of the coefficient estimate is upper-bounded as

| x i1 -x * i1 | ≤ O(t β ) + (R + 1)kη x µt √ n max i |x (0) i -x * i | + |x (0) i1 -x * i1 | δ R , = O(t β )

where

t β = O( √ k t ), δ R := (1 -η x + η x µt √ n ) R , δ (t) T = 2m exp(-C 2 O * ( 2 t ) ), δ (t) β = 2k exp(-1 O( t) )

, and µ t is the incoherence between the columns of A (t) ; see Claim 1. This result allows us to show that if the column-wise error in the dictionary decreases at each iteration t, then the corresponding estimates of the coefficients also improve.

• Step II.B: Developing an expression for the coefficient estimate-Next, we derive the expression for the coefficient estimate in the following lemma. This expression is used to analyze the dictionary update.

Lemma 4 (Expression for the coefficient estimate at the end of R-th IHT iteration).

With probability at least (1 -δ

(t) T -δ (t)

β ) the i 1 -th element of the coefficient estimate, for each i 1 ∈ supp(x * ), is given by

x i1 := x (R) i1 = x * i1 (1 -λ (t) i1 ) + ϑ (R) i1 .

Here,

EQUATION

T = 2m exp(-C 2 O * ( 2 t ) ) and δ (t) β = 2k exp(-1 O( t)

). We again observe that the error in the coefficient estimate depends on the error in the dictionary via λ 

STEP III: ANALYZING THE GRADIENT FOR DICTIONARY UPDATE:

Given the coefficient estimate we now show that the choice of the gradient as shown in ( 5) makes progress at each step. To this end, we analyze the gradient vector corresponding to each dictionary element to see if it satisfies the local descent condition of Def. 5. Our analysis of the gradient is motivated from Arora et al. (2015) . However, as opposed to the simple HT-based coefficient update step used by Arora et al. (2015) , our IHT-based coefficient estimate adds to significant overhead in terms of analysis. Notwithstanding the complexity of the analysis, we show that this allows us to remove the bias in the gradient estimate.

To this end, we first develop an expression for each expected gradient vector, show that the empirical gradient estimate concentrates around its mean, and finally show that the empirical gradient vector is (Ω(k/m), Ω(m/k), 0)-correlated with the descent direction, i.e. has no bias.

• Step III.A: Develop an expression for the expected gradient vector corresponding to each dictionary element-The expression for the expected gradient vector g

(t)

j corresponding to j-th dictionary element is given by the following lemma.

Lemma 5 (Expression for the expected gradient vector). Suppose that A (t) is ( t , 2)near to A * . Then, the dictionary update step in Algorithm 1 amounts to the following for the j-th dictionary element

E[A (t+1) j ] = A (t) j + η A g (t) j ,

where g

(t)

j is given by g

EQUATION

∆ (t) j := E[A (t) S ϑ (R) S sign(x * j )],

where

∆ (t) j = O( √ mq i,j p j t A (t)

).

•

Step III.B: Show that the empirical gradient vector concentrates around its expectation-Since we only have access to the empirical gradient vectors, we show that these concentrate around their expected value via the following lemma.

Lemma 6 (Concentration of the empirical gradient vector). Given p = Ω(mk 2 ) samples, the empirical gradient vector estimate corresponding to the i-th dictionary element, g

i concentrates around its expectation, i.e., g

(t) i -g (t) i ≤ o( k m t ). with probability at least (1 -δ (t) gi -δ (t) β -δ (t) T -δ (t)

HW ), where δ (t) gi = exp(-Ω(k)).

• Step III.C: Show that the empirical gradient vector is correlated with the descent direction-Next, in the following lemma we show that the empirical gradient vector g

(t)

j is correlated with the descent direction. This is the main result which enables the progress in the dictionary (and coefficients) at each iteration t. Lemma 7 (Empirical gradient vector is correlated with the descent direction). Suppose

A (t) is ( t , 2)-near to A * , k = O( √ n) and η A = O(m/k). Then, with prob- ability at least (1 -δ (t) T -δ (t) β -δ (t) HW -δ (t) gi ) the empirical gradient vector g (t) j is (Ω(k/m), Ω(m/k), 0)-correlated with (A (t) j -A * j )

, and for any t ∈ [T ],

A (t+1) j -A * j 2 ≤ (1 -ρ η A ) A (t) j -A * j 2 .

This result ensures for at any t ∈ [T ], the gradient descent-based updates made via (5) gets the columns of the dictionary estimate closer to the true dictionary, i.e., t+1 ≤ t . Moreover, this step requires closeness between the dictionary estimate A (t) and A * , in the spectral norm-sense, as per Def 1.

STEP IV: SHOW THAT THE DICTIONARY MAINTAINS THE CLOSENESS PROPERTY:

As discussed above, the closeness property (Def 1) is crucial to show that the gradient vector is correlated with the descent direction. Therefore, we now ensure that the updated dictionary A (t+1) maintains this closeness property. Lemma 7 already ensures that t+1 ≤ t . As a result, we show that A (t+1) maintains closeness in the spectral norm-sense as required by our algorithm, i.e., that it is still ( t+1 , 2)-close to the true dictionary. Also, since we use the gradient matrix in this analysis, we show that the empirical gradient matrix concentrates around its mean.

Step IV.A: The empirical gradient matrix concentrates around its expectation: We first show that the empirical gradient matrix concentrates as formalized by the following lemma.

Lemma 8 (Concentration of the empirical gradient matrix). With probability at least

(1 -δ (t) β -δ (t) T -δ (t) HW -δ (t) g ), g (t) -g (t) is upper-bounded by O * ( k m A * ), where δ (t) g = (n + m) exp(-Ω(m log(n)). •

Step IV.B: The "closeness" property is maintained after the updates made using the empirical gradient estimate: Next, the following lemma shows that the updated dictionary A (t+1) maintains the closeness property. Lemma 9 (A (t+1) maintains closeness). Suppose

A (t) is ( t , 2) near to A * with t = O * (1/ log(n))

, and number of samples used in step t is p = Ω(mk 2 ), then with probability at least (1 -δ

(t) T -δ (t) β -δ (t) HW -δ (t) g ), A (t+1) satisfies A (t+1) -A * ≤ 2 A * .

STEP V: COMBINE RESULTS TO SHOW THE MAIN RESULT:

Proof of Theorem 1. From Lemma 7 we have that with probability at least (1 -δ

(t) T -δ (t) β -δ (t) HW - δ (t) gi ), g (t)

j is (Ω(k/m), Ω(m/k), 0)-correlated with A * j . Further, Lemma 9 ensures that each iterate maintains the closeness property. Now, applying Lemma 15 we have that, for η A ≤ Θ(m/k), with probability at least (1 -δ

(t) alg ) any t ∈ [T ] satisfies A (t) j -A * j 2 ≤ (1 -ω) t A (0) j -A * j 2 ≤ (1 -ω) t 2

0 . where for 0 < ω < 1/2 with ω = Ω(k/m)η A . That is, the updates converge geometrically to A * . Further, from Lemma 3, we have that the result on the error incurred by the coefficients. Here,

δ (t) alg = δ (t) T + δ (t) β + δ (t) HW + δ (t) gi + δ (t)

g ). That is, the updates converge geometrically to A * . Further, from Lemma 3, we have that the error in the coefficients only depends on the error in the dictionary, which leads us to our result on the error incurred by the coefficients. This completes the proof of our main result.

C APPENDIX: PROOF OF LEMMAS:

We present the proofs of the Lemmas used to establish our main result. Also, see Table 3 for a map of dependence between the results, and Appendix D for proofs of intermediate results. Proof of Lemma 1. Let y ∈ R n be general sample generated as y = A * x * , where x * ∈ R m is a sparse random vector with support S = supp(x * ) distributed according to D.4.

The initial decoding step at the t-th iteration (shown in Algorithm 1) involves evaluating the innerproduct between the estimate of the dictionary A (t) , and y. The i-th element of the resulting vector can be written as

A (t) i , y = A (t) i , A * i x * i + w i ,

where

w i = A (t) i , A * -i x * -i . Now, since A * i -A (t) i 2 ≤ t and

A * i -A (t) i 2 2 = A * i 2 + A (t) i 2 -2 A (t) i , A * i = 2 -2 A (t) i , A * i , we have | A (t) i , A * i | ≥ 1 -2 t /2. Therefore, the term | A (t) i , A * i x * i | ≥ (1 - 2 t 2 )C , if i ∈ S, = 0

, otherwise. Now, we focus on the w i and show that it is small. By the definition of w i we have

w i = A (t) i , A * -i x * -i = =i A (t) i , A * x * = ∈S\{i} A (t) i , A * x * .

Here, since var(x * ) = 1, w i is a zero-mean random variable with variance

var(w i ) = ∈S\{i} A (t) i , A * 2 .

Now, each term in this sum can be bounded as,

A (t) i , A * 2 = ( A (t) i -A * i , A * + A * i , A * ) 2 ≤ 2( A (t) i -A * i , A * 2 + A * i , A * 2 ) ≤ 2( A (t) i -A * i , A * 2 + µ 2 n ). Next, =i A (t) i -A * i , A * 2 can be upper-bounded as ∈S\{i} A (t) i -A * i , A * 2 ≤ A * S\{i} 2 2 t .

Therefore, we have the following as per our assumptions on µ and k, (Gershgorin, 1931) . Therefore, we have ∈S\{i}

A * S\{i} 2 ≤ (1 + k µ √ n ) ≤ 2, using Gershgorin Circle Theorem

A (t) i -A * i , A * 2 ≤ 2 2 t .

Finally, we have that ∈S\{i}

A (t) i , A * 2 ≤ 2(2 2 t + k µ 2 n ) = O * ( 2 t ).

Now, we apply the Chernoff bound for sub-Gaussian random variables w i (shown in Lemma 12) to conclude that

Pr[|w i | ≥ C/4] ≤ 2 exp(-C 2 O * ( 2 t )

). Further, w i corresponding to each m should follow this bound, applying union bound we conclude that

Pr[max i |w i | ≥ C/4] ≤ 2m exp(-C 2 O * ( 2 t ) ) := δ (t) T .

Proof of Lemma 2. Consider the (r + 1)-th iterate x (r+1) for the t-th dictionary iterate, where

A (t) i -A * i ≤ t for all i ∈ [1, m]

evaluated as the following by the update step described in Algorithm 1,

EQUATION

where η

(1)

x < 1 is the learning rate or the step-size parameter. Now, using Lemma 1 we know that x (0) (3) has the correct signed-support with probability at least (1 -δ (t) T ). Further, since A (t) A * can be written as

A (t) A * = (A (t) -A * ) A * + A * A * ,

we can write the (r + 1)-th iterate of the coefficient update step using (7) as

x (r+1) = (I -η (r+1) x A (t) A (t) )x (r) -η (r+1) x (A (t) -A * ) A * x * + η (r+1) x A * A * x * .

Further, the j-th entry of this vector is given by

EQUATION

)

We now develop an expression for the j-th element of each of the term in (8) as follows. First, we can write the first term as

EQUATION

i .

Next, the second term in (8) can be expressed as

η (r+1) x ((A (t) -A * ) A * ) (j,:) x * = η (r+1) x i A (t) j -A * j , A * i x * i = η (r+1) x A (t) j -A * j , A * j x * j + η (r+1) x i =j A (t) j -A * j , A * i x * i .

Finally, we have the following expression for the third term,

η (r+1) x (A * A * ) (j,:) x * = η (r+1) x x * j + η (r+1) x i =j A * j , A * i x * i .

Now using our definition of λ

(t) j = | A (t) j -A * j , A * j | ≤ 2 t

2 , combining all the results for (8), and using the fact that since A (t) is close to A * , vectors A (t) j -A * j and A * j enclose an obtuse angle, we have the following for the j-th entry of the (r + 1)-th iterate, x (r+1) is given by

EQUATION

)

Here ξ (r+1) j is defined as

EQUATION

, we can write ξ

EQUATION

where β (t) j is defined as

EQUATION

Note that β (t) j does not change for each iteration r of the coefficient update step. Further, by Claim 2 we show that |β

EQUATION

)

where ξ

(r+1) j ≤ ξ (r+1) j

. Further, using Claim 1,

EQUATION

since x (r-1) -x * 1 = O(k). Therefore, for the (r + 1)-th iteration, we choose the threshold to be

EQUATION

and the step-size by setting the "noise" component of ( 9) to be smaller than the "signal" part, specifically, half the signal component, i.e., η

x

ξ (r+1) max ≤ (1-η (r+1) x ) 2 x (r) min + η (r+1) x 2 (1 - 2 t

2 )C, Also, since we choose the threshold as τ (r) := η

(r) x ξ (r) max , x (r) min = η (r) x ξ (r)

max , where x (0) min = C/2, we have the following for the (r + 1)-th iteration,

η (r+1) x ξ (r+1) max ≤ (1-η (r+1) x ) 2 η (r) x ξ (r) max + η (r+1) x 2 (1 - 2 t 2 )C.

Therefore, for this step we choose η

EQUATION

Therefore, η

x can be chosen as

η (r+1) x ≤ c (r+1) ( t , µ, k, n),

for a small constant c (r+1) ( t , µ, k, n), η

. In addition, if we set all η (r)

x = η x , we have that η x = Ω( k √ n ) and therefore τ (r) = τ = Ω( k 2 n ). Further, since we initialize with the hard-thresholding step, the entries in |x (0) | ≥ C/2. Here, we define ξ 

i1 = β (t) i1 + i2 =i1 | A (t) i1 , A (t) i2 ||x * i2 -x ( -1) i2 |.

From Claim 2 we have that |β

(t) i1 | ≤ t β with probability at least (1 -δ (t)

β ). Further, using Claim 1 , and letting C

( ) i := |x * i -x ( ) i | = |x ( ) i -x * i |, ξ ( )

i1 can be upper-bounded as

EQUATION

)

Rearranging the expression for (r + 1)-th update (9), and using ( 16) we have the following upperbound

C (r+1) i1 ≤ (1 -η (r+1) x )C (r) i1 + η (r+1) x λ (t) i1 |x * i1 | + η (r+1) x ξ (r+1) i1

.

Next, recursively substituting in for C (r)

i1 , where we define q= (1 -η

(q+1) x ) = 1, C (r+1) i1 ≤ C (0) i1 r q=0 (1 -η (q+1) x ) + λ (t) i1 |x * i1 | r+1 =1 η ( ) x r+1 q= (1 -η (q+1) x ) + r+1 =1 η ( ) x ξ ( ) i1 r+1 q= (1 -η (q+1) x

).

Substituting for the upper-bound of ξ ( ) i1 from ( 16),

EQUATION

Here, α (r+1) i1

is defined as

α (r+1) i1 = C (0) i1 r q=0 (1 -η (q+1) x ) + (λ (t) i1 |x * i1 | + β (t) i1 ) r+1 =1 η ( ) x r+1 q= (1 -η (q+1) x

).

(18)

Our aim now will be to express C ( )

i1 for > 0 in terms of C (0) i2 . Let each α ( ) j ≤ α ( ) i

where

j = i 1 , i 2 , . . . , i k . Similarly, let C (0) j ≤ C (0) i for j = i 1 , i 2 , .

. . , i k , and all η ( )

x = η x . Then, using Claim 3 we have the following expression for

C (R+1) i1 , C (R+1) i1 ≤ α (R+1) i1 + (k -1)η x µt √ n R =1 α ( ) max 1 -η x +η x µt √ n R- + (k -1)η x µt √ n C (0) max 1 -η x + η x µt √ n R .

Here, (1

-η x ) R ≤ (1 -η x + η x µt √ n ) R ≤ δ R .

Next from Claim 4 we have that with probability at least (1 -δ

(t) β ), R =1 α ( ) max 1 -η x + η x µt √ n R-≤ C (0) max Rδ R + 1 ηx(1- µt √ n ) ( 2 t 2 |x * max | + t β ). Therefore, for c x = µt √ n /(1 -µt √ n ) C (R+1) i1 ≤ α (R+1) i1 + (k -1)c x ( 2 t 2 |x * max | + t β ) + (R + 1)(k -1)η x µt √ n C (0) max δ R .

Now, using the definition of α (R+1) i1

, and using the result on sum of geometric series, we have

α (R+1) i1 = C (0) i1 (1 -η x ) R+1 + (λ (t) i1 |x * i1 | + β (t) i1 ) R+1 s=1 η x (1 -η x ) R-s+1 , = C (0) i1 δ R + λ (t) i1 |x * i1 | + β (t) i1 ≤ C (0) i1 δ R+1 + 2 t 2 |x * max | + t β . Therefore, C (R)

i1 is upper-bounded as

C (R) i1 ≤ (c x k + 1)( 2 t 2 |x * max | + t β ) + (R + 1)kη x µt √ n C (0) max δ R + C (0) i1 δ R . Further, since k = O( √ n/µ log(n)), kc x < 1, therefore, we have C (R) i1 ≤ O(t β ) + (R + 1)kη x µt √ n C (0) max δ R + C (0) i1 δ R , with probability at least (1 -δ (t) β ). Here, (R + 1)kη x µt √ n C (0) max δ R + C (0)

i1 δ R 0 for an appropriately large R. Therefore, the error in each non-zero coefficient is

C (R) i1 = O(t β ). with probability at least (1 -δ (t) β ).

Proof of Lemma 4. Using the expression for x (R) i1 as defined in (9), and recursively substituting for x (r) i1 we have

x (R) i1 = (1 -η x ) R x (0) j + x * i1 R r=1 η x (1 -λ (t) i1 )(1 -η x ) R-r + R r=1 η x ξ (r) i1 (1 -η x ) R-r ,

where we set all η r

x to be η x . Further, on defining

EQUATION

where γ (R)

EQUATION

)), we have

EQUATION

Note that γ (R)

i1 can be made appropriately small by choice of R. Further, by Claim 5 we have

|ϑ (R) i1 | ≤ O(t β ).

with probability at least (1 -δ

(t) β ), where t β = O( √ k t ).

Proof of Lemma 5. From Lemma 4 we have that for each j ∈ S,

x S := x (R) S = (I -Λ (t) S )x * S + ϑ (R) S ,

with probability at least (1 -δ

(t) T -δ (t)

β ). Further, let F x * be the event that sign(x * ) = sign( x),

and let 1 F x * denote the indicator function corresponding to this event. As we show in Lemma 2, this event occurs with probability at least (1 -δ

EQUATION

T ). Using this, we can write the expected gradient vector corresponding to the j-th sample as 1

F x * g (t) j = E[(A (t) x -y)sign(x * j )1 F x * ] + E[(A (t) x -y)sign(x * j )1 F x * ], = E[(A (t) x -y)sign(x * j )1 F x * ] ± γ.

Here, γ := E[(A (t) xy)sign(x * j )1 F x * ] is small and depends on δ

(t)

T and δ

(t)

β , which in turn drops with t . Therefore, γ diminishes with t . Further, since 1

F x * + 1 F x * = 1, and Pr[F x * ] = (1 -δ (t) β -δ (t)

T ), is very large, t) xy)sign(x * j )] ± γ.

g (t) j = E[(A (t) x -y)sign(x * j )(1 -1 F x * )] ± γ, = E[(A (

Therefore, we can write g

EQUATION

S ] can be made very small by choice of R, we absorb this term in γ. Therefore,

g (t) j = E[A (t) S (I -Λ (t) S )x * S + A (t) S ϑ (R) S -A * S x * S )sign(x * j )] ± γ.

Writing the expectation by sub-conditioning on the support,

g (t) j = E S [E x * S [A (t) S (I -Λ (t) S )x * S sign(x * j ) -A * S x * S sign(x * j ) + A (t) S ϑ (R) S sign(x * j )|S]] ± γ, = E S [A (t) S (I -Λ (t) S )E x * S [x * S sign(x * j )|S] -A * S E x * S [x * S sign(x * j )|S]] + E[A (t) S ϑ (R) S sign(x * j )] ± γ, = E S [p j (1 -λ (t) j )A (t) j -p j A * j ] + ∆ (t)

j ± γ, where we have used the fact that E x * S [sign(x * j )] = 0 and introduced

∆ (t) j = E[A (t) S ϑ (R) S sign(x * j )]. Next, since p j = E x * S [x * j sign(x * j )|j ∈ S], therefore, g (t) j = E S [p j (1 -λ (t) j )A (t) j -p j A * j ] + ∆ (t) j ± γ.

Further, since q j = Pr[j ∈ S] = O(k/m),

g (t) j = q j p j (1 -λ (t) j )A (t) j -A * j + 1 qj pj ∆ (t) j ± γ .

Further, by Claim 7 we have that

∆ (t) j = O( √ mq i,j p j t A (t) )].

This completes the proof.

Proof of Lemma 6. Let W = {j : i ∈ supp(x * (j) )} and then we have that

g (t) i = |W | p 1 |W | j (y (j) -A (t) x (j) )sign( x (j) (i)),

where x (j) (i) denotes the i-th element of the coefficient estimate corresponding to the (j)-th sample.

Here, for = |W | the summation

j 1 (y (j) -A (t) x (j) )sign( x (j) (i)),

has the same distribution as Σ j=1 z j , where each z j belongs to a distribution as

z := 1 (y -A (t) x)sign( x i )|i ∈ S. Also, E[(y -A (t) x)sign( x i )] = q i E[z],

where

EQUATION

)

Let w j = z j -E[z], we will now apply the vector Bernstein result shown in Lemma 11. For this, we require bounds on two parameters for these -L := w j and σ 2 := Σ j E[ w j 2 ] . Note that, since the quantity of interest is a function of x * i , which are sub-Gaussian, they are only bounded almost surely. To this end, we will employ Lemma 14 (Lemma 45 in (Arora et al., 2015) ) to get a handle on the concentration.

Bound on the norm w : This bound is evaluated in Claim 8, which states that with probability at least (1 -δ

(t) β -δ (t) T -δ (t) HW ), L := w = z -E[z] = 2 (y -A (t) x)sign( x i )|i ∈ S ≤ 2 (y -A (t) x) = O( kt β ).

Bound on variance parameter E[ w 2 ]: Using Claim 9, we have

E[ z 2 ] = O(k 2 t ) + O(kt 2 β )

. Therefore, the bound on the variance parameter σ 2 is given by

σ 2 := Σ j E[ w j 2 ] ≤ Σ j E[ z j 2 ] ≤ O( k 2 t ) + O( kt 2 β ).

From Claim 2 we have that with probability at least (1 -δ

(t) β ), t β = O( √ k t ).

Applying vector Bernstein inequality shown in Lemma 11 and using Lemma 14 (Lemma 45 in (Arora et al., 2015) ), choosing = Ω(k 3 ), we conclude

j=1 z j -E[z] = O(L) + O(σ) = o( t ), with probability at least (1 -δ (t)

gi ), where δ (t) gi = exp(-Ω(k)). Finally, substituting in (21) we have

g (t) i -g (t) i = O( k m )o( t ).

with probability at least (1 -δ

(t) gi -δ (t) β -δ (t) T -δ (t) HW ).

Proof of Lemma 7. Since we only have access to the empirical estimate of the gradient g (t)

i , we will show that this estimate is correlated with (A (t) j -A * j ). To this end, first from Lemma 6 we have that the empirical gradient vector concentrates around its mean, specifically,

g (t) i -g (t) i ≤ o( k m t ),

with probability at least (1 -δ

(t) gi -δ (t) β -δ (t) T -δ (t)

HW ). From Lemma 5, we have the following expression for the expected gradient vector

g (t) j = p j q j (A (t) j -A * j ) + p j q j (-λ (t) j A (t) j + 1 pj qj ∆ (t) j ± γ).

Let g

(t) j = 4ρ (A (t) j -A * j ) + v

, where 4ρ = p j q j and v is defined as

EQUATION

)

Then, g

i can be written as

EQUATION

)

where v = v + ( g (t) i -g (t) i ). Let v ≤ ρ A (t) i -A * i .

Using the definition of v as shown in ( 22) we have

v ≤ q j p j λ (t) j A (t) j + ∆ (t) j + o( k m t ) ± γ.

Now for the first term, since A (t) j

= 1, we have λ

(t) j = | A (t) j -A * j , A * j | = 1 2 A (t) j -A * j 2 , therefore q j p j λ (t) j A (t) j = q j p j 1 2 A (t) j -A * j 2 ,

Further, using Claim 7

∆ (t) j = O( √ mq i,j p i1 t A (t) ).

Now, since A (t) -A * ≤ 2 A * (the closeness property (Def.1) is maintained at every step using Lemma 9), and further since A * = O( m/n), we have that

A (t) ≤ A (t) -A * + A * = O( m n ).

Therefore, we have

∆ (t) j + o( k m t ) ± γ = O( √ mq i,j p i1 t A (t) ).

Here, we use the fact that γ drops with decreasing t as argued in Lemma 5. Next, using ( 23), we have

EQUATION

we have that, for

k = O( √ n) v ≤ q i p i A (t) j -A * j .

Substituting for v , this implies that g (t)

j 2 ≤ 25ρ 2 A (t) j -A * j 2

. Further, we also have the following lower-bound

g (t) j , A (t) j -A * j ≥ 4ρ A (t) j -A * j 2 -v A (t) j -A * j .

Here, we use the fact that R.H.S. can be minimized only if v is directed opposite to the direction of A (t) j -A * j . Now, we show that this gradient is (ρ , 1/100ρ , 0) correlated,

EQUATION

EQUATION

A (t) i -A * i 2 -2ρ A (t) i -A * i 2 - 25ρ 2 A (t) i -A * i 2 100ρ , ≥ ρ A (t) i -A * i 2 ≥ 0.

Therefore, for this choice of k, i.e. k = O( √ n), there is no bias in dictionary estimation in comparison to Arora et al. (2015) . This gain can be attributed to estimating the coefficients simultaneously with the dictionary. Further, since we choose 4ρ = p j q j , we have that ρ = Θ(k/m), as a result ρ + = 1/100ρ = Ω(m/k). Applying Lemma 15 we have

A (t+1) j -A * j 2 ≤ (1 -ρ η A ) A (t) j -A * j 2 , for η A = O(m/k) with probability at least (1 -δ (t) T -δ (t) β -δ (t) gi ).

Proof of Lemma 8. Here, we will prove that g (t) defined as

g (t) = j (y (j) -A (t) x (j) )sign( x (j) ) ,

concentrates around its mean. Notice that each summand (y (j) -A (t) x (j) )sign( x (j) ) is a random matrix of the form (y -A (t) x)sign( x) . Also, we have g (t) defined as

g (t) = E[(y -A (t) x)sign( x) ].

To bound g (t) -g (t) , we are interested in p j=1 W j , where each matrix W j is given by

W j = 1 p (y (j) -A (t) x (j) )sign( x (j) ) -1 p E[(y -A (t) x)sign( x) ].

Noting that E[W j ] = 0, we will employ the matrix Bernstein result (Lemma 10) to bound g (t)g (t) . To this end, we will bound W j and the variance proxy

v(W j ) = max{ p j=1 E[W j W j ] , p j=1 E[W j W j ] }.

Bound on W j -First, we can bound both terms in the expression for W j by triangle inequality as

W j ≤ 1 p (y (j) -A (t) x (j) )sign( x (j) ) + 1 p E[(y -A (t) x)sign( x) , ≤ 2 p (y -A (t)

x)sign( x) . Here, we use Jensen's inequality for the second term, followed by upper-bounding the expected value of the argument by (y -A (t) x)sign( x) .

Next, using Claim 8 we have that with probability at least (1 -δ

(t) β -δ (t) T -δ (t) HW ), y -A (t) x is O(kt β ), and the fact that sign(x) T = √ k, W j ≤ 2 p √ k (y -A (t) x) = O( k √ k p t β ).

Bound on the variance statistic v(W j )-For the variance statistic, we first look at

E[W j W j ] , E[W j W j ] = 1 p 2 E[(y (j) -A (t) x (j) )sign( x (j) ) -E[(y -A (t) x)sign( x) ] × [sign( x (j) )(y (j) -A (t) x (j) ) -(E[(y -A (t) x)sign( x) ) ]. Since E[(y -A (t) x)sign( x) ]E[(y -A (t) x)sign( x) ] is positive semidefinite, E[W j W j ] 1 p 2 E[(y (j) -A (t) x (j) )sign( x (j) ) sign( x (j) )(y (j) -A (t) x (j) ) ].

Now, since each x (j) has k non-zeros, sign( x (j) ) sign( x (j) ) = k, and using Claim 10, with probability at least (1 -δ

(t) T -δ (t) β ) E[W j W j ] ≤ k p E[(y (j) -A (t) x (j) )(y (j) -A (t) x (j) ) ] ,

= O(

k 3 t 2 β pm ) A * 2 .

Similarly, expanding E[W j W j ], and using the fact that E[(y - t) x)sign( x) ] is positive semi-definite. Now, using Claim 8 and the fact that entries of E[(sign( x (j) )sign( x (j) ) ] are q i on the diagonal and zero elsewhere, where

EQUATION

q i = O(k/m), E[W j W j ] 1 p E[(sign( x (j) )(y (j) -A (t) x (j) ) (y (j) -A (t) x (j) )sign(x (R) (j) ) ] , ≤ 1 p E[(sign( x (j) )sign( x (j) ) ] y (j) -A (t) x (j) 2 , ≤ O( k mp ) O(k 2 t 2 β ) = O( k 3 t 2 β mp

). Now, we are ready to apply the matrix Bernstein result. Since, m = O(n) the variance statistic comes out to be O(

k 3 t 2 β pm ) A * 2 ,

then as long as we choose p = Ω(mk 2 ) (using the bound on t β ), with probability at least (1 -δ

(t) β -δ (t) T -δ (t) HW -δ (t) g ) g (t) -g (t) ≤ O( k √ k p t β ) + A * O( k 3 t 2 β pm ), = O * ( k m A * ).

where δ

(t) g = (n + m) exp(-Ω(m log(n)).

Proof of Lemma 9. This lemma ensures that the dictionary iterates maintain the closeness property (Def.1) and satisfies the prerequisites for Lemma 7.

The update step for the i-th dictionary element at the s + 1 iteration can be written as

A (t+1) i -A * i = A (t) i -A * i -η A g (t) i , = A (t) i -A * i -η A g (t) i -η A ( g (t) i -g (t) i ).

Here, g

i is given by the following as per Lemma 5 with probability at least (1 -δ

(t) T -δ (t) β ) g (t) i = q i p i (A (t) i -A * i ) + q i p i (-λ (t) i A (t) i + 1 qipi ∆ (t) i ± γ).

Substituting the expression for g

(t)

i in the dictionary update step,

A (t+1) i -A * i = (1 -η A p i q i )(A (t) i -A * i ) -η A p i q i λ (t) i A (t) i -η A ∆ (t) i -η A ( g (t) i -g (t) i ) ± γ,

where

∆ (t) j = E[A (t) ϑ (R) sign(x * j )] j .

Therefore, the update step for the dictionary (matrix) can be written as

EQUATION

where, U = A (t) diag(p i q i λ (t) i ) and V = A (t) Q, with the matrix Q given by,

Q i,j = q i,j E x * S [ϑ (R)

i sign(x * j )|S], and using the following intermediate result shown in Claim 7,

E x * S [ϑ (R) i sign(x * j )|S] ≤ γ (R) i , for i = j, = O(p j t ), for i = j, we have Q i = O( √ mq i,j p i t )

. Hence, we have

Q F ≤ O(mq i,j p i t ).

Therefore,

V ≤ A (t) Q ≤ A (t) Q F = O(mq i,j p i t A * ) = O( k 2 m log(n) ) A * .

We will now proceed to bound each term in (25). Starting with (A (t) -A * )diag(1 -η A p i q i ), and using the fact that p i = O(1), q i = O(k/m), and A (t) -A * ≤ 2 A * , we have

(A (t) -A * )diag(1 -η A p i q i ) ≤ (1 -min i η A p i q i ) (A (t) -A * ) ≤ 2(1 -Ω(η A k/m)) A * .

Next, since A (t) j

= 1, we have λ

(t) j = | A (t) j -A * j , A * j | = 1 2 A (t) j -A * j 2

, and λ

(t) i ≤ 2 t /2, therefore U = A (t) diag(p i q i λ (t) i ) ≤ max i p i q i 2 t 2 A (t) -A * + A * ≤ o(k/m) A * .

Using the results derived above, and the the result derived in Lemma 8 which states that with probability at least (1 -δ

(t) β -δ (t) T -δ (t) HW -δ (t) g ), g (t) -g (t) = O * ( k m A * )) we have A (t+1) -A * = (A (t) -A * )D (1-η A piqi) + η A U + η A V + η A ( g (t) -g (t) ) ± γ, ≤ 2(1 -Ω(η A k m ) A * + o(η A k m ) A * + O(η A k 2 m log(n) ) A * + o(η A k m A * ) ± γ, ≤ 2 A * . D APPENDIX: PROOFS OF INTERMEDIATE RESULTS Claim 1 (Incoherence of A (t) ). If A * ∈ R n×m is µ-incoherent and A * i -A (t) i ≤ t holds for each i ∈ [1 . . . m], then A (t) ∈ R n×m is µ t -incoherent, where µ t = µ + 2 √ n t .

Proof of Claim 1. We start by looking at the incoherence between the columns of A * , for j = i,

EQUATION

Since

A * i , A * j ≤ µ √ n , | A (t) i , A (t) j | ≤ A * i , A * j -A * i -A (t) i , A * j -A (t) i , A * j -A (t) j , ≤ µ √ n + 2 t . Claim 2 (Bound on β (t) j : the noise component in coefficient estimate that depends on t ). With probability (1 -δ (t) β ), |β (t) j | is upper-bounded by t β = O( √ k t ), where δ (t) β = 2k exp(-1 O( t) ).

Proof of Claim 2. We have the following definition for β

EQUATION

EQUATION

A (t) j -A * j , A * i )x * i .

Here, since x * i are independent sub-Gaussian random variables, β

j is a sub-Gaussian random variable with the variance parameter evaluated as shown below

var[β (t) j ] = i =j ( A * j , A (t) i -A * i + A (t) j -A * j , A (t) i + A (t) j -A * j , A * i ) 2 ≤ 9k 2 t .

Therefore, by Lemma 12

Pr[|β (t) j | > t β ] ≤ 2exp(- t 2 β 18k 2 t ).

Now, we need this for each β (t) j for j ∈ supp(x * ), union bounding over k coefficients

Pr[max |β (t) j | > t β ] ≤ δ (t) β ,

where δ

(t) β = 2k exp(- t 2 β 18k 2 t ). Choosing t β = O( √ k t ), we have that δ (t) β = 2k exp(-1 O( t) ).

Claim 3 (Error in coefficient estimation for a general iterate (r + 1)). The error in a general iterate r of the coefficient estimation is upper-bounded as

C (r+1) i1 ≤ α (r+1) i1 + (k -1)η x µt √ n r =1 α ( ) max 1 -η x + η x µt √ n r- + (k -1)η x µt √ n C (0) max 1 -η x + η x µt √ n r .

Proof of Claim 3 . From (17) we have the following expression for

C (r+1) i1 C (r+1) i1 ≤ α (r+1) i1 + µt √ n r+1 =1 η ( ) x i2 =i1 C ( -1) i2 r+1 q= (1 -η (q+1) x

).

Our aim will be to recursively substitute for C

( -1) i1

to develop an expression for C (r+1) i1

as a function of C 0 max . To this end, we start by analyzing the iterates C

(1)

EQUATION

i1 , and so on to develop an expression for C (r+1) i1 as follows.

Expression for C:

(1)

EQUATION

Expression for C

(2)

EQUATION

i1 is given by

C (2) i1 ≤ α (2) i1 + η x µt √ n 2 =1 i2 =i1 C ( -1) i2 2 q= (1 -η x ), ≤ α (2) i1 + η x µt √ n i2 =i1 C (1) i2 + i2 =i1 C (0) i2 (1 -η x ) .

Further, we know from (26) we have

C (1) i2 = α (1) i2 + η x µt √ n i3 =i2 C (0) i3 .

Therefore, since

EQUATION

i1 ≤ α

(2)

EQUATION

Expression for C

(3)

i1 -Next, we writing C

(3)

i1 + η x µt √ n 3 =1 i2 =i1 C ( -1) i2 (1 -η x ) 3-, = α (3) i1 + η x µt √ n i2 =i1 C (0) i2 (1 -η x ) 2 + C (1) i2 (1 -η x ) + C (2) i2 , ≤ α (3) i1 + η x µt √ n i2 =i1 C (0) i2 (1 -η x ) 2 + α (1) i2 + η x µt √ n i3 =i2 C (0) i3 (1 -η x ) + C (2)

i2 .

Here, using (27) we have the following expression for C

(2) i2

C:

(2)

i2 ≤ α (2) i2 + η x µt √ n i3 =i2 α (1) i3 + η x µt √ n η x µt √ n i4 =i3,i2 C (0) i4 + i3 =i2 C (0) i3 (1 -η x ) .

Substituting for C:

(2) i2 in the expression for C

(3) i1 , and rearranging the terms in the expression for C

(3) i1 , we have

EQUATION

i4 .

(28)

Expression for C (4) i1 -Now, consider C (4) i1 C (4) i1 ≤ α (4) i1 + η x µt √ n 4 =1 i2 =i1 C ( -1) i2 (1 -η x ) 4-, ≤ α (4) i1 + η x µt √ n i2 =i1 C (0) i2 (1 -η x ) 3 + i2 =i1 C (1) i2 (1 -η x ) 2 + i2 =i1 C (2) i2 (1 -η x ) 1 + i2 =i1 C (3) i2 (1 -η x ) 0 .

(3) i2 from (28), C

(2)

EQUATION

i2 using (26), and rearranging,

C (4) i1 ≤ α (4) i1 + η x µt √ n i2 =i1 α (3) i2 + (1 -η x ) 1 i2 =i1 α (2) i2 + η x µt √ n i3 =i2,i1 α (2) i3 + i2 =i1 α (1) i2 (1 -η x ) 2 + 2η x µt √ n (1 -η x ) i3 =i2,i1 α (1) i3 + (η x µt √ n ) 2 i4 =i3,i2,i1 α (1) i4 + η x µt √ n i2 =i1 C (0) i2 (1 -η x ) 3 + 3η x µt √ n (1 -η x ) 2 i3 =i2,i1 C (0) i3 + 3(η x µt √ n ) 2 (1 -η x ) 1 i4 =i3,i2,i1 C (0) i4 + (η x µt √ n ) 3 i5 =i4,i3,i2,i1 C (0) i5 .

Notice that the terms have a binomial series like form. To reveal this structure, let each α

EQUATION

max for j = i 1 , i 2 , . . . , i k . Therefore, we have

C (4) i1 ≤ α (4) i1 + η x µt √ n (k -1)α (3) i + α (2) i (1 -η x ) 1 (k -1) + η x µt √ n (k -2) + α (1) i (k -1)(1 -η x ) 2 + 2(k -2)η x µt √ n (1 -η x ) + (k -3)(η x µt √ n ) 2 + η x µt √ n C (0) i (k -1)(1 -η x ) 3 + 3(k -2)η x µt √ n (1 -η x ) 2 + 3(k -3)(η x µt √ n ) 2 (1 -η x ) 1 + (k -4)(η x µt √ n ) 3 .

Further upper-bounding the expression, we have

C (4) i1 ≤ α (4) i1 + (k -1)η x µt √ n α (3) i + α (2) i (1 -η x ) + η x µt √ n + α (1) i (1 -η x ) 2 + 2η x µt √ n (1 -η x ) + (η x µt √ n ) 2 + (k -1)η x µt √ n C (0) i (1 -η x ) 3 + 3η x µt √ n (1 -η x ) 2 + 3(η x µt √ n ) 2 (1 -η x ) + (η x µt √ n ) 3 .

Therefore,

C (4) i1 ≤ α (4) i1 + (k -1)η x µt √ n α (3) i +α (2) i 1 -η x + η x µt √ n 1 + α (1) i 1 -η x + η x µt √ n 2 + (k -1)η x µt √ n C (0) i 1 -η x + η x µt √ n 3 . (29) Expression for C (r+1) i1

-With this, we are ready to write the general term,

C (r+1) i1 ≤ α (r+1) i1 + (k -1)η x µt √ n r =1 α ( ) max 1 -η x + η x µt √ n r- + (k -1)η x µt √ n C (0) max 1 -η x + η x µt √ n r .

Claim 4 (An intermediate result for bounding the error in coefficient calculations). With probability (1 -δ

(t) T -δ (t) β ), R =1 α ( ) max 1 -η x + η x µt √ n R-≤ C (0) i Rδ R + 1 ηx(1- µt √ n ) ( 2 t 2 |x * max | + t β ).

Proof of Claim 4 . Using (18), the quantity α

( ) i is defined as α ( ) i = C (0) i (1 -η x ) + (λ (t) i |x * i | + β (t) i ) s=1 η x (1 -η x ) -s+1 .

Published as a conference paper at ICLR 2019 Therefore, we are interested in

R =1 C (0) i (1 -η x ) 1 -η x +η x µt √ n R- + (λ (t) i |x * i | + β (t) i ) R =1 1 -η x + η x µt √ n R- s=1 η x (1 -η x ) -s+1 .

Consider the first term which depends on

C (0) i . Since (1 -η x ) ≤ (1 -η x + η x µt √ n , we have C (0) i R =1 (1 -η x ) 1 -η x + η x µt √ n R-≤ C (0) i R 1 -η x + η x µt √ n R ≤ C (0) i Rδ R ,

where δ R is a small constant, and a parameter which determines the number of iterations R required for the coefficient update step. Now, coming back to the quantity of interest

R =1 α ( ) i 1 -η x + η x µt √ n R-≤ C (0) i Rδ R + (λ (t) i |x * i | + β (t) i ) R =1 1 -η x + η x µt √ n R- s=1 η x (1 -η x ) -s+1 .

Now, using sum of geometric series result, we have that s=1 η x (1 -η x ) -s+1 , and

R =1 1 -η x + η x µt √ n R-= 1-1-ηx+ηx µt √ n R ηx-ηx µt √ n ≤ 1 ηx(1- µt √ n )

.

Therefore, with probability at least (1 -δ

(t) β ), R =1 α ( ) max 1 -η x + η x µt √ n R-≤ C (0) i Rδ R + 1 ηx(1- µt √ n ) ( 2 t 2 |x * max | + t β ),

where λ Claim 5 (Bound on the noise term in the estimation of a coefficient element in the support).

With probability (1 -δ

(t) β ), each entry ϑ (R) i1 of ϑ (R) is upper-bounded as |ϑ (R) i1 | ≤ O(t β ). Proof of Claim 5. From (19) ϑ (R)

i1 is defined as

ϑ (R) i1 := R r=1 η x ξ (r) i1 (1 -η x ) R-r + γ (R) i1 ,

where γ (R)

EQUATION

i1 is as defined in (10),

ξ (r) i1 = β (t) i1 + i2 =i1 | A (t) i1 , A (t) i2 |sign( A (t) i1 , A (t) i2 )C (r-1) i2 sign(x * i2 -x (r) i2 ).

Therefore, we have the following expression for ϑ

EQUATION

) Now ϑ (R)

i1 can be upper-bounded as

ϑ (R) i1 ≤ β (t) i1 R r=1 η x (1 -η x ) R-r + η x µt √ n R r=1 i2 =i1 C (r-1) i2 (1 -η x ) R-r + γ (R) i1 , ≤ β (t) i1 + (k -1)η x µt √ n R r=1 C (r-1) i2 (1 -η x ) R-r + γ (R) i1 .

Since from Claim 6 we have

C (r-1) i2 (1 -η x ) R-r ≤ (λ (t) max |x * max | + β (t) max ) r-1 s=1 η x (1 -η x ) R-s + kc x (1 -η x ) R-r + kη x µt √ n C (0) max δ R-2 . Further, since 1 -(1 -η x ) r-1 ≤ 1, we have that R r=1 r-1 s=1 η x (1 -η x ) R-s = R r=1 η x (1 -η x ) R-r+1 1-(1-ηx) r-1 ηx ≤ R r=1 (1 -η x ) R-r+1 ≤ 1 ηx .

Therefore,

|ϑ (R) i1 | ≤ |β (t) i1 | + (k -1) µt √ n (λ (t) max |x * max | + |β (t) max |)(1 + kc x ) + kη x µt √ n 2 RC (0) max δ R-2 + γ (R) i1 . Now, since each |β (t) i | = t β with probability at least (1 -δ (t)

β ) for the t-th iterate, and

k = O * ( √ n µ log(n) ), therefore kc x < 1, we have that |ϑ (R) i1 | ≤ O(t β ).

with probability at least (1 -δ

(t) β ).

Claim 6 (An intermediate result for

ϑ (R) i1 calculations). For c x = µt √ n /(1 -µt √ n ), we have C (r-1) i2 (1 -η x ) R-r ≤ (λ (t) max |x * max | + β (t) max ) r-1 s=1 η x (1 -η x ) R-s + kc x (1 -η x ) R-r + kη x µt √ n C (0) max δ R-2 .

Proof of Claim 6. Here, from Claim 3 we have that for any i 1 ,

C (r+1) i1 ≤ α (r+1) i1 + kη x µt √ n r =1 α ( ) max 1 -η x + η x µt √ n r-+ kη x µt √ n C (0) max 1 -η x + η x µt √ n r . therefore C (r-1) i2

is given by

C (r-1) i2 ≤ α (r-1) i2 + kη x µt √ n r-2 =1 α ( ) max 1 -η x + η x µt √ n r--2 +kη x µt √ n C (0) max 1 -η x + η x µt √ n r-2 .

Further, the term of interest C

(r-1) i2

(1 -η x ) R-r can be upper-bounded by

C (r-1) i2 (1 -η x ) R-r ≤ α (r-1) i2 (1 -η x ) R-r +(1 -η x ) R-r kη x µt √ n r-2 =1 α ( ) max 1 -η x + η x µt √ n r--2 + kη x µt √ n C (0) max 1 -η x + η x µt √ n r-2 (1 -η x ) R-r .

From the definition of α can be written as

α (r-1) i2 = C (0) max (1 -η x ) r-1 + (λ (t) max |x * max | + β (t) max ) r-1 s=1 η x (1 -η x ) r-s .

Therefore, we have

α (r-1) i2 (1 -η x ) R-r = C (0) max (1 -η x ) R-1 + (λ (t) max |x * max | + β (t) max ) r-1 s=1 η x (1 -η x ) R-s .

Next, to get a handle on α

( ) max 1 -η x + η x µt √ n

r--2 , consider the following using the definition of α ( ) i from (18), where η

(i) x = η x for all i, r =1 α ( ) max 1-η x + η x µt √ n r-= r =1 C (0) max (1 -η x ) 1 -η x + η x µt √ n r- + (λ (t) max |x * max | + β (t) max ) r =1 1 -η x + η x µt √ n r- s=1 η x (1 -η x ) -s+1 , ≤ r =1 C (0) max 1 -η x + η x µt √ n r + (λ (t) max |x * max | + β (t) max ) r =1 1 -η x + η x µt √ n r-.

Therefore,

(1 -η x ) R-r r-2 =1 α ( ) max 1 -η x + η x µt √ n r--2 ≤ r-2 =1 C (0) max 1 -η x + η x µt √ n r-2 (1 -η x ) R-r + (λ (t) max |x * max | + β (t) max )(1 -η x ) R-r r-2 =1 1 -η x + η x µt √ n r--2 , ≤ (R -2)C (0) max 1 -η x + η x µt √ n R-2 + (λ (t) max |x * max | + β (t) max ) (1-ηx) R-r ηx(1- µt √ n )

.

Therefore,

(1 -η x ) R-r r-2 =1 α ( ) max 1 -η x + η x µt √ n r--2 ≤ (r -2)C (0) max 1 -η x + η x µt √ n R-2 + (λ (t) max |x * max | + β (t) max ) (1-ηx) R-r ηx(1- µt √ n ) . Further, since (1 -η x ) ≤ (1 -η x + η x µt √ n ), kη x µt √ n C (0) max 1 -η x + η x µt √ n r-2 (1 -η x ) R-r ≤ kη x µt √ n C (0) max δ R-2 .

Therefore, combining all the results we have that, for a constant

c x = µt √ n /(1 -µt √ n ), C (r-1) i2 (1 -η x ) R-r ≤ (λ (t) max |x * max | + β (t) max ) r-1 s=1 η x (1 -η x ) R-s + kc x (1 -η x ) R-r + kη x µt √ n C (0) max δ R-2 .

Claim 7 (Bound on the noise term in expected gradient vector estimate). ∆ (t) j

where ∆ (t)

j := E[A (t) ϑ (R) sign(x * j )] is upper-bounded as, ∆ (t) j = O( √ mq i,j p j t A (t) )].

Proof of Claim 7.

∆ (t) j = E[A (t) ϑ (R) sign(x * j )] = E S [A (t) S E x * S [ϑ (R) S sign(x * j )|S]]

From (30) we have the following definition for ϑ

(R) j ϑ (R) j = β (t) j + R r=1 η x i =j | A (t) j , A (t) i |sign( A (t) j , A (t) i )C (r-1) i sign(x * i -x (r) i )(1 -η x ) R-r + γ (R) j ,

where β (t) j is defined as the following ( 11)

β (t) j = i =j ( A * j , A * i -A (t) i + A * j -A (t) j , A (t) i )x * i + i =j A (t) j -A * j , A * i x * i . Consider E x * S [ϑ (R) S sign(x * j )|S], where ϑ (R)

S is a vector with each element as defined in (30). Therefore, the elements of the vector

E x * S [ϑ (R)

S sign(x * j )|S]] are given by

E x * S [ϑ (R) i sign(x * j )|S] = E x * S [ϑ (R) i sign(x * j )|S], for i = j, E x * S [ϑ (R) j sign(x * j )|S], for i = j. Consider the general term of interest E x * S [ϑ (R) i sign(x * j )|S] ≤ R r=1 η x (1 -η x ) R-r E x * S [β i sign(x * j )|S] ♣ + µt √ n R r=1 η x (1 -η x ) R-r s =i E x * S [C (r-1) s sign(x * s -x (r) s )sign(x * j )|S] ♠ +γ (R) i . Further, since E x * S [x * i sign(x * j )|S] = 0, for i = j, p j , for i = j,

we have that

EQUATION

Further, for

♠ s := E x * S [C (r-1) s sign(x * s -x (r) s )sign(x * j )|S] we have that ♠ s = E x * S [C (r-1) j (x * j -x (r-1) j )sign(x * j )|S] ≤ C (r-1) j

, for s = j 0, for s = j.

In addition, for s =i ♠ s we have that

EQUATION

Therefore, using the results for ♣ and

s =i ♠ s , we have that E x * S [ϑ (R) j sign(x * j )|S] = γ (R) i for i = j,

and for i = j we have

EQUATION

Here, from Claim 6, for c

x = µt √ n /(1 -µt √ n ) we have C (r-1) j (1 -η x ) R-r ≤ (λ (t) max |x * max | + β (t) max ) r-1 s=1 η x (1 -η x ) R-s + kc x (1 -η x ) R-r + kη x µt √ n C (0) max δ R-2 .

Further, due to our assumptions on sparsity, kc x ≤ 1; in addition by Claim 2, and with probability at least (1 -δ

(t) β ) we have |β (t) max | ≤ t β , substituting, R r=1 C (r-1) j η x (1 -η x ) R-r ≤ (λ (t) max |x * max | + β (t) max ) R r=1 η x r-1 s=1 η x (1 -η x ) R-s + kc x R r=1 η x (1 -η x ) R-r , ≤ (λ (t) max |x * max | + t β )(1 + kc x ), = O(t β ),

with probability at least (1 -δ (t) β ). Combining results from (31), (32) and substituting for the terms in (33) using the analysis above,

E x * S [ϑ (R) i sign(x * j )|S] ≤ γ (R) i , for i = j, ≤ 3p j t + µ √ n t β + γ (R) i = O(p j t ), for i = j. Note that since γ (R) i := (1 -η x ) R (x (0) i -x * i (1 -λ (t) i )

) can be made small by choice of R. Also, since Pr[i, j ∈ S] = q i,j , we have

∆ (t) j = E S [A (t) S E x * S [ϑ (R) S sign(x * j )|S]] , ≤ O( √ mq i,j p j t A (t) ).

Claim 8 (An intermediate result for concentration results). With probability (1 -δ

(t) β -δ (t) T - δ (t) HW ) y -A (t) x is upper-bounded by O(kt β ) .

Proof of Claim 8. First, using Lemma 4 we have

x i1 := x (R) i1 = x * i1 (1 -λ (t) i1 ) + ϑ (R)

i1 . Therefore, the vector x S , for S ∈ supp(x * ) can be written as

EQUATION

where x has the correct signed-support with probability at least (1 -δ T ) using Lemma 2. Using this result, we can write y -A (t) x as

y -A (t) x = A * S x * S -A (t) S (I -Λ (t) S )x * S -A (t) S ϑ (R) S . Now, since Λ (t) ii ≤ 2 t 2 we have y -A (t) x ≤ A * S x * S -(1 - 2 t 2 )A (t) S x * S -A (t) S ϑ (R) S , = ((1 - 2 t 2 )(A * S -A (t) S ) + 2 t 2 A * S ) ♣ x * S -A (t) S ϑ (R) S ♠ .

With x * S being independent and sub-Gaussian, using Lemma 13, which is a result based on the Hanson-Wright result (Hanson and Wright, 1971) for sub-Gaussian random variables, and since

A (t) S -A * S ≤ A (t) S -A * S F ≤ √ k t , we have that with probability at least (1 -δ (t) HW ) ♣x * S = ((1 - 2 t 2 )(A * S -A (t) S ) + 2 t 2 A * S )x * S ≤ O( (1 - 2 t 2 )(A * S -A (t) S ) + 2 t 2 A * S F ),

where δ

(t) HW = exp(-1 O( t) ). Now, consider the ♣ F , since A (t) S -A * S F ≤ √ k t ♣ F := (1 - 2 t 2 )(A * S -A (t) S ) + 2 t 2 A * S F ≤ (1 - 2 t 2 ) (A * S -A (t) S ) F + 2 t 2 A * S F , ≤ √ k(1 - 2 t 2 ) t + 2 t 2 A * S F .

Consider the ♠ term. Using Claim 5, each ϑ (R) j is bounded by O(t β ). with probability at least

(1 -δ (t) β ) Therefore, ♠ = A (t) S ϑ (R) S ≤ A (t) S ϑ (R) S = A (t) S √ kO(t β ). Again, since A (t) S -A * S ≤ A (t) S -A * S F ≤ √ k t , A (t) S ≤ A (t) S -A * S + A * S ≤ A (t) S -A * S + A * S ≤ √ k t + 2.

Finally, combining all the results and using the fact that A

* S F ≤ √ k A * S ≤ 2 √ k, , y -A (t) x = O( √ k(1 - 2 t 2 ) t + 2 t √ k) + A (t) S √ kO(t β ), = O(kt β ).

Claim 9 (Bound on variance parameter for concentration of gradient vector). For z := (y -

A (t) x)sign( x i )|i ∈ S the variance parameter E[ z 2 ] is bounded as E[ z 2 ] = O(k 2 t ) + O(kt 2 β ) with probability at least (1 -δ (t) β -δ (t) T ).

Proof of Claim 9. For the variance E[ z2 ], we focus on the following,

E[ z 2 ] = E[ (y -A (t) x)sign( x i ) 2 |i ∈ S].

Here, x S is given by

EQUATION

We will now consider each term in (35) separately. We start with ♥. Since x * S s are conditionally independent of S, E[x * S x * S ] = I. Therefore, we can simplify this expression as

♥ := E[ (A * S -A (t) S (I -Λ (t) S ))x * S 2 |i ∈ S] = E[ A * S -A (t) S (I -Λ (t) S ) 2 F |i ∈ S].

Rearranging the terms we have the following for ♥,

EQUATION

) For ♥ 1 , since A (t) S ≤ √ k t + 2, we have ♥ 1 := E[ A * S Λ (t) S 2 F |i ∈ S] ≤ E[ A * S Λ (t) S 2 F |i ∈ S] ≤ A * S j∈S (λ (t) j ) 2 ≤ k( √ k t + 2) 4 t 4 . Next, since (1 -λ (t) j ) ≤ 1, we have the following bound for ♥ 2 ♥ 2 := E[ (A * S -A (t) S )(I -Λ (t) S ) 2 F |i ∈ S] ≤ E[ A * S -A (t) S 2 F |i ∈ S] ≤ A * S -A (t) S 2 F ≤ k 2 t .

Further, ♥ 3 can be upper-bounded by using bounds for ♥ 1 and ♥ 2 . Combining the results of upperbounding ♥ 1 , ♥ 2 , and ♥ 3 we have the following for (36)

♥ ≤ E[ (A * S -A (t) S (I -Λ (t) S ))x * S Next, by Claim 5, ϑ (R) j is upper-bounded as |ϑ (R) j | ≤ O(t β ). with probability (1 -δ (t)

β ). Therefore, the term ♦, the second term of (35), can be bounded as

♦ ≤ A (t) S ϑ (R) S sign( x i ) 2 ≤ ( √ k t + 2) 2 kO(t β ) 2 = O(kt 2 β ).

Finally, combining all the results, the term of interest in (35) has the following form

E[ (y -A (t) x)sign( x i ) 2 |i ∈ S] = O(k 2 t ) + O(kt 2 β ).

Claim 10 (Bound on variance parameter for concentration of gradient matrix). With probability (1 -δ

(t) T -δ (t) β ), the variance parameter E[(y -A (t) x)(y -A x) ] is upper-bounded by O( k 2 t 2 β m ) A * 2 .

Proof of Claim 10. Let F x * be the event that sign(x * ) = sign( x), and let 1 F x * denote the indicator function corresponding to this event. As we show in Lemma 2, this event occurs with probability at least (1 -δ

(t) β -δ (t) T ), therefore, E[(y -A (t) x)(y -A (t) x) ] = E[(y -A (t) x)(y -A (t) x) 1 F x * ] + E[(y -A (t) x)(y -A (t) x) 1 Fx * ], = E[(y -A (t) x)(y -A (t) x) 1 F x * ] ± γ.

Here, γ is small. Under the event F x * , x has the correct signed-support. Again, since 1 S . Therefore we have

F x * = 1 -1 Fx * , E[(y -A (t) x)(y -A (t) x) ] = E[(y -A (t) x)(y -A (t) x) (1 -1 Fx * )] ± γ, = E[(y -A (t) x)(y -A (t) x) ] ± γ.

y -A (t) x = (A * S -A (t) S (I -Λ (t) S ))x * S -A (t) S ϑ (R) S .

Using the expression above E[(y -A (t) x)(y -A (t) x) ] can be written as ] in (38) which can also be bounded similarly as

EQUATION

♣ 4 = O( k 2 m ) A * 2 .

Therefore, we have the following for ♣ in (37)

EQUATION

)

Consider ♠ in (37). Letting M = E x * S [ϑ (R) S x * S |S], and using the analysis similar to that shown in 7, we have that elements of M ∈ R k×k are given by

M i,j = E x * S [ϑ (R) i x * j |S] ≤ O(γ (R) i

), for i = j, = O( t ), for i = j.

We have the following, .

EQUATION

Again, using the result on |ϑ (R)

i1 |, we have

♦ := E S [A (t) S NA (t) S ] = m max i,j |ϑ (R) i ||ϑ (R) j |O( k 2 m 2 ) A (t) A (t) = O( k 2 t 2 β m ) A * 2 .

Combining all the results for ♣, ♠, ♥ and ♦, we have,

E[(y -A (t) x)(y-A (t) x) ] = O( k 2 m ) A * 2 + O( k 2 m ) t A * 2 + O( k 2 m ) t A * 2 + O( k 2 t 2 β m ) A * 2 , = O( k 2 t 2 β m ) A * 2 .

E ADDITIONAL EXPERIMENTAL RESULTS:

We now present some additional results to highlight the features of NOODL. Specifically, we compare the performance of NOODL (for both dictionary and coefficient recovery) with the state-of-theart provable techniques for DL presented in Arora et al. (2015) (when the coefficients are recovered via a sparse approximation step after DL)3 . We also compare the performance of NOODL with the popular online DL algorithm in Mairal et al. (2009) , denoted by Mairal '09. Here, the authors show that alternating between a 1 -based sparse approximation and dictionary update based on block co-ordinate descent converges to a stationary point, as compared to the true factors in case of NOODL.

Data Generation: We generate a (n = 1000) × (m = 1500) matrix, with entries drawn from N (0, 1), and normalize its columns to form the ground-truth dictionary A * . Next, we perturb A * with random Gaussian noise, such that the unit-norm columns of the resulting matrix, A (0) are 2/ log(n) away from A * , in 2 -norm sense, i.e., 0 = 2/ log(n); this satisfies the initialization assumptions in A.4. At each iteration, we generate p = 5000 samples Y ∈ R 1000×5000 as Y = A * X * , where X * ∈ R m×p has at most k = 10, 20, 50, and 100, entries per column, drawn from the Radamacher distribution. We report the results in terms of relative Frobenius error for all the experiments, i.e., for a recovered matrix M, we report M -M * F / M * F . To form the coefficient estimate for Mairal '09 via Lasso (Tibshirani, 1996) we use the FISTA (Beck and Teboulle, 2009) algorithm by searching across 10 values of the regularization parameter at each iteration. Note that, although our phase transition analysis for NOODL shows that p = m suffices, we use p = 5000 in our convergence analysis for a fair comparison with related techniques. 2015)) as XHT. For the techniques presented in Arora et al. (2015) , we scan across 50 values of the regularization parameter for coefficient estimation using Lasso after learning the dictionary (A), and report the optimal estimation error for the coefficients (XLasso), while for Mairal '09, at each step the coefficients estimate is chosen by scanning across 10 values of the regularization parameters. For k = 100, the algorithms of Arora et al. (2015) do not converge (shown as N/A).

Technique:

Recovered As shown in Table 4 , in comparison to NOODL the techniques described in Arora et al. (2015) still incur a large error at convergence, while the popular online DL algorithm of Mairal et al. (2009) exhibits very slow convergence rate. Combined with the convergence results shown in Fig. 2 , we observe that due to NOODL's superior convergence properties, it is overall faster and also geometrically converges to the true factors. This again highlights the applicability of NOODL in practical applications, while guaranteeing convergence to the true factors.

F APPENDIX: STANDARD RESULTS:

Definition 6 (sub-Gaussian Random variable). Let x ∼ subGaussian(σ 2 ). Then, for any t > 0, it holds that Pr[|x| > t] ≤ 2 exp t 2 2σ 2 .

F.1 CONCENTRATION RESULTS:

Lemma 10 (Matrix Bernstein (Tropp, 2015) ). Consider a finite sequence W k ∈ R n×m of independent, random, centered matrices with dimension n. Assume that each random matrix satisfies E[W k ] = 0 and W k ≤ R almost surely. Then, for all t ≥ 0,

Pr k W k ≥ t ≤ (n + m)exp -t 2 /2 σ 2 +Rt/3 ,

where

σ 2 := max{ k E[W k W k ] , k E[W k W k ] }. Furthermore, E[ k W k ] ≤ 2σ 2 log(n + m) + 1 3 R log(n + m).

:

The associated code is made available at https://github.com/srambhatla/NOODL.

µt √ n ) at each iteration r. Therefore, for a decay parameter δ R , we can choose the number of IHT iterations R, to make this component arbitrarily small. Therefore, the error in the coefficients only depends on the per column error in the dictionary, formalized by the following result.

|i ∈ S] = O(k 2 t ).

The associated code is made available at https://github.com/srambhatla/NOODL.

We use the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA)(Beck and Teboulle, 2009), which is among the most efficient algorithms for solving the 1-regularized problems. Note that, in our experiments we fix the step-size for FISTA as 1/L, where L is the estimate of the Lipschitz constant (since A is not known exactly).

Note that, although scanning across 50 values of the regularization parameter for this case would have led to better coefficient estimates and dictionary recovery, we choose 10 values for this case since it is very expensive to scan across 50 of regularization parameter at each step. This also highlights why Mairal '09 may be prohibitive for large scale applications.

When the dictionary is not known exactly, the guarantees may exist on coefficient recovery only in terms of closeness in 2-norm sense, due to the error-in-variables (EIV) model for the dictionary(Fuller, 2009;Wainwright, 2009).

Figures and Tables:

The coefficients are drawn from the distribution class D, as per Def. 4; A.3 The sparsity k satisfies k = O * ( √ n/µ log(n)); A.4 A (0) is ( 0 , 2)-close to A * as per Def. 1, and 0 = O * (1/ log(n)); A.5 The step-size for dictionary update satisfies η A = Θ(m/k); A.6 The step-size and threshold for coefficient estimation satisfies η

Figure 1: A neural implementation of NOODL. Panel (a) shows the neural architecture, which consists of three layers: an input layer, a weighted residual evaluation layer (evaluates ηx y (j) -A (t) x (r) (j) ), and an output layer. Panel (b) shows the operation of the neural architecture in panel (a). The update of x (r+1) (j)

= 10, η A = 30 k = 20, η A = 30 k = 50, η A = 15 k = 100, η A Comparative analysis of convergence properties. Panels (a-i), (b-i), (c-i),

t β with probability at least (1 -δ (t) β ) using Claim 2.

t) S E x * S [ϑ (R) S x * S |S](A * S -(I -Λ Therefore, since E S [1 i,j∈S |S] = O( k 2 m 2 ),and1 m×m = m, ♠ := E S [A (t) S M(A * S -(I -Λ

Frequently used symbols

Table content: <table><tr><td>Dictionary Related</td></tr><tr><td>Symbol Definition</td></tr><tr><td>(t) A i</td></tr></table>

Proof map: dependence of results.

Table content: <table><tr><td colspan="2">Lemmas Result</td><td/><td colspan="2">Dependence</td></tr><tr><td colspan="3">Lemma Signed-support recovery by co-efficient initialization step</td><td>-</td><td/><td>Claim 10</td></tr><tr><td colspan="3">Lemma IHT update step preserves the</td><td>Claim</td><td>1,</td></tr><tr><td/><td colspan="2">correct signed-support</td><td>Lemma</td><td>1,</td></tr><tr><td/><td/><td/><td>and Claim 2</td><td/></tr><tr><td colspan="3">Lemma Upper-bound on the error in co-</td><td>Claim</td><td>1,</td><td>Claim 9</td></tr><tr><td/><td colspan="2">efficient estimation</td><td>Claim</td><td>2,</td></tr><tr><td/><td/><td/><td colspan="2">Claim 3, and</td></tr><tr><td/><td/><td/><td>Claim 4</td><td/></tr><tr><td colspan="3">Lemma Expression for the coefficient</td><td>Claim 5</td><td/><td>Claim 8</td><td>Lemma 9</td></tr><tr><td/><td colspan="2">estimate at the end of R-th IHT</td><td/><td/></tr><tr><td/><td>iteration</td><td/><td/><td/></tr><tr><td colspan="3">Lemma Expression for the expected gra-</td><td colspan="2">Lemma 4 and</td></tr><tr><td colspan="3">dient vector Lemma Concentration of the empirical</td><td colspan="2">Claim 7 Claim 8 and</td><td>Claim 7</td><td>Lemma 8</td></tr><tr><td/><td colspan="2">gradient vector</td><td>Claim 9</td><td/></tr><tr><td colspan="3">Lemma Empirical gradient vector is cor-</td><td>Lemma</td><td>5,</td></tr><tr><td/><td colspan="2">related with the descent direc-</td><td colspan="2">Claim 7 and</td></tr><tr><td/><td>tion</td><td/><td>Lemma 6</td><td/><td>Claim 6</td><td>Lemma 7</td></tr><tr><td colspan="3">Lemma Concentration of the empirical</td><td colspan="2">Claim 8 and</td></tr><tr><td/><td colspan="2">gradient matrix</td><td>Claim 10</td><td/></tr><tr><td colspan="3">Lemma A (t+1) maintains closeness</td><td>Lemma</td><td>5,</td></tr><tr><td/><td/><td/><td colspan="2">Claim 7 and</td><td>Claim 5</td><td>Lemma 6</td><td>Theorem 1</td></tr><tr><td/><td/><td/><td>Lemma 8</td><td/></tr><tr><td>Claims</td><td>Result</td><td/><td colspan="2">Dependence</td></tr><tr><td colspan="3">Claim 1 Incoherence of A (t)</td><td>-</td><td/></tr><tr><td colspan="2">Claim 2 Bound on β</td><td>(t) j : the noise com-</td><td>-</td><td/><td>Claim 4</td><td>Lemma 5</td></tr><tr><td/><td colspan="2">ponent in coefficient estimate</td><td/><td/></tr><tr><td/><td colspan="2">that depends on t</td><td/><td/></tr><tr><td colspan="3">Claim 3 Error in coefficient estimation</td><td>-</td><td/></tr><tr><td colspan="3">for a general iterate (r + 1) Claim 4 An intermediate result for</td><td>Claim 2</td><td/><td>Claim 3</td><td>Lemma 4</td></tr><tr><td/><td colspan="2">bounding the error in coefficient</td><td/><td/></tr><tr><td/><td>calculations</td><td/><td/><td/></tr><tr><td colspan="3">Claim 5 Bound on the noise term in the</td><td>Claim 6</td><td/></tr><tr><td/><td colspan="2">estimation of a coefficient ele-</td><td/><td/><td>Claim 2</td><td>Lemma 3</td></tr><tr><td/><td colspan="2">ment in the support</td><td/><td/></tr><tr><td colspan="3">Claim 6 An intermediate result for ϑ (R) i1</td><td>Claim 3</td><td/></tr><tr><td/><td>calculations</td><td/><td/><td/></tr><tr><td colspan="3">Claim 7 Bound on the noise term in ex-</td><td colspan="2">Claim 6 and</td><td>Claim 1</td></tr><tr><td/><td colspan="2">pected gradient vector estimate</td><td>Claim 2</td><td/></tr><tr><td colspan="3">Claim 8 An intermediate result for con-</td><td>Lemma</td><td>2</td></tr><tr><td/><td colspan="2">centration results</td><td colspan="2">,Lemma 4 and</td></tr><tr><td/><td/><td/><td>Claim 5</td><td/><td>Lemma 1</td></tr><tr><td colspan="3">Claim 9 Bound on variance parameter</td><td>Claim 5</td><td/></tr><tr><td/><td colspan="2">for concentration of gradient</td><td/><td/></tr><tr><td/><td>vector</td><td/><td/><td/></tr><tr><td colspan="3">Claim Bound on variance parameter for concentration of gradient</td><td colspan="2">Lemma 2 , Lemma 4 and</td><td>Lemma 2</td></tr><tr><td/><td>matrix</td><td/><td>Claim 5</td><td/></tr></table>

Similarly, ♣ 3 (38) is also O( k 2 m ) A * 2 . Next, we consider ♣ 4 := E S [A

Table content: <table><tr><td/><td/><td/><td/><td/><td/><td colspan="3">(t) S E x  *  S [ϑ</td><td>(R) S x  *  S |S](A  *  S -(I -Λ (t) S )A</td><td>(t) S</td><td>)]</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>♠</td></tr><tr><td/><td/><td/><td colspan="6">+ E S [(A  *  S -A</td><td>(t) S (I -Λ</td><td>(t) S ))E x  *  S [x  *  S (ϑ</td><td>(R) S ) |S]A</td><td>(t) S</td><td>]</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>♥</td></tr><tr><td/><td/><td/><td colspan="3">+ E S [A</td><td colspan="3">(t) S E x  *  S [ϑ</td><td>(R) S (ϑ</td><td>(R) S ) |S]A</td><td>(t) S</td><td>]</td><td>.</td><td>(37)</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>♦</td></tr><tr><td colspan="9">Let's start with the first term (♣) of (37), which can be written as</td></tr><tr><td>♣ :≤ E S [A  *  S A  *  S ]</td><td colspan="7">+ E S [A  *  S (I -Λ (t) S )A</td><td>(t) S</td><td>]</td><td>+ E S [A</td><td>(t) S (I -Λ</td><td>(t) S )A  *  S ]</td></tr><tr><td>♣1</td><td/><td/><td/><td colspan="2">♣2</td><td/><td/><td>♣3</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>+ E S [A</td><td>(t) S (I -Λ</td><td>(t) S ) 2 A</td><td>(t) S</td><td>]</td><td>.</td><td>(38)</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>♣4</td></tr><tr><td colspan="9">Now consider each term of equation (38). First, since</td></tr><tr><td colspan="4">E S [A  *  S A  *  S ] = E S [</td><td colspan="5">A  *  i A  *  j 1 i,j∈S ] =</td><td>m</td><td>A  *  i A  *</td></tr><tr><td/><td/><td/><td colspan="2">i,j∈S</td><td/><td/><td/><td>i,j=1</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>(t) S )A</td><td>(t) S</td><td>]</td></tr><tr><td colspan="4">in (38), which can be simplified as</td><td/><td/><td/><td/></tr><tr><td>♣ 2 ≤</td><td>m i,j=1</td><td colspan="2">A  *  i A</td><td>(t) j</td><td/><td/><td/></tr><tr><td/><td/><td/><td/><td/><td colspan="4">(t) S )A</td><td>(t) S</td><td>] = O( k 2</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>(t) S (I -Λ</td><td>(t) S ) 2 A</td><td>(t) S</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>t) S</td><td>)]</td></tr><tr><td colspan="9">-E S [(A  *  S -A (t) S (I -Λ</td><td>(t) S ))E x  *  S [x  *  S (ϑ</td><td>(R) S ) |S]A</td><td>(t) S</td><td>]</td></tr><tr><td colspan="3">+ E S [A</td><td colspan="2">(t) S E x  *  S [ϑ</td><td colspan="2">(R) S (ϑ</td><td colspan="2">(R) S ) |S]A</td><td>(t) S</td><td>].</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>(t) S )A</td><td>(t) S</td><td>)]</td></tr></table>

)E S [1 i,j∈S |S] , = O( t )O( k 2 m 2 )( A (t) 1 m×m A * + A (t) 1 m×m A (t) ), = O( t )O( k 2 m ) A * 2 . O( k 2 m ) t A * 2 . Similarly, ♥ in (37) is also bounded as ♠. Next, we consider ♦ in (37). In this case, lettingE x * S [ϑ |S] = N,where N ∈ R k×k is a matrix whose each entry N i,j ≤ |ϑ

Table content: <table><tr><td/><td/><td/><td/><td/><td/><td/><td/><td/><td colspan="2">M i,j A</td><td colspan="2">(t) i (A  *  j -(1 -λ</td><td>(t) j )A</td><td>(t) j</td><td>)E S [1 i,j∈S |S] ,</td></tr><tr><td/><td/><td/><td/><td/><td/><td colspan="3">= O( t )</td><td>m i,j=1</td><td colspan="2">A</td><td>(t) i (A  *  j -(1 -λ</td><td>(t) j )A</td><td>(t) j</td></tr><tr><td>Therefore,</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="13">♠ := E S [A )] = (R) (t) S M(A  *  S -(I -Λ (t) S )A (t) S S (ϑ (R)</td><td>(R) i ||ϑ</td><td>(R) j |.</td></tr><tr><td colspan="6">Further, by Claim 5, each element ϑ</td><td>(R) j</td><td colspan="6">is upper-bounded as</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td>|ϑ</td><td colspan="5">(R) j | ≤ O(t β ).</td></tr><tr><td colspan="9">with probability at least (1 -δ β ). Therefore, (t)</td><td/><td/><td/></tr><tr><td>♦ =</td><td>m i,j=1</td><td>N i,j A</td><td>(t) i A</td><td>(t) j</td><td colspan="8">E S [1 i,j∈S |S] = max i,j |ϑ</td><td>(R) i ||ϑ</td><td>(R) j |O( k 2 m 2 )</td><td>m i,j=1</td><td>A</td><td>(t) i A</td><td>(t) j</td></tr></table>

Final error in recovery of the factors by various techniques and the computation time taken per iteration (in seconds) corresponding to Fig.2across techniques. We report the coefficient estimate after the HT step (inArora et al. (

Table content: <table/>

10 -11 8.82 × 10 -11 9.70 × 10 -11 7.33 × 10 -11 X 1.14 × 10 -11 1.76 × 10 -11 3.58 × 10 -11 4.74 × 10 -11

Table content: <table><tr><td/><td>Factor and Timing</td><td>k = 10</td><td>k = 20</td><td>k = 50</td><td>k = 100</td></tr><tr><td>NOODL</td><td colspan="2">A 9.44 × Avg. Time/iteration 46.500 sec</td><td>53.303 sec</td><td>64.800 sec</td><td>96.195 sec</td></tr><tr><td/><td>A</td><td>0.013</td><td>0.031</td><td>0.137</td><td>N/A</td></tr><tr><td/><td>X HT</td><td>0.077</td><td>0.120</td><td>0.308</td><td>N/A</td></tr><tr><td>Arora15 (''biased'')</td><td>X Lasso Avg. Time/iteration (Accounting for one</td><td>0.006 39.390 sec</td><td>0.018 39.371 sec</td><td>0.097 39.434 sec</td><td>N/A 40.063 sec</td></tr><tr><td/><td>Lasso update)</td><td/><td/><td/><td/></tr><tr><td/><td>Avg. Time/iteration (Overall Lasso search)</td><td>389.368 sec</td><td>388.886 sec</td><td>389.566 sec</td><td>395.137 sec</td></tr><tr><td/><td>A</td><td>0.011</td><td>0.027</td><td>0.148</td><td>N/A</td></tr><tr><td/><td>X HT</td><td>0.078</td><td>0.122</td><td>0.371</td><td>N/A</td></tr><tr><td>Arora15 (''unbiased'')</td><td>X Lasso Avg. Time/iteration (Accounting for one</td><td>0.005 567.830 sec</td><td>0.015 597.543 sec</td><td>0.0921 592.081 sec</td><td>N/A 686.694 sec</td></tr><tr><td/><td>Lasso update)</td><td/><td/><td/><td/></tr><tr><td/><td>Avg. Time/iteration (Overall Lasso search)</td><td>917.809 sec</td><td>947.059 sec</td><td>942.212 sec</td><td>1041.767 sec</td></tr><tr><td/><td>A</td><td>0.009</td><td>0.015</td><td>0.021</td><td>0.037</td></tr><tr><td/><td>X Lasso</td><td>0.183</td><td>0.209</td><td>0.275</td><td>0.353</td></tr><tr><td>Mairal '09</td><td>Avg. Time/iteration (Accounting for one</td><td>39.110 sec</td><td>39.077 sec</td><td>39.163 sec</td><td>39.672 sec</td></tr><tr><td/><td>Lasso update)</td><td/><td/><td/><td/></tr><tr><td/><td>Avg. Time/iteration (Overall Lasso search)</td><td>388.978 sec</td><td>388.614 sec</td><td>389.512 sec</td><td>394.566 sec</td></tr></table>

Peer review (iclr2019):
 The paper deals with the problem of recovering an exact solution for both the dictionary and the activation coefficients. As other works, the solution is based on a proper initialization of the dictionary. The authors suggest using Aurora 2015 as a possible initialization. The contribution improves Arora 2015 in that it converges linearly and recovers both the dictionary and the coefficients with no bias. The main contribution is the use of a IHT-based strategy to update the coefficients, with a gradient-based update for the dictionary (NOODL algorithm). The authors show that, combined with a proper initialization, this has exact recovery guaranties. Interestingly, their experiments show that NOODL converges linearly in number of iterations, while Arora gets stuck after some iterations. I think the paper is relevant and proposes an interesting contribution. The paper is well written and the key elements are in the body. However, there is a lot of important material in the Appendix, which I think may be relevant to the readers. It would be nice to have some more intuitive explanations at least of Theorem 1. Also, it is clear in the experiments the superiority with respect to Arora in terms of iterations (and error), but what about computational time?