Gradient Descent Dynamics of Rank-One Matrix Denoising
Authors: Zeyan Zhuang
Shenghui Song
Venue: ICLR 2026 Poster (2026)
Keywords: Random Matrix Theory
High-Dimensional Statistics
Matrix Denoising
Gradient Flow
Primary area: learning theory
TL;DR: Derive a closed-form solution for the learning dynamics of GD-based rank-one matrix denoising and reveal the BBP transition in the large-time limit.

Abstract: Matrix denoising is a crucial component in machine learning, offering valuable insights into the behavior of learning algorithms  (Bishop and Nasrabadi, 2006). This paper focuses on the rectangular matrix denoising problem, which involves estimating the left and right singular vectors of a rank-one matrix that is corrupted by additive noise. Traditional algorithms for this problem often exhibit high computational complexity, leading to the widespread use of gradient descent (GD)-based estimation methods with a quadratic cost function. However, the learning dynamics of these GD-based methods, particularly the analytical solutions that describe their exact trajectories, have been largely overlooked in existing literature. To fill this gap, we investigate the learning dynamics in detail, providing convergence proofs and asymptotic analysis. By leveraging tools from large random matrix theory, we derive a closed-form solution for the learning dynamics, characterized by the inner products of the estimates and the ground truth vectors. We rigorously prove the almost sure convergence of these dynamics as the signal dimensions tend to infinity. Additionally, we analyze the asymptotic behavior of the learning dynamics in the large-time limit, which aligns with the well-known Baik-Ben Arous-Péchée phase transition phenomenon n (Baik et al., 2005). Experimental results support our theoretical findings, demonstrating that when the signal-to-noise ratio (SNR) surpasses a critical threshold, learning converges rapidly from an initial value close to the stationary point. In contrast, estimation becomes infeasible when the ratio of the inner products between the initial left and right vectors and their corresponding ground truth vectors reaches a specific value, which depends on both the SNR and the data dimensions.

Published as a conference paper at ICLR 2026
GRADIENT DESCENT DYNAMICS OF RANK-ONE MA- TRIX DENOISING
Zeyan Zhuang & Shenghui Song
Department of Electronic and Computer Engineering The Hong Kong University of Science and Technology 2 ezhuangaceconnect wtrpoonnnnttttttt... uust.k, eeshsong@ust. ust. hk
ABSTRACT
Matrix denoising is crucial component in machine learning, offering valuable in- sights into the behavior of learning algorithms (Bishop & Nasrabadi] (206). This paper focuses the rectangular matrix denoising problem, which involves esti- mating the left and right singular vectors of a rank-one matrix that is corrupted by additive noise. Traditional algorithms for this problem often exhibit high compu- tational complexity, leading to the widespread use of gradient descent estimation methods with a quadratic cost function. However, the learning dynam- ics of these GD-based methods, particularly the analytical solutions that describe their exact trajectories, have been largely overlooked existing literature. To fill this gap, we investigate the learning dynamics in detail, providing convergence proofs and asymptotic analysis. By leveraging tools from large random matrix theory, derive closed-form solution for the learning dynamics, characterized by the inner products of the estimates and the ground truth vectors. We rigorously prove the almost sure convergence of these dynamics the signal dimensions tend to infinity. Additionally, we analyze the asymptotic behavior of the learn- ing dynamics the large-time limit, which aligns with the well-known Baik-Ben Arous-Péchhe phase transition phenomenon (Baik et al. 2005). Experimental re- sults support our theoretical findings, demonstrating that when the signal-to-nois to-sooiee ratio (SNR) surpasses critical threshold, learning converges rapidly from an ini- tial value close the stationary point. In contrast, estimation becomes infeasible when the ratio of the inner products between the initial left and right vectors and their corresponding ground truth vectors reaches a specific value, which depends on both the SNR and the data dimensions.
INTRODUCTION
Matrix denoising involves recovering signal matrix P RRxx from noisy observation X
P + Z, which a fundamental challenge in statistics with broad applications in image processing (Pedersen et al. 2009 Gorderr-Grrand
(219)) genomics (Leek) 2011). wireless communica- tions (Couillet & Hachem 2013), and other fields. This model is typically referred to as Johnstone's spiked model (Johnstone Paul| 2018), when the noise Z random matrix whose dimensions p, are large and comparable. while P is deterministic matrix with rank & min (p. n2). Extensive research (Ding& Yang 202T, Couillet Liao, 2022; [ao al., 2022 Liu al., (025a) has shown that in higs-dim settings, the left and right singular vectors of X corresponding to its top singular values could be utilized as the estimate for P.
2013)
However, when the dimensions p and are very large, the singular value decomposition (SVD) of X suffers from intolerable storage and computational overhead. To address this issue, iterative optimization and learning methods have been developed (Björck et al. 2015). among which gradient descent (GD) plays a crucial role not only in matrix denoising but also many machine learning problems. In particular, understanding the dynamics of GD is essential in explaining the remarkable performance of today's 's deep neural networks. Inspired by this, we aim to analyze the learning dynamics of the GDD-baedd rank-one signal estimation algorithm in this work.

---

Published as a conference paper at ICLR 2026
1.1 RELATED WORKS
Research on matrix denoising typically studies low-rank deformed random matrix models. Two fun- damental and widely studied variants are the spiked Wigner model "Benachh-George Nadakuditi 2011 and the spiked Wishart model Sennaych-GGorger & Nadakuditi] 2012). The seminal work by Baik, Ben Arous, and Péché (BBP) (Baik et 2005)) revealed that the extreme eigenvalues and as- sociated eigenvectors of X undergo BBP phase transition as the signal-to- to-ioe ratio (SNR) varies.
This pioneering result has since been generalized to a wide range of statistical models for signal esti-
mation, including those with correlated noise (Zhang & Mondelli 2024] Ding & Yang 2022 Gavish et al. 2023). random features (Liao & Couillet 2018). puncturing (Couillet et al., 2021). random projections Yang et al. 2021), among others. A key commonality across these studies is the use of extreme eigenvalues and eigenvectors for signal recovery.
The aforementioned studies primarily focus on the "static" or asymptotic performance of matrix denoising. Iterative algorithms for such low-rank signal recovery problems have also been estab- lished. Under the information-theereti framework (Korada & Macris 2009; Lelarge & Miolane 2017]. it can be shown that approximate message passing (AMP) can achieve minimal mean-square- error, and the corresponding SNR boundary has been determined (Barbier al. 2016 Lesieur al.) 2017 However, AMP generally requires specific prior distributions and is often tailored to partic- ular problem settings.
Another effective iterative algorithm is GD. Unlike AMP, GD does not rely on prior information and can be applied to a broad range problems (Pretorius et al. 2018). The dynamics of GD are also crucial for understanding machine learning processes. Nevertheless, results on the learning dynamics GD-based estimation are very limited. The most relevant work in the literature is (Bodin& & Macris 2021). where the authors studied the rank-one matrix denoising problem for the deformed Wigner model. In particular, a closed-fone form solution for the learning dynamics is obtained for the case where P and Z are symmetric. However, the symmetric structure for the noise may not be available in practice, where the observed data typically consists of a sequence {xx }, and the data matrix X [æ1,.... ... Fr,] is rectangular (the deformed Wishart model). Unfortunately, the overall learning dynamics for this deformed Wishart model with general structured P and Z have not been fully understood in the literature. This work aims to fill this research gap.
From a technical perspective, obtaining analytical solutions for the dynamical behavior of the Wishart model is significantly more challenging than that for the Wigner model, due its struc- tural asymmetry and the higher order of the governing differential equations. To overcome these difficulties, we employ methods from large random matrix theory (RMT) construct approximate solutions. For a comprehensive treatment of the relevant techniques, we refer the readers to (Bai
al., 2010 Pastur & Shcherbina 2011 Yao al., 2015) Erdōs & Yau 2017] Vershynin 2018 Couillet & Liao] 2022; Tao 2023). Our approach was specifically motivated by the almost sure
boundedness of extreme eigenvalues of random matrices (Yin et al. 1988; Bai & Yin 1993) and
strong convergence results for the resolvents (Bai & Silverstein
1998
1999.
1.2
CONTRIBUTIONS
The main contributions of this work are summarized as follows.
We obtain the deterministic approximations for the evolution of the inner products between the learned vectors and the ground truth, which are in closed-form. Moreover, we prove that, as the dimensions of the observation matrix approach infinity, the empirical evolution processes will converge almost surely to the deterministic approximations.
The learning dynamics are described and analyzed through a set of differential equations that involve complex variables and contour integral conditions. We investigate the analyt- ical properties of these equations. Specifically, we demonstrate that the solution admits integral representations with respect to certain class of transition kernels. Furthermore, we establish the existence and uniqueness of the solutions for the concerned equations.
We investigate the asymptotic behavior of the learning dynamics in the large-time limit. It is observed that there exists a critical threshold, below which the estimation becomes quite challenging. This phenomenon is analogous to the well-known BBP phenomenon. How- ever, the BBP phenomenon concerns the square of the inner product, which is unsigned, but we derive a signed result and reveal its relationship with the initial conditions.
2

---

Published as a conference paper at ICLR 2026
1.3 NOTATIONS AND ORGANIZATION
Throughout the paper, lowercase and uppercase boldface letters represent vectors and matrices, respectively. In represents the identity matrix of size n. For vectors and b, their inner product is defined as (a, b) = a b. Let 1 and R, denote the sets 2 = {ir : x 2 0} and I. {i z > 0}, respectively. The notation C(AA) represents the set of continuous functions defined on A. The norms 2- Il- lllleca) and Illlly denote the fn-norm for vectors, Frobenius norm for matrices, supremum norm of continuous functions on A, and the total variation of signed measures, respectively. The indicator function denoted as 1 [ and the Dirac measure at a point x is denoted as 6(x). Notation (f * (t) S f(t - a)g(a)da represents the convolution of f(t) and g().. The function (x)+ = max(x, 0). The sign function is denoted as Sgn(z) - 1 [ 0} - 1fx 0} and R.S. represents almost sure convergence.
Notations:
Organization: The rest of the paper is organized as follows. Section we introduce the problem. In Section 3 we state the main results, including the closed-form deterministic approximation for the learning dynamics and the asymptotic behavior with long-term learning. The implications and potential applications of the main results are also discussed. Experiment results are given in Section and Section concludes the paper.
PROBLEM SETUP
We consider the following rank-one Jonstone's 's spiked model
X = Z x----uu E R
(I)
where X is the observation and Z = (Zt,) denotes the random noise matrix that contains indepen- dent and identically distributed (i.i..) elements with mean zero and variance n-1, The unit vectors " Rº and v² E R represent the directions of the left and right ground truth signals, respec- tively. The parameter E (0, x)) denotes the SNR (Zhang & Mondelli, 2024). The main target is
to estimate the signal components (u', v²) from the observation matrix X in high dimensions. To this end, we consider the following cost function (Nadakuditi 2014)
H(u,, =
X ue
According to the Eckart- Young-Mirss (EYM) theorem (Eckart & Young 1936), the optimal solu- tion to the optimization problem
(ueym Veym)
arg min H(u, v), ||
(3)
given by (Maym Noym) = (ū,, 11), where ūn and 0 denote the left and right singular vectors of X corresponding to the largest singular value.
When both the dimensions p and n are very large, iterative methods such as power iteration and GD-based optimization are utilized to compute the top singular vectors (Martinsson & Tropp 2020
While power iteration exhibits rapid convergence under high SNR conditions (Wu & Zhou 2024 its effectiveness is often problem-specifis and relies favorable spectral gaps. In contrast, GD offers wider applicability and greater flexibility across diverse problem structures. Given the initial point (uo, vo) with ||æal]2 || || 1, when the step size is small, the discrete process of GD can be approximated by the gradient flow (Li et al. 2017 Liu 2017), given by
[e di vt
[proju, (Vu.uu. ux)) proje, (V (VoHluu, 0x))
= -gggggaa(a,,, vx))
(4)
where grad(-) represents the Riemannian gradient operator (Gess et al. 2024 and proj (v) =
(I xx xx= )y denotes the projection onto the tangent space. We note that (4) enforces the unit norm constraint on u, and U.. and the detailed proof is given in Appendix C During the learning process, the inner products between the ground truth and the estimates, i.e., 4r (1) (u", u.) and
(t) = (v", vz), are of interest. In particular, by ||=" x|} = 2 (xx x²), the inner product is equivalent to the distance.
In this work, we aim to derive the deterministic approximations for (t) and qo (t). To facilitate the analysis, make the following assumptions.
3

---

Published as a conference paper at ICLR 2026
Assumption 1. As - co, the ratio p/n + 0.
This assumption is common in the study of high-dimensional random matrices (Bai et al. 2010; Cui & Zdeborová 2023 Ham et al. 2025). Here, p = p(n) can be viewed as sequence indexed by n. In the following, we use p, co to denote this asymptotic regime.
Assumption 2. The random matrix Z has i..d. entries such that E[y, = 0, BllvZzj = 1, and BllvnZigl cc.
More rigorously, the elements Zi n 1/2 are sampled from a double array { (i,)
where X,y8 standardized i..d. random variables. This assumption is general and distribution- independent. The fourth-order moment condition is to ensure that the largest singular value of Z is almost surely bounded (Yin et al. 1988 Bai & Silverstein 1998).
MAIN RESULTS
3.1 LEARNING DYNAMICS
Our goal is obtain closed-form solution for the inner products 9., and g. in high dimensions. We note that qu and gu are random processes, and will demonstrate that these two random processes will almost surely converge two deterministic processes, which have closed-form expressions. These deterministic processes can be described by the integration of certain "basis" functions with respect the famous Marčenko-Pastte (MP) measure (Marčenko & Pastur 1967)
1+ (1 20) [(1 vo)2 x]]
m(dx) (1- - c +
dr.
(5)
2TTC
To simplify the notation, define the co-MP measure u(dx) -11-c) cc((dx) and the basis functions
((,) cosh(VEl), (2) = sinhi(2
2//
(a,=(1) Vx sinh(Vall, (,,((1) = cosh(2
(6)
Theorem 1. (Deterministic Approximations for qu and qu) Let un E Rp and vo E R2 be the initial vectors with unit norms and define Oa. q. (mo, ut) and ax gu (0) (vo, v*). Under
Assumptions and 2 we have. any T
sup OCTET
1qu (u((t)
a.s.
0,
(t) (t) T
sup
a.s. 0, pmn
p.n
where qu(t) Tu (t)/ / p() and (t) (ott)/v/())
Here, the deterministic functions (t),,,
(t), and p(t) are defined as
(t) = X
x (1s,x (s,0,) (t))(dx) - a.Als, (t)
JR
O c1 + 22)
[ [1e1,0 (lsz) + auli,0 (t),
(8)
X2
ane [
(1)
a (X2 + +c) /(lnee (oul
(9)
X²
p(t) = 1 + f La)da
f. da u(dx) (2au2(0o (a) + 2x² ( (a))
(10)
4

---

Published as a conference paper at ICLR 2026
/ x12 (1) al(rr)
(t,
+ 2
x { Mauulle (1x)) * (tt + 12 ( * f1,=)] (2)x(()) m(dx)
R
x (Aa, 61 (() + l= ( (I l2 [t)) p(dz)
+
(AaullT (3.x) * l (1) + 12 (7. * [3,0)] (1)) p(dx)
R
+
I (>arl(au l3,1) l(t) + >2 (á. la,ي) lt,) u(dx),
(11)
RR
with 0) = (1+ (2 x²)/2?.
Proof: The proof of Theorem is given in Appendix D
The expression for the evolution dynamics qu (t) and (t) involves a large number of integrals and convolution operations. However, since the MP measures μ and p have bounded support and density functions, the terms (t). (t). and p(t) can be computed numerically for finite t. Specifically, by equally partitioning the intervals [0, t] and [(1 /2)2, (1 + / into N₂ and N₂ grid cells respectively and applying the standard rectangle method approximate the integration, the overall computational complexity is at most O D(N,N=), Crucially, this complexity is polynomial and de- pends only the discretization parameters N, and N,, which not prohibitively large for the moderate values t.
An interesting observation is that when both a, and a, are zero, the inner products will be approx- imately zero throughout the learning process. This phenomenon can be intuitively explained by the gradient flow When us and u², as well vs and v", are orthogonal, their gradients are given by
du = -I,- u.u, Zu.,
du: = (I, 1a Zz u
(12)
dt
dt
This indicates that vectors ut and 2, will only along the directions spanned by projection of the noise and fail to learn ground truth. In Section 3.2, we will provide a more general condition under which learning fails in the large-time limit.
Theorem D has several potential applications, as discussed in following remarks. Remark 1. (Early Stopping Strategy Design) Consider the task of estimating the hime--dimensional signal directions ui and v". Assume we know the SNR X and the observation matrix X. Addition- ally, there is a performance requirement
min {Lue, ut"), 00),00)) > 7
(13)
for a given threshold y E (0, Let 2 = p/n be the estimator for c. According to Theorem ! the deterministic approximations only depend on (c, X, an, av). Thus, the following stopping time estimator can constructed
(.), = inf (1 0: min (l.())..)))) 2 77).
(14)
We note that f(,,, can be computed by line search method after obtaining que and qu. As result, the computational complexity remains dominated by O(NEN,). Furthermore, if we know the prior information about (u' v²)) (Aubin et al. 2021), can design the statistically optimal stopping time. For example, set tt Eon (tay [w)]] We note that in practice, knowing the exact SNR may be too restrictive. Fortunately, this method can also be extended to situations where E I for some known interval I. From a robustness perspective, we can set 1 = supxer T to determine the stopping time.
Remark 2. (SNR Estimation) Consider the task of estimating the SNR X. Assume the signal direc- tions (u' v²) are known. Applying gradient flow, we obtain the following estimator
A(u", ,("",v*) = arg min (ur, u²) Tu((t)llc(oorT) + +1 v²) (1)(00)()))))))
(15)
A>0
for given T > 0.
5

---

Kurchhn Published as a conference paper at ICLR 2026
We also have the following technical remarks.
Remark 3. The expressions for qu and To can be further simplified. For instance, in ( the term (fg.= ,3 (t) in the integration can be computed directly since the basis functions linear combinations of exponential finctions. Similarly, the term (Tu * (t) on the right-hand side (RHS) of (T) can also be simplified by interchanging the order of the MP integral and convolution according to Fubini's 's theorem. These simplifications can significantly reduce the complexity of the numerical evaluation.
Remark 4. (On the Role of Ju,w) The function Fa.. can be used to characterize the correlation between the estimates (ut, Ve) and the random noise Z. In particular, it can be verified by the proof in Appendix D that
sup (t, (ass, Zoril UEIET
u.s. O,
(16)
p.n-+00
where Iu, (L) (1) -v. (L)
Remark 5. In the field of statistical physics, the
(CHSCK) equation is utilized to describe the behavior of Langevin dynamics in the spherical p- spin glass model (Sarao Mannelli et al.) 2019; (020). which accounts for more complex mixed matrix-tensor structure and gradient noise. The loss function considered in the concerned rank-one matrix recovery problem can likewise be interpreted as the Hamiltonian of a corresponding spin glass system. However, the CHSCK equation is implicit and its solution is purely numerical, as it involves full-time correlations and requires solving functional fixed-poin equations for the Lagrange multipliers. In contrast, Theorem provides a closed-form expression for the dynamics. offering more tractable framework for understanding the learning behavior. Furthermore, Theorem [ has the potential to be extended to the case of stochastic gradient descent within the framework of stochastic calculus on manifolds (Du Prato & Zabczyk 2014: Gess et al. 2024), yielding more general results, which are left future investigation.
Remark 6. (Comparison with (Bodin & Macris (02T)) The proof strategy of Theorem 7 is similar to that in Bodin & Macris 2021). Specifically, we first construct system of differential equations for the characteristic functions using the resolvent, and then solve it via the Laplace transform. However, the anulytical properties and the existence and uniqueness of the solutions to the charac- teristic equation were not investigated in (Bodin & Macris [02I]. We note that these properties are important (Hachem et al. (207]). Specifically, the governing system (45) describes the dynamics of the inner products, so the uniqueness of the solution to (4) does not guarantee that of the solution to the constructed system. We show that this solution admits an integral representation via transition kernel and is unique under the given representation (cf. Theorem 3). Substituting (u,, v2) from (4) into the constructed system automatically satisfies the integral representation, which provides a rig- orous theoretical guarantee. This framework can also be applied to study the learning dynamics of other problems (Bordelon et al., 2020 Loureiro et al. 2021 Paquette et al. 2024 2021 Further- more, in contrast to the Wigner model, the rectangular model involves correlation terms between the random matrix and its resolvent. The almost sure convergence of the bilinear form is established in 36) which is not yet available in the literature. Additionall, the poles of the Laplace transforms are of higher order and are implicitly defined for several terms, which makes the problem more challenging.
Theorem iis derived for the continuous gradient flow. However, in practical applications, GD with finite learning rate is widely employed. A natural question is whether this framework can be extended to the discrete-time setting of GD. The following remark shows that the gap between the discrete and continuous processes is small. The detailed proof is given in Appendix
Remark 7. (Continuous to Discrete Process) Consider the initial point (ūo, vo) such that üolla =
||00/|
1. Assume the following canonical Riemannian GD update
úx+1
= relr(aa. 5x) (-1) gradHH(us ox)) = Trerru (-7 -pojja V VuHi(us Tx))
(17)
0k+1
retrs, (-1) proja, V V.H(ūn mx))
where the retraction mapping (utr( is given by (u))(u.))( yu) (retr, u (I vectr((,))
( |utt 1. [u+gulla)) and 1j represents the learning rate. Then, for any given T > and rj > 0,
utt. w+y.
there exists a constant C > 0 such that the following bound
max {l|ăx ukyll 1/yy- [knll2) C(n+ [
(18)
6

---

Published as a conference paper at ICLR 2026
holds with probability for large p and n. Here, (e)) denotes the solution of (4) under the same initialization, i.e., (uo, vo) = (uo, vo)- The above result demonstrates that in high-dimensten dimensional settings, the continuous gradient flow serves good approximation to the discrete GD over a finite time interval [0, . Therefore, we can characterize the behavior of GD by studying the dynamics of the system
As discussed in Section 2 the loss function H v) reaches the global minimum at (ū), U1 by the EYM theorem. We have the following corollary regarding the learning dynamics of the gap between the loss function and its optimal value.
Corollary 1. (Deterministic Approximation for the Loss Function) With the same settings as Theo- rem 7
we denote
H(u, v) = 2
x ww TooT
x
(19)
2
Further define X max(c!/4, A) and Ux.c (1+12)(c
12)
Then, we have, for any T
0,
sup 0. T
H(ue, v²) - -F((1)
a.c. م 0,
(20)
p,n-+00
where P((t) Vxx - (t) (t).
Proof: By algebra, we have
H(u, vz) Tr(X u.a )x v.u.)) 1 Tr(X ū, 0 (X 0,
-ū1, X (mt, X Xv,) 01 -(u,, Zor) Aqu (1).))))
(21)
where σi denotes the largest singular value of X. According to (Liu et al.
2025b Theorem 2).
the largest eigenvalue of X X converges to 2x,c almost surely, as p. n2 - co.
This implies Ø1
/ V0x. almost surely, by the continuous mapping theorem Van der Vaart 2000 Using Theorem
Theorem 2.3).
I and (16). 20 is proved.
We note that the loss functions H. in ( and H. in 19 are equivalent since the second term the RHS of (19) can be viewed as a constant. From Theorem [ and Corollary [ it can be observed that the learning dynamics is rotational invariant. In particular, for any Ohe O2, E [-1, all initial points (uo vo) satisfying (un, u²") = T. and (Do, v*) = Tr exhibit asymptotically same learning dynamics, i.c., the same evolution of the inner products and the loss function.
3.2 LARGE-TIME LIMIT
The following theorem demonstrates the asymptotic behavior of the deterministic approximations Tu and To. as 30.
Theorem 2. (Asymptotics of the Learning Dynamics) Define I Sgn t + 0o, we have
av
As
VI+X"
vx+ch
- 1(1 > c/4
lim (t) = Tu =I.
(22)
100
+ c
-x-1
lim (t) = Toe 1000
1+ 1(x> c/4}.
(23)
Proof:
The proof of Theorem 2 is given in Appendix
observe that when À < c1/4 the gradient flow estimation is asymptotically trivial Tn, qo + 0. When À exceeds c²/4, phase transition occurs, indicating that c2/4 is the critical threshold for estimation. This phenomenon is consistent with the well-known BBP phase transition. In the following remarks, compare Theorem 2 with the BBP phenomenon and discuss the potential applications. We note that although Theorem 2 is established for a continuous process, the BBP transition phenomenon also be extended to GD with a small learning rate by (18).
From Theorem 2

---

Published as a conference paper at ICLR 2026
Remark 8. (Comparison with BBP Phase Transition) Recall ū1 and 01 the left and right singu- lar vectors of X corresponding to the largest singular value. According to the BBP phase transition Baik al. 2003). we have
c 1
(ū u 1
0.8.
1fx c1/4),
(24)
P.n++
c
1 cx
1001 v²y
1fx> cl/43.
(25)
1+
The major difference between [2
25 and 22
23 lies in the sign and the case
V/- Vx=+c
This indicates that if we have prior information about the initial points on and ao, it is possible to determine the confidence intervals for the signs, thereby estimating the directions of the ground
truth with gradient flow. In the experiments, it will be shown that when Oco estimation becomes difficult.
Cxy 0, the
VI+X
Vx-40
Remark 9. (Estimation for X & c2/4) Even when Tu. Tu - 0, there is still a chance to estimate the signals. Assuming have prior information about Oa and aee, similar to Remark 7
we can set
Eaa. av) = arg max (T + (t)(+)))
(26)
130
and design the learning time according to i(a av)) We note that when the prior distribution of the ground truth "u","") known, the critical threshold for iterative algorithms like AMP may be lower than c4/4 (Lelarge Miolane 2017.
We also have the following phase transition phenomenon regarding the asymptotic behavior of the loss function.
Corollary 2. Define J 1f
ax 0}. As t + oG, we have
vx++1 vx+c
lim H(t)= (Vas -V- vdl t1000
(27)
Proof:
The proof of Corollary similar that of Theorem and thus omitted.
Recalling da,A = (1 + ( / it can be observed that when the SNR is smaller than the critical value, de, = (1+ (2) and the loss tends to as the learning time increases. With random initialization, often have PPJ
1, indicating that GD reach the global optimum and avoid saddle points in spherical low-rank signal estimation. In the case of tensors, GD could be more efficient, as the number of saddles on the tensor manifold grows exponentially with dimensionality (Subag & Zeitouni 2017 Auffinger et al. 2013). Furthermore, Corollary 2 implies that if J f
and À > c/4 we have
1 À
X +
X ++
lim Tu. (t)
(28)
x
x2
Note that (t) approximates the correlation between the estimates (ut v.) and the noise, i.e., Ju.e = (u,, Zo,). As A - oo, the RHS of the above equation tends to indicating that higher SNR leads to reduced dependence between (u,, vz) and the noise.
EXPERIMENTS
In this section, we verify the accuracy of the theoretical results. The simulation results are generated using GD, where the update rule follows I7 with a learning rate of r7 = 0.005.
Accuracy of the Deterministic Approximations: Figure I shows the loss function and inner prod- uct dynamics compared with their deterministic approximations. The parameters are set as = 900 and n 1200. The ground truth vectors are follows: the components of u²" and v" are in- dependently drawn from a standard Gaussian and Bernoulli distribution with p = 0.3, respectively. Both vectors are then normalized. The initial points are set as uo (10u* + 21)/ ||0u" + z l and vo = (5x" + /2)/ 50* + 22ll2- 2. where 2 E Re and 22 E R" are random vectors with i.i.d. standardized Gaussian elements. Figures Ia Ic the SNR is set A 0.3, which below the critical threshold c1/4, In contrast, in Figures 1d If the SNR is set À 1.5 1/1 By comparing
8

---

strictly positive for all Published as a conference paper at ICLR 2026
64
20% 6.11
Ses 4:11
So Ry
6.2
Tins K!!!
0.12
4.2
3 0.10
0.2
040
(a) 0.3
(b) 0.3
(c) = 0.3
0.0
Smm 10 fm
077
04-
02-
6:0
04
6.71
02-
04
05
(d) 1.5
(e) 1.5
(1) A = 1.5
Figure 1: Theoretical and Empirical Learning Dynamics.
0.6
0.5
V0 65-0
Visc R20
1.3-
6.5
10P
0.4
02
0.2
Figure 2: Learning Dynamics with Different SNR Values X.
the dark-colored lines and the light-colored lines, can be observed that the deterministic approx- imations qu. To. and H are accurate, which validates the accuracy Theorem | and Corollary Additionally, it can be seen that the deterministic approximations become more accurate as the SNR increases. This because, when the SNR is lower, random noise Z dominates, thus qu: qu have more randomness. By comparing the asymptotic values (dash-dotted lines) with the deterministic approximations, the accuracy Theorem
2 is verified.
Impact of the SNR: Figure 2 illustrates the theoretical learning curves with different SNR values. The parameters are set as c 0.5, a, = 0.211, and a, = 10.121. can be observed that when - oc. However, when
& c//4 0.841, Tu and fu approach to Th and T20 respectively. which are positive for the concerned case. This validates the accuracy of Theorem terms of the sign of the limits, because 0, + Vo+xT
c//4, they converge
VI+x²
X E [0.1,2]. 2].
can be observed that as the SNR increases, the GD-based algorithm learns the hidden information faster. From the third sub-figure in Figure 2 we note that although the loss may be larger during the early learning stage with higher SNR, it eventually decreases to lower level during the training. This indicates that learning a global process and rapid initial decrease in loss does not necessarily imply better final performance. Conversely, for data with high potential (correspond- ingly, high SNR), the loss may not decrease significantly in the initial stages of learning.
Impact of the Initial Points: Figure S illustrates the learning curves of qa with different initial values Oh.. The parameters are set
0.5 and a, 0.4. In Figure [al the SNR is set to A = 1.5 d/4 while in Figure 3b x 0.2 < d/4 The solid red curve corresponds to the case
where VT+X²* + Vatx² 0, while the dashed red curve represents the line = V
Vc.
The trend of the loss functions validates the accuracy of Corollary 2

---

Published as a conference paper at ICLR 2026
2.5
-0.14
0.5
0.5
Value of o.20)
= -0.14
1.5
0.0
Value
1.0
-1.0
0.0
-0.9
(a) À = 1.5
1.0
0.9
1.5
Value of 4,30)
Value of Ri
0.5
0.5
-0.5
-1.0
-0.9
10
(b) 0.2
Figure 3: Learning Dynamics with Different Initial Values Ou.
It can be observed that the initial value affects the speed of convergence for À > c//4, When at is close to -tt
x++c the gradient flow converges more slowly. whereas when the initial value ag is X=+1
near T the gradient flow converges faster. This also provides insights on the learning algorithms. Given the complex geometric structure of the loss surface, if the initial point is unfortunately located near some "bad" regions, it may take significant efforts for the model to learn the hidden information from the data, even though the loss decreases with GD.
CONCLUSIONS
In this work, studied the GD dynamics of the rank-one matrix denoising problem. In particular, we obtained closed-form deterministic approximations for inner products between the learned vectors and the ground truth, and proved that the random learning curves converge to the approxi- mations almost surely. Additionally, we derived the asymptotic behavior of the learning dynamics, which is consistent with the BBP phase transition phenomenon. Simulations validated the accuracy of the theoretical results. Furthermore, it was observed that the gradient flow converges faster when the SNR higher and when the initial point is close the stationary points. In contrast, there exist "bad points" where learning is disrupted. If the initial value is near the bad points, it takes longer time learn the hidden information even with high SNR.
The main results of this work, i.e., Theorems D anddd] have many potential applications and can be further extended. Given the prior information about the initial points and the ground truth, these results can be utilized to design early stopping strategies and develop SNR estimation algorithms. In this work, we have established the almost sure convergence of the gradient flow dynamies, which corresponds to the "law of large numbers." The "central limit theorem," i.e., the asymptotic fluc- tuation of qn (t) around ğ (t). h E {u, v). remains an interesting open problem. Such fluctuations capture higher-order information about the test errors of learning algorithms. Moreover, the data model considered in this work is rank-one and linear. Extending this analysis to multi-rank mod- els and nonlinear data structures, such random feature models, represents challenging direction. From a practical standpoint, investigating these aspects would yield deeper insights into the dynamics of learning.
general and
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENT
This work was supported in part by a grant from the NSFC/RGC Joint Research Scheme spon- sored by the Research Grants Council of the Hong Kong Special Administrative Region, China and National Natural Science Foundation of China (Project No. N HKUSS6656/22).
REFERENCES
Benjamin Aubin, Bruno Loureiro, Antoine Maillard, Florent Krzakala, and Lenka Zdeborová. The spiked matrix model with generative priors. IEEE Transactions on Information Theory, 67(2):
1156-1181, 2021.
Antonio Auffinger, Gérard Ben Arous, and Jifí Černy. Random matrices and complexity of spin glasses. Communications on Pure and Applied Mathemutics, 66(2):165-20) 2013.
Zhidong Bai and Jack W Silverstein. No eigenvalues outside the support of the limiting spectral distribution of dinenssinnni sample covariance matrices. The Annals of Probability, 26(1): 316-345, Jan. 1998.
Zhidong Bai and Jack W Silverstein. Exact separation of eigenvalues of large dimensional sample covariance matrices. The Annals of Probubility, 27(3):1536 1555, 1999.
Zhidong Bai and Yongquan Yin. Limit of the smallest eigenvalue of a large dimensional sample covariance matrix. The Annals Probability, 2163):1255-123 1993.
Zhidong Bai, Jack William Silverstein, et al. Spectral analysis of large dimensional random matri- ces. Springer, 2010.
Jinho Baik, Gérard Ben Arous, and Sandrine Péché. Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices. The Annals of Probability, 33(5):1643 - 1697, 2005.
Zhigang Bao, Xiucai Ding, Jingming Wang, and Ke Wang. Statistical inference for principal com- ponents of spiked covariance matrices. The Annals of Statistics, 50(2): 44441111, 2022.
Jean Barbier, Mohamad Dia. Nicolas Macris, Florent Krzakala, Thibault Lesieur, and Lenka Zde- borová. Mutual information for symmetric rank-one matrix estimation: A proof of the replica formula. In Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016.
Florent Senaych-Georger and Raj Rao Nadakuditi. The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices. Advances in Mathematics, 227(1):494-511, 2011.
Florent Senaych-Georgee and Raj Rao Nadakuditi. The singular values and vectors of low rank perturbations of large rectangular random matrices. Journal of Multivariate Analysis, 111:120- 135, 2012.
Christopher M Bishop and Nasser M Nasrabadi. Pattern recognition and machine learning, vol- ume 4. Springer, 2006.
Åke Björck et al. Numerical methods in matrix computations, volume 59. Springer, 2015.
Antoine Bodin and Nicolas Macris. Rank-one matrix estimation: analytic time evolution of gradient descent dynamics. In Conference on Learning Theory, pp. 635-678. PMIR, 2021.
Blake Bordelon, Abdulkadir Canatar, and Cengiz. Pehlevan. Spectrum dependent learning curves in kernel regression and wide neural networks. In International Conference on Machine Learning, pp. 1024-1034, 2020.
Lucilio Cordero-Grande, Daan Christiaens. Jana Hutter, Anthony N. Price, and Jo V. Hajnal. Com- plex biffuhtionweighteed image estimation via matrix recovery under general noise models. Neu- rolmage, 200:391-4000 2019. ISSN 1053-8119.

---

Published as a conference paper at ICLR 2026
Romain Couillet and Walid Hachem. Fluctuations of spiked random matrix models and failure diagnosis in sensor networks. IEEE Transactions on Information Theory, 59(1):00- 2013.
Romain Couillet and Zhenyu Liao. Random matrix methods for machine learning. Cambridge,
U.K.: Cambridge Univ. Press, 2022.
Romain Couillet, Florent Chatelain, and Nicolas Le Bihan. Two-way kernel matrix puncturing: towards resource-efficient pea and spectral clustering. In International Conference on Machine Learning. volume 139, pp. 2156-2166 2021.
Hugo Cui and Lenka 7deborová. High-dimensional asymptotics of denoising autoencoders. In Advances in Neural Information Processing Systems, pp. 11850-11890, 2023.
Giuseppe Da Prato and Jerzy Zabezyk. Stochastic equations in infinite dimensions, volume 152.
Cambridge university press, 2014.
Xiucaí Ding and Fan Yang. Spiked separable covariance matrices and principal components. The Annals of Statistics, 49(2):1113 1138, 2021.
Carl Eckart and Gale Young. The approximation of one matrix by another of lower rank. Psychome- trika, 1(3):211-218, 1936.
László Erdós and Horng-Tze Yau. A Dynamical Approach to Random Matrix Theory, volume 28. American Mathematical Soc., Providence, RI, USA, 2017.
Matan Gavish, William Leeb, and Elad Romanov. Matrix denoising with partial noise statisties: optimal singular value shrinkage of spiked f-matrices. Information and Inference: A Journal of the IMA, 222002220- 2023.
Benjamin Gess, Sebastian Kassing, and Nimit Rana. Stochastic modified flows for riemannian stochastic gradient descent. SIAM Journal Control and Optimization, 62(6):3288-33 -314,,344444 2024.
Walid Hachem, Philippe Loubaton, and Jamal Najim. Deterministic equivalents for certain func- tionals of large random matrices. The Annals of Applied Probability, (3):875-930) 2007.
Jonghyun Ham, Maximilian Fleissner, and Debarghya Ghoshdastidar. Impact of bottleneck lay- ers and skip connections on the generalization of linear denoising autoencoders. arXiv preprint arXiv:2505.24666; 2025.
Philip Hartman. Ordinary differential equations. SIAM, 2002.
lain M. Johnstone and Debashis Paul. PCA in high dimensions: An orientation. Proceedings of the IEEE, 177(6):1277-1222, Aug. 2018.
Satish Babu Korada and Nicolas Macris. Exact solution the gauge symmetric p-spin glass model complete graph. Journal of Statistical Physics, 136(2):205-230, 2009.
Jean-Françon Le Gall. Brownian motion, martingales, and stochastic calculus. Springer, 2016.
Jeffrey T Leek. Asymptotic conditional singular value decomposition for high-dimensional genomic
data. Biometrics, 674))3444-35,, 2011.
Marc Lelarge and Léo Miolane. Fundamental limits of symmetric low-rank matrix estimation. In Conference on Learning Theory, volume 65, pp. 1297-1301, 07-10 Jul 2017.
Thibault Lesieur, Florent Krzakala, and Lenka Zdeborová. Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications. Journal of Statistical Mechan- ics: Theory and Experiment, 2017(7):073403, 2017
Qianxiao Li, Cheng Tai, and Weinan F. Stochastic modified equations and adaptive stochastic gradient algorithms. In International Conference Machine Learning, volume 70, pp. 2101- 2110. PMLR, 2017.
Zhenyu Liao and Romain Couillet. On the spectrum of random features maps of high dimensional data. In International Conference Machine Learning, pp. 3063-3071, 2018.
12

---

Published as a conference paper at ICLR 2026
Qiang Liu. Stein variational gradient descent gradient flow. In Advances in Neural Information Processing Systems, 2017.
Xiaoyu Liu, Yiming Liu, Guangming Pan, Lingyue Zhang, and Zhixiang Zhang. Asymptotic lim- its of spiked eigenvalues and eigenvectors of signal-plus- noise matrices with weak signals and heteroskedastic noise. Bernoulli, 31(3):2351 2376, 2025a.
Xiaoyu Liu, Yiming Liu, Guangming Pan, Lingyue Zhang, and Zhixiang Zhang. Asymptotic lim- its of spiked eigenvalues and eigenvectors of signal-plus-nosse matrices with weak signals and heteroskedastic noise. Bernoulli, 31 31(3):2351-2377 2025b.
Bruno Loureiro, Cedric Gerbelot, Hugo Cui, Sebastian Goldt, Florent Krzakala, Mare Mezard, and Lenka Zdeborová. Learning curves of generic features maps for realistic datasets with teacher- student model. In Advances in Neural Information Processing Systems, volume 34, pp. 18137- 18151, 2021.
Per-Gunnar Martinsson and Joel A Tropp. Randomized numerical linear algebra: Foundations and algorithms. Acta Numerica, 29:403-572, 2020.
VA Marčenko and Leonid A Pastur. Distribution of eigenvalues for some sets of random matrices. Mathematics of the USSR-Sbornik, 72(114):4, 1967.
Raj Rao Nadakuditi. Optshrink: An algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage. IEEE Transactions on Information Theory, 60(5): 3002-3018, 2014.
Courtney Paquette, Kiwon Lee, Fabian Pedregosa, and Elliot Paquette. Sgd in the large: Average- case analysis, asymptotics, and stepsize eriticality. In Conference on Learning Theory, pp. 3548- 3626, 2021.
Elliot Paquette, Courtney Paquette, Lechao Xiao, and Jeffrey Pennington. 4+3 phases of compute- optimal neural scaling laws. In Advances in Neural Information Processing Systems, pp. 16459- 16537, 2024.
Leonid Andreevich Pastur and Mariya Shcherbina. Eigenvalue Distribution of Large Random Ma- Irices. American Mathematical Society, Providence, RI, USA, 2011.
Henrik Pedersen, Sebastian Kozerke, Steffen Ringgaard, Kay Nehrke, and Won Yong Kim. k-t pea: temporally constrained k-t blast reconstruction using principal component analysis. Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, 62(3):706-716, 2009.
Arnu Pretorius, Steve Kroon, and Herman Kamper. Learning dynamics of linear denoising autoen- coders. In International Conference on Machine Learning, pp. 4141-4150, 2018.
Stefano Sarao Mannelli, Florent Krzakala, Pierfrancesco Urbani, and Lenka Zdeborova. Passed & amp; spurious: Descent algorithms and local minima in spiked matrix-tense models. In Interna- tional Conference on Machine Learning, volume 97, pp. 4333-4342, 09-15 Jun 2019.
Stefano Sarao Mannelli, Giulio Biroli, Chiara Cammarota, Florent Krzakala, Pierfrancesco Urbani, and Lenka Zdeborová. Marvels and pitfalls of the langevin algorithm in noisy nigh-dimensional inference. Phys. Rev. X, 10:01 1057, Mar 2020.
Eliran Subag and Ofer Zeitouni. The extremal process of critical points of the pure p-spin spherical spin glass model. Probability theory and related fields, 168(3):773-B)))) 820, 2017.
Terence Tao. Topics random matrix theory, volume 132. American Mathematical Society, 2023.
Aad W Van der Vaart. Asymptotic statistics, volume Cambridge university press, 2000.
Roman Vershynin. High-dimensional probability: An introduction with applications in data science, volume 47. Cambridge, U.K.: Cambridge Uniy. Press, 2018.
Yuchen Wu and Kangjie Zhou. Sharp analysis of power iteration for tensor pca. Journal of Machine Learning Research, 25(195):1-42. 2024.
13

---

Published as a conference paper at ICLR 2026
Fan Yang, Sifan Liu, Edgar Dobriban, and David P. Woodruff. How to reduce dimension with PCA and random projections? IEEE Transactions on Information Theory, 67(12))8155-8188 Dec. 2021.
Jianfeng Yao, Shurong Zheng, and Zhidong Bai. Sample covariance matrices and high-dimensional data analysis. Cambridge UP, New York, 2015.
Yongquan Yin, Zhidong Bai, and Pathak R Krishnaiah. On the limit of the largest eigenvalue of the large dimensional sample covariance matrix. Probability Theory and Related Fields, 78:509-521, 1988.
Yihan Zhang and Marco Mondelli. Matrix denoising with doubly heteroscedastic noise: Fundamen- tal limits and optimal spectral methods. Advances in Neural Information Processing Systems, 37: 2024.
Zeyan Zhuang, Xin Zhang, Dongfang Xu, and Shenghui Song. No eigenvalues outside the limit- ing support of generally correlated and noncentral sample covariance matrices. arXiv preprint arXiv:2507.03356 2025.
14

---

Published as a conference paper at ICLR 2026
APPENDIX
A THE USE OF LARGE LANGUAGE MODELS (LLMS)
In this work, utilized LLMs for grammar checking and text polishing to enhance the readability.
B MATHEMATICAL TOOLS
In this section, introduce the mathematical tools used in our analysis along with related discus- sions.
B.I RESULTS ON CALCULUS OF FUNCTIONS
Lemma 1. (Grönwall's 's Lemma) Let T on [0, I] such that, for every E [0, T),
and h a nonnegative bounded measurable function
h(t) < < f h(t)dt,
(29)
for constants > 0 and b 0. Then, we have for every t [0,
h(()) & exxp(b(].
(30)
Lemma 2. (Final Value Theorem) Let f E C(TR.) be such that lim f()) = exists for some a > O and A E R. Moreover, assume that g E C(R+) satisfies Jo e- (g(t) dt < co. Then, have
lim e (J+ () = A. l00 e g(()dl. t-000
(31)
Lemma 3. Let f C((a,b) b)) be such that f(b) + 0 and 1. Then,
Stt f(x) expl Tit wit
- 2)
lim
22 a-1 1,
(32)
f(6) expl Vbtji
where I(2) = 1000 2--1 expl- ttlt denotes the Gamma function.
Proof:
The proof of Lemma is given in Appendix G
B.2 RESULTS ON RANDOM MATRICES
The deterministic approximation of the gradient flow relies the convergence of the resolvent for the covariance matrix Z, which defined as
Q(2) (ZZ 2Ip)
(33)
where z E such that Q(z) is well-defined. The co-resolvent is defined as Q(x) (Z
=I,) The convergence of the resolvent is a fundamental topic in RMT Bai et al. 2010 as characterizes the spectral distribution of the random matrix Z. Recall c limp.n- p/n. n. Define
E_ (1 V E (1 Vo)? , and the set M {0) U [E_, E+]. The following lemma
describes the asymptotic behavior of the resolvents.
Lemma 4. (Convergence of the Resolvents) Assume Assumptions 7 andd Phold and let u E Rp and
E R be deterministic vectors with ||u||| ||l||2 = 1. Let I be compact set in C such that ( > 0. Then, we have,
sup Q(z)( m(=)
a.s.
0,
(34)
=er
p.n-+00
sup ser
(TQ(()e m())
O,
(35)
p.
15

---

Published as a conference paper at ICLR 2026
sup =eT
u T Zu
0.6. 0.
(36)
p,0-+00
where m(2) is the unique solution
zem² - -1 (1 + +1
(37)
such that m(2) E C+ {2: 3()) 0} for C+, and m(2) (2) (c 1)/2.
Proof: When I is singleton, the proofs for (34) and (35) be found in (Bai al. 2010 Couillet
Liao 2022). The extension to a general compact set P follows from an e-net argument and a
similar discussion in (Bai Silverstein 1998 Eq. 3.23). The proof for (36) is given in Appendix
H
In fact, it can be shown that m(=) the Stieltjes transform BBai et 2010) of the MP distribution. Specifically, we have the integral representations
m(z) / ( m(dx))
m(2) f u(dr))
(38)
where
u(dx) - c1) 8(0) + V(x - E_) ( x)+
dz,
2TE
u(dr)= (1 - 8)) + Ve E_)+(E+ xx)+
dx.
(39)
2T
C PROOF OF THE UNIT NORM CONSTRAINT IN GRADIENT FLOW
By the evolution in (4) and taking the derivative, we can obtain
Us 12 2 ( uc dt
du. -2 (14 proja,. (V.7l(m,. mz)))
dt
2 uu, (Ip - unu. Vu?((uu vz)) 2(1 1|2 1)u, V Hull, v.).
(40)
For ease of notations, we define (t) |ur12 -1 and (t(() uu, VuH(uu v.). Thus, have ft(() fo fu (a)))a())do which implies
fu(t) C exp 2 22// gu (a((a)d)))
(41)
for some constant C. Since (o) = 0, we have C' = 0 and ||ue|| 2 1, for any / > 0. A similar argument on || yields the same result. Therefore, we have proved that the unit norm constraint is satisfied throughout the flow.
D PROOF OF THEOREM
In this section, we will prove Theorem 1 The proof consists of four parts.
1. We first construct system of differential equations for the characteristic functions (Stieltjes transforms) associated with 4 which are more tractable.
2. We establish the existence and uniqueness of the solution with specific integral represen- tation.
3. Using Lemma we develop a deterministic approximation for the original system and obtain closed-form solution.
4. Leveraging the integral representation and the uniqueness property from Step 2, show that the solution to the original system converges almost surely to that of the approximate system.
16

---

Published as a conference paper at ICLR 2026
D.1 DERIVATION OF TRACTABLE SYSTEM OF DIFFERENTIAL EQUATIONS
We start by expanding the gradient flow in (4) to
U
X u Xv:
dt
1/
X u, 010; X {uu
ut
Zu + \u v" 21 Uelb Zo; bua T u² Vi
Z u, + Ao u uu - 1,22 - z
\ v" u uu.
Zue + go'go(t) - fu)o,(tt)u
(42)
2 u, + \v (1() - [1)m;
where fu, (t) + (q.(t)go and (t.v(t) (ue, Zor). Directly solving the above system challenging due its highly coupled nature. To this end. we define the inner products (he, (v) and the resolvent-related functions as follows
Cu(z) = (u', Q(()u²), Cx(2) = ( Q(())*), Du,o(2) = ((=))""""" Gu(t,z) (u, Q(z)u)), Go(t,z)= = (v, (()))())))) B,(t, = (u, Zv)
=) = Q(=) (x) Ta(t,2)=
= (u,, Q())u)) (.((,2) -m: Q())m))))))))))))))))))))))))().
(t,
(uu. (m: Q(=)Zz))
(43)
The quantities defined in (43) can be viewed as the characteristic functions (or Stieltjes transforms)
associated with the inner products. To illustrate this, consider G (1, 2) example. By selecting a positively oriented contour I that encloses all the eigenvalues of ZZ' and applying Cauchy's
integral formula, we obtain
2TJ f. Gu(t, = =.(t,)))q((1))
(44)
To simplify notation, we define the symbolic variables set Sode (Gu, Eu. T., Y e Using (4) we can derive the differential equations for each E Sode Formally, we introduce the following system of differential equations for the complex-valued functions
Guy G. E. Eu I1 T,, T. R4 CC
0 G. 2) (t, 2) + (q.(t)C.(z) (
(45a)
0
Go :) (t, 2) + (1)(() (unaat)GG((t.)),,
(45b)
0
s. 2) (1) + =G,(t,)) (1),,,, -u.
(45c)
0
E. (t, 2) = (t) + (G, (t) (t.,
(45d)
0t
11 (1, 2) + =1,r,(4, (1, N1. (0
02
2f,
(45e)
0 (1, 2) = TT, (), (1)G (oo
(45f)
2 0t
1 0
2 Ot T. (t, = II,, (t. 2) (t)
(45g)
subjected the initial conditions
=(0 Bs(z), VS E Sode
(46)
There exists a bounded interval I C. IR such that for any given f E R, and S E Sode 2 St,,2)
analytical on L. Additionally, the following constraints hold for all f R4 and any positively oriented contour I enclosing I
(t) = (t) + (t)
(t( 1 ff (t. (x(2) 1 $ ru(t,2)dz = 2 fr
2n f. Go(t,2)dz,
2n fr
1 $ (t,
f r.(t,2)d== 1. (47)
(1,(()) 2n Jr
2m jr
17

---

Published as a conference paper at ICLR 2026
D.2 EXISTENCE AND UNIQUENESS OF THE SOLUTION
A fundamental question is the structure of the solution to the system of differential equations (45
In the following theorem, show that if the initial conditions (=(=) for S Sades together with the functions (=), Cz(z), and Pu,r(z), admit certain integral representations (which include (43)
as special case), then the solution (45) exists and unique within a certain family.
Before delving into the details, we introduce the definition of a "good transition kernel". Specifically, a function J : R ((K) + called good transition kernel if and only if
For each E R4. p(t, -) is signed measure on 14+ with bounded total variation.
The support is uniformly bounded: there exist E' < oo such that
supsup(supp(u((, ))) < 120
(48)
For every Borel set A E B(R). the mapping (u(: A) is continuous in the sense that
lim [ll(t,,) -) plt + 6, )lzv = 0.
(49)
Define the family of functions Sic as
Sx - (. {
/ p(t,dx)
1b is good transition kernel
(50)
JR
Theorem 3. Assume that functions Cu, Cv, and Du, admis the following integral representations
Cu(z)= - / (cx(dx)) Cz(z) / (dx)
J (.D...
x- Du,u(z)=
(51)
where a.c, PC.vv and P.D.vv are non-negative finite measures supported over R... Moreover, their support is bounded, ie, supp(uc.). ) supp(pc.), and supp(u C [0. E. for some E. < co and E CV E). Further, assume the initial conditions Bs, for Sode also have the integral representations of the form
vs(dz)
Bs(=) = x x-
(52)
where Vss are finite signed measures supported [0, and satisfy Vr (R4) VT, (R.) 1. Then, the system of equations (45) admits unique solution (G., G. E,, Tu, T., T.)
E 8%
Proof: We first establish the existence of local solutions via Picard iteration by inductively structing the corresponding integral representations and proving the convergence of the underlying measures. Choose the initial functions 5 Sx kernels u's,, such that
for all S E Sode represented by good transition
r T-2 - AS E Sode:
1 (g((1, dx)
(1, =
(53)
with v4. (t, 1 for all t a (0,)) = us(-) for all E Sode Suppose the induction hypothesis holds for A. We then show that 1) Sk+1 Sx: 2) the support of the underline measures u f are bounded by E, and 3) 1 (t, PT. (t, R4) = 1 under Picard's 'ss iteration scheme. In particular, the iterate G2 is given by
k+
Gl = BG. (z) + [ (a, 2) + (a) 1 (a) 2)da.
(54)
By the induction hypothesis, we can obtain
Gl+ (1,2) { VG_, (dr)
+ X S'' de (()dda f(_ PC.. (cx( (dx) S.' { ph (a)) (a, dx)
(.¹) JaL x x-
W (a, dr)da T --2
JR.
I. -
(55)
z
18

---

Published as a conference paper at ICLR 2026
where (t) f (t, R), (t) = p/c. (t, (t) pi.... (t, and 2 = - +
(a) (a). By Fubini's theorem, we can construct the following iterative scheme for J.G..
R+1 PGG (t, A) = vG. (A)+ ( 1 (a, A) + (a)
(a, A)da,
(56)
(t, (ottttt, can be bounded as for any Borel set A E (TR). The total variation for jaG..
k+1
R+1 (t, 'lrv & IPGallTV L. [I (a, (ar + (a)alllllcllll 1
+, Imē (a. (a.-) da.
(57)
which finite. The continuity condition 49 follows directly from the definition of the iteration scheme. Since u -) and VGG.() (-) are supported over [0, E], the measure PC., k++ (t,, .) also supported over [0, Fl. Therefore, MG.. k+1 is good transition kemnel and G E Sk. By construction, the relation 5*+1 E Sx. for the remaining terms E Sade can be shown in a similar manner, and we directly give the iteration scheme for the transition kernels as follows
[ - (a. (at(a))o. (A) (a)pa, (a, A)da,
HG.. k+1 (t, A) vG.. (A)
(58a)
NG. (t, A) (A) + / (a, >g (A) (a)) (a, A)da,
(58b)
k+1 (1, A) = 1=. + d'da spt. (a, + (a) (dx) - 1 x (a, dx);
"E
(58c)
k+1 (t, A) = (A)+ + [ apc dr) dr) + X((())D, - ( (a, dx),
"E
(58d)
k+1 Pll,,o (t, A) = VI (A) + S. da zurr (a, dr) + xph dx) + S(da
+ +N (a, (a, 21th (a)ttt) dx),
(58e)
1t, A) VT + !'¹ da S pt.... (a, dx) (a) (a, dx) ( (a, dr),
(58f)
*+ (1, A) = UT. (A)+ +2 S da da/ / Mie. (a, + Xo S (a)) (a, dr).
(58g)
Obviously, for any u1 k+1 S E Sode, the support of the measure M9 k+1 (t., :) is contained in [0, E].
k+1 k+1 is analogous), Here, to show that ps is good transition kernel (the argument for p. and PII.. we use the fact & E and
|ap (4, -Ivv sup / splè. (l, < Essul A
A / N/. (1, (1, = E 1 (4.))
(59)
By the induction hypothesis u (,, = u (1, 1 = 1, we have
pr. A+1 (t, vY. VY. (R+) +2 f pi... (a, (c, + (a) (q. fh xa)p44 (a, )d
=1+2 L (a) (a((a)a - 1. (a) =
(60)
The relation 1 (t, R 1 for any follows analogously. This completes the inductive step, and have thus constructed sequences (5*) k2o in Sx-
19

---

Published as a conference paper at ICLR 2026
Next, will prove the existence of the local solutions. By the continuity of the kernels, choose a sufficiently small positive number T such that for any k 3 the following holds
sup max (t, 2sll (t, 11 & U, VS E Sodes
0STET
(t)) E OSTET
sup
(61)
where U > is a constant. In fact, the bounds 61 proved by induction. Assume it holds for some h > 0. Then for + 1 and any Borel set A R. we have
|et (T, vG. (A)l & S/' In A) + (a))lecc..
|1 (a) L (a, A) va. (A)) + | (a) lua, (A) da < T (U + llue. llv + U² U lrv)
(62)
Now, choose such that 7 & (1 + Iluc, + U ||G llrv)= ) (note that since (t) )) good transition kernels, such T exists). Then, for this choice of T, obtain u (T, A) (A)) U for any Borel set A. The bounds for the other measures follow similarly and the details are omitted.
We now prove the convergence of the measures on the interval [0, T}. To this end, define
M&(1)=
max SESade
-
(63)
Hence, for each and E [0, T), the following estimates hold
(a) - ( = -|lii... (a, / (a, R) < Ma-1((a),
(a) 1 1 (a) & (2AU + 1)
(64)
Consider the difference of successive iterates for JCG.
/// -) - A)
k+1 (1, A) p/c. (4, A)) MG.
+ (a) (A) - (a)
-Fa (a)m (a, - (a)) (a, A)) da
S. [1+|ue (A) U(AAU + + U] M Ma-1(a).
(65)
Taking the supreme over all Borel sets A, obtain for all <T
A+ (t, pá. (t, Urr & 1// Mx-1(a),
(66)
where the constant C' = + ||ec, Il + U2AU + 1) + U. By analyzing the remaining terms S Sode, can get similar inequalities with difffernnt constants C. Combining the estimates for all S E Sode and using the definition of M2 (t), we obtain
M&(t) & Ō S( Ma-1(a).
(67)
where the constant Ō is independent of k. Since Mo(a) & 2U, an induction argument shows that for every & >
sup Ms(t) 2U(T+)*
(68)
0E
20

---

Published as a conference paper at ICLR 2026
In particular, we have Exs1 SUPOSIET Mx(t) c0. Therefore, for each [0, T) the sequence
of measures {uš (t, }} converges total variation to limit measure μ (t, :), for every Sode- Passing to the limit k + in the iterative scheme (58) and using the integral representation, one can verify that the functions defined by
Stt, 2) = 1 (t,
tE 10 2 EVO,
(69)
satisfy the system 455) for all S Sode. The convergence of total variation implies that each uS (f, is finite and supported within [0, E]. Moreover, since the convergence is uniform 6, i.e. suPteT (Ia)((,)) p3 (1. TV + 0 k + oc, follows that n (L,-) continuous the sense that lims-+0 (t.-) - (++ ó, ,lyy = 0 for all [0, and E Sode- This implies ps. S E Sode are good transition kernels.
To prove uniqueness, we argue by contradiction. Suppose there exist two distinct solutions S and S2 belong Syc with underlying measures jAs.1 and u4s,2 for each S E Sode- Denote
M())
sup (I1ss,((,,)) 1s,,(1,, )lrv). SESude
(70)
By the inversion formula of the Stieltjes transform (Hachem al. 2007), for i 1,2 and any interval (a, b] C R. we have
ps,(t, b]) + us. (t, (a,b))
lim 3.0
(t,(tt+++8)) da.
(71)
2
By considering the difference between the two equations associated with S and S2, and using the inversion formula 71 and together with a standard яーÀ argument (Le Gall 2016), we conclude that (t)) satisfies
M(t) < C S( M(adda, Vt E [0,
(72)
for constant C O. Since I|Ms, (ll ; 1,2, Sode bounded, it follows from Grönwall's 's lemma (Lemma I that (t) = 0 for every t = TJ. which establishes local unique- ness. We have thus demonstrated the existence and uniqueness of the solution on the interval [0, J] Using a standard extension argument for solutions of differential equations (Hartman) 2002). the solution can be extended to all > O, and we omit the details. Therefore, we have completed the proof.
can be verified that the solution to the system of differential equations (42) exists uniquely. Sub- stituting this solution (uz, v.) into (43) yields S(1,2) for each S E Sode- Moreover, the functions Bs, Cuy G. and Du,. satisfy the conditions of Theorem 3 with probability one. Consequently, the
solution {S(t, 2). E } coincides with the solution { 2). E Sode } given in Theorem ES U. Then, we have
For example, let the SVD of Z be given by Z = UAV
U UTu
(dx)
=z)
- (())))))(((
Σ
(73)
where the distribution function is given by Fc, (x) Σn representations for the other terms S Sode can be verified similarly and details omitted for brevity.
UT 11 < x). The integral
We note that the agreement between these two solutions indicates that the solution of 45 within Sx accurately captures the learning dynamics of matrix denoising.
D.3 CLOSED-FORM SOLUTION TO THE APPROXIMATED SYSTEM
Although Theorem 3 characterizes the solution structure of 45 obtaining closed-form solution still infeasible. This is due to the fact that the functions C.,, C,,, and D as well as the initial conditions Bs, random and lack explicit expressions
21

---

Published as a conference paper at ICLR 2026
However, in the high-dimensional setting, Lemma 4 implies these quantities will converge to deter- ministic limits. Specifically, as p. n + oc, for any E C,M, we have
C=(z) pin-00
m(z), Co(z) pa
mn(z, Du,((2) p.n
O,
(74)
and
[Ba (z), Ba, (z), Br, (2), Br, (2), Bz, (z), Bz, (2), BT)..
0 (a amm 2), m(=) 0,,,,,
(75)
where au = (uo, u), ag = (v0. and ml and m(z) are defined in Lemma 4 Hcuristically, may substitute the random functions Bs. C.. Cv. and Du, with their determimistic limits and consider the following asymptotic system of differential equations
Gu((,2) ((,(,)) + x4.(0)((=)
Smna(()Gu((,)),
(76a)
02
(t, 2) = (t, =) + Aqu(t)m(=) - (t)
(76b)
0t
0
E. (t, 2) = (.()) + =G,(t,2) fu.o)at)
(76c)
00
0 3. (4, =) = (1) + =G.(l,,) 2) Su.a z),
(76d)
0
lt, 2) 2 =|T,(t, 2) (t,,, (lq, (()()))))))))))))((((
(t(t)).((t,z
Ot
2.u. (t)lo,))) 2),
(76e)
1 0
fu,u, 2t, = 11, 2) + Ago fuu (t)
(76f)
2 00
1 0
2 0t (l, (1, (1)(())))(, Jun, (1) (1) (1,
(76g)
with initial conditions given in (75) and subject to constraints (47) We note that, according to Lemma 4 and (75). the initial conditions Bs and the functions Cu. C., and Du,e= (where Consequently, the system ad- mits unique solution within the family Sk. The remainder of this section is devoted to deriving the closed-form expressions for the components in (76). Note that our final objective is to solve for qu and qu since they characterize the learning dynamics.
J.D... (t,-) 0) satisfy the assumptions required by Theorem
We begin by simplifying the system of differential equations using the method of variation of pa- rameters. To this end, we define the auxiliary function Fu,e, = fo fu)x and
Gt,t, = Gt, exp exp (-Fu,m,att),, Tx T (t, 2) exp (- (-2Fu,,, 2t, = [t, 2) exp (- (t)) (t) = (t) exp (-Fu,v, h E {u, (t, = f (t, 2) exp 2F (tt),, (t) Iu. (t) (t),
(77)
To solve for 9u and qu it suffices to analyze qu. To. and Fu,v.A. Substituting into (7)), we obtain the following system of differential equations
(t, (q.(()((2)+ + Ett, 01
(78a)
0 (l, = (7)(()(()))) + +E,((,)) ()),
(78b)
01
0
(t, 2) = (t) + 2.,(t (). 2),
(78c)
0t
(t, 2) = (.(()) + (G,(tt)
(78d)
01
0 P,((,))) 2) = 2π, (1, (1,,),
(78e)
0 ft, 2) = (t, (t) (t.
(78f)
0t
22

---

Published as a conference paper at ICLR 2026
① d1
e - (fenall) (l,2) 1 (t,
+ (t), + (t) (t)
(78g)
Next, we analyze the equations 778) in the Laplace domain. Here, for function f. if there exists constant a with suptzo [f()) being finite, its Laplace transform is defined [s)
So [l for R(s) > a. Treating as a constant and taking the Laplace transform with
respect to / on both sides of (78) the following linear system of equations is obtained
sC[ 2) a.m(=)) (m))) + (E.
(79a)
I z) a.)((z) = (
(79b)
(E,,((0 =) LIū. (s) + =<)
(79c)
* .=,,)(s, z) [() + =
(79d)
s <s,,)(s, m(=) = 2C ,., [( (,,,
(79e)
s ([T.(((, m==) - 20 ol (s, z) CA.
(79f)
=L( Lle 2le ((5) =6 T. + =[T
+ x[ (+, 2) + (s,,
(79g)
We first handle 19ap 79d). Solving 2.C,1 and C[G.), we have
[1, az)m(() + \sm(() (s)
(80)
sm(2) + ( + (s)
[[,,1s,,)
(81)
s² -
Choose such that s² & M. Let I be a contour that encloses s² and does not intersect M. Then, by performing a contour integral with respect to z on both sides of (80) and (8T), we obtain ar sm (s²) + Asm( Asm (s)
(s),
(82)
a am sm (s²) + Asmn( (+) + (s).
(83)
Solving the above system of equations with respect [( and L(āu]. we have
(s] agds² mm( s²) ag (m
(84)
x2 ml
m2 (s² m)²)
[(Gu((s) Or am 222 mm 82) (s²))
(85)
1 22 m (s² ms²)
We solve (t) first. According to the fundamental equation (37 have 2m(=) + 2m(2)m(z)
-1 and =m((2) + zcmm(2mm) = -. These identities yield
(s) ands mnddsmm(s?) (s²))
ms²2)
22
+1 1 +1++2m((3)
(a) + X2 (a ( au (cm (s2) + X²
1+x²
X²
(1+ x ) + 22) - s2
(1 + X2)
3cmm 52) ms{)
+a, m2m (s2)
anc(1 x2) sm()
0,,8
X
02 5²
02 5²
12
22 -2
0x 5²
a
(s22)
an a, + X²)
s?
0n
a.(1 + +22)
X(Xx 82)
X²
Dn 82
02
52
X(12 - s2)
= Qu,1 + Qu, + Qu, + Qu.4.
(86)
where step (a) follows from the identity (2m(2) (1++)++)++++0 {) + - (1) T ) Next, we calculate the inverse Laplace transform for each term Qu, separately. Beginning with Qu,I, we have
mm) =
s2 2 (Ux ()) 22
C- [-,1Q = -
Qu C-1
ax c-1
01
8
X
23

---

Published as a conference paper at ICLR 2026
(a)
Qvy
C-1
agy
/ c-1
(82 02) x) u(drr
x
82
02
A
Oyy
J x(63,0x
(1)(xd))
(87)
where step (a) follows from Fubini's theorem. The inverse Laplace transform of Qu,2 is given by
l- 1uu2l = +uc(++19))
sm((s)
and(1 + ax
U
sp(dx)
c-l
c-1
(0A- - 8²) -
A2
01 - 5²
X2
Rod(+ + / (41.00 E3.x) (f)u(d))
(88)
X2
The inverse Laplace transforms of Q2,a and Q3 can be obtained directly from standard tables: [-*[Quaa] 0, 00, (4) and [-,,4 - --((+t>² (1). Consequently, summing all the inverse transform results yields (8. The derivation for (1) is analogous, and omit it for brevity. Note that by analyzing (7a) (9d)). explicit expressions for G, and Ex, E {u, o can also be derived. However, these intermediate quantities are not directly relevant to the objects. Instead of solving them explicitly, we use their Laplace transforms.
79g to derive the explicit expression for (F,,,,,(0) Solving for In the following, we solve (9e)
el we have Clī,,,, (s, 2)
= 111 + 11 + IIs + II,, where
s. [(s)
1l2 2m(() + 2m(2)
111
8² 42
z2 - 4z
2)x CT C.T.
(s, 2) + L Lla C.) 2)
11. =
8² - 4z
As Xs|C (s. 2) Ll E., (8, 2)
174 =
(89)
82 4z
Let both s² s²/4. E M and choose contour I that encloses the MP support M but excludes
the points 5² 24 and s². We then analyze the contour integral of each term in 88 Clearly, 2wj fr dz = 0 since II, is analytic inside and on l. For term lI2, we have
f
2m(2) + 2m(2)
2
Iu + (dr)
2T fr I2zdz = 2TJ Jr
d2
da
5² 4z
2T Jr 52 - 4z J --
J (u (dx) 32 4x
(90)
We now evaluate the integral of 173. By the definition of the Laplace transform and Fubini's theorem, we have
CT +
1 f NTddz
1 ff dz
2n Jr.
2 fr
82
2Xz (t) 2) + qut
f d=
c st
2T Jr
JO
82 4z
/
| (n(()TT,,v(()e
2A go (I)TT
stdl + 2A
ddl
(91)
where (t) - fr dz =. (t.s) and 31. (t) 2=
2nj fo dz (G.)).)). Next, solve II (t),
2nj
h E {u, v}. Due to their similarity, we will only focus on II It is more convenient to handle its Laplace transform, which is given by
1 ff GGa 2)
1
(r) = 2 f =TGGal(1. dz
3II =
dz
2TJ
32
4z
3 4z
(
2
ax)m(2) + >rm(z) ((. + G( dz.
2 J
5²
4z
22
24

---

Published as a conference paper at ICLR 2026
aur (r)
2p(dx) j d2 / (52 42) 2) - 2)
2
xlaur (r)r)(())
(92)
(s² - 4r)
where we enforce ,2 M and ,² lies outside the contour P. In step (a), the identity (80) Taking the inverse Laplace transform (with respect to the variable r)) and applying Fubini's theorem, we get
31. =
On
r
u(dx)
A 1
[lT y2
u(dx)
52
4x
12
32 4x
(a)
Aooao (1,))(1) u(dx),
u(dx)
(93)
82 4x
82 4x
where step (a) follows from the convolution property of the Laplace transform ((f*g)=((()))) Similarly, we can obtain
L arIl1,x(t) p(dx)
f
Ax(au (t)) p(dz). 22-4r 4r
I 3
(94)
s² 4x
Substituting (93) and 94 into (11) we obtain
au I (t) (t)),.)) + Arq (t) (t)
2 Jr llzd= 22
p(dz)e dt
52 4x
100
[ (. (0.7((1),, (1)6, (1) + (t) * (1) ((dr)e-*tt d
22
$2 4r
22
antll (1.m + ArC (ā. ( (1,a
A= ( * (s)
/ aantLl(T l1,z) (8) + Ar (10 ( * [I,) m (s)u(dx)
2A
JR
2A
actLl(lu (1,2) (s) + XxL {lo (?u (1,x)] * 22 ) (s)u(dr))
(95)
where in step (a) apply the convolution property. The contour integral for II, can be evaluated using a similar method that for II3, utilizing the Laplace transforms of E., and f... We provide the result directly as follows
1 2 f. 11d2 - / aaxLl(a lz,x) * (s) ArC (āo l3 * & (s)
I a.. l3,x) * (ss XeLlleu @ 4x,=)1 l (s) (96)
Since 1. frr II, 2)d2 gu. we have 2 L I - By taking the 25r inverse Laplace transform and using 90
95
and (96).
obtained.
Finally, we establish the relation between Fu.v and Fa.u. Choose I' that encloses M and let |s| be large enough such that 2², 82/ /4 lies outside Using (79e). we obtain
&
1
f 22C,7,G.]((s,2ddd
xL SLle 2f 2ng f m(z)dz 2C + 2n
(97)
which yields
2F.. = 1+ S (a)d
+ 1.' da ((dx) (a) + 2x2To(a)( l.,)
(98)
Thus, have derived the closed-form expressions for ei, and Ia.e- Substituting these results into (79c)-(79h) the closed-form expressions for (f,, î., f.,..) can be obtained. By verifying that these solutions satisfy (76a (76g). conclude that the system of equations 776) has been solved.
25

---

Published as a conference paper at ICLR 2026
D.4 ASYMPTOTIC EQUIVALENCE OF THE SOLUTIONS
Let {SO), Sodet denote the solution to (45) where the initial conditions Bs and functions Cu. Co. Du, are given in (43) with (u,, 14) = Ug, vo). Denote the underlying measures for 5 as Mss
(1) for every S E Sode- Similarly, define { S Sode the solution to (76) with S Syc
for all S E Sode with underlying measures A Let I be positive contour enclosing and all eigenvalues ZZ and M. According to (Bai & Silverstein 1998. Theorem 1.1), the largest and smallest non-zero eigenvalues of Z converge almost surely to E.+ and E__, respectively. There- fore, the deterministic contour I is well-defined with probability one. We note the identification qn (t) qn (t), qh (t) Tt (t) for h E {u, v). Define
(()) = sup max zer Sde
(99)
L sup max =er
[Bc.(2) (2,)mm()))))) |Bc, (2) a.m(()), (Br m()))
Br, (z) -m(2)). Bz, (2) B=. (3) (BT., (2)1
[C.(z) - (C) -m(=), IDu.,
(100)
In what follows, we show that Mt) + 0 almost surely. By definition, we have 1g (t) d0 (t))
I (t) - q:2h 2 (t)). (t) - (2) (t) < MT) where [T| denotes the length of r. Further-
more.
[t)(t) (a)
(1 (1+ F 2T
(1
A 0
1 + 1lo (2) ))
M().
(101)
where the variables l are omitted. In step (a). the bound 14! < follows from the existence and uniqueness of the solution (42) and its agreement to the integral representation in Theorem 3 as
discussed in Section D.2 The inequality 1q 1 < can be derived similarly. Now, by considering the difference of the equations related to E.,, obtain
< (z)1 + / (.) - |
[t, - - (t, (t)
(a, - Ga! z)
-
+ > a (a)
D (2) (=)l + (a) (c)) (a,
1l (a) ((2 (a) (a,
(102)
Due to the continuity of the good transition kernels and the finite total variation of the measures, for any fixed > O, there exists constant U such that supogier (fa)) & U. Furthermore, for every j E (1,2). S Sode and z E I, we have
Ius
(1, dx)
sup (t,) o
sup OE
sup d=1 UE
1g (t,
< U,
(103)
TV
almost surely, where d denotes the distance between the (common) support of the measures and P.
Taking the supremum over on both sides of (102), we get, for every t E [0, T).
< L+ 1' 2T
(1 (A) + v) M(a) - XLda
sup zer
- (1, -
sup
|
SET
2x
C. C2 / M(aada,
(104)
2 almost surely, where C₁ (1 + XT) and C2 - co. In contrast, the constant C2 depends only on
supzer |2| (1+2A)UJE U. By Lemma
4, it follows that C₁ + 0 almost surely as p.n
the fixed deterministic contour I and the bound U, and thus remains upper bounded. This argument extends all other components 5(1), S(2) E Soaee leading to the composite bound
M(t) < Ō E. J( M(a)da,
(105)
26

---

Published as a conference paper at ICLR 2026
almost surely. where T1 upper bounded by polynomials in L and C2 is bounded constant.
According to Grönwall's lemma (Lemmal D. we obtain
sup M(1) & Ci exp(C2T). os
(106)
Again by Lemma we have SUPOEIET Ml)) - almost surely, as p, n + 0o. This implies, in particular, that for every h E {u, v).
sup e
(t)
(2) 9h
(1) & sup M(1) O,
(107)
051ET
almost surely. This completes the proof.
E PROOF OF THE BOUND (18)
In this section, we analyze the gap between GD and gradient flow. Define
{| {|uuu L 0r ग 3
A
max
(108)
Using the fact that Ut is orthogonal to (Ip - Uu u
Xxx for any h. we have
(1 + n²² [1 tu X xxxl12) -
[u + (Ip
-t 2
X 0k
1 - k", [[[p uo u6 ) xxrel1
(109)
where G uniformly bounded term. In particular, we have |. & supzoot .vT+=
-TFT-1
}. By
109). we can approximate us++ by
d + II 26 UR Lu + IIp 2 U
0k++1
un nII, -ū,, 26 X Üs + +1 8 7,ke
(110)
where the norm of the vector 8 is uniformly bounded in k i.c., sup 8 Mulla & C 11 +
7 X l12) for some constant C > 0. A similar argument on UR yields
Uk+1 Us (I, 0g 0g X
üs + 17 67 n.k-
(111)
where the norm 8 is also bounded by C + 1 (1)) uniformly in k. Therefore, by the gradient flow and identity (110), we have
us + 172 1/28% 4.
U (k+11)
uk++1
Uky
(k+I)m
[(I, un un Xua (Ip u ) ) Ub da.
As result, by the inequality
(Lp u u. Xxi (Ip - uuzu )X 11 < || 02)112
+ [|(u) U2) u Xvi ||| + || U22 (uT -u X ||2
+ | (w) - 02) 0|
< 2 x|||X|, [lli]] - u2l12 + ||001 02112l.
(112)
for any unit norm vectors u; and v. E {1, 2]. we have
u(++100 Wot1/la 1 < | - u 1l2 + 11² 15% 1%/,0112
(k+1))
+ + X l
|/ uellz + ||
3z da
(k+1)9
< Akn + C |
+1
+ +
Aada.
(113)
km
27

---

Published as a conference paper at ICLR 2026
The same argument || (1++j Bk-1ll2 together with (13) yields the inequality (M+1) & Doy +
x
(1 +
X )rp
+
X
Jkn 1(++1)n A,d which immediately gives
ky kky + || /
km
C +
A,da,
(114)
for any non negative integers Then, for any and 1, we have us üg ue ukn + Ug and
Unn
|| uu - uxll2 & Dnn
[[ - un
Xv - (Ip uruf X da
kn
+ (t kn)
(Ip - u X 0
< 4 | A.da - C | 2 (1 + 17 122 + x ||-
(115)
Using the same method as in 115 to bound the term ||0
Th 2- we can obtain
A < | 2 /
Aada + C | (1
|X ko) + |
(116)
By Grönwall's Lemma, we have
sup A. < CII 1 Ty + x| exp [4| T).
(117)
DE
Therefore, we have derived the bound for A, with given T and 17 > 0.
Next, we study the asymptotic case. In fact. by YYn et al.
1988 Lemma 3.1), we have
nim 7 l (1 + vC)² almost surely, which implies [
2 < (1l="((*) " ll2 + | <
+ (1 Ve)² almost surely, p,
Therefore, we have P(lm sup,, SuPo A &
C'n + n)²)) 1 for constant C 0 which only depends on T, X, and c. Therefore, have
completed the proof of (88
F PROOF OF THEOREM 2
CASE I: >> e4/1
We begin by analyzing the asymptotic order of (t( t - We express (8) in the form (t) = (t) + ZZz()) (t(t) and evaluate each term separately. Through direct calculation, the following obtained
Vose
eoxe
(1) la,0x (1)
200x(U> - x) 2 200x(1> xx)
c Vaa
evac
2/z(1 - Ox) + 2 -
(118)
Since X2 > Ve, follows that U2 > (1+ vo8²= = E.. Consequently, we have
On evvae
Vievi!
(1) = À 2V
[(dx)
Oh
2X JR. 02 u(dx) + o(1)
JR 0x - x
evost
av EE Voeval
V(x E_
(EE x)
dxr o(1)
2/0213
2X
I. i2
2xx
evorla.
au E. E. (E+)))
eVE
E
(119)
200213
2Av2 E+)
13/2
where step (a) uses properties the Stieltjes transform of the MP law Couillet Liao
2022
0xm(Us) = - max(cX-2 2, X2)
(120)
28

---

Published as a conference paper at ICLR 2026
and step (6) follows from Lemma
For the term (2) we have
event
andevose
Zz(t) and
2002 +0(1).
(121)
2/0x
Similar to the evaluation of (i))), and using the identity
Txt
evext
C e-vat Vat
eVar
(3,
(122)
2(0x> - x)
2(0x -x)
2(1xx - z)
2(1x> - x)
we can get
a + 12)
/ (((x))
E+ ev V(F+ x) F
Zs(t)
dx
+ +0(1
212
J 02
(0x - 2TTC
E__
Cnc(A12 + 1)evoR
a + X²) x2) E+--- E_ eVE+1
20 2T))2((E) f0/2 +(eVV++) (123)
2(1² + cc)2²
By definition, obtained
Zt(t) a.evost
12 o(1). Combining these results, the following estimation is
(X1 - c)
(t)
au
au
2A2 / +C
X2
LVT+X VX
+c
v E+ E_
(Et(E)) au(1+ (1+)2)
eVE4
f3/2 + ooee v EEfe 3/2
-B2 //2-3//2)
2v2m(U>> - -E]]
A
12(E4)1/4
= A5 + BBeVEtity 3/2 + ole v -012)
(124)
Similar to the evaluation for (1), we get
(1) = Age + BBtt 3/2 + ole v
(125)
where
A, (X* c)
agy
Ou
(126)
2X2 X22 +
v1 X2 "aal (aa(E+)
VX2 ce
E. E_
co(c+22)
B, = 2 2 (0x 100--E44 E+)
(127)
CX
Next, we will analyze the order of the denominator p(t) for the following two cases.
F.1.1
#0
Under this condition, the solutions have the asymptotic forms (t) = AgeVext + o(evext) and (t) Ave vxt). To determine the asymptotic order of p(t), we first evaluate (t).
Writing 1II as (1) - G)()), we focus on the evaluations for G₁ and G2 as the remaining terms follow similar estimations.
Given 8* > EL, we observe - - contribute the leading order asymptotics. We then amalyze e e=2/0x+2((t) By definition, we have
- G1 + 0 as t + co. Therefore, the term G1(t) does not
20 2V (t( = 200
e 2VTx (1x) * (2xx)
JR
+ 2X
e e-2200'zz (ā. * (1,x) * (2,x)(d)) G2 + G2.(()).
(128)
RR
To analyze the limits, we apply the Final Value Theorem (Lemma 2). Since e
--vovl (1,x)()))))))))))))))))))) - 0 and [[ = 1/4(0x - I), Lemma implies G2,(()) + as + do. Applying Lemma again yields
A
lim c -o (t) O2 r
(129)
29

---

Published as a conference paper at ICLR 2026
which implies
VVXA2
lim 2v0xx (10o (āu ( +1,=]],]] x (2,))(t) 4(0x x))²
(130)
By the dominated convergence theorem and the properties of the Stieltjes transform, have
lim (2,((1) = 222. VuxAg t-00
02 2( c)
(0x (0
(131)
The asymptotic orders of G3. G4, and G5 can be evaluated similarly using Lemma the results gives
Gathering all
lim e 2-20 (t, = X2 XA ccA.))
(132)
7000
-
As result, we can get
AA + +cA
A² + 1)A7
lim 2(1x 1-000
X1
0xA2
X+
2
Oyy
au
(133)
4x2
VX +
by the observation fo f(a)da (f 1)) Therefore, provided that Au + 0, the limit values of (t) and (t) are given by lim;440 (t) = I respectively.
1-ek-
and lim (t) = I v 1-cx
Itcx
1+X
F.1.2 A, 0
According to Lemma
we have
+
eB+ eVat V(++ xx E-)
Gt
dx
8T
JE_
VI
(1+ (1+c W EE E
+ o F11g-113))
(134)
16
f3/2
Thus, there holds
V Est
((t) > 2 LL Gi(a)da > K 13/2
(135)
for sufficiently large and K > 0 a constant. Therefore, we have
(t) (Bu + o(1))e V(1))) E
(136)
VKev STevE+t4--0 E+44-3/4 which implies lim;- (1) 0. The relation lim To (1) be obtained similarly.
F.2 CASE II: À < 1/4
By applying the same method used in 124 and 125 obtain (1) olev (EF1)-3/22 h E v). This result aligns with the case discussed in Section F.1.2 which implies lim Th (1) 0, h (u,
B
F.3 CASE III: A 1/4
In this case, have 8x E.. which precisely the right endpoint of the support of the MP distribution. Fortunately, the Stieltjes transform m(x) is right continuous at x = E_ (Couillet &
Liao 2022). Using Lemma Bl and employing estimates similar to those in (24 and (125). we can obtain
Th (t) = (Ch + o(1))e V h (u,
(137)
where C. and C, constants. However, need to make precise estimate for p(t) in the denominator, as the upper bound in (136) increases to infinity. To this end, the following lemma is useful.
30

---

Published as a conference paper at ICLR 2026
Lemma 5, Assume b a a>>> 0 and f E C((a,b) b)) such that f(b) 0. Then, for any sufficiently large t, it holds that
evTy ]²
K K & [ (lo)vb-xeeue
vaydy
dx & Kye (138)
V
where Ky 2 Kr. > 0 are two constants.
Proof: It suffices to consider f(x:) and restrict the integration to _b - E, b] for small enough since evat concentrates near = 6. For given 6 0, we choose 0 such that supueju.e) [Jul & 6, where 8, is defined (166). By changing the integration variables via v = b x. a ut///. and y/t, we obtain
/ evby yy
2
e
2/bby112
Vay
dy
dr
b-8
Vue 1/2
Atbylal
No y
(1 20b
dy
2
du
b
Vae (1+80.
[ B
te
a dlo
d3
2
da
1
(139)
We next estimate the inner integral to' 8 8 atogaaa integration interval at some E (0, 1), we have
d and derive bounds for I. By splitting the
B t alod
(T.
(1 a
d3
3
d3
2e (
(144 Jar
(140)
Vra(1 80,1)
Therefore, by [8 < 147 and 140)
we derive an upper bound
2c 33-3
) Vare (- 1)
I 144+4/ where I, is some absolute constant. can be verified that ly < oo for any T E (0,
Ta3 8)2 a1/2(1 8) da lus
(141)
For the lower bound, again using 16a,2 < 8, can obtain
1 1> 1 + Voe (1+5)a
3 be
(1-)aa
d3 da
(1-51a
e-
> 1 + Vae (1+8)a
da > 1
(1- 3/2 da.
(142)
8)a
Taking & + 0o, we have lim inf +00 completed the proof for Lemma 5
B 1 for some constant 1r > 0. Therefore, we have
To determine the asymptotic order for p(t), we first give a lower bound for (t). In particular, will evaluate the term fglğ. (ÿ, {1,=)] * (xx)) Noting that the orders of l1,z and &z are dominated by evēt/p /2 and 2²v e2V4v5//4, respectively, it suffices to study the simplified integral Ja VI qu (g ev 2VE ((dz). Rewriting this expression yields
10 [ Vzlg (7r (e 2VV p(dx)
= [ V I (a(()) ( To (b)e 2 2VT(L-a)
2
/ Vaevan (f (a(e Van
(x E_ )EE x)
dr
2xer
31

---

Published as a conference paper at ICLR 2026
( [ (o(a da 1
VE+ zdæ,
(143)
E__
where f(x) vx E- [-))x((ATCC) According to 137 we know that there exists constant A > 0 such that inf 17 (t)( E+4
(e.(t)/eeVEx) E+11-11211 > Co/2. Thus, we have
e t L0.0
2v 1 f f(x)e2v3)
Sr (eVV--vai
Vaja 7
da
Vx+ rdæ
E
Va
(a) t 1/t F+ [(0)
1o0 e(VEE Vaja
2
da
-E+-adar + o(1),
(144)
Va
E+ VV) VA contributes only lower where step (a) is due So² TTTT- da 20
der termss By applying Lemma we have lim inf e / constant Kr. symmetrie argument establishes the corresponding upper bound, yielding Io.
Ie,v > Ky for some
lim s 2V 'tà Le, & Ky for some constant Ku 0.
At this stage, we can give the upper bound for (t). h E v.. Combining with (17)) we obtain
(t)
(Cn Fo (a) CneVRiLy h E (u,
(1) & pt)
(145)
V22² So Lundt
eVEsty-1/4
where Cr a constant and step (a) follows from the fact that the ratio So de (c2 (c2VE=1 converges to positive value as - co. Therefore, we have completed the proof for Theorem 2
D
G PROOF OF LEMMA 3
According to Taylor's 's lemma, we have
vb-u vb 1 (1 + 8m)u
(146)
2b
for small u. By the continuity of f. for any given 8 > O, we can choose sufficiently small such that
sup
1f(x) f(v)l < 8, sup
lóul & 6.
(147)
Define M sup supzela.t 2
[xx)
A(t)
Sh f(x) explVit](6 - z) and B(t)
f(b) exppe/tpliieeccceeeeeeeeeepbbeee bbbl + 1) 2vb a+1
By calculating |A(t)/ B(t) 1/, we have
A(t) B(t)
1|=
1 b-m
f(b - u) expl utjuu du
B(t) -
f(b) lo expl lbllbb- utju du
8 f expl df
- +
B(I)
B(I)
M expl vb Et
1a-a du
= Xi + X2(t) Xz(t).
(148)
Btt)
It is clear that lim, 00 (t) O. Next, evaluate (t). By plugging (16) into 148 and
using I(a 1) Io expl-vly" ulyddd dy. obtain
f(b) expl / bt) 10 exp
(1+64)at 2vb
du
(a) 235 exp [- + 8, y
B(t)
Jo exp exp[-ylyad yy
32

---

Published as a conference paper at ICLR 2026
So cxp 1- Jo uypl
5) dy
$ 1
(149)
8)a+13
where in step (a) we change the variable utt (2vb.. Similarly, we can derive the lower bound
explV/t] lo exp
(1+6, Jut
u du
2vb
inf
(150)
(())
(1+ 8)
Combining this with the corresponding upper bound implies
lim (1( max {11 5)-a-1 1, (1 5)-a-11
(151)
t-o0
An analogous analysis yields
5
lim sup X1(1)= &
(152)
[(bb) 8)
t+000
Since 8 > 0 can be chosen arbitrarily small, follows from (15T) and (122) that
A(t)
lim t-x0
0,
(153)
B(t)
which completes the proof.
H PROOF OF LEMMA 4
To simplify notation, let C denote a constant and C, denote a constant that depends on a. Their values may depend the context. By standard truncation argument (Bai & Silverstein 1998). suffices to assume that the elements of the random matrix Z are bounded, i.e., Ivn2.gl & C, for some absolute constant C, this truncation does not alter the asymptotic locations of the eigenval- ues.
We first prove that T (())Zo E (()) Z)) converges to almost surely for given 2 E C+.
To this end, we control the moment by the martingale difference argument and the Burkholder inequality (Bai Silverstein 1998 Lemma 2.1). Write Z [Zi,..., pn). v = [EI...., Un) T Z, = zn 0j [uI...., -- U5+1 v.]]" and Q2(z) = -Z;
[zj...., æj-1, 21+1,, Zj+1..... Define the quantities
- xIp)
1
a²(2) = Tr (2),), n)(2) 1+))(2)) 8,(2) = + 2] ()())
71
() - -z]0((z)2;
-
(154)
Let E,(.) = E(1|21...., zj) and Eo(.) = IE()) By taking the difference between u." T T(2)) Zv and
Eu Q(2) Zu, we have
n
Q(2) - TE > E, Ey-1lu
T Q(2)
Eg Eg-1 Q(z)Zu-u u - Q(=)),u,) QQ)) u (),,())Z,u,)
W1 + W2.
(155)
We first evaluate W.. By the Woodbury matrix identity, we have
Q Σ IE - E, 1+ 22(2)zz
W-
Σ
E, Ez-1Ju a Q(2)2,0;
j=l
1=1
=
Σ EB - Eg-133,n.u TQ, + EEnju (,(2)),11 = Σ
W1,15 + W1 (156)
1=1
j=1
33

---

Published as a conference paper at ICLR 2026
since the identity 8, 73 + B,ng and the independence between 2j and Qj- By apply- ing the inequality 1 - - & - the Cauchy-Schware inequality Bylayl & E3 |e| lul?, and the trace lemma Bai & Silverstein, 1998 Lemma 2.7), we have
2 (x(2)
C|2
E
1 u T
Q2 (2)
'E
(157)
(2)
where we use the bounds Zhuang et al. (z). and the moment bound E 1
2025
Hachem et al.
2007 bn (=)l |=1/3(z), 18,(z) <
Ct
On the other hand, we have
n
a Ck|=/1k Σ (2)
7
C 1k n2)3sh(=))
1.11 28
oj
2
10 u
24
2
(158)
j=1
5=1
The evaluation for
E3 Bj- W1,2 can be given by
1,25
2
oz
Σ B 1
B u Q,(3)2z}
(159)
(2(2)
Similar to 1588. we can get Σj E. W2,2 2 & Ck|213k (3:4(()))) Hence, according to the Burkholder incquality, we have
C2(12 4h + 2M)) no (38h + (4k(2))
E Wil 2h Ck
Σ 1=1,2
E, 1 W1.15 ?
+
W1.0 2k
(160)
For k
the RHS of (16)) is summable, which implies W1 - 0 almost surely by the -entelli Lemma.
We then evaluate W2, which can be given by
W2=
[EB - E5-1,8,u (,(5))zz0Q,(()333)
j=1
>
[R (5) /, 1,8,A, [=)())===]0
[ (=)
j=1
W2,1 + W2.2y-
(161)
The term W2,15 can be written as
1 W2,1 Eyn, u (2)222)2 Qz(2)
(z)
(162)
By the trace lemma and || || 1, we have
C
E W215
B| 2k
(163)
(2)
n n2kg
which implies W2.15 + 0 almost surely. The relation W2 + 0 can be shown in a similar manner, and we omit the details.
Next, we prove that Eu." T ([)2) = 0 (.)) We can show that
u Q(=) Σ
u T Q()));)
Σ 3,())u (2)))))
(164)
j=1
34

---

Published as a conference paper at ICLR 2026
Since (2(2) and xy are independent, we have
Eu
T Q(=)
E
By (u T Q)(2)z;) (8, - nz)
Σ ==1
2g Eu Q2 (2)
(165)
Using the Cauchy-Schwarc inequality E|zyl
< Ed | 2 /01
we have
|x| 1
Eu Q(z)
E 2 (z(2)z)
2 E |
(2(2)
C|2
C|2l2 (z)
My
(166)
(=)
This completes the proof for the case E C.. The same approach can be extended to analyze the E.. Furthermore, for any 0 and k > 0, it can be shown
cases where 312) 0, 0, and that,
((u ((uTQQ())) Zz))Zol <) Co. nn
(167)
for all d sufficiently large, where Cd,= is constant depending on d and dist(z,
Next, we prove the convergence for general compact set I. Since I is compact, we can choose an n -ett Nn. C I such that for every E r there exists z' E NN with |= - z'1 & n 4, and for any two distinct points x, y E Nn satisfy l vl > n-- 1. Then, have
Z 2z 2'
n Q Zu - u 2 Q(z)Zzul
(168)
dist 2) dist(A, 21)
where A {Aihessp are the eigenvalues of ZZ
Inequality 168) gives
sup =eT
Qzz)Zz
sup zET
sup s'EN u Q(2)Zu u TT Q(2')Ze
Tu Q(z')Ze
1 (A, dist(A, NEN,
sup
|u Q(')Zoul = X1 X2.
(169)
According to the no-eigenvalue property (Bai & Silverstein 1998 Theorem 1.1), have dist(A,,) r' + dist(M, A) almost surely as p. + oo and Il1, is almost surely bounded. Hence, X1 almost surely, as p, cc. By (167). for any and k 2, can choose sufficiently large such that
P(22 c)
U u Q((')' >E
Σ P ( T2()) 2c)
'EN,
ZEN,
card(Nn)Car
(a) 4(supzer |2| edn
1)
(170)
ednn
where Ca, is a constant depends on and dist(T,MM), Step (a) follows from the fact that, by the definition of n à -net, the disks of radius 1/(2//) centered at points in Nr, are mutually disjoint. Since 2, the RHS of the above inequality summable, which implies X₂ 0 almost surely by the Borel-Cantelli Lemma. Therefore, we have completed the proof.
I ADDITIONAL EXPERIMENTS
In this section, we provide supplemental experiments. The simulation settings, specifically the GD algorithm, consistent with those described in Section
Learning Dynamics with Additional (p, n, A) Settings: In Figure we plot the learning dynamics with broader range of triple (p. n, A). The settings here are consistent with those in Figure] 1 and
35

---

Published as a conference paper at ICLR 2026
077
Vile of 92161
1.5
V4lae
0.3
1.0-
08.6
0.2
633
0.0
(a) p 8, 11
20, À
1.3
0.
0.7
1.0-
Vf 3,00
2000
0.0
VA of
0.5
10MM
0.5
10
(b) p
40, n 100, À
1.3
2.0
0.6
0.4
0.4
0.10
Vokue R102
0.2
1.0
-42
02
0.0-
(c) p = 100,
200, À 1/4 + 0.05 = 0.891
6:00
0.4
0.2
of 3050
622
1.0-
0.0
0.5
-00
0.0
(d) p 100, 200, 4/4 0.05
0.791
6200
Vilue of 6:00
Vilce of 3250
1.0-
0.5-
(c)p 1000, 2000, c1/4
0.05
0.891
Figure 4: Theoretical and Empirical Learning Dynamics with Different Values of p, A.
the legends are omitted for brevity. From Figures 4a and 4b we observe that as the dimensions increase, the deterministic approximations qu and do become more accurate. Furthermore, it can be
36

---

Published as a conference paper at ICLR 2026
0.30
01200
0.15
0.1
003
100
100
130
300
(a) 0.5
(b) 3
Figure 5: Approximation Error Low Dimensions.
noted that, with high SNR, Theorem I accurate in low dimensions. However, the convergence of the deterministic approximation with low SNR is notably slower than that with high SNR.
In Figures 4e the parameters Ous Og, and e theoretical learning dynamics, i.e., the solid curves, are identical. is shown that the convergence speed for the deterministic approximation very slow, even for SNRs above the critical threshold. A comparison between Figure |4c and Figure 4d reveals that, despite the systems being in distinct phases, the learning dynamies curves closely resemble each other when their SNRs are similar. This suggests that signal detection may still be feasible even when the SNR is below the critical threshold.
be the same. In Figure 4c and Figure 4c the
Approximation Error: Figure [S] illustrates the error of the deterministic approximation in low dimensions. The initial points are set to satisfy
02, 0.447. In Figure [5a the ratio is set = K 0.5, and in Figure (bb] the ratio is set to e = & = The error is defined as E su T. (4) Yu (4)l, To (z) (o))). Here, the y-axis represents the average error computed over 50 independent trials. can be observed that the error decreases as the dimensions p and increase the same pace. Additionally, the error is smaller when the SNR is high, which is consistent with the observation in Figures [
37


Reviews:
[{"review_id": "RhRexnAuEh", "reviewer": "Reviewer_DPd9", "summary": "The paper analyzes the statistical properties of gradient flow (as a proxy for gradient descent) for rank-one matrix denoising under a deformed Wishart model (rectangular matrices) with noise that has i.i.d. entries and . It derives a deterministic term for the limit of the inner products between limiting singular vector estimates and ground truth singular vectors in the asymptotic limit of matrix dimensions becoming infinite with fixed relation fraction. As an implication of these results, the authors are able to relatively accurately predict the behavior of gradient flow depending on the problem's signal-to-noise ratio (SNR) threshold akin to the BBP [Baik, Ben Arous, Péché 2005] phase transition. The results also can be used to quantify the dependence of the dynamics on initialization (value of $\\alpha_u$/$\\alpha_v$) and reasonable stopping times can be theoretically derived (see Remark 1 and Remark 2). From a technical perspective, the results lean on the analysis of [Bodin & Macris 2021], who have showed similar results for the symmetric case. Simulations are presented that substantiate the qualitative accuracy of the asymptotic analysis in the finite sample / matrix dimension case of $p$ and $n$ fixed.", "questions": "1. In lines 141-145, it is mentioned that $\\operatorname{grad](\\cdot)$ is a \"Riemannian gradient operator, which enforces the unit norm constraint\". However, I do not see that the update equation of (4) enforces such a constraint. In some sense, this is statement is incompatible with the framework of Riemannian optimization as the Riemannian gradient lives in the tangent space onto Riemannian manifold and requires a retraction back onto the manifold (here, enforcing unit-norm vectors) to enforce the constraints.\nCan you clarify or correct this discussion? In particular, how does your studied gradient flow algorithm enforce the unit norm constraints throughout its flow?\n\n2. What are the limitations of the presented analysis for higher-rank ground-truths? Where does your current analysis fail to go through?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The analysis presented in the paper seems to be new and studies a foundational problem in high-dimensional statistics / linear algebra, the behavior of singular value decomposition under the influence of noise in the case of rectangular matrices. The noise model is rather general, which is positive. It is of interest that the gradient flow dynamics more or less matches information theoretical phase transitions that are intrinsic to the problem. \nWhile carefully checking many proofs in the appendix was beyond my abilities in the allocated time-frame as a reviewer, the results are plausible from a perspective of a reviewer who is familiar with tools for analyzing non-asymptotic high-dimensional problems.\nBeyond covering the asymmetric case, some assumptions are weaker than in the related paper [Bodin & Macris 2021], such as the finite fourth moment assumption (as a opposed to assuming existence of all moments).", "weaknesses": "A fundamental weakness of the work is that it applies only in the high-dimensional limit of $\\lim_{p, n \\to \\infty} p/n = c$, which is in contrast to many analyses of iterative algorithms in machine learning. Related to this issue, it can be pointed out that the title containing \"Gradient Descent\" is to a certain extent a misnomer as gradient flow, which is less relevant than gradient descent in practice in machine learning, is being analyzed. Thus, a lack of treatment of the discrete-time gradient descent method is a weakness of the paper given the framing of the paper.\nA more unified discussion pointing out the differences and similarities between a power method algorithm for computing the leading singular vector pair and the presented algorithm would also have been insightful - I somewhat disagree with the framing that \"SVD is intractable\" as it is clear that a reasonable algorithm for the problem would involve a partial SVD implemented via randomized techniques [see, e.g., Martinsson, Tropp 2020].\nFinally, it can be be pointed out that, while the asymmetric case being more challenging, the analyses / simulations presented are relatively close aligned to the ones of [Bodin, Macris 2021].", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "UEPHWT1bHp", "reviewer": "Reviewer_YymG", "summary": "This paper studies the gradient-flow dynamics of the rectangular rank-one spiked matrix denoising problem in the high-dimensional limit.\nUsing tools from random matrix theory and Laplace transforms over the Marčenko–Pastur law, the authors derive explicit deterministic equations describing the evolution of the overlaps between the estimated and true directions.\nThey prove convergence to a deterministic limit, identify a BBP-type phase transition at the expected signal-to-noise threshold, and show that the limiting alignment carries the same sign as the initialization (“signed BBP”).\nA kernel-based argument establishes existence and uniqueness of the solution.\nExperiments nicely confirm the theoretical predictions.\nOverall, the work extends the analysis of Bodin and Macris (2021), which treated the symmetric Wigner case, to the rectangular Wishart setting.", "questions": "I) Clarify the exact novelty relative to previous analyses.\nThe paper’s results are very close in spirit to Bodin & Macris (2021, arXiv:2105.12257), which already provided deterministic gradient-flow equations and asymptotic limits in the symmetric spiked Wigner setting.\nIt would help to spell out precisely which technical steps differ in the rectangular case and which parts of the proof had to be redone — for instance, changes in the resolvent structure, contour integration, or Laplace-transform kernel.\nAre there specific mathematical obstacles that make the rectangular case substantially more difficult, or is it mainly a matter of replacing the semicircular law with the Marčenko–Pastur one?\nA short paragraph clearly highlighting these differences would make the contribution much clearer.\n\nII) Connection to the broader literature surveyed by Macris.\nThe “Related Work” section of Bodin & Macris (2021) already offers a remarkably complete overview of the theoretical ecosystem surrounding gradient descent, AMP, and high-dimensional inference. Many of those references are directly relevant here.\nThe Bayesian analyses of the spiked Wigner and tensor models (Korada & Macris 2009; Barbier et al. 2016; Lelarge & Miolane 2018; Lesieur et al. 2017; Perry et al. 2020) provide precise information-theoretic benchmarks in the form of mutual information and MMSE.\nThe dynamical behavior of AMP and the existence of computational-to-statistical gaps (Barbier et al. 2016; Lesieur et al. 2017) are also well understood.\nGiven this context, could the authors explain what new qualitative insight is gained from their explicit time-evolution formulas?\nFor instance, does the analytic expression for the transient trajectories reveal any phenomenon that is not already implicit in the AMP state evolution or in the energy-landscape picture?\n\nIII) Relation to the matrix–tensor and Langevin dynamics literature.\nMacris notes that recent works (Sarao Mannelli et al., 2019; 2020) analyzed the optimization of mixed matrix–tensor inference problems using integro-differential Cugliandolo–Kurchan (CSHCK) equations — a fully dynamic, spin-glass-inspired formalism.\nIt would be interesting for the authors to comment on how their much simpler kernel/ODE formulation compares conceptually to those dynamical equations.\nDoes the present approach capture similar information about the saddle structure or convergence rates, but in a more tractable regime?\nOr is it strictly a deterministic “mean-field” limit without the stochastic thermal components appearing in the CSHCK-type equations?\n\nIV) Discrete gradient descent versus continuous gradient flow.\nSince the work focuses entirely on continuous-time dynamics, a natural question is whether these results extend (even approximately) to discrete gradient descent with a finite learning rate.\nPrevious works, such as Lee et al. (2016) and Ge et al. (2017), established convergence for discrete GD under the “strict saddle” property, while Saxe et al. (2013) and Mei & Montanari (2019) studied learning-rate effects in linear and nonlinear models.\nCould the authors discuss whether the signed-BBP phenomenon and transient behavior persist under discrete updates?\nEven a conjectural statement or some preliminary numerical evidence would be welcome.\n\nV) Connections to the energy landscape and spin-glass literature.\nMacris also situates the work within the broader context of non-convex optimization in random energy landscapes (Subag & Zeitouni 2017; Ros et al. 2019; Auffinger et al. 2013).\nIt would be valuable to comment on whether the present gradient flow can be interpreted as traversing a spin-glass-like energy surface with a small number of global minima and exponentially many saddles.\nDoes the deterministic flow derived here correspond to the typical trajectory that avoids these saddles in the large-n limit?\nSuch a discussion would help bridge the current mathematical analysis with the well-developed physical intuition from statistical mechanics.\n\nOverall, the paper is technically clean and the results are credible, but the authors could significantly increase its impact by situating it more deeply within the broad theoretical lineage summarized by Macris — including AMP and Bayesian limits, non-convex low-rank recovery, spin-glass Langevin dynamics, and deterministic gradient-flow analyses.\nClarifying what the present framework adds to that landscape, and where it could go next, would make the work more compelling for the ICLR audience.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The results are technically sound and contribute to the understanding of optimization dynamics in high dimension.Extending the analysis from the symmetric to the rectangular setting is nontrivial and closes a natural gap in the literature. Conceretly:\n- Provides explicit, analytic trajectories for gradient flow in the rectangular spiked model.\n- The results are rigorous and match empirical observations very well.\n- The “signed BBP” effect offers a clear dynamical interpretation of initialization dependence.\n- Writing and figures are clear; the paper is enjoyable to read.\n- Strengthens the theoretical link between random matrix theory and learning dynamics", "weaknesses": "- The work feels somewhat incremental relative to Bodin & Macris (2021); the main novelty lies in adapting the approach to the rectangular case.\n- The “signed BBP” result, while nicely explained, is largely an expected property of continuous gradient flow.\n- Experiments are limited to confirming the theory in the simplest setting; discrete-time or noisy gradient dynamics are not discussed.\n- The discussion of related literature could be more complete and precise.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "aWEbRDTASc", "reviewer": "Reviewer_RZPg", "summary": "The authors address the problem of rank-one matrix denoising, focusing on the gradient descent dynamics underlying this process. The main contribution of the paper lies in two theoretical results. The first theorem establishes a deterministic approximation for q_u and q_v, two widely used metrics that measure the alignment between the ground truth and the estimated components. The second theorem characterizes the asymptotic behavior of these deterministic approximations.", "questions": "N/A", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "The paper is well written, well structured, and clearly explained. The theoretical analysis is rigorous, and the overall presentation is easy to follow. I find the work interesting and relevant to the study of optimization dynamics in low-rank estimation problems.", "weaknesses": "My only minor concern relates to the clarity of the experimental results. In Figure 2, which illustrates the effect of the critical SNR threshold, it is visually difficult to distinguish the different curves, especially in the middle subfigure. I understand that it is challenging to convey a large amount of information within a single figure, but since this result serves as an important validation of the theoretical findings, it should be presented more clearly to enhance its impact.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "4pBfautkKM", "reviewer": "Reviewer_he9F", "summary": "This paper studies the rank-one Jonstone's spiked model in matrix denoising problem within random matrix theorem (RMT) scenerio. The main contributions are two theorems. In theorem 1, they derive a closed-form deterministic approximation for the inner products between the learned vectors and the ground truth when ratio $\\frac{p}{n}\\to c>$. In theorem 2, they show that $c^{1/4}$ is the critical threshold for the gradient flow estimation.", "questions": "In Remark 1, the assumption of knowing the SNR $\\lambda$ is practical in reality?In experiment, what will the results change if $p$ and $n$ are not that large, like $p=20$, i.e. beyond the RMT scenario.", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "- Use random matrix theory (RMT) derive a closed-form solution for the learning dynamics of matrix denoising problem. By RMT, naturally extend the matrix denoising problem to high dimensional scenario.\n- The theorems and derivations are solid, and for example, the complex but precise expression in Theorem 1 is derived.\n- The problem and assumptions are stated clearly.\n- Provide a more comprehensive understanding of the dynamics of gradient-based learning in high-dimensional matrix problems", "weaknesses": "- Lack of explanation of the application of matrix denoising in the random matrix scenario, i.e. $\\frac{p}{n}\\to c$ with $p, n \\to \\infty$.\n- The statement about computational complexity is too vague, like \"the  complexity is affordable\" on line 220 and \"We note that $\\hat{t}(\\alpha_u,\\alpha_v)$ can be efficiently computed by standard numerical methods.\" on line 244.\n- On line 244, need more details of the so called \"standard numerical methods\".\n- On line 42, \"Extensive research has shown that\" lacks reference.\n- The theoretical work on which the article is based is very classic, and there is a lack of reference from recent new theoretical work.\n- Lacks explanation of Riemannian gradient operator on line 144.\n- In experiments, only one set of $(p,n)$ is tested, need more testing sets.\n- In experients, the critical threshold is about 0.93, but the SNR is set to 0.3 and 1.5, which are far away from the threshold. More attempt on SNR need to be tested, especially the SNR near the threshold.- The ground truth is too trival.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)