Algorithmic Regularization in Over-parameterized Matrix Sensing and Neural Networks with Quadratic Activations
Yuanzhi Li Tengyu Ma Hongyang Zhang
cs.LG cs.DS math.OC stat.ML








Algorithmic Regularization in Over-parameterized Matrix Sensing and Neural Networks with Quadratic Activations
    Wei Zhang
     December 30, 2023
==============================================================================================================




    We show that the gradient descent algorithm provides an implicit regularization effect in the learning of over-parameterized matrix factorization models and one-hidden-layer neural networks with quadratic activations. 
    
Concretely, we show that given Õ(dr^2) random linear measurements of a rank r positive semidefinite matrix X^⋆, we can recover X^⋆ by parameterizing it by UU^⊤ with U∈ℝ^d× d and minimizing the squared loss, even if r ≪ d. We prove that starting from a small initialization, gradient descent recovers X^⋆ in Õ(√(r)) iterations approximately. The results solve the conjecture of Gunasekar et al. <cit.> under the restricted isometry property. 
    
The technique can be applied to analyzing neural networks with one-hidden-layer quadratic activations with some technical modifications. 





§ INTRODUCTION


Over-parameterized models are crucial in deep learning, but their workings are far from understood. Over-parameterization — the technique of using more parameters than statistically necessary — apparently improves the training: theoretical and empirical results have suggested that it can enhance the geometric properties of the optimization landscape in simplified settings <cit.> and thus make it easier to train over-parameterized models.

On the other hand,  over-parameterization often doesn't hurt the test performance, even if the number of parameters is much larger than the number of examples. Large neural networks used in practice have enough expressiveness to fit any labels of the training datasets <cit.>. The training objective function may have multiple global minima with almost zero training error, some of which generalize better than the others <cit.>. However, local improvement algorithms such as stochastic gradient descent, starting with proper initialization, may prefer some generalizable local minima to the others and thus provide an implicit effect of regularization <cit.>. 
Such regularization seems to depend on the algorithmic choice, the initialization scheme, and certain intrinsic properties of the data.

The phenomenon and intuition above can be theoretically fleshed out in the context of linear models <cit.>, whereas less is known for non-linear models whose training objectives are usually non-convex. The very important work of Gunasekar et al. <cit.> initiates the study of low-rank matrix factorization models with over-parameterization and conjectures that gradient descent prefers small trace norm solution in over-parameterized models with thorough empirical evidences.

This paper resolves the conjecture for the matrix sensing problem — recovering a low-rank matrix from linear measurements — under the restricted isometry property (RIP). We show that with a full-rank factorized parameterization, gradient descent on the squared loss with finite step size, starting with a small initialization, converges to the true low-rank matrix (which is also the minimum trace norm solution.) One advantage of the over-parameterized approach is that without knowing/guessing the correct rank,  the algorithms can automatically pick up the minimum rank or trace norm solution that fits the data. 

The analysis can be extended to learning one-hidden-layer neural networks with quadratic activations. We hope such theoretical analysis of algorithmic regularization in the non-convex setting may shed light on other more complicated models where over-parameterization is crucial (if not necessary) for efficient training. 







 §.§ Setup and Main Results

Let  be an unknown rank-r symmetric positive semidefinite (PSD) matrix in ^d× d that we aim to recover. Let _1, ⋯, _m ∈ℝ^d × d be m given symmetric measurement matrices.[Given that the matrix  is symmetric, we can assume that A_i's are symmetric without loss of generality: Because A_i,  = 1/2(A_i+A_i^⊤), for any symmetric matrix , we can always replace A_i by 1/2(A_i+A_i^⊤).] We assume that the label vector y∈^m is generated by linear measurements 

    y_i = ⟨_i, ⟩.

Here A, B = (A^⊤ B) denotes the inner product of two matrices. Our goal is to recover the matrix . [Our analysis can naturally handle a small amount of Gaussian noise in the label vector y, but for simplicity we only work with the noiseless case. ]

Without loss of generality, we assume that  has spectral norm 1.  Let σ_r() denote the r-th singular value of a matrix , and let κ = 1/σ_r() be the condition number of . We focus on the regime where r≪ d and m ≈ d·(r log d) ≪ d^2.

Let U∈ℝ^d× d be a matrix variable. We consider the following mean squared loss objective function with over-parameterization:

    min_U ∈ℝ^d × d f() = 1/2m∑_i = 1^m (y_i- ⟨_i, ^⊤⟩)^2

Since the label is generated by y_i = A_i,, any matrix U satisfying UU^⊤ = is a local minimum of f with zero training error. These are the ideal local minima that we are shooting for.  However, because the number of parameters d^2 is much larger than the number of observation m, there exist other choices of U satisfying f(U) = 0 but UU^⊤≠. 

A priori, such over-parameterization will cause over-fitting. However, we will show that the following gradient descent algorithm with small initialization converges to a desired local minimum, instead of other non-generalizable local minima:

    U_0 = α B, 
        U_t+1 = U_t - η∇ f(U_t)

The following theorem assumes the measurements matrices A_1,…, A_m satisfy restricted isometry property (RIP),  which is formally defined in Section <ref>. Casual readers may simply assume that the entries of A_i's are drawn i.i.d from standard normal distribution[Or equivalently, as discussed in the previous footnote, causal readers may assume A_i = 1/2(Q_i+Q_i^⊤) where Q_i is from standard normal distribution. Such symmetrization doesn't change the model since Q_i,  = 1/2(Q_i+Q_i^⊤), ] and the number of observations m ≲ dr^2log^3 d: it's known <cit.> that in this case A_1,…, A_m meet the requirement of the following theorem, that is, they satisfy (4r,δ)-RIP with δ≲ 1/(√(r)log d) with high probability.
[Technically, to get such RIP parameters that depends on r, one would need to slightly modify the proof of <cit.> at the end to get the dependency of m on δ. ] 


	Let c be a sufficiently small absolute constant. Assume that the set of measurement matrices (A_1,…, A_m) satisfies (4r,δ)-restricted isometry property (defined in Section <ref> formally) with δ≤ c/(κ^3√(r)log^2 d). Suppose the initialization and learning rate satisfy 0< α≤ cmin{δ√(r)κ, 1/d} and η≤ cδ. Then for every (κlog (d/α))/η≲ T ≲  1/(η√(dκα)), we have
	
    U_TU_T^⊤-_F^2 ≲α√(d) / κ^2.



Note that the recovery error U_TU_T^⊤-_F^2 can be viewed as the test error (defined in Equation (<ref>) formally) — it's the expectation of the test error on a fresh measurement A_j drawn from the standard normal distribution. The theorem above shows that gradient descent can provide an algorithmic regularization so that the generalization error depends on the size of the initialization α, instead of the number of parameters. Because the convergence is not very sensitive to the initialization, we can choose small enough α (e.g., 1/d^5) to get approximately zero generalization error. Moreover, when α is small, gradient descent can run for a long period of time without overfitting the data. We show in Section <ref> that empirically indeed the generalization error depends on the size of the initialization and gradient descent is indeed very stable. 

The analysis also applies to stochastic gradient descent, as long as each batch of the measurement matrices satisfies RIP.[Smaller batch size should also work when the learning rate is sufficiently small, although its analysis seems to require more involved techniques and is left for future work.] We also remark that our theory suggests that early stopping for this problem is not necessary when the initialization is small enough — the generalization error bounds apply until 1/(η√(dκα)) iterations. We corroborate this with empirical results in Section <ref>. 

We remark that we achieve a good iteration complexity bound 1/η≈ 1/δ≈√(r) for the gradient descent algorithm, which was not known in previous work even for low-rank parameterization, nor for the case with infinite samples (which is the PCA problem).   Part of the technical challenges is to allow finite step size η and inverse-poly initialization α (instead of exponentially small initialization). The dependency of δ on κ and r in the theorem is possibly not tight. We conjecture that δ only needs to be smaller than an absolute constant, which is left for future work. 



  
Insights of the analysis: Interestingly, our analysis ideas seem to be different from other previous work in a conceptual sense. The analysis of the logistic regression case <cit.> relies on that the iterate eventually moves to infinity. The folklore analysis of the algorithmic regularization of SGD for least squares and the analysis in <cit.> for the matrix regression with commutable measurements both follow the two-step plan: a) the iterates always stays on a low-rank manifold that only depends on the inputs (the measurement matrices) but not on the label vector y; b) generalization follows from the low complexity of the low-rank manifold. Such input-dependent but label-independent manifold doesn't seem to exist in the setting when A_i's are random. 

Instead, we show that the iterates stay in the set of matrices with approximate rank smaller or equal to the minimal possible rank that can fit the data, which is a set that depends on the labels y but not on the inputs A_i's.  
We implicitly exploit the fact that gradient descent on the population risk with small initialization only searches through the space of solutions with a lower rank than that of the true matrix .  The population risk is close to the empirical risk on matrices with rank smaller than or equal to the true rank.  Hence, we can expect the learning dynamics of the empirical risk  to be similar to that of the population risk, and therefore the iterates of GD on the empirical risk remain approximately low-rank as well. Generalization then follows straightforwardly from the low-rankness of the iterates. See Section <ref> for more high-level discussions.

We note that the factorized parameterization also plays an important role here. The intuition above would still apply if we replace UU^⊤ with a single variable X and run gradient descent in the space of X with small enough initialization. However, it will converge to a solution that doesn't generalize.  The discrepancy comes from another crucial property of the factorized parameterization: it provides certain denoising effect that encourages the empirical gradient to have a smaller eigenvalue tail. This ensures the eigenvalues tails of the iterates to grow sufficiently slowly. This point will be more precise in Section <ref> once we apply the RIP property. In section <ref>, we also empirically demonstrate that GD in the original space of  with projection to the PSD cone doesn't provide as good generalization performance as GD in the factorized space.  

Finally, we remark that the cases with rank r> 1 are technically much more challenging than the rank-1 case. For the rank-1 case, we show that the spectrum of U_t remains small in a fixed rank-(d-1) subspace, which is exactly the complement of the column span of . Hence the iterates are approximately rank one. By contrast, for the rank-r case, a direct extension of this proof strategy only gives a much weaker result compared to Theorem <ref>. Instead, we identify an adaptive  rank-(d-r) subspace in which U_t remains small. Clearly, the best choice of this adaptive subspace is the subspace of the least (d-r) left singular vectors of U_t. However, we use a more convenient surrogate. We refer the reader to Section <ref> for detailed descriptions.





 §.§ Extensions to Neural Networks with Quadratic Activations


Our results can be applied to learning one-hidden-layer neural networks with quadratic activations. We setup the notations and state results below and defer the details to Section <ref>. 

Let x∈^d be the input and U^⋆∈^d× r be the first layer weight matrix. We assume that the weight on the second layer is simply the all one's vector 1∈^r.  Formally, the label y is assumed to be generated by 

    y = 1^⊤ q(U^⋆^⊤ x)

where q(·) is the element-wise quadratic function. For simplicity, we assume that x comes from standard normal distribution 𝒩(0,_d× d). It's not hard to see that the representational power of the hypothesis class with r=d is the same as those with r > d. Thus we only focus on the case when r ≤ d.
For the purpose of this paper, the most interesting regime is the scenario when r≪ d. 

We use an over-parameterized model with a variable U∈^d× d. The prediction ŷ is parameterized by ŷ = 1^⊤ q(U^⊤ x), 
and we use the mean squared error (y-ŷ)^2 as the loss function. We use a variant of stochastic gradient descent (or gradient descent) on the mean squared loss. 

The following theorem shows that the learned model will generalize with Õ(dr^5κ^6) examples, despite that the number of parameters d^2 can be much larger than the number of samples (when d ≫ r or r is considered as a constant).[The dependency on r here is likely to be loose. Although we note that this is the first bound of this kind for this problem that shows over-parameterized models can still generalize. ] We will start with an initialization U_0 in the same way as in equation (<ref>), and denote U_1,…, U_T as the iterates. Let κ be the condition number of ^⊤. 




Given Õ(dr^5κ^6) examples, 	a variant of gradient descent (Algorithm <ref> in Section <ref>) with initialization α≲min{1/d, 1/(r^2κ^4log^2 d)} and learning rate η≲1/κ^3 r^1.5log^2 d returns a solution with generalization error at most O(dκα) at any iteration t such that  (κlog (d/α))/η≲ t ≲ 1/(η√(dκα)). 


The same analysis also applies to stochastic gradient descent as long as the batch size is at least ≳ dr^5κ^6. The analysis exploits the connection (pointed out by  <cit.>) between neural networks with quadratic activations and matrix sensing with rank-1 measurements <cit.>: one can view xx^⊤ as the measurement matrix in matrix sensing. However, these measurements don't satisfy the RIP property.  We will modify the learning algorithm slightly to cope with it. See Section <ref> for details. 


Organization: 
The rest of this paper is organized as follows:
In Section <ref>, we define notations and present a review of the restricted isometry property. 
In Section <ref>, we lay out the key theoretical insights towards proving Theorem <ref> and give the analysis for the rank-1 case as a warm-up. 
In Section <ref>, we outline the main steps for proving Theorem <ref> and Section <ref> completes the proofs of these steps. Section <ref> and Section <ref> give the proof of Theorem <ref>. Section <ref> contains numeric simulations. 
Finally, Section <ref> provide the proofs of concentration properties we have used.


  
Notations:
Let _U denotes the projection to the column span of U, and let  denotes the identity matrix.  Let U^+ denote the Moore-Penrose pseudo-inverse of the matrix U. Let · denotes the Euclidean norm of a vector and spectral norm of a matrix.  Let ·_F denote the Frobenius norm of a matrix. 
Suppose A∈^m× n, then σ_max(A) denote its largest singular value and σ_min(A) denotes its min{m,n}-th largest singular value. Alternatively, we have σ_min(A) = min_x:x=1Ax. Let A, B = (A^⊤ B) denote the inner product of two matrices. We use sin(A,B) to denote the sine of the principal angles between the columns spaces of A and B. 

Unless explicitly stated otherwise, O(·)-notation hides absolute multiplicative constants.
Concretely, every occurrence of O(x) is a placeholder for some function f(x) that satisfies ∀ x∈,  |f(x)|≤ C|x| for some absolute constant C>0. Similarly, a≲ b means that there exists an absolute constant C> 0 such that a ≲ Cb. We use the notation (n) as an abbreviation for n^O(1).  



§ PRELIMINARIES AND RELATED WORK

Recall that we assume  is rank-r and positive semidefinite. Let = ^⊤ be the eigen-decomposition of , where ∈ℝ^d × r is an orthonormal matrix and ∈ℝ^r × r is a diagonal matrix. The assumptions that =1 and σ_r() = 1/κ translate to that ∀ i ∈ [r],  1/κ≤Σ^⋆_ii≤ 1. 
Under the above notations, we see that the target solution for the variable U is equal to U = Σ^⋆^1/2R where R can be arbitrary orthonormal matrix. 
For convenience, we define the matrix M_t as

    M_t = 1/m∑_i = 1^m_i, _t_t^⊤ - _i

Then the update rule can be rewritten as 

    U_t+1 = (- η M_t)U_t

where  is the identity matrix. One of the key challenges is to understand how the matrix - η M_t transforms U_t, so that U_0 converges the target solution Σ^⋆^1/2R quickly.

Suppose that A_1,…, A_m are drawn from Gaussian distribution and  optimistically suppose that they are independent with U_t. Then, we have that M_t ≈ U_tU_t^⊤ - , since the expectation of M_t with respect to the randomness of A_i's is equal to U_tU_t^⊤ -. However, they are two fundamental issues with this wishful thinking: a) obviously U_t depends on A_i's heavily for t> 1, since in every update step A_i's are used; b) even if A_i's are independently with U_t, there are not enough A_i's to guarantee M_t concentrates around its mean U_tU_t^⊤ - in Euclidean norm. To have such concentration, we need m > d^2, whereas we only have m = d ×(r log d) samples. 




  
Restricted isometry propety: The restricted isometry property (RIP) allows us to partially circumvent both the technical issues a) and b) above. It says that using the set of linear measurement matrices A_1,…, A_m, we can preserve the Frobenius norm of any rank-r matrices approximately.


(Restricted isometry property <cit.>)  A set of linear measurement matrices  A_1,…, A_m in ℝ^d× d satisfies (r,δ)-restricted isometry property (RIP) if for any d× d matrix X with rank at most r, we have
	
    (1-δ)X_F^2 ≤1/m∑_i = 1^m A_i, X^2 ≤ (1+δ)X_F^2


The crucial consequence of RIP that we exploit in this paper is the meta statement as follows: 

    


We will state several lemmas below that reflect the principle above. The following lemma says that ⟨ℳ(X), Y⟩ behaves like ⟨ X, Y⟩ for low-rank matrices X and Y. 

<cit.>
	Let { A_i }_i=1^m be a family of matrices in ^d × d
	that satisfy (r, δ)-restricted isometry property.
	Then for any matrices X, Y ∈^d × d with rank at most r,
	we have:
	
    1/m∑_i=1^m A_iXA_iY    - XY≤δXY



The following lemma says that ℳ(X) behaves like X when multiplied by a matrix R with small operator norm.  

	Let {A_i}_i=1^m be a family of matrices in ^d × d that
	satisfy (r, δ)-restricted isometry property.
	Then for any matrix X ∈ℝ^d × d of rank at most r,
	and any matrix R ∈ℝ^d × d', where d' can be any positive integer,
	we have:
	
    1/m∑_i = 1^m ⟨_i , ⟩_i   - ≤δX·R.


Lemma <ref> is proved in Section <ref>[We suspect that Lemma <ref> is already known, however we haven't been able to find a reference.]. We can also extend Lemma <ref> to the cases when X has a higher rank (see Lemma <ref> and Lemma <ref>). The bounds are not as strong as above (which is inevitable because we only have m measurements), but are useful when X itself is relatively small. 





 §.§ Related Work




  
Generalization theory beyond uniform convergence:  
This work builds upon the remarkable work of Gunasekar et al. <cit.>, which raises the conjecture of the implicit regularization in matrix factorization models and provides theoretical evidence for the simplified setting where the measurements matrices are commutable.  Implicit regularization of gradient descent is studied in the logistic regression setting by Soudry et al. <cit.>.

Recently, the work of Hardt et al. <cit.> studies the implicit regularization provided by stochastic gradient descent through uniform stability <cit.>. Since the analysis therein is independent of the training labels and therefore it may give pessimistic bounds <cit.>. Brutzkus et al. <cit.> use a compression bound to show network-size independent generalization bounds of one-hidden-layer neural networks on linearly separable data. 

Bartlett et al. <cit.>, Neyshabur et al. <cit.>, and Cisse et al. <cit.> recently prove spectrally-normalized margin-based generalization bounds for neural networks. Dziugaite and Roy <cit.> provide non-vacuous generalization bounds for neural networks from PCA-Bayes bounds. As pointed out by Bartlett et al. <cit.>, it's still unclear why SGD without explicit regularization can return a large margin solution. This paper makes progress on explaining the regularization power of gradient descent, though on much simpler non-linear models. 




  
Matrix factorization problems: 
Early works on matrix sensing and matrix factorization problems use convex relaxation  (nuclear norm minimization) approaches and obtain tight sample complexity bounds <cit.>. Tu et al. <cit.> and Zheng and Lafferty <cit.> analyze the convergence of non-convex optimization algorithms from spectral initialization. The recent work of Ge et al. <cit.> and Bhojanapalli et al. <cit.>
shows that the non-convex objectives on matrix completion and matrix sensing with low-rank parameterization don't have any spurious local minima, and stochastic gradient descent algorithm on them converges to the global minimum. 
Such a phenomenon was already known for the PCA problem and recently shown for phase retrieval, robust PCA,
and random tensor decomposition as well (e.g., see <cit.> and references therein). 
Soltanolkotabi et al. <cit.> analyzes the optimization landscape of over-parameterized one-hidden-layer neural networks with quadratic activations. 
Empirically, Jose et al. <cit.> show that factorized parameterizations of recurrent neural networks provide additional regularization effect. 






§ PROOF OVERVIEW AND RANK-1 CASE



In this section, we demonstrate the key ideas of our proofs and give an analysis of the rank-1 case as a warm-up. 
The main intuition is that the iterate U_t stays approximately low-rank in the sense that:

	
  (a) The (r+1)-th singular value σ_r+1(U_t)  for any t≥ 0;
	
  (b) The top r singular vectors and singular values of U_tU_t^⊤ converge to those of  in logarithmic number of iterations. 

Propositions (a) and (b) can be clearly seen when the number of observations m approaches infinity and A_1,…, A_m are Gaussian measurements. Let's define the population risk f̅ as 

    f̅(U_t) = _(A_i)_kℓ∼ N(0,1)[f(U_t)] = U_tU_t^⊤- _F^2

In this case, the matrix M_t (defined in  (<ref>)) corresponds to U_tU_t^⊤ -, and therefore the update rule for U_t can be simply rewritten as 

    U_t+1    = U_t - η∇f̅(U_t) = U_t - η(U_tU_t^⊤ - )U_t 
        = U_t( - η U_t^⊤ U_t)+ η U_t

Observe that the term η U_t encourages the column span of U_t+1 to move towards the column span of , which causes the phenomenon in Proposition (b). On the other hand, the term U_t( - η U_tU_t^⊤) is performing a contraction of all the singular values of U_t, and therefore encourages them to remain small. 
As a result, U_t decreases in those directions that are orthogonal to the span of , because there is no positive force to push up those directions.

So far, we have described intuitively that the iterates of GD on the population risk remains approximately low-rank. 
Recall that the difficulty was that the empirical risk f doesn't uniformly concentrate well around the population risk f̅.[Namely, we don't have uniform convergence results in the sense that |f(U)-f̅(U)| is small for all matrices U.  (For examples, for many matrices we can have f(U) = 0 but f̅(U) ≫ 0 because we have more variables than observations.)] However, the uniform convergence can occur, at least to some extent, in the restricted set of approximately low-rank matrices! In other words, since the gradient descent algorithm only searches a limited part of the whole space, we only require restricted uniform convergence theorem such as restricted isometry property. 
Motivated by the observations above, a natural meta proof plan is that:

	
  1. The trajectory of the iterates of GD on the population risk stays in the set of approximately low-rank matrices.
	
  2. The trajectory of the empirical risk behaves similarly to that of the population risk in the set of approximately low-rank matrices.

It turns out that implementing the plan above quantitatively can be technically challenging: the distance from the iterate to the set of low-rank matrices can accumulate linearly in the number of steps. Therefore we have to augment the plan with a strong result about the rate of convergence: 

	
  3. The iterates converge to the true solution  fast enough before its effective rank increases.

For the rest of this section, we demonstrate a short proof of the rank-1 case to implement the intuitions described above.
We note that the results of the rank-r case in Section <ref> is significantly tighter than the results presented in this section.
The analysis involves more delicate techniques to control the growth of the top r eigenvectors, and requires a much stronger convergence analysis.




 §.§ Warm-up: Rank-1 Case


In this subsection, we assume that = ^⊤ for ∈^d× 1 and that = 1. 
We decompose the iterates U_t into the subspace of  and its complement:

    U_t     = _U_t +   U_t
        :=  r_t^⊤ + E_t

where we denote by r_t:= U_t^⊤ and  E_t :=  U_t.
[Observe that we have restricted the column subspace of the signal term R_t =  r_t^⊤, so that R_tR_t^⊤ is always a multiple of .
In section <ref>, we will introduce an adaptive subspace instead to decompose U_t into the signal and the error terms.]

In light of the meta proof plan discussed above,  we will show that the spectral norm and Frobenius norm of the “error term” E_t remains small throughout the iterations, whereas the “signal” term u^⋆ r_t^⊤ grows exponentially fast (in the sense that the norm of r_t grows to 1.) Note that any solution with r_t=1 and E_t =0 will give exact recovery, and for the purpose of this section we will show that r_t will converges approximately to 1 and E_t stays small. 






Under the representation (<ref>), from the original update rule (<ref>), we derive the update for E_t:
	
    E_t+1   = ( - _) · ( - η M_t) U_t 
       = E_t - η· M_t U_t


Throughout this section, we assume that r=1 and the set of measurement matrices (A_1,…, A_m) satisfies (4,δ)-RIP with δ≤ c where c is a sufficiently small absolute constant (e.g., c=0.01 suffices). 
 In the setting of this subsection, suppose α≤δ√(1/dlog1/δ) and η≲  cδ^2log^-1 (1/δα).
			Then after T = Θ(log1/αδ / η) iterations, we have:
	
    U_T U_T^⊤ - ≲δlog1/δ



As we already mentioned, Theorem <ref> is weaker than Theorem <ref> even for the case with r=1. In Theorem <ref> (or Theorem <ref>), the final error depends linearly on the initialization, whereas the error here depends on the RIP parameter. Improving Theorem <ref> would involve finer inductions, and we refer the readers to Section <ref> for the stronger results.

The following lemma gives the growth rate of E_t in spectral norm and Euclidean norm, in a typical situation when E_t and r_t are bounded above in Euclidean norm.

	In the setting of Theorem <ref>. Suppose that E_t≤ 1/2 and r_t^2 ≤ 3/2.
	Then E_t+1 can be bounded by
	
    E_t+1^2    ≤E_t^2 + 2ηδE_t r_t + 1.5 ηδE_t^2 + 9η^2 .
    E_t+1   ≤ (1 + 2ηδ) E_t + 2ηδr_t.



A recurring technique in this section, as alluded before, is to establish the approximation

    U_t+1 = U_t - η M_t U_t ≈ U_t - η (U_tU_t-)U_t

As we discussed in Section <ref>, if U_tU_t- is low-rank, then the approximation above can be established by Lemma <ref> and Lemma <ref>. However, U_tU_t- is only approximately low-rank, and we therefore we will decompose it into

    U_t U_t^⊤ -  = (U_t U_t^⊤ -  - E_t E_t^⊤)_≤ 4 + E_tE_t^⊤_

Note that U_t U_t^⊤ -  - E_t E_t^⊤ = r_t^2 u^⋆u^⋆^⊤ + u^⋆ r_t^⊤ E_t^⊤ + E_tr_tu^⋆^⊤ - has rank at most 4, and therefore we can apply Lemma <ref> and Lemma <ref>. For the term E_tE_t^⊤, we can afford to use other looser bounds (Lemma <ref> and <ref>) because E_t itself is small. 

Using the update rule (<ref>) for E_t, we have that 

    E_t+1^2
    	   = E_t^2 - 2η·E_t M_t U_t
    	+ η^2  M_t U_t^2

	When η is sufficiently small and M_t, U_t are bounded from above, the third term on the RHS is negligible compared to the second term. Therefore, we focus on the second term first.

    E_t M_t U_t
       = 1/m∑_i=1^m A_iU_t U_t^⊤ -  A_i E_t U_t^⊤

	where in the last line we rearrange the terms and use the fact that  is symmetric.
	Now we use Lemma <ref> to show that equation (<ref>) is close to 
	U_t U_t^⊤ -  E_t U_t^⊤, which is its expectation w.r.t the randomness of A_i's if A_i's were chosen from spherical Gaussian distribution.
	If U_t U_t^⊤ - was a rank-1 matrix, then this would follow from Lemma <ref> directly. However, U_tU_t^⊤ is approximately low-rank. Thus, we decompose it into a low-rank part and an error part with small trace norm as in equation (<ref>). Since U_t U_t^⊤ -  - E_t E_t^⊤ has rank at most 4,  we can apply Lemma <ref> to control the effect of A_i's, 
	
    1/m∑_i=1^m A_i(U_t U_t^⊤ -  - E_t E_t^⊤)A_iE_t U_t^⊤
       ≥U_t U_t^⊤ -  - E_t E_t^⊤, E_tU_t^⊤ - 	δU_t U_t^⊤ -  - E_t E_t^⊤_F E_t U_t^⊤
       ≥U_t U_t^⊤ -  - E_t E_t^⊤, E_tU_t^⊤ - 	1.5δE_t U_t^⊤

	where the last inequality uses that 	U_t U_t^⊤ -  - E_t E_t^⊤^2 = (1 - r_t^2)^2 + 2E_t r_t^2 ≤ 1 + E_t^2 r_t^2 ≤ 11/8.
	
	For the E_tE_t^⊤ term in the decomposition (<ref>), we have that
 	
    1/m∑_i=1^m A_iE_t E_t^⊤A_iE_t U_t^⊤   ≥1/m∑_i=1^m A_iE_t E_t^⊤A_iE_t r_t u^⋆^⊤clearly A_iE_t E_t^⊤^2 ≥ 0
       ≥E_t E_t^⊤E_t r_t u^⋆^⊤ -δE_t E_t^⊤_⋆E_t r_t u^⋆^⊤by Lemma <ref>
       ≥E_t E_t^⊤E_t r_t u^⋆^⊤ - 0.5δE_t r_t.

  Combining equation (<ref>),  (<ref>) and (<ref>), we conclude that 
 	
    E_t M_t U_t   ≥U_t U_t^⊤ - E_t U_t^⊤ - E_t E_t^⊤E_t E_t^⊤ - 2δE_t r_t - 1.5 δE_t^2,

	where we have used that E_t U_t^⊤≤E_t r_t + E_t E_t^⊤ = E_t r_t + E_t^2.
	Note that u^⋆^⊤ E_t = 0, which implies that E_t = 0 and U_t^⊤ E_t = E_t^⊤ E_t. Therefore, 
	
    U_t U_t^⊤ - E_t U_t^⊤     = U_t U_t^⊤E_t U_t^⊤
        = U_t^⊤U_T^⊤ E_t U_t^⊤ = U_t^⊤E_T^⊤ E_t U_t^⊤
        = E_tU_t^⊤ E_t U_t^⊤ = E_t U_t^⊤_F^2 ≥E_t E_t^⊤^2 ≥ 0,

because E_t U_t^⊤^2 = E_t E_t^⊤^2 + E_t r_t u^⋆^⊤^2.
We can also control the third term in RHS of equation (<ref>) by 
η^2  M_t U_t^2 ≤ 9η^2.
Since the bound here is less important (because one can always choose small enough η to make this term dominated by the first order term), we left the details to the reader. Combining the equation above with (<ref>) and  (<ref>), we conclude the proof of equation (<ref>). 

	Towards bounding the spectral norm of E_t+1, we use a similar technique to control the difference between M_t U_t and (U_t U_t^⊤ - ) U_t in spectral norm. By decomposing U_t U_t^⊤ - as in equation (<ref>) and applying Lemma <ref> 
	and Lemma <ref> respectively, we obtain that 
	
    M_t U_t -  (U_t U_t^⊤ - ) U_t   ≤ 
    		4δ(U_t U_t^⊤ -  - E_t E_t^⊤ + E_tE_t^⊤_⋆) U_t^⊤
       ≤ 8 δU_t≤ 8δ (r_t + E_t) by the assumptions that E_t_F≤ 1/2, r_t≤ 3/2

	Observing that (U_t U_t^⊤ - ) U_t = E_tU_t^⊤ U_t. Plugging the equation above into equation (<ref>), we conclude that
	
    E_t+1   ≤E_t ( - η U_t^⊤U_t) + 2ηδ (r_t + E_t) 
       ≤E_t + 2ηδ (r_t + E_t)


The next Proposition shows that the signal term grows very fast, when the signal itself is not close to norm 1 and the error term E_t is small. 

In the same setting of Proposition <ref>, we have, 
	
    r_t+1 - (1 + η(1 - r_t^2)) r_t≤ηE_t^2 r_t +  2ηδ(E_t + r_t).




By the update rule (<ref>), we have that 
	
    r_t+1   = U_t+1^⊤ = U_t^⊤ ( -η M_t^⊤) 
       = r_t - η U_t^⊤ M_t^⊤.

	By decomposing M_t = 1/m∑_i=1^m A_iU_t U_t^⊤ -  A_i as in equation (<ref>),
	and then Lemma <ref> and <ref>, we obtain that
	
    r_t+1 - (r_t - η U_t^⊤ (U_t U_t^⊤ - ) )≤δ(U_t U_t^⊤ -  - E_t E_t^⊤ + E_t E_t^⊤) U_t.

	Observe that U_t^⊤ (U_t U_t^⊤ - )  = U_t^⊤ U_t r_t - r_t = (r_t r_t^⊤ + E_t^⊤ E_t) r_t - r_t
	= (r_t^2 - 1) r_t - E_t^⊤ E_t r_t.
	Also note that
	U_t U_t^⊤ -  - E_t E_t^⊤^2 ≤ 11/8 and E_t^2 ≤ 1/4, we have that
	
    r_t+1 - (1 + η(1 - r_t^2)) r_t≤ηE_t^⊤E_t r_t +
    		2ηδU_t

	Since U_t≤r_t + E_t,
	we obtain the conclusion.

The following proposition shows that r_t converges to 1 approximately and E_t remains small by inductively using the two propositions above. 

In the setting of Theorem <ref>, after T = Θ(log(1/αδ)) / η) iterations,
	
    r_T = 1± O(δ) 
       E_T^2 ≲δ^2log(1/δ)




	We will analyze the dynamics of r_t and E_t in two stages.
	The first stage consists of all the steps such that r_t≤ 1/2,
	and the second stage consists of the rest of steps.  We will show that
	
		
  a) Stage 1 has at most O(log(1/α)/η) steps.
		Throughout Stage 1, we have
		
    E_t≤ 9δ
       r_t+1≥ (1+η/3)r_t

		
  b) In Stage 2, we have that 
		
    E_t^2_F ≲δ^2log1/δ
       r_t≤ 1+O(δ).

		And after at most O(log(1/δ)/η) steps in Stage 2, we have r_t≥ 1-O(δ). 
	

	We use induction with Proposition <ref> and  <ref>. For t=0, we have that E_0 = r_0 = α because U_0 = α B, where B is an orthonormal matrix. Suppose equation (<ref>) is true for some t, then we can prove equation (<ref>) holds by Proposition <ref>:
	
    r_t+1   ≥   (1 + η (1 - r_t^2 - 2δ - O(δ^2))) r_t
       ≥ (1+η/3)r_tby δ≤ 0.01 and r_t≤ 1/2

	Suppose equation (<ref>) holds, we can prove equation (<ref>) holds by Proposition <ref>.
	We first observe that t ≤ O(log1/α / η),
	since r_t grows by a rate of 1 + η/3.
	Denote by λ = 1 + 2ηδ. We have
	
    E_t+1≤λE_t + 2ηδr_t⇒E_t+1/λ^t+1≤E_t/λ^t + 2ηδ×(r_t/λ^t+1) 
    ⇒   E_t+1≤λ^t+1×α +
    		2ηδ×∑_i=0^t λ^i r_t/(1+η/3)^iby r_t+1≥r_t (1 + η/3).≤ 9δ,

	where the last inequality uses that
	
    λ^t+1α≤α×exp(2ηδ× O(log(1/α)) / η) = α^1 - O(δ)≤ o(δ), 
       ∑_i=1^t λ^i r_t/(1+η/3)^i≤1 + η/3/2 (η/3 - 2ηδ)


For claim b), we first apply the bound obtained in claim a) on E_t and the fact that U_t≤ 3/2
(this follows trivially from our induction, so we omit the details).
We have already proved that E_t≤ 9δ in the first stage.
For the second stage, as long as the number of iterations is less than O(log(1/δ) / η),
we can still obtain via Proposition <ref> that:

    E_t≤ 9δ + 4ηδ× O(log(1/δ / η))
    = O(δlog(1/δ)) since E_t≤ 1/2 and r_t≤ 3/2

To summarize, we bound E_t as follows:
	
    E_t^2    ≤E_0^2 + 9η^2 t + O(η t δ^2 log(1/δ)) ≲δ^2log(1/δ),

	because E_0^2 = α^2 d ≤δ^2 log1/δ,
	η t ≤ O(1),
	and	η^2 t ≤ O(η) ≤ O(δ^2 log(1/δ)).
		For the bound on r_t, we note that since E_t≤ O(δlog(1/δ)), we can simplify Proposition <ref> by 

    r_t+1 -  (1 + η (1 - r_t^2))r_t≤ 2ηδr_t + O(ηδ^2log^2(1/δ)).

	The proof that r_t = 1 ± O(δ) after at most O(log(1/δ) / η) steps follows similarly by induction.
	The details are left to the readers.

Theorem <ref> follows from Proposition <ref> straightforwardly. 

Using the conclusions of Proposition <ref>, we have
	
    U_T U_T^⊤ - ^2
    	   = (1 - r_T^2)^2 + 2 E_T r_T^2 + E_T E_T^⊤^2 
       ≤ (1 - r_T^2)^2 + 2 E_T^2 r_T^2 + E_T^4 
       ≤ O(δ^2) + O(δ^2log^21/δ) +
    	O(δ^4 log^2(1/δ))
    	= O(δ^2log^21/δ)






§ PROOF OUTLINE OF RANK-R CASE


In this section we outline the proof of Theorem <ref>. The proof is significantly more sophisticated than the rank-1 case (Theorem  <ref>), because the top r eigenvalues of the iterates grow at different speed.
Hence we need to align the signal and error term in the right way so that the signal term grows monotonically.
Concretely, we will decompose the iterates into a signal and an error term according to a dynamic subspace, as we outline below. Moreover, the generalization error analysis here is also tighter than Theorem <ref>.  
We first state a slightly stronger version of Theorem <ref>:


	There exists a sufficiently small absolute constant c > 0 such that the following is true.
	For every α∈ (0, c/d), assume that the set of measurement matrices (A_1,…, A_m) satisfies (r,δ)-restricted isometry property with δ≤ c/(κ^3√(r)log^2 d/α), η≤ cδ, and T_0 =  max{κlog(d/α)/η, 1/η√(dκα)}.
	For every t≲ T_0,
	
    U_tU_t^⊤-_F^2 ≤ (1 - η /(8 κ))^t - T_0 +O(α√(d) / κ^2).

	As a consequence, for T_1 =Θ((κlog (d/α))/η), we already have
	
    U_T_1U_T_1^⊤-_F^2 ≲α√(d) / κ^2.



When the condition number κ and rank r are both constant,  this theorem says that if we shoot for a final error ε, then we should pick our initialization _0 = α B  with α = O(ε/d). As long as the RIP-parameter δ = O(1/logd/), after O( logd/) iterations we will have that U_tU_t^⊤-_F^2 ≤.


Towards proving the theorem above, we suppose the eigen-decomposition of  can be written as = ^⊤ where ∈ℝ^d × r is an orthonormal matrix ∈ℝ^r × r is a diagonal matrix. We maintain the following decomposition of U_t throughout the iterations: 

    U_t = _S_tU_t_:=Z_t + ( - _S_t) U_t_:=E_t

Here S_t is r dimensional subspace that is recursively defined by 

    S_0    = ()  
    
    S_t    = (-η M_t-1)· S_t-1, ∀  t ≥ 1.

Here (-η M_t)· S_t-1 denotes the subspace {(-η M_t)v : v∈ S_t-1}. Note that (S_0)= () = r, and thus by induction we will have that for every t ≥ 0, 

    (Z_t)⊂(S_t), 
       (Z_t)  ≤(S_t) ≤ r.

Note that by comparison, in the analysis of rank-1 case, the subspace S_t is chosen to be () for every t, but here it starts off as () but changes throughout the iterations. We will show that S_t stays close to (). Moreover, we will show that the error term E_t — though growing exponentially fast — always remains much smaller than the signal term Z_t, which grows exponentially with a faster rate. Recall that sin(A,B) denotes the principal angles between the column span of matrices A,B. We summarize the intuitions above in the following theorem. 
	
	There exists a sufficiently small absolute constant c > 0 such that the following is true: For every α∈ (0, c/d), assume that the set of measurement matrices (A_1,…, A_m) satisfies (r,δ)-RIP with δ≤ c^4/(κ^3√(r)log^2 d/α). ρ = O(√(r)δ/κlog(d/α)), η≤ cδ.
	Then for t ≤ T_1 = Θ(κ/ηlog(d/α)), we have that
	
    sin(Z_t, ) ≲ηρ t 
       E_t≤ (1+ O(η^2 ρ t))^tE_0≤ 4E_0≤ 1/d
       σ_min(U^⋆^⊤ Z_t)≥E_t
       σ_min(^⊤ Z_t)≥min{(1+η/8κ)^tσ_min(^⊤  Z_0), 1/2√(κ)}
       Z_t≤ 5

	It follows from equation (<ref>) that after Θ(κ/ηlog(d/α)) steps, we have
	σ_min(Z_t) ≥1/2√(κ).

Note that the theorem above only shows that the least singular value of Z_t goes above 1/(2√(κ)). The following proposition completes the story by showing that once the signal is large enough, U_tU_t^⊤ converges with a linear rate to the desired solution  (up to some small error.) 



	In the setting of Theorem <ref>, suppose Z_t ≤ 5, sin(Z_t, ) ≤ 1/3, and σ_min(Z_t) ≥1/2√(κ), then we have:
	
    U_t+1U_t+1^⊤- _F^2 ≤ (1 - η/8 κ) U_tU_t^⊤-_F^2   +  O(η√(dr)E_t).



We defer the proof of Proposition <ref> to Section <ref>, which leverages the fact that function f satisfies the Polyak-Lojasiewicz condition <cit.> when U_tU_t^⊤ is well-conditioned. The rest of the section is dedicated to the proof outline of Theorem <ref>. We decompose it into the following propositions, from which Theorem <ref> follows by induction. 

The following proposition gives the base case for the induction, which straightforwardly follows from the definition _0 = α B where B is an orthonormal matrix. 
 
	In the setting of Theorem <ref>, we have 
	
    E_0 = α =  σ_min(Z_0) ≤ 1/d,  ^⊤ E_0 = 0,
    	  and  sin(Z_0, )   = 0.


The following Proposition bounds the growth rate of the spectral norm of the error E_t. 

	In the setting of Theorem <ref>,  suppose E_t≤ 1/d. Then, 
	
    E_t+1≤(1+ η O(δ√(r) + sin(Z_t, )  ) )×E_t.

	When U_tU_t^⊤ - is small, the growth of E_t
	becomes slower:
	
    E_t+1≤(1 + η O(U_tU_t^⊤ - _F + √(r)E_t)) ×E_t.


The following Proposition shows that effectively we can almost pretend that Z_t+1 is obtained from applying one gradient descent step to Z_t, up to some some error terms. 

	In the setting of Theorem <ref>,  suppose for some t we have 
	E_t≤ 1/d and Z_t≤ 5.	Then, 
	
    Z_t+1 - (Z_t - η(Z_t) - η E_t Z_t^⊤ Z_t -2 η_S_t M_t E_t )≤ητ_t,

	where τ_t ≲δ√(r)E_t.

The following proposition shows that the angle between the span of Z_t to the span of  is growing at mostly linearly in the number of steps. 

	In the setting of Theorem <ref>, assuming equation (<ref>) holds for  some t with τ_t ≤ρσ_r(U_t), Z_t≤ 5 and σ_min(Z_t) ≥E_t / 2. Then, as long as sin(Z_t, )  ≤√(ρ) we have:
	
    sin(Z_t+1, ) ≤sin(Z_t, ) + O(ηρ + ηE_t)    and  Z_t + 1≤ 5.


Then we show that the projection of the signal term Z_t to the subspace of  increases at an exponential rate (until it goes above 1/(2√(k))). Note that  is sufficiently close to the span of Z_t and therefore it implies that the least singular value of Z_t also grows. 



	In the setting of Theorem <ref>, suppose equation (<ref>) holds for some t, and Z_t≤ 5 and σ_min(Z_t) ≥E_t / 2, we have that:
	
    σ_min(^⊤ Z_t + 1)≥min{(1+η/8κ)σ_min(^⊤  Z_t), 1/2√(κ)}.


The proofs of the above propositions are deferred to Section <ref>. Now we are ready to prove Theorem <ref>. 

	When t = 0, the base case follows from Proposition <ref>.
	Assume that equations (<ref>),  (<ref>),  (<ref>), (<ref>), and (<ref>) are true before or at step t,
	we prove the conclusion for step t+1.
	By Proposition <ref> we have that equation <ref> are true with τ_t ≲δ√(r).
	We have that σ_r(U_t) ≥σ_r(Z_t), because the column subspace of Z_t and E_t are orthogonal.
	By induction hypothesis, σ_r(Z_t) ≥E_t / 2, and σ_r(U_t) ≥E_t / 2.
	Set the ρ in Proposition <ref> to O(δ√(r)/κlog(d/α)).
	When t ≤ T_1 = O(κ/ρlog(d/α)),
	we have that sin(Z_t, ) ≤ηρ t ≤√(ρ).
	Hence the assumptions of Proposition <ref> are verified.	Therefore,
	
    sin(Z_t+1, ) ≲ηρ t + O(ηρ) +  O(ηE_t) ≲ηρ(t+1)

	because E_t≤ 1/d ≪ρ, and Z_t≤ 5.
	Next,
	
    E_t+1   ≤(1 + η O(δ√(r) + sin(Z_t, )  ) )·E_t
       ≤ (1+O(η^2ρ T_1))^T_1E_0by equations  (<ref>),  (<ref>), and δ√(r)≤ O(ηρ T_1)
       ≤ 4E_0≤ 1/d since η^2ρ T_1^2 ≤ O(δ√(r)κlog(d/α)) ≤ O(1)

	Therefore we can apply Proposition <ref>
	to obtain that σ_min(^⊤ Z_t) grows by a rate of at least 1 + η/8κ.
	On the other hand by Proposition <ref>,
	E_t grows by a rate of at most 1 + η O(δ√(r) + √(ρ)) ≤ 1 + η/8κ.
	Hence we obtain Equation <ref>.


Finally we prove Theorem <ref>. 

		First of all, for t ≤ T_1 = Θ(κlogd/α / η), using Theorem <ref>, we know that the requirements in  Proposition <ref> is satisfied. Applying Proposition <ref>, we prove the theorem for t ≤ T_1. 

Then, we inductively show that the error is bounded by O(α√(d) / κ^2) from the T_1-th iteration  until the T_0-th iteration.
Suppose at iteration t, we have _t _t^⊤ - _F^2 ≲α√(d) / κ^2.
Thus, using Proposition <ref>, we know that E_t grows by a rate of at most 1 + η O(√(α) d^1/4 / κ) for this t.
This implies that for every t ∈ [T_1, T_0], we have

    E_t  + 1   ≤( 1 + η O(√(α) d^1/4 / κ) ) E_t
       ≤E_T_1(1 + η O(√(α) d^1/4 /κ))^T_0≤ 4 E_T_1≤ 16 E_0.

By inductive hypothesis we recall _t _t^⊤ - _F^2 ≲κ√(d) / κ^2, which implies by elementary calculation that σ_min(_t) ≥1/2 √(κ), Z_t ≤ 5 and sin(Z_t, ) ≤1/3 (by using _t = _t + _t and E_t≤ 16 E_0). Thus, the requirements in Proposition <ref> hold, and applying it we obtain that U_t+1U_t+1^⊤- _F^2≲α√(d) / κ^2. This completes the induction. 





 

§ NEURAL NETWORKS WITH QUADRATIC ACTIVATIONS

 
In this section, we state the algorithms and generalization bounds for learning over-parameterized neural nets with quadratic activations, and give the key lemma for the analysis. 

Let (x_1, y_1),…, (x_m, y_m) be n examples where x_i's are from distribution 𝒩(0, _d× d) and y_i's are generated according to equation (<ref>). 
Let  ŷ = 1^⊤ q(U^⊤ x)
be our prediction. We use mean squared loss as the empirical loss. 
For technical reasons, we will optimize a truncated version of the empirical risk as

    f̃(U) = 1/m∑_i=1^n (ŷ_i-y_i)^21_U^⊤ x^2≤ R

for some parameter R that will be logarithmic in dimension later.  We design a variant of gradient descent as stated in Algorithm <ref>.
We remark mostly driven by the analysis, our algorithm has an explicit re-scaling step. It resembles the technique of weight decay <cit.>, which has similar effect to that of an ℓ_2 regularization. In the noiseless setting, the issue with vanilla weight decay or ℓ_2 regularizer is that the recovery guarantees will depend on the strength of the regularizer and thus cannot achieve zero. An alternative is to use a truncated ℓ_2 regularizer  that only penalizes when U_t has norm bigger than a threshold. 
Our scaling that we are using is dynamically decided, and in contrast to this truncated regularizer, it scales down the iterate when the norm of U_t is small, and it scales up the iterate when U_t has norm bigger than the norm of .[But note that such scaling up is unlikely to occur because the iterate stays low-rank] Analyzing standard gradient descent is left for future work. 

 We note that one caveat here is that for technical reason, we assume that we know the Frobenius norm of the true parameter .  It can be estimated by taking the average of the prediction 1/m∑_i  = 1^m y_i since [y]  = _F^2, and the algorithm is likely to be robust to the estimation error of _F^2. However, for simplicity, we leave such a robustness analysis for future work. 

  




As alluded before, one-hidden-layer neural nets with quadratic activation closely connects to matrix sensing because we can treat write the neural network prediction by: 

    1^⊤  q(^⊤ x) = xx^⊤, ^⊤

Therefore, the i-th example x_i corresponds to  the i-th measurement  matrix in the matrix sensing via A_i = x_ix_i^⊤.  Assume { x_1x_1^⊤, …,  x_mx_m^⊤} satisfies RIP, then we can re-use all the proofs for matrix sensing. However, this set of rank-1 measurement matrices doesn't satisfy RIP with high probability. The key observation is that if we truncated the observations properly, then we can make the truncated set of these rank-1 measurements  satisfy RIP property again. Mathematically, we prove the following Lemma. 



Let (A_1,…, A_m) = { x_1x_1^⊤, …,  x_mx_m^⊤} where x_i's are  i.i.d. from ∼𝒩(0, ). Let  R = log(1/δ). Then, for every q, δ∈ [ 0, 0.01] and m ≳ d   log^4 d/q δ/δ^2, we have:  with probability at least 1 - q, for every symmetric matrix X,

    1/m∑_i = 1^m ⟨ A_i, X ⟩ A_i 1_|⟨ A_i, X ⟩| ≤ R - 2 X - (X) ≤δ_⋆


Suppose X has rank at most  r matrices and spectral norm at most 1, we have, 

    1/m∑_i = 1^m ⟨ A_i, X ⟩ A_i 1_|⟨ A_i, X ⟩| ≤ R - 2 X - (X) ≤ rδ


We can see that equation (<ref>) implies Lemma <ref> with a simple change of parameters (by setting δ to be a factor of r smaller). The proof uses standard technique from supreme of random process, and is deferred to Section <ref>. 

Theorem <ref> follows straightforwardly from replacing the RIP property by the Lemma above. We provide a proof sketch below. 


The basic ideas is to re-use the proof of Theorem <ref> at every iteration.  First of all, we will replace all the RIP properties [We only require RIP on symmetric matrices.] in (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>)  by Lemma <ref>.  The only difference is that we will let the δ when applying Lemma <ref> to be 1/r smaller than the δ in Theorem <ref>.

We note that in Lemma <ref>, there is an additional scaling of the identity term compared to Lemma <ref>. This is the reason why we have to change our update rule.
We note that the update rule in Algorithm <ref> undo the effect of this identity term and is identical to the update rule for matrix sensing problem:
Let c_t = (U_t U_t^⊤) - ().
Denote by M_t' = M_t - c_t.
The update in Algorithm <ref> can be re-written as:

    _t + 1   = ( - η M'_t - η c_t  )_t
    
       = (1 - η c_t) ( - η/1 - η c_t M'_t ) _t

Hence we still have _t + 1 = 1/1 - η c_t_t + 1 = ( - η/1 - η c_t M'_t ) _t. Thus the update rule here corresponds to the update rule for the matrix sensing case, and the rest of the proof follows from the proof of Theorem <ref>.



















§ SIMULATIONS


In this section, we present simulations to complement our theoretical results.
In the first experiment, we show that the generalization performance of gradient descent depends on the choice of initialization.
In particular, smaller initializations enjoy better generalization performance than larger initializations.
In the second experiment, we demonstrate that gradient descent can run for a large number of  iterations and the test error keeps decreasing, which suggests early stopping is not necessary. 
In the third experiment, we show that a natural projected gradient descent procedure works poorly compared to gradient descent on the factorized model, which suggests the power of using a factorized model.
In the last experiment, we report results for running stochastic gradient descent on the quadratic neural network setting, starting from a large initialization.

We generate the true matrix by sampling each entry of  independently from a standard Gaussian distribution and let = ^⊤. 
Each column of  is normalized to have unit norm, so that the spectral norm of  is close to one.
For every sensing matrix A_i, for i = 1,…,m, we sample the entries of A_i independently from a standard Gaussian distribution.
Then we observe b_i = A_i.
When an algorithm returns a solution X̂, we measure training error by:

    √(∑_i=1^m (A_iX̂ - b_i)^2/∑_i=1^m b_i^2.)

We measure test error by:

    X̂ - /.

For the same , we repeat the experiment three times, by resampling the set of sensing matrices {A_i}_i=1^m.
We report the mean and the error bar.



  
Choice of initialization.
Let U_0 = α.
We use m = 5dr samples, where rank r = 5.
We plot the training and test error for different values of α.
Figure <ref> shows that the gap between the training and test error
narrows down as α decreases.
We use step size 0.0025 and run gradient descent for 10^4 iterations.



In Figure <ref>, we run for longer iterations to further compare the generalization performance of initialization U_0 = α for α = 1.0, 10^-3.
We report the mean values at each iteration over three runs.
When α = 1.0, despite the training error decreases below 10^-4,
the test error remains to be on the order of 10^-1.







  
Accuracy.
In this experiment, we fix the initialization to be U_0 = 0.01,
and apply the same set of parameters as the first experiment.
We keep gradient descent running for 10^5 iterations, to see if test error
keeps decreasing or diverges at some point.
Figure <ref> confirms that test error goes down gradually.







  
Projected gradient descent.
In this experiment, we consider the following natural projected gradient descent (PDG) procedure.
Let f(X) = 1/m∑_i=1^m (A_iX - b_i )^2.
At every iteration, we first take a gradient step over f(X),
then project back to the PSD cone.
We consider the sample complexity of PGD by varying the number of sensing
matrices m from 5d to 35d.
Here the rank of  is 1.
We found that the performance of PGD is much worse compared to gradient descent on the factorized model.
For both procedures, we use step size equal to 0.0025.
We stop when the training error is less than 0.001, or when the number of iterations reaches 10^4.
Figure <ref> shows that gradient descent on the factorized model
consistently recovers  accurately.
On the other hand, the performance of projected gradient descent gets even worse as d increases from 100 to 150.





  
Stochastic gradient descent.
In this experiment, we complement our theoretical results by running stochastic gradient descent from large initializations.
We generate m = 5dr random samples and compute their true labels as the training dataset.
For stochastic gradient descent, at every iteration we pick a training data point uniformly at random from the training dataset.
We run one gradient descent step using the training data point.
We initialize with U_0 = and use step size 8 × 10^-5.
Figure <ref> shows that despite the training error decreases to 10^-7, the test error remains large.
We also report the results of running gradient descent on the same instance for comparison.
As we have already seen, Figure <ref> shows that gradient descent also gets stuck at a point with large test error, despite the training error being less than 10^-7.







§ CONCLUSIONS


The generalization performance of over-parameterized non-linear models, especially neural networks, has been a very intriguing research topic. This paper theoretically explains the regularization power of the optimization algorithms in learning matrix factorization models and one-hidden-layer neural nets with quadratic activations. In these cases, the gradient descent algorithm prioritizes to searching through the low complexity models. 

It's an very interesting open question to establish similar results for deeper neural networks with other activations (e.g., ReLU) and loss functions (e.g., logistic loss). We remark that likely such a result will require not only a better understanding of statistics, but also a deep grasp of the behavior of the optimization algorithms for non-linear models, which in turns is another fascinating open question.  



 §.§ Acknowledgments: 


We thank Yuxin Chen, Yann Dauphin, Jason D. Lee, Nati Srebro, and Rachel A. Ward for insightful discussions at various stages of the work. 
Hongyang Zhang is supported by NSF grant 1447697.



§ PROOF OF MAIN PROPOSITIONS




 §.§ Proof of Proposition <ref>

We start off with a straightforward triangle inequality for bounding E_t+1 given E_t.  

	By update rule (equation (<ref>)) and the definition of E_t (equation (<ref>)), we have that 
	
    E_t+1 = (-_S_t+1)( - η M_t)E_t

	It follows that 
	
    E_t+1≤( - η M_t)E_t≤E_t + ηM_tE_t


Therefore, next we will bound the norm of M_tE_t. The key idea is to use the restricted isometry property to control the effect of M_t (see the meta claim equation (<ref>) in Section <ref> for more intuitions.)

	In the setting of Proposition <ref>, we have that 
	
    M_tE_t≤E_t(δZ_tZ_t^⊤ -  + 2δZ_tE_t^⊤ + δE_tE_t^⊤_⋆ + Z_t E_t^⊤ + ( - _S_t))

	As a direct consequence, using the assumption E_t≤ 1/d,
	
    E_t + 1≤E_t( 1 + O (ηδ√(r)  + η( - _S_t)) )




	We first note that E_t = ( - _S_t)E_t by the update rule and definition <ref>. It follows that M_tE_t = M_t ( - _S_t)E_t. 
	Next, since U_tU_t^⊤ - -E_tE_t^⊤ has rank at most 4r, by Lemma <ref> we have that 
	
    1/m∑_i = 1^m A_i, U_tU_t^⊤ - -E_tE_t^⊤A_i( - _S_t)E_t - (U_tU_t^⊤ - -E_tE_t^⊤)( - _S_t)E_t
    ≤    δ×U_tU_t^⊤ - -E_tE_t^⊤( - _S_t)E_t
    
    	=    δ×(Z_tZ_t^⊤ -  + 2Z_tE_t^⊤_F )E_t

	Note that (U_tU_t^⊤ - -E_tE_t^⊤)( - _S_t)E_t = Z_tE_t^⊤ E_t - ( - _S_t)E_t. It follows that 
	
    1/m∑_i = 1^m A_i, U_tU_t^⊤ - -E_tE_t^⊤A_i( - _S_t)E_t
       ≤( - _S_t)E_t + Z_tE_t^⊤ E_t + δ(Z_tZ_t^⊤ -  + 2Z_tE_t^⊤_F)E_t
       ≤E_t(δZ_tZ_t^⊤ -  + 2δZ_tE_t^⊤ + Z_tE_t^⊤ + ( - _S_t))

	By Lemma <ref>, we have that 
	
    1/m∑_i = 1^m A_i, E_tE_t^⊤A_i( - _S_t)E_t - E_tE_t^⊤ ( - _S_t)E_t≤δE_tE_t^⊤_⋆E_t

	Combining equation (<ref>) and (<ref>) we complete the proof of equation (<ref>). To prove equation (<ref>), we will use equation (<ref>) and that Z_tZ_t^⊤ - ≲√(r) (Corollary <ref>),
	Z_t^2 ≤√(r)×Z_t Z_t^⊤≲ r.


When U_t U_t^⊤ - is small, the growth of E_t becomes slower.

	In the setting of Proposition <ref>, we have that
	
    M_t≤U_tU_t^⊤ -  + δU_tU_t^⊤ -  + 3√(r)E_t.

	As a direct consequence, M_t≤ O(1). And
	
    M_tE_t≲E_t(U_tU_t^⊤ - _F + √(r)E_t)




	By the definition of M_t, we have that
	
    M_t    = 1/m∑_i=1^m A_iU_tU_t^⊤ - A_i
        = 1/m∑_i=1^m A_iU_tU_t^⊤ -  - E_tE_t^⊤A_i + 1/m∑_i=1^m A_iE_tE_t^⊤A_i
       ≤U_tU_t^⊤ -  + δ×(U_tU_t^⊤ -  - E_tE_t^⊤ + E_t E_t^⊤) 
       ≤ (1 + δ) U_tU_t^⊤ -  + 3√(r)E_t

	where the second to last line is because of Lemma <ref> and Lemma <ref>.
	For the last line, we use
	
    U_tU_t^⊤ -  - E_tE_t^⊤ = Z_tZ_t^⊤ -  + 2 Z_tE_t^⊤≲Z_tZ_t^⊤ -  + √(r)E_t,

	because of Corollary <ref>.
	And Z_t Z_t^⊤- _F ≤U_t U_t^⊤ - _F  + O(√(r)E_t), because
	
    U_t U_t^⊤ - _F    = (Z_t  +E_t) (Z_t + E_t)^⊤ - _F 
       ≥Z_t Z_t^⊤- _F - 2Z_t_F E_t - E_tE_t_F



Finally we complete the proof of Proposition <ref>.

	Using the fact that  has spectral norm less than 1, wee can bound the term ( - _S_t) in equation (<ref>) by
	
    ( - _S_t)≤^⊤ ( - _S_t) = ( - _S_t)   = sin(S_t, )

	Since S_t is the column span of Z_t, using the equation above and equation (<ref>) we conclude the proof. 





 §.§ Proof of Proposition <ref>


We first present a simpler helper lemma that will be used in the proof. 
One can view X in the following lemma as a perturbation. The lemma bounds the effect of the perturbation to the left hand side of equation (<ref>). 

	Let S ∈^d × r be a column orthonormal matrix and
	S^ =  - SS^⊤ be its orthogonal complement.
	Let X ∈^d× d be any matrix where X < 1/3.
	We have:
	
    _( - X) S ( - X) S^ + 2_S X S^≲X^2.
 	



	By definition, we know that:
	
    _( - X) S   = ( - X) S ( S^⊤( - X)^⊤ ( - X) S )^-1  S^⊤( - X)^⊤
       = ( - X) S ( - S^⊤ (X^⊤  + X) S  + S^⊤ X^⊤ X S )^-1 S^⊤( - X)^⊤

	Denote by Y = S^⊤ (X^⊤ + X) S - S^⊤ X^⊤ X S.
	We have that Y≤ 2 X + X^2 < 7 X / 3.
	By expanding ( - Y)^-1, we obtain:
	
    ( - Y)^-1 - ( + Y)   ≤∑_i=2^∞Y^i 
       ≤Y^2/1 - Y≤ 25 X^2.

	Hence we get:
	
    ( - Y)^-1 - ( + S^⊤ (X^⊤ + X) S)≤ 26 X^2.

	Denote by
	
    A    = ( - X) S ( + S^⊤ (X^⊤ + X) S) S^⊤ ( - X^⊤) 
       = ( - X) (_S + _S (X^⊤ + X) _S) ( - X^⊤).

	We separate the terms in A which has degree 1 in X:
	
    A_1    = _S - X _S - _S X^⊤ + _S (X^⊤ + X) _S

	Consider the spectral norm of A ( - X) S^.
	We have that A_1 S^ = -_S X S^.
	For - A_1 X S^, the only term which has degree 1 in X is -_S X S^.
	To summarize, we obtain by triangle inequality that:
	
    _( - X) S ( - X) S^⊤ + 2_S X S^   ≤A + 2_S X S^ + 26 X^2 
       ≲A_1 + 2_S X S^ + X^2 ≲X^2



Now we are ready present the proof of Proposition <ref>.


	We first consider the term _t Z_t:
	
    _t _t    = 1/m∑_i = 1^m ⟨_i, _t _t^⊤ - ⟩_i _t 
        = 1/m∑_i = 1^m ⟨_i, Z_t Z_t^⊤ -⟩_i _t 	
        + 1/m∑_i = 1^m ⟨_i, _t Z_t^⊤ + Z_t _t^⊤⟩_iZ_t +  1/m∑_i = 1^m ⟨_i,  _t _t^⊤⟩_i _t

	Note that ∇ f(Z_t) =  1/m∑_i = 1^m ⟨_i, Z_t Z_t^⊤ -⟩_i _t. 
	Using  Lemma <ref> and <ref> on the two terms in line (<ref>) with the fact that Z_t^⊤ E_t = 0, we have that
	
    M_tZ_t - E_t Z_t^⊤ Z_t - (Z_t)≤ 
    	2 √(r)δ_t Z_t^2 +  δ d _t ^2 Z_t

We decompose Z_t+1 by: 
	
    Z_t + 1   = _S_t + 1 U_t + 1 = U_t + 1 - ( - _S_t + 1) U_t + 1
        = ( - η M_t ) Z_t + ( - η M_t ) _t  -  ( - _S_t + 1)  ( - η M_t ) _t  
       =  Z_t - η M_t  Z_t  + _S_t + 1 (  - η M_t ) E_t,

	where in the second equation we use the fact that
	( - _S_t+1) ( - η M_t) Z_t = 0,
	since the column subspace of Z_t is S_t.
	Setting S = _S_t and X = η M_t in Lemma <ref>, 
	we have that _S_t+1 = _( - X) S. Applying Lemma <ref>, we conclude:
	
    _S_t+1 ( - η M_t) S^ + 2_S ( - η M_t) S^≤η^2 M_t^2,

	which implies that:
	
    _S_t+1 ( - η M_t) E_t + 2_S ( - η M_t) E_t≤η^2 E_tM_t^2.

	Since M_t≤ O(1) by Lemma <ref>,
	the conclusion follows by combining Equation (<ref>) and (<ref>) with the equation above.





 §.§ Proof of Proposition <ref>


Towards proving Proposition <ref>, we further decompose Z_t into 

    Z_t = ( + F_t)R_t

where R_t∈^r× d, F_t∈^d× r are defined as

    R_t = ^⊤ Z_t,  F_t = (- _)Z_t R_t^+

Recall that X^+ denotes the pseudo-inverse of X. We first relate the the spectral norm of F_t with our target in Proposition <ref>, the principal angle between Z_t and .  Up to third order term, we show that F_t is effective equal to the principle angle, and this is pretty much the motivation to decompose Z_t in equation (<ref>). 




	Let F_t be defined as in equation (<ref>). Then, if F_t < 1/3, we have that 
	
    F_t - F_t^3 ≤sin(Z_t, ) ≤F_t



	Let S_t = ( + F_t)( + F_t^⊤ F_t)^-1/2. By the fact that ^⊤ F_t = 0, we have S_t^⊤ S_t = and S_t has the same column span as Z_t. Therefore, the columns of S_t form an orthonormal basis of Z_t, and we have that 
	
    sin(Z_t, ) = ( - _)S_t = F_t( + F_t^⊤ F_t)^-1/2

	Suppose F_t has singular value σ_j, j=1,…, r, then it's straightforward to show that F_t( + F_t^⊤ F_t)^-1/2 has singular values σ_j/√(1+σ_j^2), j=1,…, r. The conclusion then follows basic calculus and the fact that maxσ_j ≤ 1/3. 


Therefore, it suffices to bound the spectral norm of F_t.  However, the update rules of F_t or Z_t are difficult to reason about. Therefore, we introduce the following intermediate term Z̃_t that bridges Z_t+1, F_t+1 with Z_t and F+t. We define Z̃_t as:

    Z̃_t = ( -  η E_t Z_t^⊤)Z_t (  -  2 η Z_t^+_S_t M_t E_t )


The motivation of defining such Z_t is that Z_t+1 depends on Z_t via relatively simple formula as the lemma below shows:

	In the setting of Proposition <ref>,

    Z_t + 1 - (Z̃_t - η∇ f(Z_t) )= O(ητ_t)



The proof of Proposition <ref> is deferred to the later part of this section. We also decompose Z̃_t same as Z_t to Z̃_t = ( + F̃_t) R̃_t, where 

    R̃_t = ^⊤Z̃_t, F̃_t = (- _)Z̃_t R̃_t^+


We will prove that R̃_t is close to R_t and F̃_t is close to F_t in the following sense:
In the setting of Proposition <ref>,

    σ_min(R̃_t ) ≥σ_min(R_t)( 1 - η/100 κ)


In the setting of Proposition <ref>,

    F̃_t - F_t≲ηE_t  + ηρ



We will prove these propositions in the following sections. After that, we will focus on the update from Z̃_t to Z_t + 1. In particular, we will bound R_t + 1 and F_t + 1 directly using R̃_t and F̃_t.



  §.§.§ Proofs of Proposition <ref>, <ref>, and <ref>




By definition of Z̃_t, 

    Z̃_t    =  ( -  η E_t Z_t^⊤)Z_t (  -  2 η Z_t^+_S_t M_t E_t )
    
       = Z_t - η( E_t Z_t^⊤ Z_t + 2 _S_t M_t E_t ) + 2 η^2  E_t Z_t^⊤_S_t M_t E_t


Recall that by Proposition <ref>, the update rule of Z_t + 1 satisfies 

    Z_t + 1 - (Z_t  - η∇ f(Z_t)- η E_t Z_t^⊤ Z_t - 2 η_S_t M_t E_t  )≤ητ_t


Putting the above two formulas together (using the bound of M_t E_t as in Lemma <ref>) we have that

    Z_t + 1 - (Z̃_t - η∇ f(Z_t) ) = O(η^2 E_t+  ητ_t)


On the other hand, a direct calculation also shows that Z̃_t - Z_t = O(ηE_t). Moreover, Z̃_t is a rank r matrix and Z_t, Z̃_t = O(1). Therefore, let us denote by Z̃_t = Z_t + Δ_t, we have:

    ∇ f(Z̃_̃t̃)    = 1/m∑_i = 1^m ⟨ A_i, Z̃_t Z̃_t^⊤  - X^* ⟩ A_i Z̃_t 
    
        = 1/m∑_i = 1^m ⟨ A_i, Z_t Z_t^⊤  - X^* ⟩ A_i Z̃_t  +  1/m∑_i = 1^m ⟨ A_i, Z_t Δ_t^⊤ + Δ_t Z_t^⊤⟩ A_i Z̃_t  +  1/m∑_i = 1^m ⟨ A_i, Δ_t Δ_t^⊤⟩ A_i Z̃_t 
    
       = 1/m∑_i = 1^m ⟨ A_i, Z_t Z_t^⊤  - X^* ⟩ A_i Z_t+  1/m∑_i = 1^m ⟨ A_i, Z_t Z_t^⊤  - X^* ⟩ A_i  Δ_t
    
       +  1/m∑_i = 1^m ⟨ A_i, Z_t Δ_t^⊤ + Δ_t Z_t^⊤⟩ A_i Z̃_t  +  1/m∑_i = 1^m ⟨ A_i, Δ_t Δ_t^⊤⟩ A_i Z̃_t



Note that Δ_t is at most rank 2r, therefore, we can apply the RIP property of { A_i }_i = 1^m (Lemma <ref>) and Z̃_̃t̃, 
Z_t = O(1) to conclude that 

    1/m∑_i = 1^m ⟨ A_i, Z_t Z_t^⊤  - X^* ⟩ A_i  Δ_t
    +    ∑_i = 1^m ⟨ A_i, Z_t Δ_t^⊤ + Δ_t Z_t^⊤⟩ A_i Z̃_t  +    ∑_i = 1^m ⟨ A_i, Δ_t Δ_t^⊤⟩ A_i Z̃_t  = O(Δ_t)


Recall that Δ_t = O(η E_t ), the equation above and the formula for ∇ f(Z̃) implies that

    Z_t + 1 - (Z̃_t - η∇ f(Z̃_̃t̃) ) = O(η^2 E_t+  ητ_t) = O(ητ_t)




The above proposition implies that Z_t + 1 is very close to doing one step of gradient descent from Z̃_t. Thus, we will mainly focus on Z̃_t in the later section. 


By definition of R̃_t, we know that 

    R̃_t    = ^⊤Z̃_t = ^⊤( -  η E_t Z_t^⊤)Z_t (  -  2 η Z_t^+_S_t M_t E_t )
       = ^⊤Z_t (  -  2 η Z_t^+_S_t M_t E_t ) -   η^⊤ E_t Z_t^⊤ Z_t (  -  2 η Z_t^+_S_t M_t E_t )



By definition of R_t and the assumption that F_t ≤ 1/3, we have that σ_r(Z_t) ≥1/2σ_min (R_t), which implies that Z_t^+≤ 2 R_t^+. Thus,
using the bound of M_t E_t as in Lemma <ref>,

    Z_t^+_S_t M_t E_t ≲ M_t E_t/σ_min(R_t)≲E_t (δ√(r) + sin(Z_t, ))/σ_min(R_t)



Similarly, we can get: 

    ^⊤E_t  = ^⊤( - _S_t)E_t ≤E_t^⊤( - _S_t) =  E_tsin(Z_t, )

Therefore, bounding the terms in equation (<ref>) using the bounds above, we have:

    σ_min(R̃_t )    ≥σ_min(^⊤Z_t) (1 - η Z_t^+_S_t M_t E_t) - η^⊤ E_t Z_t^⊤ Z_t (  -  2 η Z_t^+_S_t M_t E_t )
       ≥σ_min(R_t)  (1 - η Z_t^+_S_t M_t E_t) - 2η^⊤ E_t
       ≥σ_min(R_t)( 1 - η O(  E_t (δ√(r) + sin(Z_t, ))/σ_min(R_t)) )



Again, using σ_r(Z_t) ≥1/2σ_min (R_t) we know that 

    E_t (δ√(r) + sin(Z_t, ))/σ_min(R_t)   ≤ 2E_t (δ√(r) + sin(Z_t, ))/σ_min(Z_t)


By assumption of Proposition <ref> and that σ_min(Z_t) ≥1/2E_t, we obtain:

    2E_t (δ√(r) + sin(Z_t, ))/σ_min(Z_t)   ≤ 2(δ√(r) + sin(Z_t, ))
       ≲√(ρ)by Assumption in proposition <ref>


Thus, by equation (<ref>) above, and our choice of parameter such that √(ρ)≲ 1/κ, we conclude that

    σ_min(R̃_t ) ≥σ_min(R_t)( 1 - η/100 κ)




Next we prove prove Proposition <ref>, which focus on F̃_t:


We know that right multiply Z_t by an invertible matrix does not change the column subspace as Z_t, thus won't change the definition of F_t, so we can just focus on ( -  η E_t Z_t^⊤)Z_t. 
We know that 

    F̃_t    = (- _) ( -  η E_t Z_t^⊤)Z_t (^⊤ ( -  η E_t Z_t^⊤)Z_t)^+
       = (- _) Z_t (^⊤ Z_t - η^⊤ E_t Z_t^⊤ Z_t )^+
        - η  (- _)E_t Z_t^⊤(^⊤ Z_t - η^⊤ E_t Z_t^⊤ Z_t )^+



Since Z_t is rank r, we can do the SVD of Z_t as Z_t =  V Σ_Z W^⊤ for column orthonormal matrices V, W ∈ℝ^d × r and diagonal matrix Σ_Z ∈ℝ^r × r. Then, we can write (^⊤ Z_t - η^⊤ E_t Z_t^⊤ Z_t )^+ as:

    (^⊤ Z_t - η^⊤ E_t Z_t^⊤ Z_t )^+   = ((^⊤ V- η^⊤ E_t Z_t^⊤V ) Σ_Z W^⊤)^+
       =  W^⊤Σ_Z^-1( ^⊤V - η^⊤_t Z_t^⊤ V)^-1



By our assumption that sin(Z_t, ) ≤1/3 we know that σ_min(^⊤ V) = Ω(1). Thus, by Z_t = O(1) and Woodbury matrix inversion formula we have:

    ( ^⊤V - η^⊤_t Z_t^⊤ V)^-1 -  (^⊤V )^-1   ≲η^⊤_t Z_t^⊤ V≲η^⊤_t





Recall that F_t  = ( - _) Z_t R_t^+ = ( - _) Z_t W^⊤Σ_Z^-1 (^⊤)^-1. Thus, using σ_r(Z_t) ≥1/2σ_min (R_t), an elementary calculation gives us:

    F̃_t - F_t   ≲ηE_t + ηF_t^⊤ E_t /σ_min(R_t)


By assumption E_t≲σ_min(Z_t) in Proposition <ref>, together with σ_r(Z_t) ≥1/2σ_min (R_t) and inequality <ref>, we have

    F_t^⊤ E_t /σ_min(R_t)≲F_tE_tsin(Z_t, )/σ_min(Z_t)≲F_t^2

By our choice of parameter, we know that F_t^2 ≤ρ, therefore, we proved that F̃_t - F_t≤η(  E_t_2 + ρ ) as desired.


Now, we can just focus on Z̃_t. One of the crucial fact about the gradient (Z̃_t) is that it can be decomposed into 

    (Z̃_t) = N_t R̃_t

where N_t is a matrix defined as 

    N_t =  1/m∑_i = 1^m_i, Z̃_tZ̃_t^⊤ - _i( + F̃_t)

Therefore, Z̃_t and (Z̃_t) share the row space and we can factorize the difference between Z̃_t and η_t(Z̃_t) as

    Z̃_tN_t - η(Z̃_t) = (F̃_t - η N_t)R̃_t


Note that the definition of N_t depends on the random matrices A_1,…, A_t. The following lemma show that for our purpose, we can essentially view N_t as its population version — the counterpart of N_t when we have infinitely number of examples. The proof uses the RIP properties of the matrices A_1,…, A_m. 




	In the setting of Proposition <ref>, let N_t be defined as in equation (<ref>). Then, 
	
    N_t - (Z̃_tZ̃_t^⊤ - X^⋆)( + F̃_t)≤ 2δZ̃_tZ̃_t^⊤ - X^⋆_F



Recalling the definition of N_t, 
by Lemma <ref>, we have that 

    N_t - (Z̃_tZ̃_t^⊤ - X^⋆)( + F̃_t)   ≤δZ̃_tZ̃_t^⊤ - X^⋆_F  + F̃_t
       ≤ 2δZ̃_tZ̃_t^⊤ - X^⋆_F by the assumption F̃_t< 1/3





For any t ≥ 0, suppose Z_t + 1 - (Z̃_t - η𝒢(Z̃_t))≤ητ, we have

    Z_t + 1≤Z̃_t ( 1 - 1/2ηZ̃_t^2 ) + 2 ηZ̃_t  X^⋆ + ητ



	By Lemma <ref>, we have: 
	
    𝒢(Z̃_t) - (Z̃_t Z̃_t^⊤ - X^⋆ )Z̃_t  ≤δZ̃_t Z̃_t^⊤ - X^⋆_FZ̃_t

	Therefore, 
	
    Z_t + 1   ≤Z̃_t  - η𝒢(Z̃_t)  + ητ
       ≤Z̃_t   - η (Z̃_t Z̃_t^⊤ - X^⋆ )Z̃_t  + ηδZ̃_t Z̃_t^⊤ - X^⋆_FZ̃_t  + ητ
       ≤( - ηZ̃_t Z̃_t^⊤ + η X^⋆)Z̃_t  + ηδ√(r)Z̃_t(Z̃_t^2 + X^⋆) + ητ
       ≤ (1 - ηZ̃_t) Z̃_t  + 1/2ηZ̃_t( Z̃_t^2 + 4X^⋆)+ ητby δ√(r)≤ 1/2
	


As a direct corollary, we can inductive control the norm of Z̃_t. 

	In the setting of Proposition <ref>, we have that 
	
    Z_t+1≤ 5, R̃_t≤ 6

	Moreover, 
	
    Z_t+1Z_t+1^⊤ - X^⋆_F≲√(r), N_t≲√(r)



	Using the assumption that =1 and the assumption that equation (<ref>) holds, then we have that Z̃_t≤Z_t(1 + O(ηE_t)) ≤ 5(1 + O(ηE_t)). Applying Lemma <ref> with τ = O(τ_t), we have that Z_t+1≤ 5. We also have that R̃_t≤Z̃_t≤ 6. Moreover, we have Z_t+1Z_t+1^⊤ - X^⋆_F ≤Z̃_t_⋆ + X^⋆_F ≲√(r)(Z̃_t+ X^⋆)≲√(r). As a consequence, N_t≤ (1 + δ)Z̃_t Z̃_t^⊤ - X^⋆_F( U^⋆ + F̃_t)≲√(r). 



We start off with a lemma that controls the changes of R̃_t relatively to R_t+1.  

	In the setting of Proposition <ref>, then we have that R̃_tR_t+1^-1 can be written as
	
    R̃_tR_t+1^+ =  + ηR̃_tR̃_t^⊤ - η + ξ_t^(R)

where ξ_t^(R)≲ηδ√(r) + ηρ + ηF̃_t + η^2. It follows that R̃_tR_t+1^+≤ 4/3 and τ≤ 2 ρσ_min(R_t+1). 


	By the definition of R̃_t and equation (<ref>), we have that
	
    ητ   ≥_t + 1 - R̃_t - η^⊤(Z̃_t) 
        = _t + 1 - ( - η^⊤ N_t)R̃_t

Form this we can first obtain a very weak bound on σ_min(R_t+1):
		
    σ_min(_t + 1)    ≥σ_min(( - η^⊤ N_t)R̃_t) - ητ
       ≥3/4·σ_min(R̃_t) - O(ηρσ_min(R̃_t)) by τ≤ O(ρσ_min(R̃_t)) and η^⊤ N_t≤η N_t≲η√(r)≤ 1/4
       ≥1/2·σ_min(R̃_t)

	Re-using equation (<ref>), we have 
		
    - ( - η^⊤ N_t)R̃_tR_t+1^+   ≤ητ·R_t+1^+ = ητ/σ_min(R_t+1)
       ≤ 2ητ/σ_min(R̃_t) ≤ O(ηρ)

where we used equation (<ref>) and τ≤ O(ρσ_min(R̃_t)). This also implies a weak bound for R̃_tR_t+1^+ that R̃_tR_t+1^+≤ 2. 
By Lemma <ref>, we have that 

    η^⊤ N_t -η^⊤ (Z̃_tZ̃_t^⊤ - X^⋆)( + F̃_t)≤ 2δηZ̃_tZ̃_t^⊤ - X^⋆_F≲δη√(r)

where we use the fact that Z̃_tZ̃_t^⊤ - X^⋆_F≲√(r).  Note that = ^⊤ and Z̃_t = ( + F̃_t)R̃_t,  we have η^⊤ N_t -η (R̃_tR̃_t^⊤-) (+F̃_t)^⊤( +F̃_t)≲δη√(r). 

Bounding the higher-order term, we have that 

    η (R̃_tR̃_t^⊤-)(+F̃_t)^⊤ ( + F̃_t)- η (R̃_tR̃_t^⊤-)≲ηF̃_tby R̃_t≤ 6

which implies that 

    η^⊤ N_t -η (R̃_tR̃_t^⊤-)≲δη√(r) + ηF̃_t

Combining the equation above with equation (<ref>) and R̃_tR_t+1^+≤ 2, we have that 
	
    - ( - ηR̃_tR̃_t^⊤ + η) R̃_tR_t+1^+   ≲ηδ√(r)  + ηρ + ηF̃_t

For η≲ 1, we know that 

    ( - ηR̃_tR̃_t^⊤ + η)  ( + ηR̃_tR̃_t^⊤ - η) - ≲η^2

This implies that 

    ( + ηR̃_tR̃_t^⊤ - η) - R̃_tR_t+1^+   ≲ηδ√(r) + ηρ + ηF̃_t + η^2

which completes the proof. 




We express F_t+1 as a function of F̃_t and other variables. 

In the setting of Proposition <ref>, let N_t be defined as in equation (<ref>). Then, 

    F_t+1 = F̃_t( -ηR̃_tR̃_t^⊤ )R̃_tR_t+1^+ + ξ_t^(F)
	
where 
ξ_t^(F)≲δη√(r) + ηρ+ ηF̃_t^2.



By equation (<ref>), we have that 

    ( - _)·(Z_t+1 - (Z̃_t - η(Z̃_t)))≤ητ_t

which, together with the decomposition (<ref>), implies  

    ητ   ≥_t + 1_t + 1 -  F̃_t R̃_t + η ( - _) (Z̃_t) 
       = 	_t + 1_t + 1 -  (F̃_t -η ( - _) N_t)R̃_t

Recall that τ_t ≤ 2ρσ_min(R_t+1) (by Lemma <ref>), we conclude

    _t + 1-  (F̃_t -η( - _)  N_t)R̃_t R_t+1^+≤ 2ηρ

Note that ( -_)X^⋆ = 0 and that ( -_)Z̃_tZ̃_t^⊤ = F̃_tR̃_tR̃_t^⊤(+F̃_t)^⊤. We obtain that 


    ( -_)N_t-  F̃_tR̃_tR̃_t^⊤(+F̃_t)^⊤( + F̃_t)   ≤ 2δZ̃_tZ̃_t^⊤ - X^⋆_F ≲δ√(r)

where we used the fact that Z̃_tZ̃_t^⊤ - X^⋆_F  ≲√(r) (by Lemma <ref>). 

Bounding the higher-order term, we have that

    F̃_tR̃_tR̃_t^⊤(+F̃_t)^⊤ ( + F̃_t)- F̃_tR̃_tR̃_t^⊤≲F̃_t^2 by R̃_t≤ 6 from Corollary <ref>

Combining equation (<ref>), (<ref>) and the equation above, and using the fact that R̃_tR_t+1^+≤ 2, we complete the proof.






Combining Lemma <ref> and Lemma <ref>, we can relate the F_t+1 with F̃_t:

	In the setting of Proposition <ref>, we have that F̃_t can be written as 
	
    F_t+1 = F̃_t( - η) + ξ_t

	where ξ_t≤ O(ηδ√(r)  + ηF̃_t^2 +  ηρ + η^2). As a consequence, 
	
    F_t+1≤F̃_t + O(ηF̃_t^2 + ηρ)



	Combine Lemma <ref> and Lemma <ref>, we have that 
	
    F_t + 1   = F̃_t( -ηR̃_tR̃_t^⊤ )R̃_tR_t+1^+ + ξ_t^(F)
       = F̃_t( -ηR̃_tR̃_t^⊤)( + ηR̃_tR̃_t^⊤ - η + ξ_t^(R)) + ξ_t^(F)

	
	Thus, with the bound on ξ_t^(R) and ξ_t^(F) from Lemma <ref> and Lemma <ref>, we know that 
	
    F_t + 1 - ( - η) F̃_t≲ηF̃_t +  F̃_tξ_t^(R)  + ξ_t^(F)   ≲ηF̃_t^2 + ηδ√(r) + ηρ + η^2
       ≲ηF̃_t^2 + ηρ



The proof of Proposition <ref> follow straightforwardly from Lemma <ref> and Lemma <ref>. 


	Using the assumption that F_t≲√(ρ) (Thus F_t^2 ≲ρ). Since we have showed that F̃_̃t̃ - F_t≲η (ρ + E_t), the proof of this proposition followings immediately from Lemma <ref>. 





 §.§ Proof of Proposition <ref>


We first prove the following technical lemma that characterizes how much the least singular value of a matrix changes when it got multiplied by matrices that are close to identity. 

Suppose Y_1∈ℝ^d× r and Σ is a PSD matrix in ℝ^r × r. For some η > 0, let 

    Y_2 =  ( + ηΣ)Y_1 ( - η Y_1^⊤  Y_1)

Then, we have:  

    σ_min (Y_2) ≥(1 + ησ_min(Σ) )(1 - ησ_min(Y_1)^2) σ_min(Y_1)




First let's consider the matrix  Y := Y_1 ( - η Y_1^⊤  Y_1). We have that that σ_min(  Y ) = (1 - ησ_min(Y_1)^2) σ_min(Y_1). Next, we bound the least singular value of ( + ηΣ)Y: 

    σ_min(( + ηΣ)Y)   ≥σ_min( + ηΣ)σ_min(Y)= (1+σ_min(Σ))σ_min(Y)

where we used the facts that σ_min(AB)≥σ_min(A)σ_min(B), and that for any symmetric PSD matrix B, σ_min( + B) = 1 + σ_min(B). 

Now we are ready to prove Proposition <ref>. Note that the least singular value of Z_t is closely related to the least singular value of R̃_t because F_t is close to 0. Using the machinery in the proof of Lemma <ref>, we can write R_t+1 as some transformation of R̃_t, and then use the lemma above to bound the least singular value of R_t+1 from below. 

Recall that in equation (<ref>) in the the proof of Lemma <ref>, we showed  that 

    - ( - ηR̃_tR̃_t^⊤ + η) R̃_tR_t+1^+   ≲ηδ√(r)  + ηρ + ηF_t

Together with R̃_t≤ 6 and R̃_t R_t + 1^+≤ 2 by Corollary <ref> and Lemma <ref> respectively, we have:

    - ( + η) R̃_t ( - ηR̃_t^⊤R̃_t )  R_t + 1^+   ≲ηδ√(r)  + ηρ + ηF_t + η^2

Denote by  ξ = - ( + η) R̃_t ( - ηR̃_t^⊤R̃_t )  R_t + 1^+, we can rewrite

    R_t + 1 - ( + η) R̃_t ( - ηR̃_t^⊤R̃_t )   = ξ R_t + 1

Without loss of generality, let us assume that σ_min(R̃_t) ≤1/1.9 √(κ).  By Lemma <ref> that

    σ_min(R_t + 1) ≥(1 + ησ_min ()) (1 - ησ_min(R̃_t)^2) σ_min(R̃_t) /1  +O( ηδ√(r)  +ηρ +ηF_t + η^2 )

Since F_t≲ηρ t by Proposition <ref>, we have:

    σ_min(R_t + 1)    ≥ (1 + η / κ) (1 - ησ_min(R̃_t)^2) σ_min(R̃_t) (1 - O ( ηδ√(r)  + ηρ + η^2ρ t ) )
       ≥(1 + η(1/3κ -  O ( δ√(r)  +  ρ + η(δρ) t ) )  )σ_min(R̃_t)  by σ_min(R̃_t) ≤1/1.9√(κ)
       ≥(1 +   η/4κ)σ_min(R̃_t)   by t ≤c/ηκρ

Since  σ_min(R̃_t ) ≥σ_min(R_t)( 1 - η/100 κ) by Proposition <ref>,
the conclusion follows.






 §.§ Proof of Proposition  <ref>

	Since U_t+1 = U_t - η(U_t), we first write down the error for step t+1:
	
    U_t+1U_t+1^⊤ - ^2
    		= (U_t - η(U_t)) (U_t^⊤ - η(U_t)^⊤) - ^2 
        = U_t U_t^⊤ - ^2 - 2η(U_t) U_t^⊤U_t U_t^⊤ - 
        + 2η^2 (U_tU_t^⊤ - )(U_t)(U_t)^⊤ + 4η^2 (U_t)U_t^⊤^2 
        - 4η^3 (U_t) U_t^⊤(U_t)(U_t)^⊤
    		+ η^4 (U_t)(U_t)^⊤^2

	Denote by
	
    Δ = U_tU_t^⊤ - E_tE_t^⊤ = E_t Z_t^⊤ + Z_t E_t^⊤ + E_tE_t^⊤.

	First we have that
	
    Δ≤ 2  E_t  Z_t  +  E_t _2^2  ≤ O(E_t).

	We first consider the degree one term of η in equation (<ref>).
	
		In the setting of this subsection, we have that:
		
    (U_t)U_t^⊤U_tU_t^⊤ - ≥ (1 - O(η)) (Z_tZ_t^⊤ - )Z_t^2
    		- O(δ) Z_tZ_t^⊤ - ^2  - O(√(dr)E_t).

	
	
		First, we have that
		
    (U_t) U_t^⊤U_tU_t^⊤ - ≥(U_t) U_t^⊤Z_tZ_t^⊤ -  - (U_t) U_t^⊤Δ

		Since U_t≤ O(1), we focus on the gradient of U_t.
		We divide (U_t) into the sum of three parts:
		
    Y_1 = 1/m∑_i=1^m A_iZ_tZ_t^⊤ -  A_i Z_t, 
        Y_2 = 1/m∑_i=1^m A_iZ_tZ_t^⊤ -  A_i E_t, 
        Y_3 = 1/m∑_i=1^m A_iΔ A_i U_t

		By Lemma <ref>, we have that:
		
    Y_1≤(Z_tZ_t^⊤ - ) Z_t + δ×Z_tZ_t^⊤ - Z_t≤ O(1),

		where Z_tZ_t^⊤ - ≤ O(√(r)) by Corollary <ref>. Since Y_1 has rank at most r, we get that Y_1≤ O(√(r)). For Y_2, we apply Lemma <ref> again:
		
    Y_2≤Z_tZ_t^⊤ - E_t + δ×Z_tZ_t^⊤ - E_t≤ O(E_t).

		By assumption E_t≤ 1 / d, hence Y_2≤ O(1).
		Finally we apply Lemma <ref> for Y_3:
		
    Y_3≤ΔU_t + δ× (2E_tZ_t^⊤ + E_tE_t^⊤) U_t≤ O(E_t)

		To summarize, we have shown that (U_t)≤ O(√(r)).
		Back to equation (<ref>), we have shown that:
		
    (U_t) U_t^⊤Δ≤ O(√(dr)) E_t.

						For the other part in equation (<ref>),
		
    (U_t)U_t^⊤Z_tZ_t^⊤ -    ≥M_t Z_tZ_t^⊤Z_tZ_t^⊤ - 
    			- M_t ΔZ_tZ_t^⊤ - 
       ≥M_t Z_tZ_t^⊤Z_tZ_t^⊤ -  - O(√(dr)) E_t,

		because M_t≤ O(1) from Lemma <ref>. Lastly, we separate out the Δ term in M_t as follows:
		
    M_t Z_tZ_t^⊤Z_tZ_t^⊤ - 
    			=    (Z_t)Z_tZ_tZ_t^⊤ -  + 
       1/m∑_i=1^m A_iE_tZ_t^⊤ + Z_tE_t^⊤A_iZ_tZ_t^⊤ (Z_tZ_t^⊤ - ) + 
       1/m∑_i=1^m A_iE_tE_t^⊤A_iZ_tZ_t^⊤ (Z_tZ_t^⊤ - )

		Since Z_t is a rank-r matrix, by Lemma <ref>, Equation (<ref>) is at least:
		
    ⟨ Z_t Z_t^⊤ -  , Z_t Z_t^⊤(Z_t Z_t^⊤ - )⟩ - δZ_t^2 Z_t Z_t^⊤- ^2.

		Equation (<ref>) is similarly bounded by O(√(dr)E_t) using Lemma <ref>, since E_tZ_t^⊤ also has rank at most r.
		Finally for equation (<ref>), while E_tE_t^⊤ may have rank d, we can still apply Lemma <ref> and afford to lose a factor of d, because E_tE_t^⊤≤ O(1 / d^2).
		To summarize, we have shown that:
		
    (U_t) U_t^⊤U_tU_t^⊤ - ≥(Z_t) Z_t^⊤Z_tZ_t^⊤ -  - O(√(dr)) E_t.

	The conclusion follows by combining the above equation with (<ref>).
	
	Next we work out the degree two term in equation (<ref>).	
	
		In the setting of this subsection, we have that both
		
    |(U_tU_t^⊤ - )(U_t)(U_t)^⊤|
    		 (U_t) U_t^⊤^2 
       ≲(Z_tZ_t^⊤ - ) Z_t^2 + δZ_tZ_t^⊤ - ^2 + r E_t

	
	
		The idea of the proof is similar to that for equation (<ref>) and (<ref>). First of all, we have
		
    U_tU_t^⊤ - (U_t) (U_t)^⊤≥Z_t Z_t - (U_t) (U_t)^⊤
    			- Δ×(U_t)(U_t)^⊤

		The second term is at most O(r) ×E_t from our proof.
		Next,
		
    Z_tZ_t^⊤ - (U_t) (U_t)^⊤
    			= Z_tZ_t^⊤ - M_t Z_tZ_t^⊤ M_t^⊤
    			+ Z_tZ_t^⊤ - M_t Δ M_t^⊤

		The second term is at most:
		
    Δ×M_t^⊤ (Z_tZ_t^⊤ - ) M_t≲ r E_t

		Lastly, we expand out M_t to obtain:
		
    Z_tZ_t^⊤ - M_t Z_tZ_t^⊤ M_t^⊤
    			= 1/m∑_i=1^m A_iU_tU_t^⊤ - Z_tZ_t^⊤ - A_iZ_tZ_t^⊤M_t^⊤
        = 1/m∑_i=1^m A_iZ_tZ_t^⊤ - A_iW + 1/m∑_i=1^m A_iΔA_iW

		where we denote by W = (Z_tZ_t^⊤ - )M_tZ_tZ_t^⊤.
		Clearly, the rank of W is at most r and it is not hard to see that the spectral norm of W is O(1).
		For the first part, we apply Lemma <ref> to obtain:
		
    |Z_tZ_t^⊤ - W| + δ×Z_tZ_t^⊤ - W≲(Z_tZ_t^⊤ - )Z_t^2 + Z_tZ_t^⊤ - ^2

		where we used that M_t≤ O(1) and δ≲ 1 / √(r).
		For the second part, we apply Lemma <ref> and Lemma <ref> together on Δ to obtain:
		
    |ΔW| + δ× (2E_tZ_t^⊤ + E_tE_t^⊤) W≤ O(r) ×E_t

		where we used the assumption that E_t≤ 1 /d.
		The proof for (U_t)U_t^⊤^2 is similar.
		The difference is that we will obtain W' = Z_tZ_t^⊤ M_t Z_tZ_t^⊤ instead.
		To bound Z_tZ_t - W', we use Lemma <ref> and Lemma <ref> to control M_t. The details are left to the readers.
	
		Finally we consider the degree three and four terms of η in equation (<ref>):
	
    4η^3 (U_t) U_t^⊤(U_t)(U_t)^⊤≤ O(η^3) ×(U_t)^2,

	because U_tU_t^⊤ - ≤ O(1) and (U_t)≤ O(1). For the gradient of U_t, we have already decomposed it to the sum of Y_1, Y_2 and Y_3.
	And our proof already implies that:
	
    Y_1≤ O(√(r)) ×Z_tZ_t^⊤ - , 
       Y_2, Y_3≤ O(√(d)) ×E_t

	Combining all results together, we get that:
	
    (U_t)^2 ≤ O(r) Z_tZ_t^⊤ - ^2 + O(d E_t^2)

	Hence equation (<ref>) is at most:
	
    O(r η^3) Z_tZ_t^⊤ - ^2 + O(E_t)

	Combining the above equation with Claim <ref> and <ref>, we have shown that:
	
    U_t+1U_t+1^⊤ - ^2
    		≤   U_tU_t^⊤ - ^2 - (η - O(η^2)) (Z_tZ_t^⊤ - )Z_t^2 + O(√(dr)) E_t
        + O(ηδ + η^2δ + rη^3) Z_tZ_t^⊤ - ^2

	Lastly we show that:
	
    (Z_tZ_t^⊤ - )Z_t^2 ≳1/κZ_tZ_t^⊤ - ^2.

	The conclusion follows since it is not hard to show that
	
    U_tU_t^⊤ - ^2 = Z_tZ_t^⊤ - ^2 ± O(√(dr)) E_t.

The rest of the proof is dedicated to equation (<ref>).
Denote by Z_t = U Σ V^⊤ its SVD.
Recall that X^⋆ = U^⋆Σ^⋆U^⋆^⊤.
We have

    (Z_t Z_t^⊤ - U^⋆) Z_t^2
         = (U Σ^2  - U^⋆Σ^⋆U^⋆^⊤ U) Σ^2 
       ≥σ_min(Σ) U Σ^2 - U^⋆Σ^⋆U^⋆^⊤ U^2 
       = σ_min(Σ) (Σ^2^2 + Σ^⋆ Y^2 - 2Σ^2Y^⊤Σ^⋆ Y),

where we denote by Y = U^⋆^⊤ U.
Next, we have

    Z_t Z_t^⊤ - X^⋆^2
      = Σ^2^2 + Σ^⋆^2 - 2Σ^2Y^⊤Σ^⋆ Y.

By assumption we have sin(Z_t, U^*) ≤1/3, which gives us σ_min(Y) ≥ 1/4.
Meanwhile, the spectral norm of Y is at most 1.
Based on the two facts, we will prove equation (<ref>) by showing:

    Σ^2^2 + Σ^⋆ Y^2 - 2Σ^2Y^⊤Σ^⋆Y≥ c (Σ^2^2 + Σ^⋆^2 - 2Σ^2Y^⊤Σ^⋆ Y).

where c ≤σ_min(Y)^2 / (1 + σ_min(Y)^2) (e.g. c = 1/17 suffices).
By Cauchy-Schwarz inequality,

    2(1-c) Σ^2Y^⊤Σ^⋆ Y≤ (1-c) Σ^2^2 + (1-c) Y^⊤Σ^⋆ Y^2.

And then

    Σ^⋆ Y^2 - c Σ^⋆^2 - (1-c) Y^⊤Σ^⋆ Y^2 
       = (1-c)×(Σ^⋆ YY^⊤Σ^⋆ ( - YY^⊤)) - c×(Σ^⋆^2 ( - YY^⊤)) 
       ≥((1-c) σ_min(Y)^2 - c) ×(Σ^⋆Σ^⋆ ( - YY^⊤)) ≥ 0,

where the last line is because - YY^⊤ is PSD since Y≤ 1.
By combining equation (<ref>) and (<ref>), we have obtained equation (<ref>).
















































By the assumption that σ_min(Z_t)^2 ≥ 1/(4 κ) = Ω(δ), we complete the proof of equation (<ref>).

 


§ MISSING PROOFS IN SECTION <REF>



Let us first consider the case when X is a rank-1 matrix.  Suppose X  = aa^⊤  with a = 1. We then have:

    1/m∑_i = 1^m ⟨ A_i, X ⟩  A_i 1_|⟨ A_i, X ⟩| ≤ R   = 1/m∑_i = 1^m ⟨ A_i, a a^⊤⟩ A_i 
       =  1/m∑_i = 1^m⟨ x_i, a ⟩^2  x_i x_i^⊤ 1_⟨ x_i, a ⟩^2 ≤ R^2

We define: 
 
    H(x_1, ⋯, x_m) := sup_u, v ∈ℝ^d,  u _2=   v_2 = 1| 1/m∑_i = 1^m ⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2 1_⟨ x_i,u  ⟩^2 ≤ R^2 - 2 ⟨ u, v ⟩^2  - 1|



It suffices to bound H because by definition, for every X  = aa^⊤ being a rank one matrix, with a = 1, we have that  

    1/m∑_i⟨ x_i, a ⟩^2  x_i x_i^⊤ 1_⟨ x_i, a ⟩^2 ≤ R^2 - 2X  - I≤  H(x_1, ⋯, x_m)

Let us further decompose H into two terms:

    H(x_1, ⋯, x_m)    ≤sup_u, v ∈ℝ^d,  u _2=   v_2 = 1| 1/m∑_i = 1^m⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2 1_⟨ x_i,u  ⟩^2 ≤ R^2  1_⟨ x_i,v  ⟩^2 ≤ R^2  - 2 ⟨ u, v ⟩^2 - 1 |
        + sup_u, v ∈ℝ^d,  u _2=   v_2 = 1| 1/m∑_i = 1^m ⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2   1_⟨ x_i,u  ⟩^2 ≤ R^2 1_⟨ x_i,v  ⟩^2 > R^2 |

 
 
Let us bound the two term separately. For the first term, for every unit vectors u, v, we define functions f_u, v: ℝ^d →ℝ as f_u, v(x) = ⟨  u, x ⟩⟨ v, x ⟩  1_⟨ x_i,u  ⟩^2 ≤ R^2 1_⟨ x_i,v  ⟩^2 ≤ R^2. We have that f_u, v(x) ≤ R^4. Thus by the symmetrization technique and the contraction principle(e.g.,  see Corollary 4.7 in <cit.>), we have:

    [ sup_u, v ∈ℝ^d,  u _2=   v_2 = 1| ∑_i  = 1^m (f_u, v(x_i)^2 - [f_u, v(x_i)^2] ) |] ≤ 8R^2[  sup_u, v ∈ℝ^d,  u _2=   v_2 = 1| ∑_i = 1^m _i f_u, v(x_i)| ]

 where {_i}_i = 1^m is a set of i.i.d.  Rademacher random variables.  We can further bound the right hand side of the inequality above by:
 
    [  sup_u, v ∈ℝ^d,  u _2=   v_2 = 1| ∑_i = 1^m _i f_u, v(x_i)| ]    ≤[  ∑_i = 1^m _i x_i x_i^⊤]

 A standard bound on the norm of Gaussian random matrices gives us: [  ∑_i = 1^m _i x_i x_i^⊤] ≲√(md). Therefore, we conclude that 
 
    [ sup_u, v ∈ℝ^d,  u _2=   v_2 = 1| ∑_i  = 1^m (f_u, v(x_i)^2 - [f_u, v(x_i)^2 ] ) |]  ≲ R^2 √(md)

 Now let us consider the expectation of f_u, v(x_i)^2, a direct calculation shows that 
 
    | [f_u, v(x_i)^2]  - 2⟨ u, v ⟩^2 - 1 |    =| [⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2 1_⟨ x_i,u  ⟩^2 ≤ R^2  1_⟨ x_i,v  ⟩^2 ≤ R^2 ] - [⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2  ]|
       ≤[⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2 1_⟨ x_i,u  ⟩^2 ≤ R^2  1_⟨ x_i,v  ⟩^2 > R^2 ]
    + [⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2 1_⟨ x_i,u  ⟩^2 > R^2 ]
       ≤ 2 [⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2 1_⟨ x_i,u  ⟩^2 > R^2  ]

An elementary calculation of Gaussian variables gives us: 

    [⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2 1_⟨ x_i,u  ⟩^2 > R^2  ] ≲ R^4 e^- R^2/2

Putting everything together, for R ≥ 1 we are able to bound the first term of equation (<ref>) by:

    [ sup_u, v ∈ℝ^d,  u _2=   v_2 = 1| 1/m∑_i = 1^m⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2 1_⟨ x_i,u  ⟩^2 ≤ R^2  1_⟨ x_i,v  ⟩^2 ≤ R^2  - 2 ⟨ u, v ⟩^2 - 1 |  ]≲ R^4 ( √(d/m)+ e^- R^2 /2)

Moreover, for every u, v ∈ℝ^d,  u =   v = 1 we know that f_u, v(x)  ≤ R^2, we can apply <cit.> to transform the bound above into a high probability bound. We have that for every s ∈ [0, 1], with probability at least 1 - e^-Ω(s^2 m / (dR^4)), 

    sup_u, v ∈ℝ^d,  u _2=   v_2 = 1| 1/m∑_i = 1^m⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2 1_⟨ x_i,u  ⟩^2 ≤ R^2  1_⟨ x_i,v  ⟩^2 ≤ R^2  - 2 ⟨ u, v ⟩^2 - 1 | ≲ R^4 ( √(d/m)+ e^- R^2 /2) + s


Picking s = (R^2 √(d)/√(m)log1/q) with R =  Θ(log(1/δ) ) we obtain that 

    sup_u, v ∈ℝ^d,  u _2=   v_2 = 1| 1/m∑_i = 1^m⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2 1_⟨ x_i,u  ⟩^2 ≤ R^2  1_⟨ x_i,v  ⟩^2 ≤ R^2  - 2 ⟨ u, v ⟩^2 - 1 | ≤δ

For the second term of equation (<ref>) , we have that 

    sup_u, v ∈ℝ^d,  u _2=   v_2 = 1| 1/m∑_i = 1^m ⟨ x_i, u ⟩^2 ⟨ x_i , v ⟩^2   1_⟨ x_i,u  ⟩^2 ≤ R^2 1_⟨ x_i,v  ⟩^2 > R^2 |
       ≤sup_ v ∈ℝ^d,  v_2 = 1| 1/m∑_i = 1^m R^2 ⟨ x_i , v ⟩^2   1_⟨ x_i,v  ⟩^2 > R^2 |

By <cit.>, we have that for every s > 0, with probability at 1 - e^-Ω(s √(d)):

    sup_ v ∈ℝ^d,  v_2 = 1| 1/m∑_i = 1^m R^2 ⟨ x_i , v ⟩^2   1_⟨ x_i,v  ⟩^2 > R^2 | ≲ R^2 s^2 ( d^2/m + d^2/m s^2 R^-2log^2 m/n)


Taking s =Ω( log1/q/√(d)), putting everything together we prove the Lemma for the case when =aa^⊤ is rank one. 

For  of general rank, the proof follows by decomposing  to a sum of rank one singular vectors and apply triangle inequality directly. 




§ RESTRICTED ISOMETRY PROPERTIES


In this section we list additional properties we need for the set of measurement matrices {A_i}_i=1^m.
Lemma <ref> follows from the definition of RIP matrices.
The rest three Lemmas are all direct implications of Lemma <ref>.

For every x ∈ℝ^d, y ∈ℝ^d' of norm at most 1,
we have: 

    1/m∑_i = 1^m ⟨_i , ⟩ x^⊤_i  y - x^⊤ y
    	   = 1/m∑_i = 1^m ⟨_i , ⟩⟨_i, x y^⊤^⊤⟩ - x^⊤ y 
       ≤⟨, xy^⊤^⊤⟩  + δ_F xy^⊤^⊤ - x^⊤ y 
       ≤δ_F _2

The first inequality uses Lemma <ref>.


The following Lemmas deal with matrices that may have rank bigger than r.
The idea is to decompose the matrix into a sum of rank one
matrices via SVD, and then apply Lemma <ref>.


	Let {A_i}_i=1^m be a family of matrices in ^d × d
	that satisfy (r, δ)-restricted isometry property.
	Then for any matrices X, Y ∈^d × d,
	where the rank of Y is at most r,
	we have:
	
    1/m∑_i=1^m A_iXA_iY 		- XY≤δXY




	Let X = U D V^⊤ be its SVD.
	We decompose D = ∑_i=1^d D_i where each D_i contains only the i-th	diagonal entry of D,
	and let X_i = U D_i V^⊤ for each i=1,…,d.
	Then we have:
	
    1/m∑_i=1^m A_iXA_iY   = ∑_j=1^d ( 1/m∑_i=1^m A_iX_jA_iY) 
       ≤∑_j=1^d (X_jY + δX_jY) = XY + δXY




	Let {A_i}_i=1^m be a family of matrices in ^d × d
	that satisfy (1, δ)-restricted isometric property.
	Then for any matrix X ∈^d × d and matrix R ∈^d × d', where d' can be any positive integer,
	we have:
	
    1/m∑_i = 1^m A_iX  A_i R - X R≤δX×R.

	The following variant is also true:
	
    1/m∑_i=1^m A_iX  U A_i R - U X R≤δX×U×R,

 	where U is any matrix in ^d × d.



	Let X = U D V^⊤ be its SVD.
	We define X_i and D_i the same as in the proof of Lemma <ref>,
	for each i=1,…,d.
	
	For every x ∈^d, y ∈^d' with norm at most one, we have:
	
    1/m∑_i=1^m A_iX  x^⊤ A_i R y
    			- x^⊤ y 
    
    		=   ∑_j=1^d ( 1/m∑_i=1^m _i_j_ix y^⊤^⊤)
    			- x^⊤ y 
    ≤   ∑_j=1^d ( _jx y^⊤^⊤ + δ_j)
    			- x^⊤ y
    		= δ.

	The variant can be proved by the same approach (details omitted).


Asymmetric sensing matrices.
Recall that when each A_i is asymmetric, we simply use (A_i + A_i^⊤) / 2 instead of A_i as our sensing matrix.
While {(A_i + A_i) / 2}_i=1^m may only ensure the restricted isometry property for symmetric matrices,
we have the same inequality when the matrix X in Lemma <ref> is symmetric, which is the case for all our applications of Lemma <ref>:
[More precisely, X corresponds to any one of U_tU_t^⊤ -, E_tE_t^⊤, Z_tZ_t^⊤, or their linear combinations.]

    1/m∑_i=1^m A_i + A_i^⊤/2X (A_i + A_i^⊤) R / 2 - XR≤δXR.

Since X is symmetric, A_iX = A_i^⊤X.
The above equation then follows by applying Lemma <ref> twice,
with {A_i}_i=1^m and {A_i^⊤}_i=1^m as sensing matrices respectively.

For the applications of Lemma <ref> in Equations (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>),
we note that either X or Y is symmetric in all applications.
Suppose that X is symmetric, then we have the following
when we use (A_i + A_i^⊤) / 2 as the i-th sensing matrix:

    1/m∑_i=1^m A_i + A_i^⊤/2XA_i + A_i^⊤/2Y - XY≤δXY + Y^⊤/2

It is straightforward to verify that our proof still holds using the above inequality instead. The details for left for the readers.


plain


