On the Crucial Role of Initialization for Matrix Factorization

## Abstract

Abstract This work revisits the classical low-rank matrix factorization problem and unveils the critical role of initialization in shaping convergence rates for such nonconvex and nonsmooth optimization. We introduce Nyström initialization, which significantly improves the global convergence of Scaled Gradient Descent (ScaledGD) in both symmetric and asymmetric matrix factorization tasks. Specifically, we prove that ScaledGD with Nyström initialization achieves quadratic convergence in cases where only linear rates were previously known.
Furthermore, we extend this initialization to low-rank adapters (LoRA) commonly used for finetuning foundation models.
Our approach, NoRA, i.e., LoRA with Nyström initialization, demonstrates superior performance across various downstream tasks and model scales, from 1B to 7B parameters, in large language and diffusion models.

## 1 Introduction

Compared with learning rates and descent directions, initialization has been a relatively overlooked aspect of optimization.
In the widely studied smooth optimization literature (Nesterov, ; Ghadimi and Lan, ), as long as a suitable (small) learning rate is chosen, most of optimization algorithms such as GD provably converge to a stationary point at the same rate, regardless of initialization.
This work goes beyond stationary points and highlights the crucial role of initialization for global optimality of Burer-Monteiro factorization (Burer and Monteiro, ) – the same algorithm can exhibit markedly different behaviors, such as linear vs. quadratic convergence, depending on initialization.

We consider matrix factorization as a canonical example, where the goal is to solve i) symmetric problems, $\min_{\mathbf{X}}\|\mathbf{X}\mathbf{X}^{\top}-\mathbf{A}\|_{\mathsf{F}}^{2}$; and ii) asymmetric ones, $\min_{\mathbf{X},\mathbf{Y}}\|\mathbf{X}\mathbf{Y}^{\top}-\mathbf{A}\|_{
\mathsf{F}}^{2}$. While these classical problems can be handled via various approaches, they are notoriously challenging for optimization, since they are nonconvex, nonsmooth (albeit differentiable), non-coercive (for asymmetric problems), and do not satisfy Polyak-Lojasiewicz (PL) condition (Chi et al., ). Let $\mathbf{A}\in\mathbb{R}^{m\times n}$ (or $\mathbf{A}\in\mathbb{R}^{m\times m}$) for asymmetric (symmetric) problems, $\mathbf{X}\in\mathbb{R}^{m\times r}$ and $\mathbf{Y}\in\mathbb{R}^{n\times r}$. Building on the relation of $\text{rank}(\mathbf{A})$ and $r$, we can categorize matrix factorization into three setups: exact-parametrized ($\text{rank}(\mathbf{A})\!=\!r$), over-parametrized ($\text{rank}(\mathbf{A})\!<\!r$), and under-parametrized ($\text{rank}(\mathbf{A})\!>\!r$).

**Table 1: Comparison of complexity for global optimality in (a)symmetric matrix factorization in various settings. Here, EP, OP, and UP are abbreviations for exact-, over- and under- parametrization. $\epsilon$ is the prescribed optimality error, and $\kappa$ denotes the condition number of $\mathbf{A}$. Note that our bounds for UP depict the complexity to near optima; see formal descriptions in Defs. 1 and 2. The “special” initialization in AltGD is still a small initialization, but with more careful designs that will be clear in Sec. 3.1. Works marked with * are designed for another related setting (hence the comparison may not be fair).**
| setting | alg. | ref. | init. | rate |  |
| --- | --- | --- | --- | --- | --- |
| Asymmetric | EP | GD | (Ye and Du, ) | small | ${\cal O}\big{(}\kappa^{4}+\kappa^{4}\log(1/\epsilon)\big{)}$ |
| AltGD | (Ward and Kolda, ) | special | ${\cal O}\big{(}\kappa^{2}\log(1/\epsilon)\big{)}$ |  |  |
| ScaledGD | (Jia et al., ) | small | ${\cal O}(\log(1/\epsilon))$ |  |  |
| ScaledGD | Theorem 3 | Nyström | ${\cal O}(1)$ |  |  |
| OP | AltGD | (Ward and Kolda, ) | special | ${\cal O}\big{(}\kappa^{2}\log(1/\epsilon)\big{)}$ |  |
| ScaledGD | Theorem 6 | Nyström | ${\cal O}(1)$ |  |  |
| UP | GD | (Du et al., ) | small | asymptotic |  |
| ScaledGD | Theorem 4 | Nyström | ${\cal O}(1)$ |  |  |
| Symmetric | EP | GD* | (Stöger and Soltanolkotabi, ) | small | ${\cal O}\big{(}\kappa^{8}+\kappa^{2}\log(1/\epsilon)\big{)}$ |
| ScaledGD | Theorem 1 | Nyström | ${\cal O}\big{(}\kappa^{3}+\log\log(1/\epsilon)\big{)}$ |  |  |
| OP | GD* | (Stöger and Soltanolkotabi, ) | small | ${\cal O}\big{(}\kappa^{8}+\kappa^{6}\log(\kappa/\epsilon)\big{)}$ |  |
| ScaledGD | Theorem 5 | Nyström | ${\cal O}\big{(}\kappa^{3}+\log\log(1/\epsilon)\big{)}$ |  |  |
| UP | ScaledGD | Theorem 2 | Nyström | ${\cal O}(1/\epsilon\cdot\log(1/\epsilon))$ |  |

The asymmetric problem ii) is thoroughly explored in the literature.
For the exact- and over- parametrized cases, global convergence has been established for GD, Alternating GD (AltGD), and ScaledGD (Du et al., ; Ye and Du, ; Ward and Kolda, ; Jia et al., ; Tong et al., ), where most of them admit a linear rate. Regarding under-parametrized settings, only asymptotic global convergence of GD is established in (Du et al., ) to the best of our knowledge.
Common to above algorithms is the small initialization with $\mathbf{X}_{0}\sim{\cal N}(0,\zeta_{x}^{2})$ and $\mathbf{Y}_{0}\sim{\cal N}(0,\zeta_{y}^{2})$ for some sufficiently small $\zeta_{x}^{2}$ and $\zeta_{y}^{2}$.
However, such initialization results in unfavorable performance both theoretically and empirically, partly because of the need of escaping from a saddle point $(\mathbf{0},\mathbf{0})$.(^1^11It is not hard to see that $(\mathbf{0},\mathbf{0})$ is a stationary point by computing its gradient. To show that it is a saddle point, assume without loss of generality that $\mathbf{A}\neq\mathbf{0}$ and $\mathbf{A}_{1,1}>0$. Let $\tilde{\mathbf{X}}$ and $\tilde{\mathbf{Y}}$ be zero matrices except for one entry, i.e., $[\tilde{\mathbf{X}}]_{1,1}=[\tilde{\mathbf{Y}}]_{1,1}=\delta$ for arbitrarily small $\delta>0$. It follows that that $\|\tilde{\mathbf{X}}\tilde{\mathbf{Y}}^{\top}-\mathbf{A}\|_{\mathsf{F}}^{2} \leq\|\mathbf{0}\mathbf{0}^{\top}-\mathbf{A}\|_{\mathsf{F}}^{2}\leq\|-\tilde{ \mathbf{X}}\tilde{\mathbf{Y}}^{\top}-\mathbf{A}\|_{\mathsf{F}}^{2}$. Hence $(\mathbf{0},\mathbf{0})$ is a saddle point.)

This work proposes Nyström initialization to effectively bypass the aforementioned saddle point. More importantly, it significantly enhances the global convergence rates when applied on top of ScaledGD. In the exact- and over-parametrized settings, Nyström initialization boosts ScaledGD to converge at a quadratic rate (i.e., ${\cal O}(\log\log(1/\epsilon))$) on symmetric problems and enables a one-step convergence for asymmetric problems.
For the more challenging case with under-parametrization, we prove that with our Nyström initialization, ScaledGD converges at a linear rate to the neighbor of a global optimum on symmetric problems, and then exhibits a sublinear rate to a more fine-grained neighboring area.
Overall, Nyström initialization enables us to improve existing rates in exact-, over-, and under-parametrized settings; see more detailed comparisons in Tab. 1.

Our results highlight that the convergence of ScaledGD is critically determined by the initialization. Taking symmetric and exact-parametrized problems as an example, our quadratic rate slows down to a linear one when adopting either small initialization or slightly perturbed Nyström initialization.

After demonstrating the theoretical merits of Nyström initialization, we further extend its applications to another scenario with Burer-Monteiro factorization, in the context of LoRA for finetuning deep neural networks (Hu et al., ). This is motivated by the fact that asymmetric matrix factorization is equivalent to LoRA applied on linear models with whitened data (Arora et al., ; Jiang et al., ), and is in line with several recent works that take insights from matrix factorization to improve LoRA (Zhang and Pilanci, ; Yaras et al., ). Compared with existing strategies for initializing LoRA (Büyükakyüz, ; Meng et al., ; Wang et al., ), our Nyström initialization for LoRA (abbreviated as NoRA) is more economical and aligns better with existing deployment pipelines. The effectiveness of NoRA is demonstrated on downstream tasks from various domains, through both diffusion and large language models (LLMs). In a nutshell, our contributions can be summarized as:

- ❖
Faster rates. Nyström initialization is provably beneficial to ScaledGD. For symmetric problems, it catalyzes not only the first quadratic rate in exact- and over- parameterized settings, but also a (sub)linear rate for under-parametrization where only asymptotic results were known. It also allows more remarkable improvement on asymmetric problems; see details in Tab. 1. Moreover, these improved rates are obtained through a unified analysis framework.
- ❖
Critical role of initialization. Our theoretical results convey an intriguing message for nonconvex (nonsmooth) optimization: the behaviors of the same algorithm, whether converging at a quadratic or linear rate, are critically determined by initialization.
- ❖
Practical implications. We further illustrate the power of Nyström initialization for finetuning diffusion and large language models (LLMs). The resultant approach, NoRA, effectively improves the performance of LoRA on several representative tasks.

Notation. Bold lowercase (capital) letters denote column vectors (matrices); $(\cdot)^{\top}$, $(\cdot)^{\dagger}$ and $\|\cdot\|_{\mathsf{F}}$ refer to transpose, pseudo inverse, and Frobenius norm of a matrix; $\|\cdot\|$ is the $\ell_{2}$ (spectrum) norm of a vector (matrix); $\sigma_{i}(\cdot)$ and $\lambda_{i}(\cdot)$ denote the $i$-th largest singular value and eigenvalue, respectively. A (scalar) Gaussian random variable is denoted as ${\cal N}(\mu,\sigma^{2})$, where $\mu$ is its mean and $\sigma^{2}$ is the variance. We also use $\mathbf{X}\sim{\cal N}(\mu,\sigma^{2})$ for simplicity, which means that $\mathbf{X}$ is a random matrix whose entries are i.i.d. Gaussian variables with specified mean and variance.

### 1.1 Related works

We only streamline results on the convergence of matrix factorization under the broad umbrella of optimization under forth order growth, where the objective is to minimize a forth-order polynomial with first order approahces. Other closely related topics, such as LoRA variants, can be found in Apdx. A.1.

Optimization under forth-order growth. Matrix factorization problems considered in this work are classical examples of forth-order growth functions. It involves a complex landscape characterized by nonconvexity, nonsmoothness, and the absence of PL condition. Similar to other works listed in Tab. 1, the goal of this work is to recap this classical problem and to unveil intriguing behaviors from an optimization perspective. Recent works have examined the convergence of several algorithms, such as GD, AltGD, and ScaledGD (Du et al., ; Ye and Du, ; Ward and Kolda, ; Jia et al., ; Jiang et al., ) in exact- and over- parametrized settings. Most of them admit linear convergence with different dependences on the condition number of the factorized matrix $\mathbf{A}$. GD for matrix square root problems is studied in (Jain et al., ).
Another closely related setting within forth-order growth is matrix sensing; see e.g., (Stöger and Soltanolkotabi, ; Zhang et al., ; Jin et al., ; Xiong et al., ; Xu et al., ). Linear rates are obtained for problems with exact- and over- parametrization, despite some of them demand early stopping. Similar to matrix factorization, not too much is known for under-parametrization.
There are other approaches to tackle general forth-order growth optimization. For example, relative smoothness is considered in (Lu et al., ); adaptive step sizes induced by fine-grained geometry are studied in (Davis et al., ). The work of (Dragomir and Nesterov, ) also copes with such problems but requires convexity of the objective.

## 2 The power of initialization for symmetric matrix factorization

### 2.1 Preliminaries

We start to examine the critical role of initialization on symmetric matrix factorization problems. Consider the following objective

$$ $\displaystyle\min_{\mathbf{X}\in\mathbb{R}^{m\times r}}\frac{1}{4}\|\mathbf{X} \mathbf{X}^{\top}-\mathbf{A}\|_{\mathsf{F}}^{2}.$ (1) $$

Within this section, we assume that $\mathbf{A}\in\mathbb{R}^{m\times m}$ is positive semidefinite (PSD), otherwise one can employ the asymmetric formulation as in later sections. Problem (1) also closely links with matrix sensing, particularly under a sufficient number of Gaussian measurements (Xiong et al., ). From an optimization perspective, problem (1) is nonconvex and has no global Lipschitz gradient. These undesirable properties pose challenges for analyzing the convergence of classical optimization approaches (Tu et al., ; Chi et al., ).

Notationally, let $r_{A}:=\text{rank}(\mathbf{A})$ and further denote the compact eigendecomposition as $\mathbf{A}=\mathbf{Q}\mathbf{\Sigma}\mathbf{Q}^{\top}$, where $\mathbf{Q}\in\mathbb{R}^{m\times r_{A}}$ and $\mathbf{\Sigma}\in\mathbb{R}^{r_{A}\times r_{A}}$. Since PSD matrices share the same eigen and singular values, we employ $\sigma_{i}(\cdot)$ to denote both in this section. Without loss of generality, we assume that the largest and smallest singular values are $\sigma_{1}(\mathbf{A})=1$ and $\sigma_{r_{A}}(\mathbf{A})=1/\kappa$ such that the condition number is $\kappa$.

ScaledGD as our optimizer. We investigate the power of initialization on ScaledGD (Tong et al., ),
a preconditioned version of GD; see detailed discussions in e.g., (Tong et al., ; Jia et al., ). Starting from $t=0$ with a learning rate $\eta>0$, the update of ScaledGD is given by

$$ $\displaystyle\mathbf{X}_{t+1}=\mathbf{X}_{t}-\eta\underbrace{(\mathbf{X}_{t} \mathbf{X}_{t}^{\top}-\mathbf{A})\mathbf{X}_{t}}_{\text{gradient}}\cdot \underbrace{(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})^{-1}}_{\text{preconditioner}}.$ (2) $$

The inversion of the $r\times r$ matrix $\mathbf{X}_{t}^{\top}\mathbf{X}_{t}$ is computationally feasible in the low-rank setting with $r\ll m$. Small initialization is widely adopted, i.e., $[\mathbf{X}_{0}]_{ij}\sim{\cal N}(0,\zeta^{2})$, where $\zeta$ is a sufficiently small positive number. Under such initialization, ScaledGD converges linearly for exact-parametrization ($r=r_{A}$),(^2^22This linear rate is indicated by our numerical results in Fig. 1 (a). While we are not aware of a direct proof for this observation, it is presumable that the analysis in the asymmetric and exact-parametrized setting (Jia et al., ) could be adapted to provide some guarantees.) yet less is known for under- and over-parametrization; see more in Tab. 1. Next, we show that a simple yet effective initialization can provoke faster convergence of ScaledGD.

### 2.2 Nyström initialization

To improve the convergence rates, it is essential to ensure that the initialization satisfies two conditions for exact- and under-parametrized problems(^3^33For the ease of presentation, the over-parametrized setting is considered in the appendix.): i) each column of $\mathbf{X}_{0}$ is in the column space of $\mathbf{A}$, and ii) $\mathbf{X}_{0}$ is full rank, i.e., $\text{rank}(\mathbf{X}_{0})=r$. The analytical rationale will be elucidated in the subsequent sections.
A straightforward means to meet these conditions is via Nyström sketch (Gittens and Mahoney, )

where $\mathbf{\Omega}\in\mathbb{R}^{m\times r}$ is a Gaussian random matrix. From this initialization, it is not difficult to see that condition i) is satisfied already. Our next lemma shows that the condition ii) holds w.h.p.

###### Lemma 1 (Initialization for exact- and under- parametrization) .

For some universal constant $\tau>0$, $\sigma_{r}(\mathbf{X}_{0})\geq\xi\tau(\sqrt{r_{A}}-\sqrt{r-1})\sigma_{r_{A}}(
\mathbf{A})$ is satisfied with high probability, i.e., $\text{rank}(\mathbf{X}_{0})=r$ w.h.p.

The detailed expression for this “high probability” in Lemma 1 can be found in Apdx. B.1.1.
Note that there will be a “w.h.p.” over the initialization in most of our results. This refers to that $\text{rank}(\mathbf{X}_{0})=r$ in needed for exact- and under-parametrized settings, and $\text{rank}(\mathbf{X}_{0})=r_{A}$ is needed when over-parametrized.

###### Lemma 1 (Initialization for exact- and under- parametrization) .

### 2.3 Nyström initialization in the exact-parametrized setting

We start with Nyström initialization for exact-parametrized problems, i.e., $r_{A}=r$. Our first result dives into the implicit regularization induced by the ScaledGD under the proposed initialization.

###### Lemma 2 .

If $\mathbf{X}_{0}$ is obtained by Nyström initialization (3) and $\text{rank}(\mathbf{X}_{0})=r$ is satisfied, ScaledGD in (2) ensures that for all $t\geq 0$

- •
every column of $\mathbf{X}_{t}$ is in the column space of $\mathbf{A}$, and $\mathbf{X}_{t}=\mathbf{Q}\mathbf{\Phi}_{t}$ for some $\mathbf{\Phi}_{t}\in\mathbb{R}^{r\times r}$; and,
- •
the smallest eigenvalue of $\mathbf{X}_{t}\mathbf{X}_{t}^{\top}$ satisfies that
$\displaystyle\sigma_{r}(\mathbf{X}_{t+1}\mathbf{X}_{t+1}^{\top})$
$\displaystyle\geq(1-\eta)^{2t+2}\sigma_{r}(\mathbf{X}_{0}\mathbf{X}_{0}^{\top}
)+(1-\eta)\sigma_{r}(\mathbf{A})-(1-\eta)^{2t+3}\sigma_{r}(\mathbf{A}).$

Lemma 2 implies the full rankness of $\mathbf{X}_{t}$ over the trajectory, i.e., $\text{rank}(\mathbf{X}_{t})=\text{rank}(\mathbf{\Phi}_{t})=r,\forall t$. This ensures an invertible preconditioner $\mathbf{X}_{t}^{\top}\mathbf{X}_{t}$. In other words, iteration (2) is well-defined. The most important implication of Lemma 2 is the alignment of $\mathbf{X}_{t}$ with the directions of eigenvectors of $\mathbf{A}$, that is, $\mathbf{X}_{t}=\mathbf{Q}\mathbf{\Phi}_{t}$. This can be equivalently understood as the elimination of the residual space, i.e., $(\mathbf{I}-\mathbf{Q}\mathbf{Q}^{\top})\mathbf{X}_{t}=\mathbf{0},\forall t$. While we will expand this discussion shortly, this alignment in directions enables us to establish a quadratic rate for ScaledGD.

###### Theorem 1 .

With Nyström initialization (3), ScaledGD in (2) exhibits a two-phase behavior w.h.p over the initialization.

- •
Phase 1 (linear convergence). Let $\eta={\cal O}(\frac{1}{\kappa^{3}\|\mathbf{A}\|_{\mathsf{F}}})$. After $T_{1}:={\cal O}(\kappa^{3}\sqrt{r}\log\kappa)$ iterations, ScaledGD ensures that $\|\mathbf{X}_{T_{1}}\mathbf{X}_{T_{1}}^{\top}-\mathbf{A}\|_{\mathsf{F}}\leq{
\cal O}(1/\kappa^{2})$; and,
- •
Phase 2 (quadratic convergence). After Phase I, ScaledGD converges quadratically with $\eta=0.5$. In particular, $\|\mathbf{X}_{T}\mathbf{X}_{T}^{\top}-\mathbf{A}\|_{\mathsf{F}}\leq\epsilon$ is achieved after $T={\cal O}\big{(}\log\log(\frac{1}{\kappa\epsilon})\big{)}$ iterations.

Theorem 1 establishes that global optimality of (1) is attained by ScaledGD within ${\cal O}(\kappa^{3}\sqrt{r}\log\kappa+\log\log\frac{1}{\kappa\epsilon})$ iterations. ScaledGD first converges to a local region satisfying $\|\mathbf{X}_{t}\mathbf{X}_{t}^{\top}-\mathbf{A}\|_{\mathsf{F}}\leq{\cal O}(
\frac{1}{\kappa^{2}})$ linearly, after which a quadratic rate can be granted. This is, to the best of our knowledge, the first quadratic rate for symmetric matrix factorization (1). Interestingly, it is achieved without requiring (exact) Hessian on a nonconvex and nonsmooth problem. Note that the high probability in Theorem 1 can be equivalently stated as the requirement on $\text{rank}(\mathbf{X}_{0})=r$.
A graphical illustration of this quadratic rate can be found in Fig. 1 (a) using synthetic data detailed in Apdx. E.1. It is observed that ScaledGD with Nyström initialization outperforms linearly converging algorithms such as GD and ScaledGD with small initialization.
Moreover, it is worth emphasizing that Theorem 1 has no requirement on the magnitude of Nyström initialization – it does not need $\xi$ in (3) to be small. Compared with a small initialization, i.e., $\mathbf{X}_{0}\approx\mathbf{0}$, this avoids escaping from the stationary point $\mathbf{0}$. The convergence of ScaledGD under various choices of $\xi$ can be found in (the solid lines of) Fig. 1(b).

The critical role of initialization.
As shown in Lemma 2, Nyström initialization aligns $\mathbf{X}_{t}$ to the directions of eigenvectors $\mathbf{Q}$, thereby eliminating the residual space, i.e., $(\mathbf{I}-\mathbf{Q}\mathbf{Q}^{\top})\mathbf{X}_{t}=\mathbf{0},\forall t$. This is in stark contrast with most of existing works (Du et al., ; Ye and Du, ; Jia et al., ), where small initialization only guarantees that $\|(\mathbf{I}-\mathbf{Q}\mathbf{Q}^{\top})\mathbf{X}_{t}\|_{\mathsf{F}}$ converges to 0 0 at a linear rate. By getting rid of the residual space via Nyström initialization, ScaledGD can achieve a quadratic rate. We graphically illustrate this point in Fig. 1 (b), where we perturb Nyström initialization slightly to inject noise into the residual space. Reflected in the dotted lines, even if the noise is so small such that the earlier convergence does not differ from Nyström initialization, only a linear convergence can be observed for the perturbed initialization.

Extensions to the case of over-parametrization. Nyström initialization is further extended to cope with over-parametrized case ($r>r_{A}$) in Apdx. B.4.
For this specific setup, we slightly modify ScaledGD by substituting the possibly non-invertible $(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})^{-1}$ in (2) with $(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})^{\dagger}$; see (26).
Unlike previous works (Xu et al., ; Zhang et al., ), our modification requires no damping parameters thanks to our Nyström initialization.
This leads to, as far as we know, the first quadratic rate for over-parametrized problems. Additional numerical experiments on over-parametrized problems are provided in Fig. 4 in appendix to validate the established quadratic rate.

Figure: Figure 1: Convergence of ScaledGD under Nyström initialization (optimality error vs. iteration) in different settings. (a) Comparison of GD, and ScaledGD with small / Nyström initialization (ours). (b) Solid lines show that our initialization is not sensitive to magnitude of $\xi$; and dotted lines illustrate that quadratic convergence cannot be obtained after perturbing the initialization, i.e., $\mathbf{X}_{0}=\mathbf{A}\mathbf{\Omega}+\mathbf{N}$, where $[\mathbf{N}]_{ij}~{}\sim{\cal N}(0,\xi_{n}^{2})$. (c) Comparison of ScaledGD under Nyström initialization with various $\eta$.
Refer to caption: x1.png

###### Lemma 2 .

###### Theorem 1 .

### 2.4 Nyström initialization in the under-parametrized setting

Next, we consider the under-parametrized case of (1), i.e., $r<r_{A}$. To the best of our knowledge, only asymptotic convergence is established for GD on such problems (Du et al., ). This is partially because that even the local PL condition is challenging to be verified. With Nyström initialization, we will show that ScaledGD converges under a slightly weaker criterion.

###### Definition 1 (Weak optimality) .

Matrix $\mathbf{X}\in\mathbb{R}^{m\times r}$ is weakly optimal to (1) if $\mathbf{X}^{\top}\mathbf{A}^{\dagger}\mathbf{X}-\mathbf{I}_{r}=\mathbf{0}$.

Our first result characterizes that all global optima are also weakly optimal. In other words, if weak optimality is ensured, this algorithm has a chance to reach a global optimum as well.

###### Lemma 3 .

All globally optimal solutions to (1) are also weakly optimal.

We then focus on the convergence of ScaledGD to weak optimality. In the case of under-parametrization, Nyström initialization also aligns $\mathbf{X}_{t}$ to the directions of eigenvectors of $\mathbf{A}$.

###### Lemma 4 .

If ScaledGD in (2) is equipped with Nyström initialization (3), one can write $\mathbf{X}_{t}=\mathbf{Q}\mathbf{\Phi}_{t},\forall t$ for some $\mathbf{\Phi}_{t}\in\mathbb{R}^{r_{A}\times r}$.

Lemma 4 shows that $(\mathbf{I}-\mathbf{Q}\mathbf{Q}^{\top})\mathbf{X}_{t}=\mathbf{0},\forall t$ also holds, namely, Nyström initialization eliminates the residual space. Building upon this, the convergence of ScaledGD can be established.

###### Theorem 2 .

The following holds w.h.p. for ScaledGD (2) with Nyström initialization (3):

i) (Linear convergence to neighborhood of weak optima). If one chooses a constant $\eta\leq 1$, ScaledGD ensures that $\|\mathbf{X}_{t}^{\top}\mathbf{A}^{\dagger}\mathbf{X}_{t}-\mathbf{I}_{r}\|_{
\mathsf{F}}\leq{\cal O}(\eta r)+\epsilon$ in ${\cal O}(\log\frac{1}{\epsilon})$ iterations; or,

ii) (Convergence to weak optima). Let $\eta={\cal O}(\epsilon/r)$, weak optimality is ensured by ScaledGD after ${\cal O}(\frac{r}{\epsilon}\log\frac{1}{\epsilon})$ iterations, i.e., $\|\mathbf{X}_{t}^{\top}\mathbf{A}^{\dagger}\mathbf{X}_{t}-\mathbf{I}_{r}\|_{
\mathsf{F}}\leq\epsilon$.

If one chooses a constant learning rate e.g, $\eta=0.1$, linear convergence can be established until reaching a neighboring area of a weakly optimal solution. The error $\|\mathbf{X}_{t}^{\top}\mathbf{A}^{\dagger}\mathbf{X}_{t}-\mathbf{I}_{r}\|_{
\mathsf{F}}={\cal O}(\eta r)$ is low, given that $r$ is typically small in practice. A graphical illustration of this linear rate can be found in Fig. 1 (c). On the other hand, if the learning rate is chosen according to the prescribed accuracy $\epsilon$, one can obtain a sublinear rate ${\cal O}(\frac{r}{\epsilon}\log\frac{1}{\epsilon})$ to exact weak optimality. These behaviors clearly indicate a step scheduling of learning rates (e.g., setting $\eta=0.1,0.01,\ldots$ every a few iterations) for both fast convergence and exact weak optimality in practice.
It is also worth mentioning that the convergence under both choices of $\eta$ has no dependence on $\kappa$. This aligns with the presumption in previous works (Tong et al., ; Jia et al., ) that ScaledGD performs well on ill-conditioned problems, providing the first rigorous justification for the under-parametrized setting.

Finally, we show that even in the worst case, ScaledGD guarantees that $\mathbf{X}_{t}$ converges to a point that is adequately close to a global solution, and the relative distance is sublinear in $r$.

###### Lemma 5 .

Let $\mathbf{Q}_{1}$ be the first $r$ column on $\mathbf{Q}$, and $\mathbf{\Sigma}_{1}$ be the top-left $r\times r$ sub-block of $\mathbf{\Sigma}$. Denote an optimal solution to (1) as $\mathbf{X}_{*}=\mathbf{Q}_{1}\mathbf{\Sigma}_{1}^{1/2}$. W.h.p. over the initialization, ScaledGD (2) with Nyström initialization (3) ensures

$$ $\displaystyle\lim_{t\rightarrow\infty}\|\mathbf{X}_{t}-\mathbf{X}_{*}\|_{ \mathsf{F}}\leq{\cal O}(r^{3/4}).$ $$

###### Definition 1 (Weak optimality) .

###### Lemma 3 .

###### Lemma 4 .

###### Theorem 2 .

###### Lemma 5 .

## 3 The power of initialization for asymmetric matrix factorization

### 3.1 Initialization and modified ScaledGD

This section demonstrates that the power of initialization is even more striking in solving asymmetric matrix factorization than symmetric ones. Given $\mathbf{A}\in\mathbb{R}^{m\times n}$, consider the following problem

$$ $\displaystyle\min_{\mathbf{X}\in\mathbb{R}^{m\times r},\mathbf{Y}\in\mathbb{R} ^{n\times r}}\frac{1}{2}\|\mathbf{X}\mathbf{Y}^{\top}-\mathbf{A}\|_{\mathsf{F} }^{2}.$ (4) $$

Denote $\text{rank}(\mathbf{A})=r_{A}$, and the compact SVD as $\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$, where $\mathbf{U}\in\mathbb{R}^{m\times r_{A}}$, $\mathbf{\Sigma}\in\mathbb{R}^{r_{A}\times r_{A}}$, and $\mathbf{V}\in\mathbb{R}^{n\times r_{A}}$. Similar to the previous section, we assume that $\sigma_{1}(\mathbf{A})=1$ and $\sigma_{r_{A}}(\mathbf{A})=1/\kappa$.

Nyström initialization. We adopt an asymmetric manner to initialize $\mathbf{X}_{0}$ and $\mathbf{Y}_{0}$ for (4), i.e.,

where $\mathbf{\Omega}$ is a Gaussian random matrix of $\mathbb{R}^{n\times r}$ with $[\mathbf{\Omega}]_{ij}\sim{\cal N}(0,\xi^{2}),\forall i,\forall j$. We can follow the same steps of Lemma 1 to show that $\mathbf{X}_{0}$ in (5) is rank $r$ w.h.p. in exact- and under-parametrized settings.
Moreover, there is no requirement on the magnitude of $\xi$, meaning that it is possible to start far from the saddle point $(\mathbf{0},\mathbf{0})$.
This asymmetry of $\mathbf{X}_{0}$ and $\mathbf{Y}_{0}$ in (5) is in contrast with small initialization which typically induces $\|\mathbf{X}_{0}\|_{\mathsf{F}}\approx\|\mathbf{Y}_{0}\|_{\mathsf{F}}$ (Du et al., ; Jia et al., ). The merits will become clear shortly.
Note that AltGD (Ward and Kolda, ) also adopts sketch at initialization, i.e., $\mathbf{X}_{0}={\cal O}(\mathbf{A}\mathbf{\Omega}_{1}/\sigma_{1}(\mathbf{A}))$ and $\mathbf{Y}_{0}={\cal O}(\sigma_{1}(\mathbf{A})\mathbf{\Omega}_{2})$, where $\mathbf{\Omega}_{1}$ and $\mathbf{\Omega}_{2}$ are Gaussian random matrices. Besides the requirement on small variance of $\mathbf{\Omega}_{1}$ and $\mathbf{\Omega}_{2}$ and the explicit need of $\sigma_{1}(\mathbf{A})$, this initialization cannot eliminate the residual space. Consequently, AltGD demands early stopping in exact- and over-parametrized problems, and little is known for under-parametrized case.

Modified ScaledGD. To adapt to the non-invertible $\mathbf{Y}_{0}^{\top}\mathbf{Y}_{0}=\mathbf{0}$ in Nyström initialization (5), we modify the first iteration of ScaledGD. More precisely, the updates are summarized below

$$ $\displaystyle\mathbf{X}_{1}=\mathbf{X}_{0},~{}\text{and}~{}\mathbf{X}_{t+1}= \mathbf{X}_{t}-\eta(\mathbf{X}_{t}\mathbf{Y}_{t}^{\top}-\mathbf{A})\mathbf{Y}_ {t}(\mathbf{Y}_{t}^{\top}\mathbf{Y}_{t})^{-1},\forall t\geq 1;$ (6a) $\displaystyle\mathbf{Y}_{t+1}=\mathbf{Y}_{t}-\eta(\mathbf{X}_{t}\mathbf{Y}_{t} ^{\top}-\mathbf{A})^{\top}\mathbf{X}_{t}(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})^ {-1},\forall t\geq 0.$ (6b) $$

### 3.2 Nyström initialization in the exact-parametrized setting

We start with the exact-parametrized case, i.e., $r_{A}=r$ in (4). The benefit of Nyström initialization (5) for iteration (6) is again the alignment of $\mathbf{X}_{t}$ and $\mathbf{Y}_{t}$ to the directions of singular vectors.

###### Lemma 6 .

The modified ScaledGD in (6) under Nyström initialization (5) guarantees that $\mathbf{X}_{t}=\mathbf{U}\mathbf{\Phi}_{t}$ and $\mathbf{Y}_{t}=\mathbf{V}\mathbf{\Psi}_{t}$, $\forall t\geq 0$ for some $\mathbf{\Phi}_{t}\in\mathbb{R}^{r\times r}$ and $\mathbf{\Psi}_{t}\in\mathbb{R}^{r\times r}$.

Similar to the symmetric problems, the implication of Lemma 6 is the elimination of residual space, i.e., $(\mathbf{I}-\mathbf{U}\mathbf{U}^{\top})\mathbf{X}_{t}=\mathbf{0}$ and $(\mathbf{I}-\mathbf{V}\mathbf{V}^{\top})\mathbf{Y}_{t}=\mathbf{0}$. This turns out to be even more beneficial for asymmetric problems, as it induces one-step convergence of ScaledGD.

###### Theorem 3 (One-step convergence) .

With $\eta=1$ and Nyström initialization (5), the modified ScaledGD (6) guarantees $\mathbf{X}_{1}\mathbf{Y}_{1}^{\top}=\mathbf{A}$ w.h.p. over the initialization. In other words, global convergence is achieved in one step.

Theorem 3 has a fundamental implication, that is, optimization is also a competitive tool for matrix factorization. This is because ScaledGD is the first optimization approach to share the same ${\cal O}(mnr)$ complexity as (compact) SVD given $r\leq\min\{m,n\}$. Comparing to symmetric matrix factorization (cf. Theorem 1), Theorem 3 suggests that problem (4) requires less iterations to be solved owing to the asymmetry of $\mathbf{X}_{0}$ and $\mathbf{Y}_{0}$ at initialization (5).
This partially agrees with results in (Xiong et al., ), which illustrate the benefit of asymmetry in Burer-Monterio factorization for matrix sensing.

Lastly, we present a result that may be of independent interest – the asymmetric and symmetric problems are interconnected under our Nyström initialization. This link is made clear in the proof of the following corollary (to Theorem 1), which states that ScaledGD admits quadratic convergence under different choices of step sizes.

###### Corollary 1 (Quadratic convergence) .

With Nyström initialization (5) and different choices of step sizes, modified ScaledGD in (6) have a similar behavior as Theorem 1, that is, a two-phase behavior w.h.p. over the initialization:

- •
Phase 1 (linear convergence). Let $\eta={\cal O}(\frac{1}{\kappa^{3}\|\mathbf{A}\|_{\mathsf{F}}})$. After $T_{1}:={\cal O}(\kappa^{3}\sqrt{r}\log\kappa)$ iterations, ScaledGD ensures that $\|\mathbf{X}_{T_{1}}\mathbf{Y}_{T_{1}}^{\top}-\mathbf{A}\|_{\mathsf{F}}\leq{
\cal O}(1/\kappa^{2})$.
- •
Phase 2 (quadratic convergence). After Phase I, ScaledGD converges quadratically with $\eta=0.5$. In particular, $\|\mathbf{X}_{T}\mathbf{Y}_{T}^{\top}-\mathbf{A}\|_{\mathsf{F}}\leq\epsilon$ is ensured after $T={\cal O}\big{(}\log\log(\frac{1}{\kappa\epsilon})\big{)}$ iterations.

Extensions to over-parametrization. One-step global convergence can also be established for over-parametrized asymmetric problems under Nyström initialization. More on this can be found in Apdx. C.3, where we provide the first convergence result on ScaledGD under such a setup.

###### Lemma 6 .

###### Theorem 3 (One-step convergence) .

###### Corollary 1 (Quadratic convergence) .

### 3.3 Nyström initialization in the under-parametrized setting

Lastly, we tackle the case of under-parametrization in the asymmetric problem (4), where $r_{A}>r$. Similar to the symmetric case in Sec.2.4, we consider a slightly weaker version of optimality.

###### Definition 2 (Generalized weak optimality) .

We say $(\mathbf{X},\mathbf{Y})$ is weakly optimal if $\mathbf{Y}^{\top}\mathbf{A}^{\dagger}\mathbf{X}-\mathbf{I}_{r}=\mathbf{0}$.

Generalized weak optimality is satisfied by any global optimum, which is proved in Lemma 13 in the appendix. With this preparation, we are ready to show that ScaledGD converges in a single step.

###### Theorem 4 .

If $\eta=1$, ScaledGD in (6) with Nyström initialization (5) ensures generalized weak optimality in one iteration w.h.p., i.e., $\mathbf{Y}_{1}^{\top}\mathbf{A}^{\dagger}\mathbf{X}_{1}-\mathbf{I}_{r}=\mathbf
{0}$.

The critical role of initialization. Through the theoretical analyses in the previous two sections, it is evident that the convergence of ScaledGD for matrix factorization is highly dependent on the initialization. Here is an intuitive, though not strictly rigorous, summary: Small initialization results in behaviors similar to first-order optimizers, i.e., linear convergence (Jia et al., ). In contrast, the proposed Nyström initialization catalyzes quadratic rates and even one-step convergence, resembling the optimization trajectory of second-order approaches such as Newton’s method (Nesterov, ).

###### Definition 2 (Generalized weak optimality) .

###### Theorem 4 .

## 4 NoRA: Nyström low rank adapters

Our theoretical results highlight the merits of suitable initialization for matrix factorization problems. One of the key insights is that the Burer-Monterio factorization benefits from good directions of $\mathbf{X}_{0}$ and $\mathbf{Y}_{0}$ at initialization; cf. Lemmas 2 and 6. We term this as directional alignment.
In this section, we extend the benefit of initialization to practical scenarios, showing that directional alignment is also beneficial for low-rank adapters (LoRA) in finetuning deep neural networks (Hu et al., ).

LoRA enhances parameter efficiency of finetuning by approximating the unknown parameter-change $\Delta\mathbf{W}\in\mathbb{R}^{m\times n}$ through Burer-Monterio factorization

$$ $\displaystyle\mathbf{W}_{0}+\Delta\mathbf{W}\approx\mathbf{W}_{0}+\mathbf{X} \mathbf{Y}^{\top}$ (7) $$

where $\mathbf{W}_{0}\in\mathbb{R}^{m\times n}$ is the pretrained weight (of a particular layer), and $\mathbf{X}\in\mathbb{R}^{m\times r}$ and $\mathbf{Y}\in\mathbb{R}^{n\times r}$ with $r\ll\min\{m,n\}$. A more detailed recap of LoRA can be found in Apdx. A.1. Directional alignment can be achieved if singular vectors for $\Delta\mathbf{W}$ are leveraged to initialize $\mathbf{X}_{0}$ and $\mathbf{Y}_{0}$. While $\Delta\mathbf{W}$ is unavailable a priori, empirical wisdom suggests that there exists a set of well-performed adapters that lie in the column (row) span of the pretrained weight matrix (Lingam et al., ), i.e., $\text{ColSpan}(\Delta\mathbf{W})\subseteq\text{ColSpan}(\mathbf{W}_{0})$ and $\text{RowSpan}(\Delta\mathbf{W})\subseteq\text{RowSpan}(\mathbf{W}_{0})$. In other words, $\mathbf{W}_{0}$ can be adopted as a suitable replacement of $\Delta\mathbf{W}$ for directional alignment.

Figure: Figure 2: Which singular values have the largest change after finetuning with LoRA of rank $r$? Orange: top-$r$ singular values; blue: other singular values. Note that here we only plot the first 64 singular values as others rarely have sufficiently large change.
Refer to caption: x4.png

Having $\text{ColSpan}(\mathbf{W}_{0})$ alone is insufficient for directional alignment, since it does not specify which directions are more crucial.
To answer this question, we examine the singular values that undergo the most significant change after LoRA finetuning on a few-shot learning task (Malladi et al., ). OPT-1.3B is chosen as the base model and LoRA is applied to its query and value matrices with $r=8$; more details can be found in Apdx. E.3.
For each LoRA layer, we count the indices of $r$ singular values that exhibit the largest changes after finetuning, and summarize their frequencies across all layers in Fig. 2.
It is observed that the top-$r$ singular values tend to have larger change, explaining the success of LoRA initialization approaches that aligns $\mathbf{X}_{0}$ with the directions corresponding to these singular values, such as PiSSA and OLoRA (Meng et al., ; Büyükakyüz, ).
However, across all tested datasets, a substantial portion of non-top-$r$ singular values also demonstrate significant variation, and the frequency is positively linked to the singular values.
In other words, the directions corresponding to larger singular values tend to be more important.
This is akin to the principle of Nyström initialization $\mathbf{X}_{0}=\mathbf{W}_{0}\mathbf{\Omega}$, evidenced by its spectrum, i.e., $\mathbb{E}[\mathbf{X}_{0}\mathbf{X}_{0}^{\top}]\propto\mathbf{W}_{0}\mathbf{W}
_{0}^{\top}$.

Building upon these observations, and considering the accelerated convergence with Nyström initialization in ScaledGD, we propose two novel variants of LoRA:

- •
Nyström LoRA (NoRA): This approach applies (5) directly on top of LoRA, that is, $\mathbf{X}_{0}=\mathbf{W}_{0}\mathbf{\Omega}$ and $\mathbf{Y}_{0}=\mathbf{0}$.
- •
Nyström preconditioned LoRA (NoRA+): This approach not only advances LoRA initialization with (5), but also leverages ScaledGD for optimization.

We note that ScaledGD has already been applied for LoRA training in (Zhang and Pilanci, ), which we refer to as LoRA-P (P for preconditioning). We will show that both LoRA and LoRA-P benefit significantly from Nyström initialization. Due to space limitation, we summarize NoRA and NoRA+ in Algs. 1 and 2, respectively in the appendix, with additional explanations in Apdx. A.3.

**Table 2: Test accuracy of NoRA and NoRA+ for few-shot learning with OPT-1.3B.**
| OPT-1.3B | SST-2 | WSC | BoolQ | CB | RTE | ReCoRD | MultiRC | SQuAD | avg ($\uparrow$) |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Prefix | 92.9±0.9 | 59.6±1.6 | 73.1±2.3 | 71.6±2.9 | 65.2±2.6 | 69.7±1.0 | 64.4±3.2 | 82.2±1.4 | 72.3 |
| LoRA | 93.1±0.2 | 59.1±2.0 | 70.6±5.2 | 72.6±3.7 | 69.1±4.7 | 70.8±1.0 | 68.0±1.4 | 81.9±1.8 | 73.2 |
| OLoRA | 92.7±0.5 | 60.0±2.3 | 70.9±3.1 | 80.3±2.7 | 69.7±1.0 | 71.3±1.2 | 66.7±0.9 | 80.0±1.4 | 74.0 |
| PiSSA | 92.7±0.6 | 60.6±3.7 | 70.4±0.7 | 78.0±7.2 | 70.4±2.8 | 70.9±1.2 | 67.9±2.1 | 82.1±0.4 | 74.1 |
| NoRA | 93.4±0.7 | 60.6±3.8 | 73.2±0.6 | 79.2±5.2 | 72.0±1.3 | 71.3±1.0 | 68.5±1.2 | 81.8±0.7 | 75.0 |
| NoRA+ | 93.2±0.5 | 61.2±0.6 | 72.9±1.3 | 79.5±5.8 | 72.4±3.6 | 71.5±0.9 | 68.4±1.2 | 82.0±0.9 | 75.1 |

Deployment efficiency. NoRA offers practical advantages over other initialization methods such as PiSSA and OLoRA. It not only bypasses the computationally expensive SVD or QR decomposition, but also avoids the need to modify to the pretrained weights. NoRA is thus an off-the-shelf solution to enhance LoRA without altering existing pipelines. We expand on this in Apdx. A.3.

**Table 3: Training loss of NoRA and NoRA+ with stable-diffusion.**
| loss($\downarrow$) | LoRA | LoRA-P | NoRA | NoRA+ |
| --- | --- | --- | --- | --- |
| avg | 0.092$\pm 0.012$ | 0.093$\pm 0.012$ | 0.084$\pm 0.017$ | 0.084$\pm 0.015$ |

## 5 Numerical results for NoRA

The efficiency of proposed NoRA and NoRA+ is demonstrated on large-scale finetuning tasks involving diffusion and LLMs. The experiments are conducted with PyTorch (Paszke et al., ) on NVIDIA H100 GPUs. Details on datasets and experimental procedures can be found in Apdx. E.

### 5.1 Few-shot learning with OPT-1.3B

Our evaluation starts with a few-shot learning task following (Malladi et al., ). The objective is to rapidly adapt a language model with a small training set. The datasets for this experiment are drawn from GLUE and SuperGLUE benchmarks (Wang et al., , ). Consistent with (Malladi et al., ), we randomly sample 1,000 data points for training and another 1,000 for testing.

We embrace OPT-1.3B as our base model (Zhang et al., ) and apply LoRA to the query and value matrices in the attention module. This aligns with common practice for models of this size. The rank of LoRA is set to 8, leading to approximately 1.5M trainable parameters, which is significantly less than the model size. We compare the proposed NoRA and NoRA+ with LoRA, prefix tuning (Li and Liang, ), OLoRA (Büyükakyüz, ), and PiSSA (Meng et al., ). Note that the latter two serve as alternative methods for initializing LoRA.

The performance of different algorithms is summarized in Tab. 2. It is evident that OLoRA, PiSSA, NoRA, and NoRA+ all outperform LoRA because their initialization strategies have provided more favorable directions for optimization. Among these initialization approaches, NoRA and NoRA+ have the best average accuracy, with absolute improvement over LoRA by 1.8 and 1.9, respectively.

### 5.2 Subject-driven image generation with stable-diffusion

Figure: Figure 3: Generated images from NoRA and NoRA+ with stable-diffusion.
Refer to caption: x7.png

**Table 4: Test accuracy of various algorithms for commonsense reasoning on LLaMA-7B. HS and WG are abbreviations for HellaSwag and WinoGrande, respectively.**
| LLaMA-7B | BoolQ | PIQA | SIQA | HS | WG | ARC-e | ARC-c | OBQA | avg ($\uparrow$) |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| LoRA | 66.42 | 80.03 | 77.84 | 82.88 | 81.85 | 79.92 | 63.40 | 77.20 | 76.19 |
| LoRA-P | 68.96 | 80.95 | 77.43 | 81.54 | 80.27 | 78.83 | 64.16 | 79.20 | 76.41 |
| NoRA | 68.20 | 80.79 | 78.40 | 85.09 | 80.27 | 79.17 | 62.80 | 78.80 | 76.69 |
| NoRA+ | 69.85 | 81.83 | 77.38 | 82.09 | 80.03 | 79.67 | 64.25 | 78.60 | 76.71 |

Next, we focus on subject-driven image generation (Ruiz et al., ). The goal of this task is to finetune a diffusion model with only a few user-specific images (typically less than $10$) so that the modal can generate the same object in various contexts.
The base model is selected as StableDiffusion v1.4 (Rombach et al., ) (0.98B parameters in total). We adhere to the default setting and finetune the U-Net with LoRA. The rank of LoRA is set as 4, amounting to 0.8M trainable parameters. The diffusion model is finetuned on a user-specific training set containing pictures of a dog labeled “a photo of $V_{\text{dog}}$,” with the aim to generate proper images under the prompt “a $V_{\text{dog}}$ eating nachos.”

To demonstrate the power of initialization, we compare NoRA and NoRA+ with LoRA and LoRA-P. The averaged training loss of considered approaches are summarized in Tab. 3. It can be seen that NoRA and NoRA+ have $9.6\%$ smaller training loss compared with LoRA and LoRA-P, demonstrating the benefits of directional alignment at initialization.
The generated images are listed in Fig. 3. Some of images generated by LoRA are not natural. For instance, the third one does not have a nice expression for nachos, and the tenth is not vivid. For LoRA-P, the dog in the third image is also not natural. NoRA and NoRA+, on the other hand, both generate high-fidelity pictures. However, there is a floating plate in the 8th image of NoRA+, but ensuring diffusion models to follow physical laws goes beyond the scope of this work. Additional results are provided in Apdx. E.5, where we finetune on images of a cat toy. The generated images from NoRA and NoRA+ have more lively facial details compare to those not using Nyström initialization.

### 5.3 Commonsense reasoning with LLaMA-7B and LLaMA2-7B

Our evaluation is further scaled to LLMs using LLaMA and LLaMA2-7B (Touvron et al., , ). We tackle commonsense reasoning tasks following the setup in (Hu et al., ). Training data are merged from 8 datasets listed in Tab. 4. The test sets remain separate for individual evaluation. These reasoning tasks are intended to push the model beyond pattern recognition, requiring commonsense and knowledge to make proper inferences. The rank of LoRA is chosen as 32.

The results on LLaMA-7B are summarized in Tab. 4. It is observed that NoRA improves the average accuracy by 0.5 over LoRA, while NoRA+ also surpasses LoRA-P. These results underscore the significance of initialization for optimizing LoRA. The numerical results on LLaMA2-7B are presented in Tab. 5. In this case, it is observed that LoRA is unstable, henceforth the results for LoRA are taken from (Liu et al., ). This instability is not observed in other approaches tested. In this experiment, the benefit of the Nyström initialization is particularly pronounced, as the absolute improvement is even greater compared to the results on LLaMA-7B.

**Table 5: Test accuracy of various algorithms for commonsense reasoning on LLaMA2-7B. The results marked with $\ddagger$ are taken from (Liu et al., ).**
| LLaMA2-7B | BoolQ | PIQA | SIQA | HS | WG | ARC-e | ARC-c | OBQA | avg ($\uparrow$) |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| LoRA^‡ | 69.8 | 79.9 | 79.5 | 83.6 | 82.6 | 79.8 | 64.7 | 81.0 | 77.6 |
| LoRA-P | 71.47 | 81.50 | 78.81 | 85.97 | 80.43 | 81.14 | 66.55 | 81.00 | 78.35 |
| NoRA | 71.16 | 83.08 | 79.53 | 85.90 | 81.85 | 80.64 | 66.13 | 81.80 | 78.76 |
| NoRA+ | 70.52 | 81.94 | 79.07 | 87.66 | 82.24 | 82.70 | 67.06 | 80.20 | 78.92 |

### 5.4 Math reasoning with Gemma-7B

Our last evaluation tackles mathematical reasoning. Gemma-7B (Gemma-team et al., ) is finetuned for 2 epochs on the MetaMathQA-100K dataset (Yu et al., ). For this task, LoRA rank is set as 32, leading to 100M trainable parameters. The performance is assessed on GSM8K (Cobbe et al., ).

The test accuracy of various approaches is summarized in Tab. 6. We also include PiSSA (Meng et al., ) into the comparison. Note that PiSSA uses LoRA rank as $64$ but is only finetuned for a single epoch. Despite this difference, the computational cost on backward passes is the same for PiSSA and NoRA. The results clearly show that NoRA (NoRA+) outperforms LoRA (LoRA-P), highlighting the effectiveness of our Nyström initialization.

**Table 6: Test accuracy of different algorithms for math reasoning tasks. The results marked with $\ddagger$ are taken from (Meng et al., ).**
| GSM8K | LoRA | PiSSA^‡ | NoRA | LoRA-P | NoRA+ |
| --- | --- | --- | --- | --- | --- |
| Gemma-7B | 76.72 | 77.94 | 78.62 | 77.03 | 78.47 |

## 6 Concluding remarks and discussions

This work characterizes how initialization can crucially determine the convergence behavior of the same optimization algorithm on matrix factorization problems. We prove that Nyström initialization can significantly improve the complexity bounds of ScaledGD under a wide spectrum of settings; see details in Tab. 1. One of the key improvements is that Nyström initialization enables a quadratic convergence for exact- and over-parametrized problems, whereas small initialization only guarantees a linear rate on ScaledGD. This performance gap calls for more careful investigation into the role of initialization in optimization. Additionally, the proposed Nyström initialization offers practical merits when applied on finetuning with LoRA, delivering deployment flexibility and promising numerical performance on large-scale problems with LLMs and diffusion models.

An alternative interpretation. Our theoretical results can also be interpreted as emphasizing the importance of identifying the correct directions when optimizing functions with fourth-order growth. In the matrix factorization and sensing literature, it is common for algorithms such as GD or ScaledGD to exhibit a two-phase behavior: an initial phase of selecting the correct direction (also known as the spectral phase), followed by a second phase of rapid (e.g., linear) local convergence. Our findings clearly indicate that the correct directions in the first phase are critical – not only for faster termination of this phase, more importantly, for exponentially impacting the convergence rate in the second phase of ScaledGD.

Future directions. Our findings present several avenues for further exploration. From a theoretical standpoint, this work focuses on the impact of initialization in canonical matrix factorization problems. We believe our results can be extended to more complex settings, such as matrix sensing and tensor factorization, which are part of our future research plans.
On the practical side, our work hints at the potential for further gains by leveraging priors embedded in pretrained weights within the context of Burer-Monteiro factorization with LoRA. Investigating how to better uncover and utilize this hidden information is another attractive direction for future research.

## Appendix A Missing details

### A.1 More on related work

Convergence of over-parametrized matrix factorization problems. Consider again the asymmetric problem as an example, i.e., $\min_{\mathbf{X},\mathbf{Y}}\|\mathbf{X}\mathbf{Y}^{\top}-\mathbf{A}\|^{2}$ with $\mathbf{A}\in\mathbb{R}^{m\times n}$, $\mathbf{X}\in\mathbb{R}^{m\times r}$ and $\mathbf{Y}\in\mathbb{R}^{n\times r}$. Over-parametrization refers to the case where $\text{rank}(\mathbf{A})\leq r$. The gradient flow on the extreme over-parametrized problems, where $r\geq\max\{m,n\}$, is studied in (Tarmoun et al., ).
There are also papers (Stöger and Soltanolkotabi, ; Zhang et al., ; Jin et al., ; Xiong et al., ) considering the matrix sensing problem, which partially relates to our problem when there are sufficient Gaussian measures. The work of (Arora et al., ) considers deeper problem (i.e., having more than $3$ layers) while assuming $\mathbf{A}$ is full rank. Our results on over-parametrization can be found in Apdx. B.4 and Apdx. C.3 for symmetric and asymmetric problems, respectively. The comparison of ScaledGD with other works on over-parametrized problems can be found in Tab. 1.

LoRA and parameter-efficient finetuning. LoRA (Hu et al., ) is a notable example of parameter-efficient finetuning (PEFT) (Dong et al., ) approaches. The goal of PEFT is to reduce the resource requirement for finetuning LLMs on downstream tasks. Other commonly adopted PEFT methods include, e.g., adapters (Houlsby et al., ), zeroth-order optimizers (Malladi et al., ; Zhang et al., , ), and prefix tuning (Li and Liang, ). There are also various efforts to further enhance LoRA via adaptivity (Zhang et al., ), chaining (Lialin et al., ; Xia et al., ), regularization (Li et al., ; Li and Giannakis, ), low-bit training (Dettmers et al., ; Li et al., ), modifications for long-sequences (Chen et al., ), weight decomposition (Liu et al., ), and combining with sparsity (Nikdan et al., ). Additionally, there are several approaches aiming at further reducing the number of trainable parameters in LoRA; examples include (Kopiczko et al., ; Lingam et al., ; Gao et al., ; Zhu et al., ; Hao et al., ; Bałazy et al., ). While originally designed for finetuning LLMs, LoRA also finds its applications in other domains, such as image generation (Gu et al., ) and continual learning (Smith et al., ).

LoRA initialization. When first proposed, LoRA initialization was largely overlooked.
The work of (Hayou et al., ) justifies that whether setting $\mathbf{X}_{0}$ or $\mathbf{Y}_{0}$ to be $\mathbf{0}$ affects performance from a stability perspective. Recent works (Büyükakyüz, ; Meng et al., ) observe a fundamental difference between initialization of LoRA and neural networks, emphasizing the availability of prior knowledge. These works experimentally demonstrate that pretrained model can serve as prior to guide the direction of adapters, and hence perform QR or SVD on the pretrained matrix and using (scaled) top-$r$ singular vectors for LoRA initialization. Follow-up study (Wang et al., ) exploits stability for further improvement. However, these initialization methods are computationally expensive and lack flexibility for deployment. The proposed NoRA initialization overcomes these limitations.

Nyström sketch. Nyström sketch has well-documented success in signal processing and machine learning for coping with large-scale matrices under memory constraints. It has been applied in various settings; see e.g., (Gittens and Mahoney, ; Tropp et al., ; Frangella et al., ). This work only employs this approach to ensure full-rankness required in certain settings, yet the properties for recovery is not explored. We believe that other sketches are also applicable once full rankness is ensured.

### A.2 LoRA for linear models as asymmetric matrix factorization

We argue that LoRA applied on linear models given a whitened dataset is equivalent to the asymmetric matrix factorization problem. The whitened dataset is widely adopted for theoretical analyses, and we refer to (Arora et al., ; Jiang et al., ; Yaras et al., ) for more details.

Assume that we have a pretrained (linear) model $\mathbf{W}_{0}\in\mathbb{R}^{m\times n}$. Applying LoRA on this layer with whitened data $\mathbf{B}$ is equivalent to solving the following problem

$$ $\displaystyle\frac{1}{2}\|(\mathbf{W}_{0}+\mathbf{X}\mathbf{Y}^{\top})-\mathbf {B}\|_{\mathsf{F}}^{2}.$ (8) $$

It is clearly that this problem (8) is the same as (4) by setting $\mathbf{A}=\mathbf{B}-\mathbf{W}_{0}$.

Unfortunately, existing works provide no theoretical support on the most widely adopted initialization approach for LoRA in practice – either $\mathbf{X}_{0}$ or $\mathbf{Y}_{0}$ is chosen as $\mathbf{0}$ to preserve $\mathbf{W}_{0}+\mathbf{X}_{0}\mathbf{Y}_{0}^{\top}=\mathbf{W}_{0}$. In this sense, our Nyström initialization in (5) is the first means of initialization that justifies one variable can be set to $\mathbf{0}$.

Additional similarities between LoRA and matrix factorization. LoRA and matrix factorization share similar mathematical properties. For example, they both have no spurious local minima (Du et al., ; Ge et al., ; Jang et al., ). There are also recent efforts using insights from matrix factorization to further improve LoRA; see e.g., (Yaras et al., ; Nikdan et al., ).

### A.3 More on NoRA and NoRA+

| Algorithm 1 NoRA for a specific LoRA layer<br>1:Initialize: $\xi$ – standard deviation of random matrix $\mathbf{\Omega}$<br>2:Set $\mathbf{X}_{0}$ and $\mathbf{Y}_{0}$ via Nyström initialization (5)<br>3:Standard training process | Algorithm 2 NoRA+ for a specific LoRA layer<br>1:Initialize: $\xi$ – standard deviation of random matrix $\mathbf{\Omega}$; $\lambda$ – numerical stability of matrix inversion<br>2:Set $\mathbf{X}_{0}$ and $\mathbf{Y}_{0}$ via Nyström initialization (5)<br>3:for $t=0,\dots,T-1$ do<br>4:     Get gradient $\mathbf{G}_{\mathbf{X}_{t}}$ and $\mathbf{G}_{\mathbf{Y}_{t}}$<br>5:     if $t>0$ then<br>6:         $\mathbf{G}_{\mathbf{X}_{t}}\leftarrow\mathbf{G}_{\mathbf{X}_{t}}(\mathbf{Y}_{t<br>}^{\top}\mathbf{Y}_{t}+\lambda\mathbf{I}_{r})^{-1}/\|(\mathbf{Y}_{t}^{\top}<br>\mathbf{Y}_{t}+\lambda\mathbf{I}_{r})^{-1}\|_{\mathsf{F}}$<br>7:     end if<br>8:     $\mathbf{G}_{\mathbf{Y}_{t}}\leftarrow\mathbf{G}_{\mathbf{Y}_{t}}(\mathbf{X}_{t<br>}^{\top}\mathbf{X}_{t}+\lambda\mathbf{I}_{r})^{-1}/\|(\mathbf{X}_{t}^{\top}<br>\mathbf{X}_{t}+\lambda\mathbf{I}_{r})^{-1}\|_{\mathsf{F}}$<br>9:     Optimizer update<br>10:end for |
| --- | --- |

As discussed in Sec. 4, LoRA can significantly benefit from the aligned directions at initialization. Besides the theoretical benefits of applying Nyström initialization on ScaledGD (NoRA+), Nyström initialization can also be used directly with Adam (or AdamW), i.e., NoRA. There are several reasons for this. First, directional alignment from initialization is beneficial to most optimizers. While our theoretical results focus on ScaledGD, we believe that the aligned directions also improve GD. Despite the improvement may be less significant as in ScaledGD, we conjecture that the linear term in (Ye and Du, , Theorem 1.1) can be removed with Nyström initialization, because it can be roughly understood as the price of searching for proper directions. In other words, the benefits of Nyström initialization extend to other optimizers as well. Second, Adam also affords an explanation of preconditioning, and the preconditioner for $\mathbf{X}_{t}$ is also closely related to $\mathbf{Y}_{t}$. In other words, Adam shares similarities with ScaledGD in (6). These two reasons prompt the proposed NoRA, as summarized in Alg. 1. For NoRA+ in Alg. 2, we modify the vanilla ScaledGD iterations in (6) with two add-ons. First, a small parameter $\lambda$ is introduced for numerical stability of matrix inversion. This is a standard practice for numerical optimizers such as Adam (Kingma and Ba, ; Loshchilov and Hutter, ). Second, the gradient is normalized by the Frobenius norm of its preconditioner. The reason is that an optimal $\lambda$ is difficult to tune as shown in (Zhang and Pilanci, ), where they use $\lambda$ from $10^{-6}$ to $100$. With this normalizer, we can set $\lambda=10^{-6}$ in all our experiments without any tuning. Moreover, this normalizer is useful to prevent the instability in earlier iterations due to the non-invertable $\mathbf{Y}_{0}=\mathbf{0}$.

Deployment efficiency of NoRA. One benefit of NoRA (as well as NoRA+) is that it can be deployed jointly with adapters trained with LoRA – and hence there is no need to modify the current pipeline for deployment. This is because both of NoRA and LoRA do not need to modify the pretrained parameters, and the finetuned model is just $\mathbf{W}_{0}+\mathbf{X}_{T}\mathbf{Y}_{T}^{\top}$, where $\mathbf{W}_{0}$ is the pretrained model, and $\mathbf{X}_{T}$ and $\mathbf{Y}_{T}$ are finetuned adapter weights. On the contrary, other initialization approaches such as PiSSA and OLoRA (Meng et al., ; Büyükakyüz, ) are less efficient for using jointly with LoRA at deployment because both approaches modify the pretrained weights, so that the finetuned model becomes $\widehat{\mathbf{W}}_{0}+\mathbf{X}_{T}\mathbf{Y}_{T}^{\top}$, where $\widehat{\mathbf{W}}_{0}=\mathbf{W}_{0}-\mathbf{X}_{0}\mathbf{Y}_{0}^{\top}$. The use of $\widehat{\mathbf{W}}_{0}$ comes from the fact that initialization in PiSSA and OLoRA does not satisfy $\mathbf{X}_{0}\mathbf{Y}_{0}^{\top}=\mathbf{0}$. Consequently, when deploying PiSSA jointly with LoRA, one needs to store both $\mathbf{W}_{0}$ (for LoRA) and $\widehat{\mathbf{W}}_{0}$ (for PiSSA), leading to reduced memory efficiency.

## Appendix B Missing proofs for symmetric settings

### B.1 Initialization of exact- and under-parametrized problems

#### B.1.1 Proof of Lemma 1

###### Proof.

Let the compact eigenvalue decomposition of $\mathbf{A}$ be $\mathbf{A}=\mathbf{Q}\mathbf{\Sigma}\mathbf{Q}^{\top}$, where $\mathbf{Q}\in\mathbb{R}^{m\times r_{A}}$ and $\mathbf{\Sigma}\in\mathbb{R}^{r_{A}\times r_{A}}$. We then have that

$$ $\displaystyle\mathbf{X}_{0}=(\mathbf{Q}\mathbf{\Sigma})(\mathbf{Q}^{\top} \mathbf{\Omega}).$ (9) $$

It is not hard to verify that the matrix $\mathbf{Q}^{\top}\mathbf{\Omega}\in\mathbb{R}^{r_{A}\times r}$ is also a Gaussian random matrix, where each entry follows ${\cal N}(0,\xi^{2})$. Applying Lemma 19 on $\mathbf{Q}^{\top}\mathbf{\Omega}$, it can be seen that

$$ $\displaystyle\mathbb{P}\Big{(}\frac{\sigma_{r}(\mathbf{Q}^{\top}\mathbf{\Omega })}{\xi}\leq\tau(\sqrt{r_{A}}-\sqrt{r-1})\Big{)}\leq(C_{1}\tau)^{r_{A}-r+1}+e^ {-C_{2}r_{A}}:=\delta$ $$

where $C_{1}$ and $C_{2}$ are universal constants independent of $r_{A}$ and $r$. This inequality shows that with probability at least $1-\delta$, $\sigma_{r}(\mathbf{Q}^{\top}\mathbf{\Omega})\geq\xi\tau(\sqrt{r_{A}}-\sqrt{r-1})$.

Note that inequality $\sigma_{\min}(\mathbf{C}\mathbf{D})\geq\sigma_{\min}(\mathbf{C})\sigma_{\min}(
\mathbf{D})$ holds given full column rank of $\mathbf{C}$; see Lemma 17. Applying it to (9), we have that

$$ $\displaystyle\sigma_{r}(\mathbf{X}_{0})$ $\displaystyle\geq\sigma_{r_{A}}(\mathbf{Q}\mathbf{\Sigma})\sigma_{r}(\mathbf{Q }^{\top}\mathbf{\Omega})=\sigma_{r_{A}}(\mathbf{A})\sigma_{r}(\mathbf{Q}^{\top }\mathbf{\Omega})$ $\displaystyle\stackrel{{\scriptstyle(a)}}{{\geq}}\xi\tau(\sqrt{r_{A}}-\sqrt{r- 1})\sigma_{r_{A}}(\mathbf{A})$ $$

where (a) holds with probability at least $1-\delta$.
∎

###### Proof.

### B.2 Missing proofs for the symmetric and exact-parametrized setting

In the exact-parametrized setting, it is convenient to define

$$ $\displaystyle\mathbf{B}_{t}:=\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}$ (10) $$

where $\mathbf{\Phi}_{t}\in\mathbb{R}^{r\times r}$ comes from Lemma 2, i.e., $\mathbf{X}_{t}=\mathbf{Q}\mathbf{\Phi}_{t}$. The notation $\mathbf{B}_{t}$ will be used frequently in this subsection. With the help of Lemma 2, $\mathbf{B}_{t}$ can be understood as the “core” part of $\mathbf{X}_{t}\mathbf{X}_{t}^{\top}$, because $\mathbf{X}_{t}\mathbf{X}_{t}^{\top}=\mathbf{Q}\mathbf{\Phi}_{t}\mathbf{\Phi}_{
t}^{\top}\mathbf{Q}^{\top}=\mathbf{Q}\mathbf{B}_{t}\mathbf{Q}^{\top}$. Once proving Lemma 2, it allows us to study dynamics using a simpler but equivalent notion $\|\mathbf{B}_{t}-\mathbf{\Sigma}\|_{\mathsf{F}}$, i.e.,

$$ $\displaystyle\|\mathbf{X}_{t}\mathbf{X}_{t}^{\top}-\mathbf{A}\|_{\mathsf{F}}= \|\mathbf{Q}(\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}-\mathbf{\Sigma})\mathbf {Q}^{\top}\|_{\mathsf{F}}=\|\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}-\mathbf{ \Sigma}\|_{\mathsf{F}}=\|\mathbf{B}_{t}-\mathbf{\Sigma}\|_{\mathsf{F}}.$ $$

#### B.2.1 Proof of Lemma 2

###### Proof.

The proof relies on $\mathbf{B}_{t}$ defined in (10).
We will prove this lemma by induction. Since $\mathbf{X}_{0}=\mathbf{A}\mathbf{\Omega}$ in Nyström initialization,
we have that $\mathbf{\Phi}_{0}=\mathbf{\Sigma}\mathbf{Q}^{\top}\mathbf{\Omega}$. Moreover, our base assumption $\sigma_{r}(\mathbf{B}_{0})>0$ is true because $\text{rank}(\mathbf{B}_{0})=\text{rank}(\mathbf{X}_{0}\mathbf{X}_{0}^{\top})=r$, which is the result of Lemma 1.

For induction, assume that $\mathbf{X}_{t}$ can be written as $\mathbf{X}_{t}=\mathbf{Q}\mathbf{\Phi}_{t}$ with a full rank $\mathbf{\Phi}_{t}\in\mathbb{R}^{r\times r}$ at iteration $t$. By the update (2), we have that

$$ $\begin{split}\mathbf{X}_{t+1}&=\mathbf{X}_{t}-\eta(\mathbf{X}_{t}\mathbf{X}_{t }^{\top}-\mathbf{A})\mathbf{X}_{t}(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})^{-1}\\ &=\mathbf{Q}\mathbf{\Phi}_{t}-\eta\mathbf{Q}(\mathbf{\Phi}_{t}\mathbf{\Phi}_{t }^{\top}-\mathbf{\Sigma})\mathbf{Q}^{\top}\mathbf{Q}\mathbf{\Phi}_{t}(\mathbf{ \Phi}_{t}^{\top}\mathbf{Q}^{\top}\mathbf{Q}\mathbf{\Phi}_{t})^{-1}\\ &\stackrel{{\scriptstyle(a)}}{{=}}\mathbf{Q}\bigg{(}\mathbf{\Phi}_{t}-\eta( \mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}-\mathbf{\Sigma})\mathbf{\Phi}_{t}( \mathbf{\Phi}_{t}^{\top}\mathbf{\Phi}_{t})^{-1}\bigg{)}\\ &\stackrel{{\scriptstyle(b)}}{{=}}\mathbf{Q}\underbrace{\bigg{(}(1-\eta) \mathbf{\Phi}_{t}+\eta\mathbf{\Sigma}\mathbf{\Phi}_{t}^{-\top}\bigg{)}}_{:= \mathbf{\Phi}_{t+1}},\end{split}$ (11) $$

where (a) uses $\mathbf{Q}^{\top}\mathbf{Q}=\mathbf{I}_{r}$; and (b) uses $\mathbf{\Phi}_{t}$ is full rank (hence invertible). Note that $\mathbf{Q}$ and $\mathbf{A}$ share the same column space. This proves the first claim i) of this lemma.

Next we show that the smallest eigenvalue of $\mathbf{B}_{t+1}$ is bounded away from 0 0, or equivalently, $\mathbf{\Phi}_{t+1}$ is full rank. To start with, we have that from the expression of $\mathbf{\Phi}_{t+1}$ in (11),

$$ $\begin{split}\mathbf{B}_{t+1}&=\mathbf{\Phi}_{t+1}\mathbf{\Phi}_{t+1}^{\top}=( 1-\eta)^{2}\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}+2\eta(1-\eta)\mathbf{ \Sigma}+\eta^{2}\mathbf{\Sigma}\mathbf{\Phi}_{t}^{-\top}\mathbf{\Phi}_{t}^{-1} \mathbf{\Sigma}\\ &=(1-\eta)^{2}\mathbf{B}_{t}+2\eta(1-\eta)\mathbf{\Sigma}+\eta^{2}\mathbf{ \Sigma}\mathbf{B}_{t}^{-1}\mathbf{\Sigma}.\end{split}$ (12) $$

Note that $\mathbf{B}_{t+1}$ is a PSD matrix by definition (hence the eigenvalues and singular values are the same).
To see the smallest eigenvalue of $\mathbf{B}_{t+1}$ is lower bounded, we will apply Lemma 15 on (12) twice, i.e.,

$$ $\begin{split}&~{}~{}~{}~{}~{}\sigma_{r}(\mathbf{B}_{t+1})\\ &\stackrel{{\scriptstyle(c)}}{{\geq}}2\eta(1-\eta)\sigma_{r}(\mathbf{\Sigma})+ \sigma_{r}\Big{(}(1-\eta)^{2}\mathbf{B}_{t}+\eta^{2}\mathbf{\Sigma}\mathbf{B}_ {t}^{-1}\mathbf{\Sigma}\Big{)}\\ &\stackrel{{\scriptstyle(d)}}{{\geq}}2\eta(1-\eta)\sigma_{r}(\mathbf{\Sigma})+ (1-\eta)^{2}\sigma_{r}\big{(}\mathbf{B}_{t}\big{)}\\ &\stackrel{{\scriptstyle(e)}}{{\geq}}(1-\eta)^{2t+2}\sigma_{r}(\mathbf{B}_{0}) +2\eta(1-\eta)\sigma_{r}(\mathbf{\Sigma})\frac{1-(1-\eta)^{2t+2}}{2\eta-\eta^{ 2}}\\ &\stackrel{{\scriptstyle(f)}}{{\geq}}(1-\eta)^{2t+2}\sigma_{r}(\mathbf{B}_{0}) +(1-\eta)\sigma_{r}(\mathbf{\Sigma})-(1-\eta)^{2t+3}\sigma_{r}(\mathbf{\Sigma} ),\end{split}$ (13) $$

where (c) and (d) are because of Lemma 15; (e) is by unrolling $\sigma_{r}(\mathbf{B}_{t})$ using (d); and (f) is by $\frac{2\eta}{2\eta-\eta^{2}}\geq 1$. Combining (11) and (13) concludes the induction.
∎

###### Proof.

#### B.2.2 Proof of Theorem 1

###### Proof.

The proof is by combining Lemmas 7 and 8.
∎

###### Lemma 7 (Phase I. Linear convergence to near optima) .

Let $\eta={\cal O}(\frac{1}{\kappa^{3}\|\mathbf{A}\|_{\mathsf{F}}})$. After ${\cal O}(\kappa^{3}\sqrt{r}\log\kappa)$ iterations, ScaledGD (2) with Nyström initialization (3) ensures that $\|\mathbf{X}_{t}\mathbf{X}_{t}^{\top}-\mathbf{A}\|_{\mathsf{F}}\leq{\cal O}(1/
\kappa^{2})$.

###### Proof.

Subtracting $\mathbf{\Sigma}$ from both sides of (12), we can obtain that

$$ $\displaystyle\mathbf{B}_{t+1}-\mathbf{\Sigma}=(1-\eta)^{2}(\mathbf{B}_{t}- \mathbf{\Sigma})-\eta^{2}\mathbf{\Sigma}+\eta^{2}\mathbf{\Sigma}\mathbf{B}_{t} ^{-1}\mathbf{\Sigma}.$ $$

This implies that

$$ $\displaystyle~{}~{}~{}~{}~{}\|\mathbf{B}_{t+1}-\mathbf{\Sigma}\|_{\mathsf{F}}$ $\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}(1-\eta)^{2}\|\mathbf{B}_{t}- \mathbf{\Sigma}\|_{\mathsf{F}}+\eta^{2}\|\mathbf{\Sigma}\|_{\mathsf{F}}+\eta^{ 2}\|\mathbf{\Sigma}\mathbf{B}_{t}^{-1}\|_{2}\|\mathbf{\Sigma}\|_{\mathsf{F}}$ $\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}(1-\eta)^{2}\|\mathbf{B}_{t}- \mathbf{\Sigma}\|_{\mathsf{F}}+\eta^{2}\|\mathbf{\Sigma}\|_{\mathsf{F}}+\eta^{ 2}\|\mathbf{\Sigma}\|_{2}\|\mathbf{B}_{t}^{-1}\|_{2}\|\mathbf{\Sigma}\|_{ \mathsf{F}}$ $\displaystyle\leq(1-\eta)\|\mathbf{B}_{t}-\mathbf{\Sigma}\|_{\mathsf{F}}+\eta^ {2}\|\mathbf{\Sigma}\|_{\mathsf{F}}+\eta^{2}\frac{\sigma_{1}(\mathbf{\Sigma}) \|\mathbf{\Sigma}\|_{\mathsf{F}}}{\sigma_{r}(\mathbf{B}_{t})}$ $$

where (a) is by $\|\mathbf{M}\mathbf{N}\|_{\mathsf{F}}\leq\|\mathbf{M}\|_{2}\|\mathbf{N}\|_{
\mathsf{F}}$; and (b) follows from the sub-multiplicity of $\|\cdot\|_{2}$.

By Lemma 2, if $\eta\leq 2/3$ and there exists $T_{1}$ such that $\sigma_{r}(\mathbf{B}_{T_{1}})\geq\sigma_{r}(\Sigma)/3$, then it holds that $\sigma_{r}(\mathbf{B}_{t})\geq\sigma_{r}(\mathbf{\Sigma})/3,\forall t\geq T_{1}$.
According to Lemma 1, we can choose $\xi$ in (3) sufficiently large such that $\sigma_{r}(\mathbf{B}_{0})\geq\sigma_{r}(\mathbf{\Sigma})/3$, i.e., $T_{1}=0$.
Alternatively, to avoid such a requirement on $\xi$, we can simply choose a constant step size, e.g., $\eta=0.5$, and run a constant number of steps, $T_{1}={\cal O}(1/\eta)$, to ensure $\sigma_{r}(\mathbf{B}_{T_{1}})\geq\sigma_{r}(\Sigma)/3$; see Lemma 2.
For simplicity of the results, our proof below goes with the first method, i.e., $T_{1}=0$.

$$ $\displaystyle~{}~{}~{}~{}~{}\|\mathbf{B}_{t+1}-\mathbf{\Sigma}\|_{\mathsf{F}}$ $\displaystyle\leq(1-\eta)\|\mathbf{B}_{t}-\mathbf{\Sigma}\|_{\mathsf{F}}+\eta^ {2}\|\mathbf{\Sigma}\|_{\mathsf{F}}+\eta^{2}\frac{\sigma_{1}(\mathbf{\Sigma}) \|\mathbf{\Sigma}\|_{\mathsf{F}}}{\sigma_{r}(\mathbf{B}_{t})}$ $\displaystyle\leq(1-\eta)\|\mathbf{B}_{t}-\mathbf{\Sigma}\|_{\mathsf{F}}+\eta^ {2}\|\mathbf{\Sigma}\|_{\mathsf{F}}+3\eta^{2}\frac{\sigma_{1}(\mathbf{\Sigma}) \|\mathbf{\Sigma}\|_{\mathsf{F}}}{\sigma_{r}(\mathbf{\Sigma})}$ $\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}\eta\|\mathbf{\Sigma}\|_{ \mathsf{F}}+3\eta\kappa\|\mathbf{\Sigma}\|_{\mathsf{F}}+(1-\eta)^{t+1-T_{1}}\| \mathbf{B}_{T_{1}}-\mathbf{\Sigma}\|_{\mathsf{F}}$ $\displaystyle=\eta\|\mathbf{A}\|_{\mathsf{F}}+3\eta\kappa\|\mathbf{A}\|_{ \mathsf{F}}+(1-\eta)^{t+1-T_{1}}\|\mathbf{B}_{T_{1}}-\mathbf{\Sigma}\|_{ \mathsf{F}}$ $$

where (c) is by Lemma 14. From this inequality it is not difficult to see that once $\eta={\cal O}(\frac{1}{\kappa^{3}\|\mathbf{A}\|_{\mathsf{F}}})$, one will have $\|\mathbf{B}_{t+1}-\mathbf{\Sigma}\|_{\mathsf{F}}\leq{\cal O}(1/\kappa^{2})$ within the stated iterations.
∎

###### Lemma 8 (Phase II. Quadratic convergence to global optima) .

If we choose $\eta=0.5$ and suppose that after $T_{2}$ iterations, $\sigma_{r}(\mathbf{B}_{T_{2}})\geq\sigma_{r}(\mathbf{\Sigma})/3$ and $\|\mathbf{B}_{T_{2}}-\mathbf{\Sigma}\|_{\mathsf{F}}\leq 2/(3\kappa^{2})$ are satisfied, ScaledGD then ensures that for any $t\geq T_{2}$,

$$ $\displaystyle\|\mathbf{X}_{t+1}\mathbf{X}_{t+1}^{\top}-\mathbf{A}\|_{\mathsf{F }}=\|\mathbf{B}_{t+1}-\mathbf{\Sigma}_{r}\|_{\mathsf{F}}\leq\frac{4}{3\kappa^{ 2}}\frac{1}{2^{2^{t+1}}}.$ $$

###### Proof.

Let $\mathbf{C}_{t}=\mathbf{\Sigma}^{-1}\mathbf{B}_{t}$. We can rewrite (12) as

$$ $\displaystyle\mathbf{C}_{t+1}=(1-\eta)^{2}\mathbf{C}_{t}+2\eta(1-\eta)\mathbf{ I}_{r}+\eta^{2}\mathbf{C}_{t}^{-1}.$ $$

Subtracting $\mathbf{I}_{r}$ and rearranging it, we arrive at

$$ $\displaystyle\mathbf{C}_{t+1}-\mathbf{I}_{r}=(1-2\eta)(\mathbf{C}_{t}-\mathbf{ I}_{r})+\eta^{2}\mathbf{C}_{t}^{-1}(\mathbf{C}_{t}-\mathbf{I}_{r})^{2}.$ $$

By choosing $\eta=0.5$, we have that

$$ $\displaystyle\mathbf{C}_{t+1}-\mathbf{I}_{r}=\frac{1}{4}\mathbf{C}_{t}^{-1}( \mathbf{C}_{t}-\mathbf{I}_{r})^{2}.$ $$

Multiplying both sides with $\mathbf{\Sigma}$, we have that

$$ $\displaystyle\mathbf{B}_{t+1}-\mathbf{\Sigma}$ $\displaystyle=\frac{1}{4}\mathbf{\Sigma}\mathbf{B}_{t}^{-1}\mathbf{\Sigma}( \mathbf{C}_{t}-\mathbf{I}_{r})(\mathbf{C}_{t}-\mathbf{I}_{r})$ $\displaystyle=\frac{1}{4}\mathbf{\Sigma}\mathbf{B}_{t}^{-1}(\mathbf{B}_{t}- \mathbf{\Sigma})\mathbf{\Sigma}^{-1}(\mathbf{B}_{t}-\mathbf{\Sigma}).$ $$

This implies that

$$ $\displaystyle\|\mathbf{B}_{t+1}-\mathbf{\Sigma}\|_{\mathsf{F}}$ $\displaystyle\leq\frac{1}{4}\|\mathbf{\Sigma}\|_{2}\|\mathbf{B}_{t}^{-1}\|_{2} \|\mathbf{B}_{t}-\mathbf{\Sigma}\|_{\mathsf{F}}\|\mathbf{\Sigma}^{-1}\|_{2}\| \mathbf{B}_{t}-\mathbf{\Sigma}\|_{\mathsf{F}}$ $\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\frac{3}{4}\frac{\sigma_{1}( \mathbf{\Sigma})}{\sigma_{r}^{2}(\mathbf{\Sigma})}\|\mathbf{B}_{t}-\mathbf{ \Sigma}\|_{\mathsf{F}}^{2}\stackrel{{\scriptstyle(b)}}{{=}}\frac{3\kappa^{2}}{ 4}\|\mathbf{B}_{t}-\mathbf{\Sigma}\|_{\mathsf{F}}^{2}$ $$

where (a) is by Lemma 2, i.e., once $\sigma_{r}(\mathbf{B}_{T_{2}})\geq\sigma_{r}(\mathbf{\Sigma})/3$, then $\sigma_{r}(\mathbf{B}_{t})\geq\sigma_{r}(\mathbf{\Sigma})/3$ holds for all $t\geq T_{2}$; and (b) is by $\sigma_{1}(\mathbf{\Sigma})=1$ and $\sigma_{r}(\mathbf{\Sigma})=1/\kappa$.

Finally, applying Lemma 16, it can be seen that a quadratic rate can be established long as $\|\mathbf{B}_{T_{2}}-\mathbf{\Sigma}\|_{\mathsf{F}}\leq\frac{2}{3\kappa^{2}}$, and this condition is satisfied from Lemma 7.
∎

###### Proof.

###### Lemma 7 (Phase I. Linear convergence to near optima) .

###### Proof.

###### Lemma 8 (Phase II. Quadratic convergence to global optima) .

###### Proof.

### B.3 Missing proofs for the symmetric and under-parametrized setting

We start with some notation that would be helpful for this subsection. Let the compact eigenvalue decomposition of $\mathbf{A}=\mathbf{Q}\mathbf{\Sigma}\mathbf{Q}^{\top}$, where $\mathbf{Q}\in\mathbb{R}^{m\times r_{A}}$, and $\mathbf{\Sigma}\in\mathbb{R}^{r_{A}\times r_{A}}$.

In Lemma 4, we will prove that $\mathbf{X}_{t}=\mathbf{Q}\mathbf{\Phi}_{t}$ always holds if we employ Nyström initialization and ScaledGD in (2), where $\mathbf{\Phi}_{t}\in\mathbb{R}^{r_{A}\times r}$. We also denote $\mathbf{\Theta}_{t}:=\mathbf{\Phi}_{t}(\mathbf{\Phi}_{t}^{\top}\mathbf{\Phi}_{
t})^{-1}$, where the invertibility of $(\mathbf{\Phi}_{t}^{\top}\mathbf{\Phi}_{t})$ will become clear in the proof.

Lastly, let $\mathbf{B}_{t}:=\mathbf{\Phi}_{t}^{\top}\mathbf{\Sigma}^{-1}\mathbf{\Phi}_{t}$. Note that $\mathbf{B}_{t}\in\mathbb{R}^{r\times r}$ and $\mathbf{B}_{t}=\mathbf{X}_{t}^{\top}\mathbf{A}^{\dagger}\mathbf{X}_{t}$.

#### B.3.1 Proof of Lemma 3

###### Proof.

We start with rewriting $\mathbf{A}$,

$$ $\displaystyle\mathbf{A}=[\mathbf{Q}_{1},\mathbf{Q}_{2}]\begin{bmatrix}\mathbf{ \Sigma}_{1}&\mathbf{0}\\ \mathbf{0}&\mathbf{\Sigma}_{2}\end{bmatrix}\begin{bmatrix}\mathbf{Q}_{1}^{\top }\\ \mathbf{Q}_{2}^{\top}\end{bmatrix}=\mathbf{Q}_{1}\mathbf{\Sigma}_{1}\mathbf{Q} _{1}^{\top}+\mathbf{Q}_{2}\mathbf{\Sigma}_{2}\mathbf{Q}_{2}^{\top}$ (14) $$

where $\mathbf{Q}_{1}\in\mathbb{R}^{m\times r}$ and $\mathbf{Q}_{2}\in\mathbb{R}^{m\times(r_{A}-r)}$ are the first $r$ and other columns of $\mathbf{Q}$, respectively; and $\mathbf{\Sigma}_{1}\in\mathbb{R}^{r\times r}$ and $\mathbf{\Sigma}_{2}\in\mathbb{R}^{(r-r_{A})\times(r-r_{A})}$ are diagonal matrices formed by the first $r$ and the rest diagonal entries of $\mathbf{\Sigma}$.

It is not difficult to see that the optimal solution of (1) is $\mathbf{X}_{*}=\mathbf{Q}_{1}\mathbf{\Sigma}_{1}^{1/2}\mathbf{U}^{\top}$, where $\mathbf{U}\in\mathbb{R}^{r\times r}$ is any unitary matrix that accounts for rotation. Note that the pseudo-inverse of $\mathbf{A}$ can be written as $\mathbf{A}^{\dagger}=\mathbf{Q}\mathbf{\Sigma}^{-1}\mathbf{Q}^{\top}$. Plugging $\mathbf{X}_{*}$ into the definition of weak optimality, we arrive at

$$ $\displaystyle\mathbf{X}_{*}^{\top}\mathbf{A}^{\dagger}\mathbf{X}_{*}$ $\displaystyle=\mathbf{U}\Sigma_{1}^{1/2}\mathbf{Q}_{1}^{\top}(\mathbf{Q}_{1} \mathbf{\Sigma}_{1}^{-1}\mathbf{Q}_{1}^{\top}+\mathbf{Q}_{2}\mathbf{\Sigma}_{2 }^{-1}\mathbf{Q}_{2}^{\top})\mathbf{Q}_{1}\mathbf{\Sigma}_{1}^{1/2}\mathbf{U}^ {\top}\stackrel{{\scriptstyle(a)}}{{=}}\mathbf{I}_{r}$ $$

where in (a) we use the facts $\mathbf{Q}_{1}^{\top}\mathbf{Q}_{1}=\mathbf{I}_{r}$ and $\mathbf{Q}_{1}^{\top}\mathbf{Q}_{2}=\mathbf{0}_{r\times(r_{A}-r)}$. This concludes the proof.
∎

###### Proof.

#### B.3.2 Proof of Lemma 4

###### Proof.

The proof is based on induction. First we have that $\mathbf{X}_{0}=\mathbf{A}\mathbf{\Omega}=\mathbf{Q}\mathbf{\Sigma}\mathbf{Q}^{
\top}\mathbf{\Omega}$. It is clear that $\mathbf{\Phi}_{0}=\mathbf{\Sigma}\mathbf{Q}^{\top}\mathbf{\Omega}$. Now suppose that one can write $\mathbf{X}_{t}=\mathbf{Q}\mathbf{\Phi}_{t}$, following the update (2), it is not hard to see that

$$ $\begin{split}\mathbf{\Phi}_{t+1}&=\mathbf{\Phi}_{t}-\eta\big{(}\mathbf{\Phi}_{ t}\mathbf{\Phi}_{t}^{\top}-\mathbf{\Sigma}\big{)}\mathbf{\Phi}_{t}(\mathbf{ \Phi}_{t}^{\top}\mathbf{\Phi}_{t})^{-1}\\ &=(1-\eta)\mathbf{\Phi}_{t}+\eta\mathbf{\Sigma}\underbrace{\mathbf{\Phi}_{t}( \mathbf{\Phi}_{t}^{\top}\mathbf{\Phi}_{t})^{-1}}_{:=\mathbf{\Theta}_{t}}.\end{split}$ (15) $$

The variable $\mathbf{\Theta}_{t}\in\mathbb{R}^{r_{A}\times r}$ can be roughly viewed as a pseudo-inverse of $\mathbf{\Phi}_{t}^{\top}$ because $\mathbf{\Phi}_{t}^{\top}\mathbf{\Theta}_{t}=\mathbf{I}_{r}$.
We note that the invertibility of $(\mathbf{\Phi}_{t}^{\top}\mathbf{\Phi}_{t})$ will become clear in Lemma 9.
∎

###### Proof.

#### B.3.3 Proof of Theorem 2

###### Proof.

Using $\mathbf{\Phi}_{t}^{\top}\mathbf{\Theta}_{t}=\mathbf{I}_{r}$, definition of $\mathbf{B}_{t}=\mathbf{\Phi}_{t}^{\top}\mathbf{\Sigma}^{-1}\mathbf{\Phi}_{t}$ (at the start of Apdx. B.3), and the update of $\mathbf{\Phi}_{t+1}$ in (15), it is not difficult to verify that

$$ $\displaystyle\mathbf{B}_{t+1}=(1-\eta)^{2}\mathbf{B}_{t}+2\eta(1-\eta)\mathbf{ I}_{r}+\eta^{2}\mathbf{\Theta}_{t}^{\top}\mathbf{\Sigma}\mathbf{\Theta}_{t}.$ (16) $$

Subtracting $\mathbf{I}_{r}$ on both sides of (16), we can get

$$ $\displaystyle\mathbf{B}_{t+1}-\mathbf{I}_{r}=(1-\eta)^{2}(\mathbf{B}_{t}- \mathbf{I}_{r})-\eta^{2}\mathbf{I}_{r}+\eta^{2}\mathbf{\Theta}_{t}^{\top} \mathbf{\Sigma}\mathbf{\Theta}_{t}.$ $$

This ensures that

$$ $\displaystyle~{}~{}~{}~{}\|\mathbf{B}_{t+1}-\mathbf{I}_{r}\|_{\mathsf{F}}$ $\displaystyle\leq(1-\eta)^{2}\|\mathbf{B}_{t}-\mathbf{I}_{r}\|_{\mathsf{F}}+ \eta^{2}\sqrt{r}+\eta^{2}\|\mathbf{\Theta}_{t}^{\top}\mathbf{\Sigma}\mathbf{ \Theta}_{t}\|_{\mathsf{F}}$ $\displaystyle\leq(1-\eta)^{2}\|\mathbf{B}_{t}-\mathbf{I}_{r}\|_{\mathsf{F}}+ \eta^{2}\sqrt{r}+\eta^{2}\frac{r}{\sigma_{r}(\mathbf{B}_{t})}$ $$

where the last inequality is because of Lemma 10. Suppose that $\eta\leq 2/3$, from Lemma 9, one can see that there exists a time $T_{1}$ such that $\sigma_{r}(\mathbf{B}_{t})\geq 1/3,\forall t\geq T_{1}$. We assume $T_{1}=0$ following the same argument (i.e., initialized large with large $\xi$) as previous proofs. With these arguments, we obtain that

$$ $\begin{split}&~{}~{}~{}~{}\|\mathbf{B}_{t+1}-\mathbf{I}_{r}\|_{\mathsf{F}}\\ &\leq(1-\eta)\|\mathbf{B}_{t}-\mathbf{I}_{r}\|_{\mathsf{F}}+\eta^{2}\sqrt{r}+3 r\eta^{2}\\ &\leq\eta\sqrt{r}+3\eta r+(1-\eta)^{t+1-T_{1}}\|\mathbf{B}_{T_{1}}-\mathbf{I}_ {r}\|_{\mathsf{F}}\\ &\leq\eta\sqrt{r}+3\eta r+(1-\eta)^{t+1-T_{1}}\|\mathbf{B}_{T_{1}}-\mathbf{I}_ {r}\|_{\mathsf{F}}.\end{split}$ (17) $$

This implies a linear rate, i.e, $\|\mathbf{B}_{t+1}-\mathbf{I}_{r}\|_{\mathsf{F}}\leq{\cal O}(\eta r)+\epsilon$ if $\eta={\cal O}(1)$ with sufficient iterations.

Inequality (17) also implies that choosing $\eta={\cal O}(\epsilon/r)$, $\|\mathbf{B}_{t+1}-\mathbf{I}_{r}\|_{\mathsf{F}}\leq\epsilon$ at a rate of ${\cal O}(\frac{r}{\epsilon}\log\frac{1}{\epsilon})$. The proof is thus completed.
∎

###### Proof.

#### B.3.4 Proof of Lemma 5

###### Proof.

We start with notation. Let

$$ $\displaystyle\mathbf{\Sigma}=\begin{bmatrix}\mathbf{\Sigma}_{1}&\mathbf{0}\\ \mathbf{0}&\mathbf{\Sigma}_{2}\end{bmatrix},~{}~{}~{}~{}~{}\mathbf{\Phi}_{t}= \begin{bmatrix}\mathbf{M}_{t}\\ \mathbf{N}_{t}\end{bmatrix},$ (18) $$

where $\mathbf{\Sigma}_{1}\in\mathbb{R}^{r\times r}$ is the learnable eigenvalues, while $\mathbf{\Sigma}_{2}\in\mathbb{R}^{(r_{A}-r)\times(r_{A}-r)}$ are the unlearnable eigenvalues, and $\mathbf{M}_{t}\in\mathbb{R}^{r\times r}$ and $\mathbf{N}_{t}\in\mathbb{R}^{(r_{A}-r)\times r}$. Ideally at global convergence, we hope that $\mathbf{M}_{t}\rightarrow\mathbf{\Sigma}_{1}^{1/2}$ up to rotation; while $\mathbf{N}_{t}\rightarrow\mathbf{0}$.

We consider a scenario with $t\rightarrow\infty$, i.e., $\epsilon\rightarrow 0$ and $\mathbf{B}_{t}=\mathbf{I}_{r}$. Using (18) to rewrite $\mathbf{B}_{t}=\mathbf{I}_{r}$, we have that

$$ $\displaystyle\mathbf{M}_{t}^{\top}\mathbf{\Sigma}_{1}^{-1}\mathbf{M}_{t}+ \mathbf{N}_{t}^{\top}\mathbf{\Sigma}_{2}^{-1}\mathbf{N}_{t}=\mathbf{I}_{r}.$ (19) $$

The above equation implies that

$$ $\displaystyle\text{Tr}(\mathbf{M}_{t}^{\top}\mathbf{\Sigma}_{1}^{-1}\mathbf{M} _{t})$ $\displaystyle=\text{Tr}(\mathbf{M}_{t}^{\top}\mathbf{\Sigma}_{1}^{-1/2}\mathbf {\Sigma}_{1}^{-1/2}\mathbf{M}_{t})$ (20) $\displaystyle=\|\mathbf{\Sigma}_{1}^{-1/2}\mathbf{M}_{t}\|_{\mathsf{F}}^{2} \stackrel{{\scriptstyle(a)}}{{\leq}}r$ $$

where (a) is by (19) and Lemma 18.

Since we hope $\mathbf{\Sigma}_{1}^{-1/2}\mathbf{M}_{t}\rightarrow\mathbf{I}_{r}$, we have that

$$ $\begin{split}&~{}~{}~{}~{}\|\mathbf{\Sigma}_{1}^{-1/2}\mathbf{M}_{t}-\mathbf{I }_{r}\|_{\mathsf{F}}^{2}\\ &=\text{Tr}\bigg{(}(\mathbf{\Sigma}_{1}^{-1/2}\mathbf{M}_{t}-\mathbf{I}_{r})^{ \top}(\mathbf{\Sigma}_{1}^{-1/2}\mathbf{M}_{t}-\mathbf{I}_{r})\bigg{)}\\ &=\text{Tr}\big{(}\mathbf{M}_{t}^{\top}\mathbf{\Sigma}_{1}^{-1/2}\mathbf{ \Sigma}_{1}^{-1/2}\mathbf{M}_{t}\big{)}+\text{Tr}(\mathbf{I}_{r})-2\text{Tr}( \mathbf{M}_{t}^{\top}\mathbf{\Sigma}_{1}^{-1/2})\\ &\stackrel{{\scriptstyle(a)}}{{\leq}}\text{Tr}\big{(}\mathbf{M}_{t}^{\top} \mathbf{\Sigma}_{1}^{-1/2}\mathbf{\Sigma}_{1}^{-1/2}\mathbf{M}_{t}\big{)}+ \text{Tr}(\mathbf{I}_{r})+2r^{3/2}\\ &\stackrel{{\scriptstyle(b)}}{{\leq}}2r+2r^{3/2},\end{split}$ (21) $$

where (a) is because that i) for any $r\times r$ matrix $\mathbf{C}$, we have that $\text{Tr}(\mathbf{C})\geq r\min_{i}\mathbf{C}_{ii}\geq-r\|\mathbf{C}\|_{
\mathsf{F}}$, ii) take $\mathbf{C}=\mathbf{M}_{t}^{\top}\mathbf{\Sigma}_{1}^{-1/2}$ and then apply (20); and (b) is by (20).

To bound $\mathbf{N}_{t}$, it can be seen that

$$ $\displaystyle\frac{1}{\sigma_{r+1}(\mathbf{A})}\text{Tr}\big{(}\mathbf{N}_{t}^ {\top}\mathbf{N}_{t}\big{)}\leq\text{Tr}\big{(}\mathbf{N}_{t}^{\top}\mathbf{ \Sigma}_{2}^{-1}\mathbf{N}_{t}\big{)}\stackrel{{\scriptstyle(c)}}{{\leq}}r$ (22) $$

where (c) is by applying Lemma 18 on (19).
This suggests that $\|\mathbf{N}_{t}\|_{\mathsf{F}}\leq\sqrt{r\sigma_{r+1}(\mathbf{A})}$.

Lastly, note that $\mathbf{X}_{*}$ can be written as $\mathbf{X}_{*}=\mathbf{Q}[\mathbf{\Sigma}_{1}^{1/2},\mathbf{0}]^{\top}$ and $\mathbf{X}_{t}=\mathbf{Q}\mathbf{\Phi}_{t}$. Using this fact and combining (21) and (22), we have that

$$ $\begin{split}\|\mathbf{X}_{t}-\mathbf{X}_{*}\|_{\mathsf{F}}^{2}&=\|\mathbf{M}_ {t}-\mathbf{\Sigma}_{1}^{1/2}\|_{\mathsf{F}}^{2}+\|\mathbf{N}_{t}\|_{\mathsf{F }}^{2}\\ &=\|\mathbf{\Sigma}_{1}^{1/2}(\mathbf{\Sigma}_{1}^{-1/2}\mathbf{M}_{t}-\mathbf {I}_{r})\|_{\mathsf{F}}^{2}+\|\mathbf{N}_{t}\|_{\mathsf{F}}^{2}\\ &\leq\sigma_{1}(\mathbf{\Sigma}_{1}^{1/2})^{2}\|\mathbf{\Sigma}_{1}^{-1/2} \mathbf{M}_{t}-\mathbf{I}_{r}\|_{\mathsf{F}}^{2}+\|\mathbf{N}_{t}\|_{\mathsf{F }}^{2}\\ &={\cal O}(r^{3/2}),\end{split}$ (23) $$

where we used $\sigma_{1}(\mathbf{\Sigma})=1$ and $\sigma_{r+1}(\mathbf{\Sigma})\leq 1$. The proof is thus completed.
∎

###### Proof.

#### B.3.5 Useful lemmas for symmetric and under-parametrized problems

It is clear that $\mathbf{B}_{t}$ is symmetric by definition, i.e., $\mathbf{B}_{t}=\mathbf{\Phi}_{t}^{\top}\mathbf{\Sigma}^{-1}\mathbf{\Phi}_{t}$. This enables us to give a lower bound on $\sigma_{r}(\mathbf{B}_{t})$ using Lemma 15.

###### Lemma 9 .

$\sigma_{r}(\mathbf{B}_{t})$ is lower bounded by

$$ $\displaystyle\sigma_{r}(\mathbf{B}_{t+1})\geq(1-\eta)-(1-\eta)^{2t+3}+(1-\eta) ^{2t+2}\sigma_{r}(\mathbf{B}_{0}).$ $$

###### Proof.

Given the definition of $\mathbf{B}_{t}$, it is not difficult to see that $\mathbf{B}_{t}$ is PSD for all $t$. We can then apply Lemma 15 on (16) to arrive at

$$ $\displaystyle~{}~{}~{}~{}~{}\sigma_{r}(\mathbf{B}_{t+1})$ $\displaystyle\geq 2\eta(1-\eta)+\sigma_{r}\big{(}(1-\eta)^{2}\mathbf{B}_{t}+ \eta^{2}\mathbf{\Theta}_{t}^{\top}\mathbf{\Sigma}\mathbf{\Theta}_{t}\big{)}$ $\displaystyle\geq 2\eta(1-\eta)+(1-\eta)^{2}\sigma_{r}\big{(}\mathbf{B}_{t} \big{)}$ $\displaystyle\stackrel{{\scriptstyle(a)}}{{\geq}}(1-\eta)^{2t+2}\sigma_{r}( \mathbf{B}_{0})+2\eta(1-\eta)\frac{1-(1-\eta)^{2t+2}}{2\eta-\eta^{2}}$ $\displaystyle\stackrel{{\scriptstyle(b)}}{{\geq}}(1-\eta)^{2t+2}\sigma_{r}( \mathbf{B}_{0})+(1-\eta)-(1-\eta)^{2t+3}$ $$

where (a) uses Lemma 14 to unroll $\sigma_{r}(\mathbf{B}_{t})$; and (b) is because $\frac{2\eta}{2\eta-\eta^{2}}\geq 1$.
∎

###### Lemma 10 .

Let $\mathbf{\Theta}_{t}$ and $\mathbf{B}_{t}$ defined the same as those in Apdx. B.3. It is guaranteed to have that

$$ $\displaystyle\|\mathbf{\Theta}_{t}^{\top}\mathbf{\Sigma}\mathbf{\Theta}_{t}\|_ {\mathsf{F}}\leq\frac{r}{\sigma_{r}(\mathbf{B}_{t})}.$ $$

###### Proof.

Using the inequality $\|\mathbf{A}^{\top}\mathbf{A}\|_{\mathsf{F}}\leq\|\mathbf{A}\|_{\mathsf{F}}^{2}$, we have that

$$ $\displaystyle\|\mathbf{\Theta}_{t}^{\top}\mathbf{\Sigma}\mathbf{\Theta}_{t}\|_ {\mathsf{F}}=\|\mathbf{\Theta}_{t}^{\top}\mathbf{\Sigma}^{1/2}\mathbf{\Sigma}^ {1/2}\mathbf{\Theta}_{t}\|_{\mathsf{F}}\leq\|\mathbf{\Sigma}^{1/2}\mathbf{ \Theta}_{t}\|_{\mathsf{F}}^{2}.$ (24) $$

Now let $\mathbf{E}_{t}:=\mathbf{\Sigma}^{1/2}\mathbf{\Theta}_{t}$ and $\mathbf{F}_{t}:=\mathbf{\Sigma}^{-1/2}\mathbf{\Phi}_{t}$. Since we have that $\mathbf{F}_{t}^{\top}\mathbf{E}_{t}=\mathbf{I}_{r}$, we have that

$$ $\displaystyle\|\mathbf{F}_{t}^{\top}\mathbf{E}_{t}\|_{\mathsf{F}}=\|\mathbf{I} _{r}\|_{\mathsf{F}}=\sqrt{r}.$ $$

Since we also have that

$$ $\displaystyle\sqrt{r}=\|\mathbf{F}_{t}^{\top}\mathbf{E}_{t}\|_{\mathsf{F}} \stackrel{{\scriptstyle(a)}}{{\geq}}\sigma_{r}(\mathbf{F}_{t})\|\mathbf{E}_{t} \|_{\mathsf{F}}\stackrel{{\scriptstyle(b)}}{{=}}\sqrt{\sigma_{r}(\mathbf{B}_{t })}\|\mathbf{E}_{t}\|_{\mathsf{F}},$ (25) $$

where (a) holds because $\mathbf{E}_{t}$ and $\mathbf{F}_{t}$ share the same column space and row space and both of them have rank $r$, which implies that $\langle Null(\mathbf{F}),[\mathbf{E}_{t}]_{i}\rangle=\mathbf{0},\forall i$ ($[\mathbf{E}_{t}]_{i}$ is the $i$th column of $\mathbf{E}_{t}$). Note that (a) does not hold true for general two matrices $\mathbf{E}_{t}$ and $\mathbf{F}_{t}$. (b) is because $\mathbf{F}_{t}^{\top}\mathbf{F}_{t}=\mathbf{B}_{t}$, which means that the singular values of $\mathbf{F}_{t}$ are just square root of eigenvalues of $\mathbf{B}_{t}$. This implies that $\|\mathbf{E}_{t}\|_{\mathsf{F}}\leq\sqrt{r}/\sqrt{\sigma_{r}(\mathbf{B}_{t})}$. Combining this inequality with (24), we have that

$$ $\displaystyle\|\mathbf{\Theta}_{t}^{\top}\mathbf{\Sigma}\mathbf{\Theta}_{t}\|_ {\mathsf{F}}\leq\|\mathbf{\Theta}_{t}^{\top}\mathbf{\Sigma}^{1/2}\|_{\mathsf{F }}^{2}=\|\mathbf{E}_{t}\|_{\mathsf{F}}^{2}\leq\frac{r}{\sigma_{r}(\mathbf{B}_{ t})}.$ $$

The proof is thus completed.
∎

###### Lemma 9 .

###### Proof.

###### Lemma 10 .

###### Proof.

### B.4 Symmetric and over-parametrized setting

Nyström initialization for over-parametrization. While the initialization still follows (3), we need to adapt Lemma 1 to the over-parameterized setting, i.e., $r>r_{A}$.

###### Lemma 11 (Initialization for over-parametrization) .

There exists a universal constant $\tau>0$ such that $\sigma_{r_{A}}(\mathbf{X}_{0})\geq\xi\tau(\sqrt{r}-\sqrt{r_{A}-1})\sigma_{r_{A
}}(\mathbf{A})$ is satisfied with high probability. In other words, $\text{rank}(\mathbf{X}_{0})=r_{A}$ w.h.p.

###### Proof.

Similar to the proof of Lemma 1, let the compact eigenvalue decomposition of $\mathbf{A}$ be $\mathbf{A}=\mathbf{Q}\mathbf{\Sigma}\mathbf{Q}^{\top}$, where $\mathbf{Q}\in\mathbb{R}^{m\times r_{A}}$ and $\mathbf{\Sigma}\in\mathbb{R}^{r_{A}\times r_{A}}$. This implies that $\mathbf{X}_{0}=(\mathbf{Q}\mathbf{\Sigma})(\mathbf{Q}^{\top}\mathbf{\Omega})$.

It is not hard to verify that the matrix $\mathbf{Q}^{\top}\mathbf{\Omega}\in\mathbb{R}^{r_{A}\times r}$ is also a Gaussian random matrix, where each entry follows ${\cal N}(0,\xi^{2})$. Applying Lemma 19 on $(\mathbf{Q}^{\top}\mathbf{\Omega})^{\top}$, and using the fact $(\mathbf{Q}^{\top}\mathbf{\Omega})^{\top}$ and $(\mathbf{Q}^{\top}\mathbf{\Omega})$ share the same singular values, it can be seen that

$$ $\displaystyle\mathbb{P}\Big{(}\frac{\sigma_{r_{A}}(\mathbf{Q}^{\top}\mathbf{ \Omega})}{\xi}\leq\tau(\sqrt{r}-\sqrt{r_{A}-1})\Big{)}\leq(C_{1}\tau)^{r-r_{A} +1}+e^{-C_{2}r}:=\delta_{2}$ $$

where $C_{1}$ and $C_{2}$ are universal constants independent of $r_{A}$ and $r$. This inequality shows that with probability at least $1-\delta_{2}$, $\sigma_{r_{A}}(\mathbf{Q}^{\top}\mathbf{\Omega})\geq\xi\tau(\sqrt{r}-\sqrt{r_{
A}-1})$.

Note that inequality $\sigma_{\min}(\mathbf{C}\mathbf{D})\geq\sigma_{\min}(\mathbf{C})\sigma_{\min}(
\mathbf{D})$ holds given full column rank of $\mathbf{C}$; see Lemma 17. Applying it to (9), we have that

$$ $\displaystyle\sigma_{r_{A}}(\mathbf{X}_{0})$ $\displaystyle\geq\sigma_{r_{A}}(\mathbf{Q}\mathbf{\Sigma})\sigma_{r_{A}}( \mathbf{Q}^{\top}\mathbf{\Omega})=\sigma_{r_{A}}(\mathbf{A})\sigma_{r_{A}}( \mathbf{Q}^{\top}\mathbf{\Omega})$ $\displaystyle\stackrel{{\scriptstyle(a)}}{{\geq}}\xi\tau(\sqrt{r}-\sqrt{r_{A}- 1})\sigma_{r_{A}}(\mathbf{A})$ $$

where (a) holds with probability at least $1-\delta_{2}$.
∎

Next, we provide additional results of Nyström initialization on over-paramtrized setting of problem (1), where we have $r_{A}<r$. For a desirable convergence rate, we need to slightly modify the ScaledGD update to

$$ $\displaystyle\mathbf{X}_{t+1}=\mathbf{X}_{t}-\eta(\mathbf{X}_{t}\mathbf{X}_{t} ^{\top}-\mathbf{A})\mathbf{X}_{t}(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})^{ \dagger}.$ (26) $$

Compared with iteration (2) for exact-parametrization, the modification is on $(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})^{\dagger}$. This pseudo-inverse is necessary because $(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})$ is not necessarily invertible in the over-parametrized setting. We note that unlike previous work (Xu et al., ) which modifies the same term to $(\mathbf{X}_{t}^{\top}\mathbf{X}_{t}+\lambda\mathbf{I})^{-1}$, (26) does not need the damping parameter $\lambda\mathbf{I}$ in the preconditioner. We will observe shortly in Fig. 4 that the quadratic rate is not achieved with the damping factor.

Let the compact eigendecomposition of $\mathbf{A}=\mathbf{Q}\mathbf{\Sigma}\mathbf{Q}^{\top}$ for $\mathbf{Q}\in\mathbb{R}^{m\times r_{A}}$, and $\mathbf{\Sigma}\in\mathbb{R}^{r_{A}\times r_{A}}$. We can also establish that $\mathbf{X}_{t}$ affords a simpler representation.

###### Lemma 12 .

Under the Nyström initialization (3) and iteration (26), the variable $\mathbf{X}_{t}$ can be written as $\mathbf{X}_{t}=\mathbf{Q}\mathbf{\Phi}_{t}$ for some $\mathbf{\Phi}_{t}\in\mathbb{R}^{r_{A}\times r}$. Moreover, we have that

$$ $\displaystyle\mathbf{\Phi}_{t+1}=(1-\eta)\mathbf{\Phi}_{t}+\eta\mathbf{\Sigma} (\mathbf{\Phi}_{t}^{\dagger})^{\top}.$ (27) $$

###### Proof.

We prove this by induction. Clearly, our initialization satisfies this because $\mathbf{X}_{0}=\mathbf{A}\mathbf{\Omega}=\mathbf{Q}\mathbf{\Sigma}\mathbf{Q}^{
\top}\mathbf{\Omega}$, i.e., $\mathbf{\Phi}_{0}:=\mathbf{\Sigma}\mathbf{Q}^{\top}\mathbf{\Omega}$. Now suppose that $\mathbf{X}_{t}=\mathbf{Q}\mathbf{\Phi}_{t}$ holds for $t$. We then show that $\mathbf{X}_{t+1}=\mathbf{Q}\mathbf{\Phi}_{t+1}$ to finish the induction. In particular, plugging $\mathbf{X}_{t}=\mathbf{Q}\mathbf{\Phi}_{t}$ into (26), we arrive at

$$ $\displaystyle\mathbf{X}_{t+1}=\mathbf{Q}\underbrace{\bigg{[}\mathbf{\Phi}_{t}- \eta(\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}-\mathbf{\Sigma})\mathbf{\Phi}_{ t}(\mathbf{\Phi}_{t}^{\top}\mathbf{\Phi}_{t})^{\dagger}\bigg{]}}_{:=\mathbf{ \Phi}_{t+1}}.$ $$

Clearly, the term inside the brackets is $\mathbf{\Phi}_{t+1}$. The induction is thus finished.

Now we proof the second part of this lemma. Let the SVD of $\mathbf{\Phi}_{t}:=\mathbf{U}_{t}\mathbf{\Sigma}_{t}\mathbf{V}_{t}^{\top}$, where $\mathbf{U}_{t}\in\mathbb{R}^{r_{A}\times r_{A}}$, $\mathbf{\Sigma}_{t}\in\mathbb{R}^{r_{A}\times r_{A}}$, and $\mathbf{V}_{t}\in\mathbb{R}^{r\times r_{A}}$. We note that $\mathbf{U}_{t}$ is unitary for this case. With the SVD, we have that $\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}=\mathbf{U}_{t}\mathbf{\Sigma}_{t}^{2
}\mathbf{U}_{t}^{\top}$, and $(\mathbf{\Phi}_{t}^{\top}\mathbf{\Phi}_{t})^{\dagger}=\mathbf{V}_{t}\mathbf{
\Sigma}_{t}^{-2}\mathbf{V}_{t}^{\top}$. Plugging these into $\mathbf{\Phi}_{t+1}$ defined earlier, we arrive at

$$ $\displaystyle\mathbf{\Phi}_{t+1}$ $\displaystyle=\mathbf{\Phi}_{t}-\eta(\mathbf{U}_{t}\mathbf{\Sigma}_{t}^{2} \mathbf{U}_{t}^{\top}-\mathbf{\Sigma})\mathbf{U}_{t}\mathbf{\Sigma}_{t}\mathbf {V}_{t}^{\top}\mathbf{V}_{t}\mathbf{\Sigma}_{t}^{-2}\mathbf{V}_{t}^{\top}$ $\displaystyle=\mathbf{\Phi}_{t}-\eta(\mathbf{U}_{t}\mathbf{\Sigma}_{t}^{2} \mathbf{U}_{t}^{\top}-\mathbf{\Sigma})\mathbf{U}_{t}\mathbf{\Sigma}_{t}^{-1} \mathbf{V}_{t}^{\top}$ $\displaystyle=\mathbf{\Phi}_{t}-\eta\mathbf{U}_{t}\mathbf{\Sigma}_{t}\mathbf{V }_{t}^{\top}+\eta\mathbf{\Sigma}\mathbf{U}_{t}\mathbf{\Sigma}_{t}^{-1}\mathbf{ V}_{t}^{\top}$ $\displaystyle=(1-\eta)\mathbf{\Phi}_{t}+\eta\mathbf{\Sigma}(\mathbf{\Phi}_{t}^ {\dagger})^{\top}.$ $$

This completes the proof.
∎

Next, let $\mathbf{B}_{t}=\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}$. With (27) we have that

$$ $\begin{split}\mathbf{B}_{t+1}&=(1-\eta)^{2}\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^ {\top}+\eta(1-\eta)\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\dagger}\mathbf{\Sigma} +\eta(1-\eta)\mathbf{\Sigma}(\mathbf{\Phi}_{t}^{\dagger})^{\top}\mathbf{\Phi}_ {t}^{\top}+\eta^{2}\mathbf{\Sigma}(\mathbf{\Phi}_{t}^{\dagger})^{\top}\mathbf{ \Phi}_{t}^{\dagger}\mathbf{\Sigma}\\ &\stackrel{{\scriptstyle(a)}}{{=}}(1-\eta)^{2}\mathbf{B}_{t}+2\eta(1-\eta) \mathbf{\Sigma}+\eta^{2}\mathbf{\Sigma}(\mathbf{\Phi}_{t}^{\dagger})^{\top} \mathbf{\Phi}_{t}^{\dagger}\mathbf{\Sigma}\\ &\stackrel{{\scriptstyle(b)}}{{=}}(1-\eta)^{2}\mathbf{B}_{t}+2\eta(1-\eta) \mathbf{\Sigma}+\eta^{2}\mathbf{\Sigma}\mathbf{B}_{t}^{-1}\mathbf{\Sigma},\end {split}$ (28) $$

where in (a) we used the SVD of $\mathbf{\Phi}_{t}:=\mathbf{U}_{t}\mathbf{\Sigma}_{t}\mathbf{V}_{t}^{\top}$, where $\mathbf{U}_{t}\in\mathbb{R}^{r_{A}\times r_{A}}$, $\mathbf{\Sigma}_{t}\in\mathbb{R}^{r_{A}\times r_{A}}$ and $\mathbf{V}_{t}\in\mathbb{R}^{r\times r_{A}}$, $\mathbf{\Phi}_{t}^{\dagger}=\mathbf{V}_{t}\mathbf{\Sigma}_{t}^{-1}\mathbf{U}_{
t}^{\top}$, and $\mathbf{U}_{t}$ is unitary; and in (b) we assume that $\mathbf{B}_{t}$ is full rank. Note that this assumption can be easily verified given $\text{rank}(\mathbf{B}_{0})=r_{A}$; and the iteration on $\mathbf{B}_{t}$ (28) is exactly the same as in exact-parametrized cases (12). The latter allows us to bound $\sigma_{r_{A}}(\mathbf{B}_{t})$ away from 0 0 in the same way as Lemma 2.

In other words, the over-parametrized case under our initialization reduces to the exact-parametrized case given the same iteration on $\mathbf{B}_{t}$ (28) (cf. (12)). This allows as to use the same argument of Theorem 1 to derive a quadratic rate for over-parametrized case.

###### Theorem 5 .

With high probability over the initialization, the behavior of update (26) under Nyström initialization (3) can be described as:

Phase 1 (linear convergence). Let $\eta={\cal O}(\frac{1}{\kappa^{3}\|\mathbf{A}\|_{\mathsf{F}}})$. After $T_{1}:={\cal O}(\kappa^{3}\sqrt{r}\log\kappa)$ iterations, ScaledGD ensures that $\|\mathbf{X}_{T_{1}}\mathbf{X}_{T_{1}}^{\top}-\mathbf{A}\|_{\mathsf{F}}\leq{
\cal O}(1/\kappa^{2})$.

Phase 2 (quadratic convergence). After Phase I, ScaledGD converges quadratically with $\eta=0.5$. In particular, $\|\mathbf{X}_{T}\mathbf{X}_{T}^{\top}-\mathbf{A}\|_{\mathsf{F}}\leq\epsilon$ is ensured after $T={\cal O}\big{(}\log\log(\frac{1}{\kappa\epsilon})\big{)}$ iterations.

###### Proof.

The proof is the same as Theorem 1 given the same iteration on $\mathbf{B}_{t}$ in (28). We omit it to avoid redundancy.
∎

Numerical illustration. A numerical illustration for ScaledGD under Nyström initialization in over-parametrized case can be found in Fig. 4. We adopt ScaledGD-($\lambda$) (Xu et al., ), the damping version of ScaledGD, as another baseline. It can be seen that only our approach achieves a quadratic rate; see Fig. 4(a). We also slightly perturb our initialization with small noise, and it can be seen that the quadratic convergence breaks down immediately. This demonstrate the critical role of initialization: i) it helps to get rid of damping using pseudo-inverse; and ii) it ensures a quadratic rate.

Figure: Figure 4: Convergence of ScaledGD under Nyström initialization (optimality error vs. iteration) on over-parametrized problems detailed in Apdx. E.1. (a) Comparison of GD, ScaledGD-$(\lambda)$ with small initialization, and ScaledGD with our initialization. (b) Solid lines show that our initialization is not sensitive to magnitude; and dotted lines illustrate that quadratic convergence cannot be obtained even with slightly perturbed initialization, i.e., $\mathbf{X}_{0}=\mathbf{A}\mathbf{\Omega}+\mathbf{N}$, where $[\mathbf{N}]_{ij}~{}\sim{\cal N}(0,\xi_{n}^{2})$.
Refer to caption: x11.png

###### Lemma 11 (Initialization for over-parametrization) .

###### Proof.

###### Lemma 12 .

###### Proof.

###### Theorem 5 .

###### Proof.

## Appendix C Missing proofs for asymmetric settings

### C.1 Missing proofs for asymmetric and exact-parametrized setting

#### C.1.1 Proof of Lemma 6

###### Proof.

The proof is finished by induction. From our Nyström initialization, one has that $\mathbf{\Psi}_{0}=\mathbf{0}$ and $\mathbf{\Phi}_{0}=\mathbf{\Sigma}\mathbf{V}^{\top}\mathbf{\Omega}$. Now assume that one can write $\mathbf{X}_{t}=\mathbf{U}\mathbf{\Phi}_{t}$ and $\mathbf{Y}_{t}=\mathbf{V}\mathbf{\Psi}_{t}$ for some iteration $t$. We will show that $\mathbf{X}_{t+1}=\mathbf{U}\mathbf{\Phi}_{t+1}$ and $\mathbf{Y}_{t+1}=\mathbf{V}\mathbf{\Psi}_{t+1}$ under iteration (6). Let us start with $\mathbf{X}_{t+1}$. Note that if $t=0$, $\mathbf{X}_{1}=\mathbf{U}\mathbf{\Phi}_{1}$ is trivial. We only focus on $t\geq 1$, where we have

$$ $\displaystyle\mathbf{X}_{t+1}$ $\displaystyle=\mathbf{X}_{t}-\eta(\mathbf{X}_{t}\mathbf{Y}_{t}^{\top}-\mathbf{ A})\mathbf{Y}_{t}(\mathbf{Y}_{t}^{\top}\mathbf{Y}_{t})^{-1}$ $\displaystyle=\mathbf{U}\mathbf{\Phi}_{t}-\eta(\mathbf{U}\mathbf{\Phi}_{t} \mathbf{\Psi}_{t}^{\top}\mathbf{V}^{\top}-\mathbf{U}\mathbf{\Sigma}\mathbf{V}^ {\top})\mathbf{V}\mathbf{\Psi}_{t}(\mathbf{\Psi}_{t}^{\top}\mathbf{V}^{\top} \mathbf{V}\mathbf{\Psi}_{t})^{-1}$ $\displaystyle=\mathbf{U}\mathbf{\Phi}_{t}-\eta\mathbf{U}(\mathbf{\Phi}_{t} \mathbf{\Psi}_{t}^{\top}-\mathbf{\Sigma})\mathbf{\Psi}_{t}(\mathbf{\Psi}_{t}^{ \top}\mathbf{\Psi}_{t})^{-1}$ $\displaystyle=\mathbf{U}\underbrace{\bigg{(}\mathbf{\Phi}_{t}-\eta(\mathbf{ \Phi}_{t}\mathbf{\Psi}_{t}^{\top}-\mathbf{\Sigma})\mathbf{\Psi}_{t}(\mathbf{ \Psi}_{t}^{\top}\mathbf{\Psi}_{t})^{-1}\bigg{)}}_{:=\mathbf{\Phi}_{t+1}}.$ $$

Note that the invertible of $(\mathbf{\Psi}_{t}^{\top}\mathbf{\Psi}_{t})$ will become clear in the proof of Corollary 1.

Using a similar argument, it is not hard to show that $\mathbf{Y}_{t}=\mathbf{V}\mathbf{\Psi}_{t}$ for all $t$. We do not repeat here.
∎

###### Proof.

#### C.1.2 Proof of Theorem 3

###### Proof.

Based on the initialization (5) and iteration (6), we can obtain that

$$ $\displaystyle\mathbf{\Phi}_{1}=\mathbf{\Phi}_{0}$ (29a) $\begin{split}\mathbf{\Psi}_{1}&=\mathbf{V}^{\top}\mathbf{Y}_{1}=\mathbf{0}- \eta\mathbf{V}^{\top}(\mathbf{0}-\mathbf{A})^{\top}\mathbf{U}\mathbf{\Phi}_{0} (\mathbf{\Phi}_{0}^{\top}\mathbf{U}^{\top}\mathbf{U}\mathbf{\Phi}_{0})^{-1}\\ &=\eta\mathbf{V}^{\top}\mathbf{V}\mathbf{\Sigma}\mathbf{U}^{\top}\mathbf{U} \mathbf{\Phi}_{0}(\mathbf{\Phi}_{0}^{\top}\mathbf{U}^{\top}\mathbf{U}\mathbf{ \Phi}_{0})^{-1}\\ &=\eta\mathbf{\Sigma}\mathbf{\Phi}_{0}(\mathbf{\Phi}_{0}^{\top}\mathbf{\Phi}_{ 0})^{-1}\\ &=\eta\mathbf{\Sigma}\mathbf{\Phi}_{0}^{-\top}.\end{split}$ (29b) $$

This ensures that

$$ $\displaystyle\mathbf{\Phi}_{1}\mathbf{\Psi}_{1}^{\top}=\eta\mathbf{\Sigma}.$ $$

Choosing $\eta=1$ completes the proof.
∎

###### Proof.

#### C.1.3 Proof of Corollary 1

###### Proof.

The corollary is proved through an asymmetric-to-symmetric reduction.

Step 1. Positive definiteness of $\mathbf{\Phi}_{t}\mathbf{\Psi}_{t}^{\top}$.
We will first show that $\mathbf{\Phi}_{t}\mathbf{\Psi}_{t}^{\top}$ is symmetric and positive definite (PD) for any $t\geq 1$. From the proof of Theorem 3, it can be seen that $\mathbf{\Phi}_{1}\mathbf{\Psi}_{1}^{\top}=\eta\mathbf{\Sigma}$ is symmetric and PD. This means that the base case of induction holds. Now suppose that $\mathbf{\Phi}_{t}\mathbf{\Psi}_{t}^{\top}$ is symmetric and PD at iteration $t$. Based on Lemma 6, we can write the iteration as

$$ $\displaystyle\mathbf{\Phi}_{t+1}$ $\displaystyle=(1-\eta)\mathbf{\Phi}_{t}+\eta\mathbf{\Sigma}\mathbf{\Psi}_{t}^{ -\top}$ (30a) $\displaystyle\mathbf{\Psi}_{t+1}$ $\displaystyle=(1-\eta)\mathbf{\Psi}_{t}+\eta\mathbf{\Sigma}\mathbf{\Phi}_{t}^{ -\top}.$ (30b) $$

This gives that

$$ $\displaystyle\mathbf{\Phi}_{t+1}\mathbf{\Psi}_{t+1}^{\top}=(1-\eta)^{2}\mathbf {\Phi}_{t}\mathbf{\Psi}_{t}^{\top}$ $\displaystyle+2\eta(1-\eta)\mathbf{\Sigma}+\eta^{2}\mathbf{\Sigma}(\mathbf{ \Phi}_{t}\mathbf{\Psi}_{t}^{\top})^{-1}\mathbf{\Sigma}.$ (31) $$

The symmetry of $\mathbf{\Phi}_{t+1}\mathbf{\Psi}_{t+1}^{\top}$ directly follows from (31). For the positive definiteness of $\mathbf{\Phi}_{t+1}\mathbf{\Psi}_{t+1}^{\top}$, we can apply Lemma 15 to get

$$ $\displaystyle\lambda_{\min}(\mathbf{\Phi}_{t+1}\mathbf{\Psi}_{t+1}^{\top})\geq (1-\eta)^{2}\lambda_{\min}(\mathbf{\Phi}_{t}\mathbf{\Psi}_{t}^{\top})$ $\displaystyle+2\eta(1-\eta)\lambda_{\min}(\mathbf{\Sigma})+\eta^{2}\lambda_{ \min}(\mathbf{\Sigma}(\mathbf{\Phi}_{t}\mathbf{\Psi}_{t}^{\top})^{-1}\mathbf{ \Sigma})>0.$ $$

This concludes the PD of $\mathbf{\Phi}_{t+1}\mathbf{\Psi}_{t+1}^{\top}$.

Step 2. Define $\mathbf{B}_{t}:=\mathbf{\Phi}_{t}\mathbf{\Psi}_{t}^{\top}$, then (31) can be rewritten as

$$ $\displaystyle\mathbf{B}_{t+1}=(1-\eta)^{2}\mathbf{B}_{t}+2\eta(1-\eta)\mathbf{ \Sigma}+\eta^{2}\mathbf{\Sigma}\mathbf{B}_{t}^{-1}\mathbf{\Sigma}$ (32) $$

which is exactly the same iteration as (12) for the symmetric exact-parametrized case. Based on the results from Step 1, that is, $\mathbf{\Phi}_{t+1}\mathbf{\Psi}_{t+1}^{\top}$ is symmetric and PD, we can apply the same analysis steps for symmetric exact-parametrized problems, i.e., Theorem 1 to get the bounds stated in this corollary. We do not repeat for conciseness.
∎

###### Proof.

### C.2 Missing proofs for asymmetric and under-parametrized setting

#### C.2.1 How good is weak optimality?

###### Lemma 13 .

Every global optimum for (4) is also weakly optimal.

###### Proof.

We start with rewriting the SVD of $\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ as

$$ $\displaystyle\mathbf{A}=[\mathbf{U}_{1},\mathbf{U}_{2}]\begin{bmatrix}\mathbf{ \Sigma}_{1}&\mathbf{0}\\ \mathbf{0}&\mathbf{\Sigma}_{2}\end{bmatrix}\begin{bmatrix}\mathbf{V}_{1}^{\top }\\ \mathbf{V}_{2}^{\top}\end{bmatrix}=\mathbf{U}_{1}\mathbf{\Sigma}_{1}\mathbf{V} _{1}^{\top}+\mathbf{U}_{2}\mathbf{\Sigma}_{2}\mathbf{V}_{2}^{\top}$ (33) $$

where $\mathbf{U}_{1}\in\mathbb{R}^{m\times r}$ and $\mathbf{U}_{2}\in\mathbb{R}^{m\times(r_{A}-r)}$ are the first $r$ and other columns of $\mathbf{U}$, respectively; $\mathbf{\Sigma}_{1}\in\mathbb{R}^{r\times r}$ and $\mathbf{\Sigma}_{2}\in\mathbb{R}^{(r-r_{A})\times(r-r_{A})}$ are diagonal matrices formed by the first $r$ and rest diagonal entries of $\mathbf{\Sigma}$; and $\mathbf{V}_{1}\in\mathbb{R}^{n\times r}$ and $\mathbf{V}_{2}\in\mathbb{R}^{n\times(r_{A}-r)}$ are the first $r$ and other columns of $\mathbf{V}$.

It is not hard to see that the optimal solutions of (1) are $\mathbf{X}_{*}=\mathbf{U}_{1}\mathbf{\Sigma}_{1}^{1/2}\mathbf{Q}$ and $\mathbf{Y}_{*}=\mathbf{V}_{1}\mathbf{\Sigma}_{1}^{1/2}\mathbf{Q}^{-\top}$, where $\mathbf{Q}\in\mathbb{R}^{r\times r}$ is any invertible matrix. Using these notation, we have that

$$ $\displaystyle\mathbf{Y}_{*}^{\top}\mathbf{A}^{\dagger}\mathbf{X}_{*}$ $\displaystyle=\mathbf{Q}^{-1}\mathbf{\Sigma}_{1}^{1/2}\mathbf{V}_{1}^{\top}( \mathbf{V}_{1}\mathbf{\Sigma}_{1}^{-1}\mathbf{U}_{1}^{\top}+\mathbf{V}_{2} \mathbf{\Sigma}_{2}^{-1}\mathbf{U}_{2}^{\top})\mathbf{U}_{1}\mathbf{\Sigma}_{1 }^{1/2}\mathbf{Q}$ $\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathbf{I}_{r}$ $$

where in (a) we use the facts $\mathbf{U}_{1}^{\top}\mathbf{U}_{1}=\mathbf{I}_{r}$ and $\mathbf{U}_{1}^{\top}\mathbf{U}_{2}=\mathbf{0}_{r\times(r_{A}-r)}$. This concludes the proof.

∎

###### Lemma 13 .

###### Proof.

#### C.2.2 Proof of Theorem 4

###### Proof.

The update in (6) ensures that

$$ $\displaystyle\mathbf{\Phi}_{1}=\mathbf{\Phi}_{0},~{}~{}~{}~{}~{}~{}$ (34a) $\begin{split}\mathbf{\Psi}_{1}&=\mathbf{V}^{\top}\mathbf{Y}_{1}=\mathbf{0}- \eta\mathbf{V}^{\top}(\mathbf{0}-\mathbf{A})^{\top}\mathbf{U}\mathbf{\Phi}_{0} (\mathbf{\Phi}_{0}^{\top}\mathbf{U}^{\top}\mathbf{U}\mathbf{\Phi}_{0})^{-1}\\ &=\eta\mathbf{V}^{\top}\mathbf{V}\mathbf{\Sigma}\mathbf{U}^{\top}\mathbf{U} \mathbf{\Phi}_{0}(\mathbf{\Phi}_{0}^{\top}\mathbf{U}^{\top}\mathbf{U}\mathbf{ \Phi}_{0})^{-1}\\ &=\eta\mathbf{\Sigma}\mathbf{\Phi}_{0}(\mathbf{\Phi}_{0}^{\top}\mathbf{\Phi}_{ 0})^{-1}\\ &\stackrel{{\scriptstyle(a)}}{{:=}}\eta\mathbf{\Sigma}\mathbf{\Theta}_{0}\end{split}$ (34b) $$

where in (a) we define $\mathbf{\Theta}_{t}:=\mathbf{\Phi}_{t}(\mathbf{\Phi}_{t}^{\top}\mathbf{\Phi}_{
t})^{-1}$.

From the Definition 2, we can see that

$$ $\displaystyle\mathbf{Y}_{1}^{\top}\mathbf{A}^{\dagger}\mathbf{X}_{1}$ $\displaystyle=\mathbf{\Psi}_{1}^{\top}\mathbf{V}^{\top}\mathbf{V}\mathbf{ \Sigma}^{-1}\mathbf{U}^{\top}\mathbf{U}\mathbf{\Phi}_{1}=\mathbf{\Psi}_{1}^{ \top}\mathbf{\Sigma}^{-1}\mathbf{\Phi}_{1}$ $\displaystyle=\eta\mathbf{\Theta}_{0}^{\top}\mathbf{\Sigma}\mathbf{\Sigma}^{-1 }\mathbf{\Phi}_{0}=\eta\mathbf{I}_{r}.$ $$

This means that when $\eta=1$, generalized weak optimality can be achieved in one step for under-parametrized problems.
∎

###### Proof.

### C.3 Asymmetric and over-parametrized setting

Next, we establish the one step convergence with Nyström initialization in the asymmetric over-parametrized setting, where $r_{A}<r$. We also need to slightly modify the ScaledGD update to

$$ $\displaystyle\mathbf{X}_{1}=\mathbf{X}_{0},~{}\text{and}~{}\mathbf{X}_{t+1}= \mathbf{X}_{t}-\eta(\mathbf{X}_{t}\mathbf{Y}_{t}^{\top}-\mathbf{A})\mathbf{Y}_ {t}(\mathbf{Y}_{t}^{\top}\mathbf{Y}_{t})^{\dagger},\forall t\geq 1$ (35a) $\displaystyle\mathbf{Y}_{t+1}=\mathbf{Y}_{t}-\eta(\mathbf{X}_{t}\mathbf{Y}_{t} ^{\top}-\mathbf{A})^{\top}\mathbf{X}_{t}(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})^ {\dagger},\forall t\geq 0.$ (35b) $$

Comparing with (6), the difference is that here we use pseudo-inverse to bypass the possible non-invertibility of $(\mathbf{X}_{t}^{\top}\mathbf{X}_{t})$ and $(\mathbf{Y}_{t}^{\top}\mathbf{Y}_{t})$ in the over-parametrized case. We also note that to the best of our knowledge, there is no previous result that establishes the convergence of ScaledGD (or its variants) for asymmetric over-parametrized problems.

###### Theorem 6 .

Under Nyström initialization (5), the modified ScaledGD iterations (35) converge globally in a single step w.h.p. over the initialization, i.e., $\mathbf{X}_{1}\mathbf{Y}_{1}^{\top}=\mathbf{A}$ if the learning rate is chosen as $\eta=1$.

###### Proof.

Let the compact eigendecomposition of $\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ for $\mathbf{U}\in\mathbb{R}^{m\times r_{A}}$, $\mathbf{\Sigma}\in\mathbb{R}^{r_{A}\times r_{A}}$, and $\mathbf{V}\in\mathbb{R}^{n\times r_{A}}$.

The Nyström initialization ensures that $\mathbf{X}_{0}=\mathbf{X}_{1}=\mathbf{U}\mathbf{\Phi}_{0}$, where $\mathbf{\Phi}_{0}\in\mathbb{R}^{r_{A}\times r}$ and clearly $\mathbf{\Phi}_{0}=\mathbf{\Sigma}\mathbf{V}^{\top}\mathbf{\Omega}$. Using the expression of $\mathbf{X}_{1}$, iteration (35) gives that

$$ $\displaystyle\mathbf{Y}_{1}=\eta\mathbf{V}\mathbf{\Sigma}\mathbf{\Phi}_{0}( \mathbf{\Phi}_{0}^{\top}\mathbf{\Phi}_{0})^{\dagger}.$ $$

Let the compact SVD of $\mathbf{\Phi}_{0}:=\mathbf{P}\mathbf{D}\mathbf{Q}^{\top}$, where $\mathbf{P}\in\mathbb{R}^{r_{A}\times r_{A}}$, $\mathbf{D}\in\mathbb{R}^{r_{A}\times r_{A}}$ and $\mathbf{Q}\in\mathbb{R}^{r\times r_{A}}$. Note that $\mathbf{P}$ is unitary. With the compact SVD of $\mathbf{\Phi}_{0}$, we have that $(\mathbf{\Phi}_{0}^{\top}\mathbf{\Phi}_{0})^{\dagger}=\mathbf{Q}\mathbf{D}^{-2
}\mathbf{Q}^{\top}$, which implies that

$$ $\displaystyle\mathbf{X}_{1}\mathbf{Y}_{1}^{\top}=\eta\mathbf{U}\mathbf{P} \mathbf{D}\mathbf{Q}^{\top}\mathbf{Q}\mathbf{D}^{-2}\mathbf{Q}^{\top}\mathbf{Q }\mathbf{D}\mathbf{P}^{\top}\mathbf{\Sigma}\mathbf{V}^{\top}\stackrel{{ \scriptstyle(a)}}{{=}}\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}=\mathbf{A}$ $$

where (a) is because $\mathbf{P}$ is unitary and the choice of $\eta=1$.
∎

###### Theorem 6 .

###### Proof.

## Appendix D Other useful lemmas

###### Lemma 14 .

Let $A_{t+1}=(1-\theta)A_{t}+\beta$ with some $\alpha\in(0,1)$ and $\beta\geq 0$, then we have

$$ $\displaystyle A_{t+1}=(1-\theta)^{t+1}A_{0}+\beta\frac{1-(1-\theta)^{t+1}}{ \theta}\leq(1-\theta)^{t+1}A_{0}+\frac{\beta}{\theta}.$ $$

###### Proof.

The proof can be completed by simply unrolling $A_{t+1}$ and using the fact $1+\alpha+\alpha^{2}+\ldots+\alpha^{t}\leq\frac{1}{1-\alpha}$. ∎

###### Lemma 15 .

If $\mathbf{A}\in\mathbb{R}^{n\times n}$ and $\mathbf{B}\in\mathbb{R}^{n\times n}$ are positive semi-definite matrices, we have $\lambda_{\min}(\mathbf{A}+\mathbf{B})\geq\lambda_{\min}(\mathbf{A})+\lambda_{
\min}(\mathbf{B})$.

###### Proof.

The smallest eigenvalue of $\mathbf{A}+\mathbf{B}$ can be expressed as

$$ $\displaystyle\lambda_{\min}(\mathbf{A}+\mathbf{B})=\min_{\mathbf{x}\neq\mathbf {0}}\frac{\mathbf{x}^{\top}(\mathbf{A}+\mathbf{B})\mathbf{x}}{\mathbf{x}^{\top }\mathbf{x}}=\min_{\mathbf{x}_{1}\neq\mathbf{0},\mathbf{x}_{1}=\mathbf{x}_{2}} \frac{\mathbf{x}_{1}^{\top}\mathbf{A}\mathbf{x}_{1}}{\mathbf{x}_{1}^{\top} \mathbf{x}_{1}}+\frac{\mathbf{x}_{2}^{\top}\mathbf{B}\mathbf{x}_{2}}{\mathbf{x }_{2}^{\top}\mathbf{x}_{2}}.$ (36) $$

On the other hand, we also have that

$$ $\displaystyle\lambda_{\min}(\mathbf{A})+\lambda_{\min}(\mathbf{B})=\min_{ \mathbf{x}_{1}\neq\mathbf{0},\mathbf{x}_{2}\neq\mathbf{0}}\frac{\mathbf{x}_{1} ^{\top}\mathbf{A}\mathbf{x}_{1}}{\mathbf{x}_{1}^{\top}\mathbf{x}_{1}}+\frac{ \mathbf{x}_{2}^{\top}\mathbf{B}\mathbf{x}_{2}}{\mathbf{x}_{2}^{\top}\mathbf{x} _{2}}.$ (37) $$

Because (36) is a constrained version of the minimization problem (37), they share the same objective, but (36) has shrinked feasible region. It is not difficult to see that $\lambda_{\min}(\mathbf{A}+\mathbf{B})\geq\lambda_{\min}(\mathbf{A})+\lambda_{
\min}(\mathbf{B})$. The proof is thus completed.
∎

###### Lemma 16 .

Consider a sequence $\{A_{t}\}_{t}$ with $A_{t}\geq 0,\forall t$. If there exists $\alpha$ such that $A_{t+1}\leq\alpha A_{t}^{2}$ and $A_{0}\leq\frac{1}{2\alpha}$, $A_{t}$ converges to 0 0 at a quadratic rate, i.e.,

$$ $\displaystyle A_{t+1}\leq\frac{1}{\alpha}\frac{1}{2^{2^{t+1}}}.$ $$

###### Proof.

Unrolling $A_{t+1}$, we get that

$$ $\displaystyle A_{t+1}\leq\alpha A_{t}^{2}\leq\alpha^{3}A_{t-1}^{4}\leq\alpha^{ 7}A_{t-2}^{8}\leq\frac{1}{\alpha}(\alpha A_{0})^{2^{t+1}}\leq\frac{1}{\alpha} \frac{1}{2^{2^{t+1}}}.$ $$

The proof is thus completed.
∎

###### Lemma 17 .

Let $\mathbf{A}\in\mathbb{R}^{m\times n}$ be a matrix with full column rank and $\mathbf{B}\in\mathbb{R}^{n\times p}$ be a non-zero matrix. Let $\sigma_{\min}(\cdot)$ be the smallest non-zero singular value. Then it holds that $\sigma_{\min}(\mathbf{A}\mathbf{B})\geq\sigma_{\min}(\mathbf{A})\sigma_{\min}(
\mathbf{B})$.

###### Proof.

Using the min-max principle for singular values,

$$ $\displaystyle\sigma_{\min}(\mathbf{A}\mathbf{B})$ $\displaystyle=\min_{\|\mathbf{x}\|=1,\mathbf{x}\in\text{ColSpan}(\mathbf{B})} \|\mathbf{A}\mathbf{B}\mathbf{x}\|$ $\displaystyle=\min_{\|\mathbf{x}\|=1,\mathbf{x}\in\text{ColSpan}(\mathbf{B})} \Big{\|}\mathbf{A}\frac{\mathbf{B}\mathbf{x}}{\|\mathbf{B}\mathbf{x}\|}\Big{\| }\cdot\|\mathbf{B}\mathbf{x}\|$ $\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\min_{\|\mathbf{x}\|=1,\|\mathbf {y}\|=1,\mathbf{x}\in\text{ColSpan}(\mathbf{B}),\mathbf{y}\in\text{ColSpan}( \mathbf{B})}\|\mathbf{A}\mathbf{y}\|\cdot\|\mathbf{B}\mathbf{x}\|$ $\displaystyle\geq\min_{\|\mathbf{y}\|=1,\mathbf{y}\in\text{ColSpan}(\mathbf{B} )}\|\mathbf{A}\mathbf{y}\|\cdot\min_{\|\mathbf{x}\|=1,\mathbf{x}\in\text{ ColSpan}(\mathbf{B})}\|\mathbf{B}\mathbf{x}\|$ $\displaystyle\geq\min_{\|\mathbf{y}\|=1}\|\mathbf{A}\mathbf{y}\|\cdot\min_{\| \mathbf{x}\|=1,\mathbf{x}\in\text{ColSpan}(\mathbf{B})}\|\mathbf{B}\mathbf{x}\|$ $\displaystyle=\sigma_{\min}(\mathbf{A})\sigma_{\min}(\mathbf{B})$ $$

where (a) is by changing of variables, i.e., $\mathbf{y}=\mathbf{B}\mathbf{x}/\|\mathbf{B}\mathbf{x}\|$.
∎

###### Lemma 18 .

For PSD matrices $\mathbf{A}$ and $\mathbf{B}$, if $\mathbf{A}+\mathbf{B}=\mathbf{I}_{r}$, then we have $\text{Tr}(\mathbf{A})\leq r$ and $\text{Tr}(\mathbf{B})\leq r$.

###### Proof.

The proof is straightforward and is omitted here. ∎

###### Lemma 19 ( Rudelson and Vershynin ( 2009 ) ) .

Let $\bf{W}$ be an $d\times r$ matrix with $d\geq r$. The entries of $\bf{W}$ are drawn independently from ${\cal N}(0,1)$. Then for every $\tau>0$, we have that

$$ $\displaystyle\mathbb{P}\big{(}\sigma_{r}({\bf W})\leq\tau(\sqrt{d}-\sqrt{r-1}) \big{)}\leq(C_{1}\tau)^{d-r+1}+e^{-C_{2}d}.$ $$

where $C_{1}$ and $C_{2}$ are universal constants independent of $d$ and $r$.

###### Lemma 14 .

###### Proof.

###### Lemma 15 .

###### Proof.

###### Lemma 16 .

###### Proof.

###### Lemma 17 .

###### Proof.

###### Lemma 18 .

###### Proof.

###### Lemma 19 ( Rudelson and Vershynin ( 2009 ) ) .

## Appendix E Missing experimental details

### E.1 Details for problems with synthetic data

This subsection contains the detailed setup for the problems with synthetic data in Figs. 1 and 4. Recall that here we focus on symmetric problems under exact-, under-, and over-parametrization.

For the exact-parametrized problem in Fig. 1 (a) and (b), we choose the PSD matrix $\mathbf{A}\in\mathbb{R}^{m\times m}$ in the following manner. We set $m=1000$ and $r=r_{A}=20$. The non-zero singular values are set as $\{1.0,0.99,0.98,\ldots,0.82,0.01\}$, where we intentionally set $\sigma_{r_{A}}=0.01$ to enlarge the condition number. We choose the step size of GD as $0.01$ to avoid divergence. The learning rate for ScaledGD is $0.5$.

For the under-parametrized problem in Fig. 1 (c), we choose PSD matrix $\mathbf{A}\in\mathbb{R}^{m\times m}$ in the following manner. We set $m=1000$ and $r_{A}=40$. The singular values of $\mathbf{A}$ are $\{1.0,0.99,0.98,\ldots,0.65,0.64,0.05,$ $0.025,0.01\}$. We choose $r=20$ to ensure the under-parametrized nature of this problem.

For the over-parametrized case in Fig. 4 (a) and (b), we choose PSD matrix $\mathbf{A}\in\mathbb{R}^{m\times m}$ in the following manner. We set $m=1000$ and $r_{A}=20$. The non-zero singular values are chosen as $\{1.0,0.99,0.98,\ldots,0.82,0.01\}$, where we intentionally set $\sigma_{r_{A}}=0.01$ to enlarge the condition number. We set $\mathbf{X}$ to be over-parametrized by letting $r=60$.
We choose the step size of GD as $0.01$. The learning rate of ScaledGD-$\lambda$ is set as 0.5, and its damping parameter $\lambda$ is chosen as $0.01$. The learning rate for ScaledGD with Nyström initialization is $0.5$.

Figure: Figure 5: The dog dataset.
Refer to caption: x13.png

Figure: Figure 6: The cat-toy dataset.
Refer to caption: x14.png

### E.2 Datasets

The evaluation of NoRA and NoRA+ is carried out on commonly adopted datasets in the literature.

GLUE benchmark. GLUE is designed to provide general-purpose evaluation of language
understanding (Wang et al., ). Those adopted in our work include SST-2 (sentiment analysis, (Socher et al., )), RTE(^4^44[https://paperswithcode.com/dataset/rte](https://paperswithcode.com/dataset/rte)) (inference). These datasets are released under different permissive licenses.

SuperGLUE benchmark. SuperGLUE (Wang et al., ) is another commonly adopted benchmark for language understanding, and it is more challenging compared with GLUE. The considered datasets include CB (inference, (De Marneffe et al., )), ReCoRD (question answering, (Zhang et al., )), WSC (coreference resolution, (Levesque et al., )), BoolQ (question answering, (Clark et al., )), and MiltiRC (question answering, (Khashabi et al., )). These datasets are released under different permissive licenses.

Commonsense reasoning. These datasets are a collection tasks that require commonsense reasoning to answer. The considered datasets include WinoGrande (Sakaguchi et al., ), PIQA (Bisk et al., ), SOCIAL-I-QA (SIQA) (Sap et al., ), HellaSwag (Zellers et al., ), ARC-easy, ARC-challenge (Chollet, ) and OpenbookQA (Mihaylov et al., ). These datasets are released under different permissive licenses.

Math. For mathematical problems, we consider GSM8K (Cobbe et al., ) dataset that consists
of high quality linguistically diverse school math problems created by human problem writers. This dataset is under MIT license. We also adopt MetaMathQA dataset (Yu et al., ), which is constructed through bootstrapping mathematical questions by rewriting the question from multiple perspectives. This dataset is under MIT license.

Additional datasets. We also use SQuAD (question answering, (Rajpurkar et al., )) in our experiments, which is released under license CC BY-SA 4.0.

Datasets for DreamBooth. The datasets (dog and cat-toy) used for Sec. 5.2 are obtained directly from Huggingface. The dog dataset(^5^55[https://huggingface.co/datasets/diffusers/dog-example](https://huggingface.co/datasets/diffusers/dog-example)) contains 5 dog images; see Fig. 5. The cat-toy(^6^66[https://huggingface.co/datasets/diffusers/cat-toy-example](https://huggingface.co/datasets/diffusers/cat-toy-example)) dataset has 4 images; see Fig. 6. Both datasets are representative examples for the purpose of DreamBooth – finetuning with only few images for personalized generalization.

### E.3 Details for Fig. 2

The experiment setting and training protocols are the same as few-shot learning with OPT-1.3B in the following subsection. Here, we are interested in the change of singular values after LoRA finetuning. For each LoRA layer, we compare the singular values of $\mathbf{W}_{0}$ and $\mathbf{W}_{0}+\mathbf{X}_{T}\mathbf{Y}_{T}^{\top}$, where $\mathbf{X}_{T},\mathbf{Y}_{T}$ are LoRA weights after training, and find out the indices of $r$ singular values that have the largest change after finetuning. We then count the indices across all LoRA layers. Fig. 2 plots indices vs. counts.

### E.4 Few-shot learning with OPT-1.3B

For this experiment, we first search for the best batchsizes for LoRA, and the same batchsize is applied for other tested algorithms as well. Then we search additionally for the best learning rate for each algorithm. This ensures that different algorithms see the same amount of data, while still having their best performed learning rate. The hyperparameters adopted are searched over values in Tab. 7. Adam is adopted for optimization.

**Table 7: Hyperparameters used for few-shot learning with OPT-1.3B.**
| Hyperparameters | Values |
| --- | --- |
| LoRA $r$ | 8 |
| LoRA $\alpha$ | 16 |
| LoRA module | q_proj, v_proj |
| # epochs | 5 |
| batchsize | 2, 4, 8 |
| learning rate | 1$\times 10^{-5}$, 5$\times 10^{-5}$, 1$\times 10^{-4}$ |
| NoRA $\xi$ | 0.05, 0.1, 0.2 |

### E.5 DreamBooth with stable-diffusion

Stable Diffusion V1.4 (Rombach et al., ) is adopted as base model, where LoRA is applied to the UNet. The text-encoder is not finetuned. We adopt the default parameter-choice from Huggingface, which is summarized in Tab. 8. We adopt AdamW as the optimizer with a weight decay of $0.01$.

**Table 8: Hyperparameters used for DreamBooth with stable-diffusion.**
| Hyperparameters | Values |
| --- | --- |
| LoRA $r$ | 4 |
| LoRA $\alpha$ | 4 |
| LoRA module | to_q, to_k, to_v, to_out |
| # iterations | 500 |
| batchsize | 1 |
| learning rate | 1$\times 10^{-4}$ |
| NoRA $\xi$ | 0.1 |

We provide additional results to further support the efficiency of NoRA by finetuning the stable-diffusion-v1.4 model using the same protocol as in Sec. 5.2. Here we adopt a dataset with $4$ toy-cat images; see Fig. 6. After finetuning $500$ steps using prompt “a photo of toy cat”, our goal is to generate images “a toy cat wearing glasses.” The generated images are shown in Fig. 7. In general, all tested algorithms do not distinguish the hands and the tail of toy cat well. However, both LoRA and LoRA-P generate images with less accurate facial details. For example, the glasses are not wearing well, or the eyes are not clear. However, the details of faces generated by NoRA and NoRA+ are quite clear.

Figure: Figure 7: Generated images from NoRA and NoRA+ with stable-diffusion.
Refer to caption: x15.png

### E.6 Commonsense reasoning with LLaMA and LLaMA2

The base models considered are LLaMA-7B and LLaMA2-7B. The experimental setup and choices of hyperparameters follow (Liu et al., ). The hyperparameters are summarized in Tab. 9.

**Table 9: Hyperparameters used for commonsense reasoning with LLaMA-7B and LLaMA2-7B.**
| Hyper-parameters | Values |
| --- | --- |
| LoRA $r$ (rank) | 32 |
| LoRA $\alpha$ | 64 |
| LoRA module | q_proj, k_proj, v_proj, up_proj, down_proj |
| epoch | 3 |
| learning rate | $3\times 10^{-4}$ |
| batchsize | 16 |
| cutoff length | 256 |
| NoRA $\xi$ | 0.02, 0.05, 0.1 |

### E.7 Math reasoning with Gemma-7B

Our last evaluation tackles mathematical reasoning. Gemma-7B (Gemma-team et al., ) is finetuned for 2 epochs on MetaMathQA-100K dataset (Yu et al., ). LoRA rank is set as 32, leading to 100M trainable parameters. The performance is assessed on GSM8K (Cobbe et al., ), and hyperparameters are summarized in Tab. 10.

**Table 10: Hyperparameters used for math reasoning with Gemma-7B.**
| Hyper-parameters | Values |
| --- | --- |
| LoRA $r$ (rank) | 32 |
| LoRA $\alpha$ | 64 |
| LoRA module | q_proj, k_proj, v_proj, o_proj, up_proj, down_proj, gate_proj |
| epoch | 2 |
| learning rate | $3\times 10^{-4}$, $4\times 10^{-4}$, $5\times 10^{-4}$ |
| batchsize | 128 |
| NoRA $\xi$ | 0.02, 0.05, 0.1 |

## Peer Reviews

confidence: 3
contribution: 4
initial_rating: 5
presentation: 4
questions: .
soundness: 4
strengths: On one hand, I feel that the quadratic convergence with the Nyström initialization is significant. I have not directly worked in the area of matrix completion, so I do not have full confidence in this judgment, but this type of quadratic convergence result seems quite important.

At the same time, however, I'm not quite sure if the quadratic convergence is coming from the scaled GD algorithm or from the Nyström initialization. Since the scaled GD algorithm is basically a Newton update, shouldn't the quadratic convergence happen asymptotically, regardless of the initialization? Is it that the Nyström initialization allows the algorithm to reach the quadratically convergent phase more quickly? Perhaps the authors can clarify this point.
summary: The paper shows that the low-rank matrix factorization problem can be significantly accelerated using a better initialization method. Specifically, the authors introduce the Nyström initialization, which uses a sketching idea to initialize the low-rank factors to be in the appropriate eigenspaces. Then, using Scaled Gradient Descent (ScaledGD), one can obtain quadratic convergence.

The authors then present Nyström-initialized LoRA (NoRA) and present experimental results that show the practical benefits of using such an initialization compared to using the standard LoRA.
weaknesses: On the other hand, I feel that the matrix completion result and the Nyström-initialized LoRA (NoRA) are strongly connected. The theory doesn't apply to the NoRA setting, but I do see that the Nyström initialization is a nice conceptual motivation for a potentially better LoRA initialization.

However, now with the flood of empirical LoRA papers, I think the community has a very high bar for recognizing work that purports to improve LoRA in a practical sense. I am not sure if the experimental validations of this work really clear that bar.

---

confidence: 5
contribution: 2
initial_rating: 3
presentation: 2
questions: See "Weaknesses"
soundness: 3
strengths: Good coverage on matrix factorization problems (symmetric/asymmetric, exact/over/under-parametrization); Writings are clear; Some connection to practical problems.
summary: This paper studies the convergence of ScaledGD, a special preconditioned GD for matrix factorization (MF) problems, on symmetric/asymmetric MF problems. The authors show that initializing the factorized model as some random sketching of the target matrix significantly achieves faster convergence. Then, the authors propose to initialize the Low-rank adaptation (LoRA) using some random sketching of the pre-trained weights.
weaknesses: The reviewer has concerns regarding the significance and practical relevance of the theoretical results and the way results are connected to LoRA. Specifically:

1. The matrix factorization problem is not a practical problem to be solved, but rather a problem through which one builds a theoretical understanding of GD and its variants in a nonconvex setting. In this regard, any procedure using the target matrix itself, in the reviewer's opinion, is not allowed: if we know the target matrix, we can directly factorize it (which is why the MF problem is not a practical problem), thus any procedure (other than computing gradient, we are studying GD after all) using some information about the true target matrix is not allowed as it simply moves the factorized model one step closer to the target thus reduce the amount of work GD needs to converge. As such, the crucial role of initialization, as the authors stated, is just taking advantage of knowing the target matrix, which deviates from the original purpose of studying the MF problem: to understand how GD minimizes non-convex objectives even with random initialization. (Disclaimer: This point only pertains to MF problems. For practical problems such as matrix sensing, and matrix completion, the initializations can be special as they are non-trivially improving the convergence of GD)

2. as the reviewer already pointed out in 1., if we know the target matrix, then we can directly factorize it, which is exactly what the one-step modified ScaledGD with Nystrom initialization is doing, so Section 3 is trivial.

3. The connection to LoRA is puzzling in two ways: First, what is proposed is not remotely close to what is theoretically studied, Section 2, 3 studies initialization with the targe matrix sketched. The proposed initialization for LoRA is with the pre-trained weights sketched. Based on this, one should get the conclusion that "pre-trained weights" $\approx$ "targe matrix", which is clearly wrong. (Don't get too serious about this argument, the reviewer is just taking the shortcut to make the point). Moreover, the theoretical claims are special initialization accelerates GD, but none of the empirical experiments are about training time. On top of that, the performance improvement from special initialization seems very little to the reviewer.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. In Figure 2b, it kind of suggests that smaller $\xi$ is preferred for quicker convergence? So does it hint that Nystrom initialization benefits from small initializaion as well, which also induces zero saddle point concerns?

2. In NoRA, say we have two pretrained weights $W_0,W_1$ of pretty different magnitudes, does it indicate corresponding NoRA initialization also scales with such magnitudes since $X_\{0,1\}=W_\{0,1\}\Omega$? Probably adding certain weight normalization would make the finetuning more stable?


minor writing issues:
line 340: "modified ScaledGD in (6) have" should be "has"
soundness: 3
strengths: The current work studies how initialization affects ScaledGD method in matrix factorization problems, and considers its extension to LoRA model. I personally find the topic valuable. The theoretic bound improvement seems huge for classic matrix factorization problem. For LoRA model, recent work has studied customized optimizers and stepsize scheduler, NoRA proposes a new initialization scheme which has its theoretic grounds. Moreover, compared to PiSSA, NoRA is naturally zero-initialized and is thus a more convenient initialization to LoRA.
summary: The current work describes a novel initialization method for matrix factorization problems which accelerates convergence of ScaledGD both theoretically and empirically. The authors also consider extension of Nystrom initialization to LoRA model and shows its superiority compared to different baselines experimentally.
weaknesses: 1. The convergence bound to under-parameterized case involves weak optimality. I'm not very familiar with this notion, the authors show all globally optimal solutions are weakly optimal, is there any intuition of this weak optimality or is it just curated from technical stuff? How larger is the space of weak optimality compared to true optimality?

2. Ideally to apply Nystrom initialization to LoRA, one should learn $\Delta W$, which suggests that one should know how much weight change is in order to initialize finetuning model, and is thus unrealistic.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: What is the role of Gaussian playing in your analysis? To put it in other words, do you think it's possible to replace Gaussian with something that can be applied faster to $A$, or Gaussian is essential for the analysis?
soundness: 3
strengths: This paper presents a simple initialization framework for matrix factorization that significantly improves convergence rates, which is to align the principle directions of $X_0$ with $A$ through a rank-$r$ Gaussian. Algorithmically, this does not incur (asymptotically) extra cost as multiplying $A$ with a small matrix is cheap. It also showcases an interesting phenomenon: in the exact parametrization setting, the proposed algorithm has a linear convergence for iterations proportion to $\kappa^3$, then rapidly switches to quadratic convergence, resembles that of second-order method without explicitly resorting to second-order information. The application to LoRA is nice and natural.
summary: This paper considers the matrix factorization problem, where one is given a matrix $A$ and the goal is to compute a rank-$r$ factorization that approximates $A$. It studies the setting when $A$ is PSD and the factorization is also PSD, and when $A$ is not symmetric. Three settings are studied: exact, over and under parametrization, meaning the rank of $A$ is equal, smaller or larger than $r$. Their algorithm is preconditioned gradient, and the key of this paper it show that a good initialization significantly improves convergence rate. The initialization is to set $X_0=A \Omega$ for a zero-mean Gaussian matrix $\Omega$. Due to rotational invariance of Gaussian, this initialization has many desired property, such as it aligns with the desired top-$r$ eigenspace. Using this technique, they obtain improved rates for both symmetric matrix factorization and asymmetric variants. They also show a major application for this theory, that is, the initialization for LoRA training. Standard LoRA setting initializes $X_0$ to random Gaussian matrix, and the theory developed in this paper suggests a better alternative is to use $X_0=W_0 \Omega$ where $W_0$ is the pre-trained weights. Experiments are performed on various few-shot learning tasks on benchmarks to show advantages over LoRA and its variants.
weaknesses: I think the paper is overall well-written with no obvious weaknesses. The empirical improvements over LoRA seem rather marginal but I'm not an expert on experiments.

Meta-review: This paper presents a Nyström initialization scheme for low-rank matrix factorization. The authors demonstrate that this initialization is critical for improving the convergence rates of the preconditioned ScaledGD algorithm in three different settings: (a) under-parameterized, (b) exactly-parameterized, and (c) over-parameterized. The paper addresses both symmetric and asymmetric matrices. Beyond deriving convergence rates for matrix factorization, the authors also explore applying the Nyström initialization in the LoRA setting.

The reviewers recognized the paper’s technical soundness, clarity of presentation, and overall contributions. During the rebuttal phase, the authors addressed most of the reviewers’ concerns, including the significance of the theoretical results in the non-symmetric case, the applicability of the results to more general problems (such as matrix sensing), and the relatively marginal improvements observed in the LoRA experiments. Although some issues remain and the paper could be further improved, I believe that (a) the overall quality of the work, (b) the insights into the importance of initialization for convergence rates, and (c) the contributions made are all clearly above the acceptance threshold. I therefore recommend its acceptance to ICLR.

Decision: Accept (Poster)