Eigenvector Overlaps of Random Covariance Matrices and their

Submatrices
Attal Elie∗, Allez Romain†

January 16, 2025

Abstract
arXiv:2501.08768v1 [math.PR] 15 Jan 2025

We consider the singular vectors of any m×n submatrix of a rectangular M ×N Gaussian

matrix and study their asymptotic overlaps with those of the full matrix, in the macroscopic
regime where N / M , m / M as well as n / N converge to fixed ratios. Our method makes use
of the dynamics of the singular vectors and of specific resolvents when the matrix coefficients
follow Brownian trajectories. We obtain explicit forms for the limiting rescaled mean squared
overlaps for right and left singular vectors in the bulk of both spectra, for any initial matrix
A . When it is null, this corresponds to the Marchenko-Pastur setup for covariance matrices,
and our formulas simplify into Cauchy-like functions.

1 Introduction
Suppose A is a deterministic M × N matrix with M ≥ N and Bt has the same dimensions and
contains independent Brownian motions. The matrix
1
Xt := A + √ Bt (1.1)
N
can be viewed as a noisy observation of A . For m < M and n < N , we are interested in
comparing Xt with X̃t , defined by
(
Xtij , if i ≤ m and j ≤ n ,
X̃tij = (1.2)
0 , otherwise,

through their singular vectors. We focus on the case n ≤ m so that X̃t has rank n almost surely.
 11 
Xt · · · Xt1n 0 · · · 0
 .. .. .. .. 
 . . . .
 
 .. .. .. .. . . .. 
 . . . . . .
 .. .. .. .. 
 
X̃t = 
 . . . . .
Xtm1 · · · Xtmn 0 · · · 0

 
 0 ··· 0 0 · · · 0
 .. .. .. . . .. 
 
..
 . . . . . .
0 ··· 0 0 ··· 0

∗
Ecole Polytechnique, CMAP, elie.attal@polytechnique.edu
†
Qube Research & Technologies, romain.allez@qube-rt.com

1
Let us introduce:

• X T the Singular Values Decomposition (SVD) of X , with singular values

pt = Ut Λt Vt p t
t
λ1 ≥ ... ≥ λtN > 0 , left singular vectors ut1 , ... , utM and right singular vectors
v1t , ... , vN
t .

p p q
• X̃t = Ũt Mt ṼtT the SVD of X̃t , with singular values µt1 ≥ ... ≥ µtn > 0 = µtn+1 =
p
... = µtN , left singular vectors ũt1 , ... , ũtM and right singular vectors ṽ1t , ... , ṽN
t . We add
T
the following condition: the null space of X̃t can be divided into two parts. The singular
vectors ũtn+1 , ... , ũtm have all their M − m last components equal to zero, representing
the fact that m ≥ n so that the first n columns do not form a free family of vectors.
Furthermore, the vectors ũtm+1 , ... , ũtM have all their first m components equal to zero,
representing the part of the null space due to the shape of X̃t and its M − m null columns.
Specifically, the latter can be seen as em+1 , ... , eM (where ei has all his coefficients null
except the i-th which equals 1). Note that this condition corresponds to taking certain
linear combinations of the vectors of the null space, and therefore does not modify the
formulas obtained for the singular vectors associated with non-zero singular values.

We are interested in the limiting behaviour of the overlaps ⟨ũti |utj ⟩ and ⟨ṽit |vjt ⟩ for any t , as
M , N , m , n → ∞ with N / M → q , as well as n / N → α and m / M → β , i.e. the macroscopic
regime. Specifically, we study these limits for singular vectors in the bulk of both spectra.
This is equivalent to studying the overlaps between the eigenvectors of the square matrices
Rt := XtT Xt , R̃t := X̃tT X̃t , Lt := Xt XtT and L̃t := X̃t X̃tT which are empirical covariance
matrices of Xt or X̃t . When A is null, one can view Xt as a dataset of M independent samples
of N independent Gaussian variables of mean zero and variance t , and X̃t is a subselection of
a macroscopic number of samples and features. Our work allows one to compare the Principal
Component Analysis (PCA) of Xt with X̃t ’s eigenvector by eigenvector, under the assumption
of independent features, which corresponds to the Marchenko-Pastur setup. Note that similarly
to [4], the time t is the variance of the noise added to A , but it is also a way to derive dynamics
that allow us to obtain our results.
As mentioned in [4], there is no trivial deterministic relation between the eigenvectors of
a symmetric matrix and those of one of its principal minors. In that context, the Random
Matrix Theory approach has proved to be a powerful tool allowing to obtain explicit asymptotic
formulas for the expectations of the squared overlaps. The case of Wishart matrices we are
considering here is no different, we expect to obtain similar results for the overlaps of left and
right singular vectors using random matrices.
Moreover, the use of random matrices accounts for the noise measured on top of a relevant
signal. The Marchenko-Pastur distribution of singular values (see [28, 34, 33]) for perturbed
data or image such as (1.1) has been widely used for denoising in many different contexts,
including MRI images [38, 37, 40], financial data [7, 23, 36] and wireless communications [5, 35].
Other results that focus on the eigenvalues of Wishart matrices (squared singular values of Xt
when A ≡ 0 in our setup) such as the BBP phase transition [6] and the Tracy-Widom law for
extreme eigenvalues [22] have found applications in various fields [32, 27]. Although these results
mainly focus on the eigenvalues of such random matrices, their eigenvectors have gained interest
over the years. The main focus is to derive estimators of the population covariance matrix while
observing a sample covariance, such as in [24, 29, 9, 26]. Additionally, minors of Wishart matrices
have been increasingly studied in recent years, with applications in conditional independence in
covariance matrices [15], compressed sensing [11, 21] and percolation theory [1, 14].
In our previous work [4], we derived explicit formulas in the context of symmetric Gaussian
matrices for the limiting rescaled mean squared overlaps between the eigenvectors of a principal
submatrix and those of the full matrix. Our approach was based on analysing the eigenvec-
tor flow under the Dyson Brownian motion and deriving the dynamics of a specific resolvent.

2
However, these findings were confined to symmetric matrices and their principal minors. In the
present article, we extend this method to the singular vectors of rectangular Gaussian matrices,
or equivalently, to the eigenvectors of Wishart matrices. By examining the singular vectors’ dy-
namics in this context, we establish analogous results for the limiting overlaps in the macroscopic
regime.
Our work therefore reaches two main domains of application. On the one hand, we study the
information contained in a subimage of a rectangular noisy image through their singular vectors.
On the other hand, we establish a link between the Principal Component Analysis (PCA) of
a sample covariance matrix with identity population covariance, and the PCA obtained when
removing a macroscopic number of features or samples. In particular, we believe our results
can bring new insights into Incremental PCA algorithms [20, 3, 39], PCA with missing data
[30, 18], Risk Management or Portfolio Optimization by financial sector [13] or time-dependent
PCA methods [25, 12].
In Section 2, we introduce the dynamics of the eigenvalues and the eigenvectors and derive
the correlation structure of the different Brownian motions in presence. We then recall some
results on the Stieltjes transforms of the spectral densities and their limiting Burgers equations.
The special case A ≡ 0 is shown to be that of the Marchenko-Pastur distribution. Section 3
contains the resolution of our problem. We introduce the quantities we want to study, and define
three resolvents that have forms similar to the one used in [4]. We prove that in the scaling limit,
they become solutions of a deterministic system of coupled differential equations (3.1) that we
are able to solve explicitly. Using an inversion formula, we obtain explicit forms for the limiting
rescaled mean squared overlaps, for a general matrix A (see (3.2)). In the case A ≡ 0 , we have
the following limits (3.3) for λ , µ > 0 ,
h
2
i (1 − α) t µ̄ + α (1 − β) t λ̄ + (1 − αβ) (α + 1q ) t2
N E ⟨ṽitN |vjtN ⟩ −→ q ,
(1 − αβ)2 t2 + q (λ̄ − µ̄) (αβ λ̄ − µ̄)

h
2
i (1 − β) t µ̄ + β (1 − α) t λ̄ + (1 − αβ) (1 + βq ) t2
N E ⟨ũtiN |utjN ⟩ −→ q ,
(1 − αβ)2 t2 + q (λ̄ − µ̄) (αβ λ̄ − µ̄)

√
(1 − αβ) t λ µ
⟨ṽitN |vjtN ⟩ ⟨ũtiN |utjN ⟩
 
NE −→ q ,
(1 − αβ)2 t2 + q (λ̄ − µ̄) (αβ λ̄ − µ̄)
N n m
, N , M → (q , α , β) as well as µtiN → µ , λtjN → λ and using the


as M , N , m , n → ∞ with M
   
notations µ̄ := µ − α + βq t and λ̄ := λ − 1 + 1q t .

2 Eigenvalue and Eigenvector Dynamics

In 1989, Bru ([8]) derived the dynamics of the eigenvalues and eigenvectors of Rt = XtT Xt . For
any 1 ≤ j ≤ N , we have
N
2 q t M 1 X λtj + λtk
dλtj = √ λj dbj (t) + dt + dt , (2.1)
N N N λt − λtk
k=1 j
k̸=j
q q
N N λt dw (t) + λtk dwkj (t)
1 X λtj + λtk t 1 X j jk
dvjt = − t t vj dt + √ vkt , (2.2)
2N (λj − λk )2
N λtj − λtk
k=1 k=1
k̸=j k̸=j

where {bj | 1 ≤ j ≤ N } and {wjk | 1 ≤ j ≤ M , 1 ≤ k ≤ N , j ̸= k} are two families of independent

Brownian motions, independent of each other. Specifically, in the proof of these dynamics given

3
by Bru, we can identify these processes as dbj (t) = ⟨utj |dBt vjt ⟩ and dwjk (t) = ⟨utj |dBt vkt ⟩ . These
dynamics are different from those obtained in the symmetric Brownian case, i.e. the Dyson
Brownian motion [17, 33, 34, 19], but we can find some similarities. First, we remark that the
eigenvalues are still subject to a repulsion force, which is not exactly inversely proportional to
their distance. Moreover, since the family of Brownian motions dw is independent of db , the
eigenvectors’ dynamics can be seen as diffusion processes in a random environment, given by the
eigenvalues trajectories, which is also the case for the Dyson Brownian motion. Note that this
is due to the fact that the Brownian motions in Bt are uncorrelated, otherwise, the coefficients
of the population covariance matrix in the v t and ut bases would appear in both dynamics.
By replacing Xt with XtT , we can obtain the dynamics of the left singular vectors utj for
1 ≤ j ≤ M . They are distinguished into two cases depending on whether j ≤ N or whether
j > N (corresponding to a vector in the null space of XtT ). For 1 ≤ j ≤ N , we have
N
1 X λtj + λtk N −M t
dutj = − ut dt + u dt
2N (λtj − λtk )2 j 2N λtj j
k=1
k̸=j

q q
1
N
X λtj dwkj (t) + λtk dwjk (t) 1
M
X dwkj (t) t
+√ utk +√ uk ,
λtj λtk
q
N − N λtj
k=1 k=N +1
k̸=j

whereas for N + 1 ≤ j ≤ M ,
N N
1 X 1 t 1 X dwjk (t) t
dutj = − u dt + √ uk .
λtk j
q
2N N λ t
k=1 k=1 k

Note that the roles of dwjk and dwkj are exchanged for the left singular vectors. Additionally,
these dynamics are identical to
q q
M t + λt M λ t dw (t) + λtk dwjk (t)
t 1 X λ j k t 1 X j kj
duj = − u dt + √ utk ,
2N (λtj − λtk )2 j N k=1 λt − λt
j k
k=1
k̸=j k̸=j

for any 1 ≤ j ≤ M , if we set λtN +1 = ... = λtM = 0 and using the convention 0 / 0 = 0 . This
form will be used throughout this paper to simplify our computations. Obviously, the notation
wjk with k > N is not properly defined, but with our convention it is always multiplied by a
null factor.
The truncated matrix R̃t = X̃tT X̃t has null coefficients outside of its n×n top left submatrix,
and has rank n almost surely. Therefore, its eigenvectors associated with non-zero eigenvalues
only have non-zero coefficients on their first n components, and inversely, its eigenvectors in the
null space have all their first n components equal to zero. Consequently, there is no interaction
with the null space and we can deal with ṽ1t , ... , ṽnt only. They behave
√ the same way the vit do
if we replace N with n and M with m (the scaling remains in 1 / N ). Similarly, the dynamics
of the µti can be derived from those of the λtj , meaning we have for any 1 ≤ i ≤ n ,
n
2 m 1 X µti + µtl
q
dµti =√ t
µi db̃i (t) + dt + dt ,
N N N µti − µtl
l=1
l̸=i
p q
n
1 X µti + µtl t 1 X
n µti dw̃il (t) + µtl dw̃li (t)
t
dṽi = − ṽ dt + √ ṽlt ,
2N (µti − µtl )2 i N l=1 µti − µtl
l=1
l̸=i l̸=i

4
 n o n o
where db̃i := ⟨ũti |dB̃t ṽit ⟩ | 1 ≤ i ≤ n and dw̃il := ⟨ũti |dB̃t ṽlt ⟩ | 1 ≤ i ≤ m , 1 ≤ l ≤ n , i ̸= l
are independent of each other. Here we defined by B̃t the truncated version of the Brownian
matrix Bt , the way we defined X̃t from Xt . For ũti , using the convention µtn+1 = ... = µtm = 0
and 0 / 0 = 0 , we get for any 1 ≤ i ≤ m ,
p q
m t t m µ t dw̃ (t) + µtl dw̃il (t)
t 1 X µi + µl t 1 X i li
dũi = − ũ dt + √ ũtl .
2N (µti − µtl )2 i N l=1 µ t − µt
i l
l=1
l̸=i l̸=i

Finally, if we want to study the overlaps between the singular vectors Xt and those of X̃t ,
we need to compute the correlations between the different Brownian motions. In Appendix A,
we prove that
⟨utj |dBt vkt ⟩ ⟨ũti |dB̃t ṽlt ⟩ = ⟨ũti |utj ⟩ ⟨ṽlt |vkt ⟩ dt ,
for any (i, l, j, k) ∈ {1 ; ... ; m} × {1 ; ... ; n} × {1 ; ... ; M } × {1 ; ... ; N } . Thus, when the following
correlations are properly defined, we have



 dwjk (t) db̃i (t) = ⟨ũti |utj ⟩ ⟨ṽit |vkt ⟩ dt ,
db (t) dw̃ (t) = ⟨ũt |ut ⟩ ⟨ṽ t |v t ⟩ dt ,

j il i j l j


 dbj (t) db̃i (t) = ⟨ũi |uj ⟩ ⟨ṽi |vjt ⟩ dt ,
t t t
dw (t) dw̃ (t) = ⟨ũt |ut ⟩ ⟨ṽ t |v t ⟩ dt .

jk il i j l k

Since our work focuses on eigenvectors in the bulk, we need to make the following assumption:
the spectrum of AT A , i.e. λ01 ≥ ... ≥ λ0N (recall that X0 = A), has an empirical distribution
converging to a continuous density ρ(·, 0) :
N
1 X
δλ0 (dλ) −→ ρ(λ, 0) dλ .
N j
j=1

Similarly, we assume
n
1 X
δµ0 (dλ) −→ ρ̃(µ, 0) dµ ,
n i
i=1

where ρ̃(·, 0) is also continuous. For any time t , we denote the (continuous) limiting density of
Rt ’s spectrum (respectively of the non-zero part of R̃t ’s spectrum) by ρ(·, t) (respectively ρ̃(·, t)).
We can therefore define the Stieltjes transforms associated with both spectra,
N n
1 X 1 1X 1
GN (z, t) := and G̃N (z̃, t) := ,
N z − λtj n z̃ − µti
j=1 i=1

and write their respective limits as N → ∞ as

Z Z
ρ(λ, t) ρ̃(µ, t)
G(z, t) := dλ and G̃(z̃, t) := dµ .
R z−λ R z̃ − µ

These functions, defined for z , z̃ ∈ C \ R , are classical tools used to study the limiting behaviour
of the spectral densities. Indeed, one can recover ρ from G using the Sokhotski-Plemelj formula

lim G(λ ± i ε , t) = v(λ, t) ∓ iπ ρ(λ, t) , (2.3)

ε→0+

ρ(λ′ ,t)
dλ′ is the Hilbert transform of ρ and P.V. denotes Cauchy’s
R
where v(λ, t) := P.V. R λ−λ′
principal value.

5
 Applying Itô’s lemma, we find, in the scaling limit, the following Burgers equation (see
Appendix B.1),  
1
∂t G = 1 − − 2z G ∂z G − G2 , (2.4)
q
which was originally found in [16] . Notice that this limiting differential equation is deterministic,
which confirms the intuition that the spectral density becomes deterministic in the scaling limit
and the eigenvalues stick to their quantiles. In the context of symmetric Gaussian matrices of
[4], we also find a Burgers equation, however it does not contain the additive G2 term. Using
the method of characteristics (see Appendix B.2), equation (2.4) can be solved, leading to an
implicit equation on G in the general case:
G (zt zt′ , 0)
G(z, t) = , (2.5)
1 + t G (zt zt′ , 0)

where zt := 1 − t G(z, t) and zt′ := z (1 − t G(z, t)) − q −1 − 1 t . In the case A ≡ 0 , we have

G(z, 0) = 1 / z and the equation gives G(z, t) as a zero of a second-order polynomial. We can
find the correct root due to the fact that G(z, t) ∼ 1 / z as |z| → ∞ . We obtain
q
z − ( 1q − 1) t − (z − (1 + √1q )2 t) (z − (1 − √1q )2 t)
G(z, t) = .
2zt
It corresponds to the Stieltjes transform of the Marchenko-Pastur distribution
p
(λ+ − λ)(λ − λ− )
ρ(λ, t) =
2πλt
√ 2
with λ± = (1 ± 1 / q) t , see [28] .
Similarly, we find for G̃ the limiting equation
 
β
∂t G̃ = α − − 2αz̃ G̃ ∂z̃ G̃ − α G̃2 , (2.6)
q
leading to the implicit equation

G̃ (z̃t z̃t′ , 0)
G̃(z, t) = , (2.7)
1 + αt G̃ (z̃t z̃t′ , 0)
 
where z̃t := 1 − αt G̃ (z̃, t) and z̃t′ := z̃ 1 − αt G̃ (z̃, t) − (β / q − α) t . When A ≡ 0 , it is the
 p 2
Stieltjes transform of the Marchenko-Pastur density with µ± = 1 ± β / αq αt .
These equations will be pivotal for the remainder of our computations.

3 Limiting Behaviour of the Overlaps

3.1 The General Case
We now define the quantities under investigation. Let us introduce the notations
2
Vij (t) := ⟨ṽit |vjt ⟩ , for 1 ≤ i ≤ n and 1 ≤ j ≤ N ,

2
Uij (t) := ⟨ũti |utj ⟩ , for 1 ≤ i ≤ m and 1 ≤ j ≤ M ,

Wij (t) := ⟨ṽit |vjt ⟩ ⟨ũti |utj ⟩ , for 1 ≤ i ≤ n and 1 ≤ j ≤ N .

6
The normalisation constraints of the orthonormal bases indicate that these objects vanish as
1 / N in the bulk, so that our goal is to compute the limits of N E [Vij (t)] , N E [Uij (t)] and
N E [Wij (t)] . More precisely, if in / n → x ∈ [0 , 1] and jN / N → y ∈ [0 , 1] , we have
N E [Vin jN (t)] → V (x, y, t) , where the limiting overlapping function V is the object we want to
explicit. Similarly, N E [Uin jN (t)] → U (x, y, t) and N E [Win jN (t)] → W (x, y, t) . For the left
singular vectors, we have three other cases:

• If n + 1 ≤ in ≤ m and jN / N → y ∈ [0 , 1] , then N E [Uin jN (t)] → U (1) (y, t) which does

not depend on i because the roles of utn+1 , ... , utm can be exchanged.

• If in / n → x ∈ [0 , 1] and N + 1 ≤ jN ≤ M , then N E [Uin jN (t)] → U (2) (x, t) .

• If n + 1 ≤ in ≤ m and N + 1 ≤ jN ≤ M , then N E [Uin jN (t)] → U (3) (t) .

We can define the quantile functions λ(·, t) and µ(·, t) of the limiting spectral densities at time
t as Z ∞ Z ∞
x= ρ(λ, t) dλ = ρ̃(µ, t) dµ .
λ(x,t) µ(x,t)

They allow us to define more suitable target functions using a change of variable. We define V̄ ,
Ū and W̄ with V̄ (µ(x, t) , λ(y, t) , t) = V (x , y , t) . The function Ū can be extended to the three
other cases with:

• Ū (0 , λ(y, t) , t) = U (1) (y, t) ,

• Ū (µ(x, t) , 0 , t) = U (2) (x, t) ,

• Ū (0 , 0 , t) = U (3) (t) .

Note that in most applications, we are only interested in the overlaps of singular vectors as-
sociated with non-zero singular values. Moreover, numerical simulations of overlaps involving
singular vectors of the null space can vary depending on the chosen vector, but their roles are
theoretically exchangeable in the scaling limit. We include these cases in our study (only for
the left singular vectors, since we treat the case M ≥ N and m ≥ n) because contrary to the
symmetric case of [4], the vectors of the null space appear in the dynamics of Section 2, and are
therefore needed to achieve the calculations.
Similarly to [4] , the dynamics of Vij , Uij and Wij (see Appendix C) are difficult to deal
with directly. In order to find an explicit expression for their limits, we introduce three complex
functions of the variables z , z̃ ∈ C \ R ,
n N
(N ) 1 XX Vij (t)
SV (z, z̃, t) := ,
N (z̃ − µti )(z − λtj )
i=1 j=1

m M
(N ) 1 XX Uij (t)
SU (z, z̃, t) := ,
N (z̃ − µti )(z − λtj )
i=1 j=1
q
1
n X
N µti λtj Wij (t)
(N )
X
SW (z, z̃, t) := .
N (z̃ − µti )(z − λtj )
i=1 j=1

They have self-averaging properties as the sum of many different random variables, and still
encode all the information of the squared overlaps as Stieltjes transforms. They play a role
similar to the resolvents used in [4], [10], [24] and [31] . We typically expect these quantities to
converge to deterministic integrals involving the goal functions V̄ , Ū and W̄ .

7
 This intuition is confirmed in Appendix D, as we show that applying Itô’s lemma to all three
resolvent gives a deterministic system of coupled partial differential equations in the scaling
limit:
  


∂ S
t V = g(z, t) ∂ S
z V + g̃(z̃, t) ∂ S
z̃ V + 2 SW − G(z, t) − α G̃(z̃, t) SV




  
 1 β
−1 −α
q q
∂t SU = g(z, t) ∂z SU + g̃(z̃, t) ∂z̃ SU + 2 SW − z − G(z, t) − z̃ − α G̃(z̃, t) SU (3.1)









∂t SW = g(z, t) ∂z SW + g̃(z̃, t) ∂z̃ SW + SW 2 + z z̃ S S ,
V U

where g(z, t) := 1 − 1q − 2z G(z, t) and g̃(z̃, t) := α − βq − 2αz̃ G̃(z̃, t) . Since the characteristics
of these equations are the same as those of (2.4) and (2.6), we can solve them using the method
of characteristics. If we introduce the notation S 0 (t) := S(zt zt′ , z̃t z̃t′ , 0) for any of our three
Stieltjes transforms, then our solutions are given by
zt z̃t SV0 (t)

 SV (z, z̃, t) = 2
(1−t SW (t)) −zt zt′ z̃t z̃t′ SV0 (t) SU0 (t) t2
0








zt′ z̃t′ SU
0 (t)


SU (z, z̃, t) =  2

 z z̃ (1−t SW (t)) −zt zt z̃t z̃t′ SV0 (t) SU
0 ′ 0 (t) t2





S 0 (t) (1−t SW 0 (t) +z z ′ z̃ z̃ ′ S 0 (t) S 0 (t) t
) t t t t V


SW (z, z̃, t) = W U
,


0 (t) 2 −z z ′ z̃ z̃ ′ S 0 (t) S 0 (t) t2
( 1−t SW ) t t t t V U

where zt , zt′ , z̃t and z̃t′ are defined in (2.5) and (2.7). For the detailed resolution, see Appendix
E. We stress that identifying this limiting differential system and solving it is the most crucial
part of our work.
(N )
Since SV converges to a deterministic limit SV , we deduce that it is also the limit of its
mean. The eigenvalues being deterministic in the scaling limit (we expect them to stick to the
quantiles of their limiting deterministic distribution), the expectation is asymptotically taken
only on the overlaps, meaning we have
Z Z
V̄ (µ, λ, t) ρ̃(µ, t) ρ(λ, t)
SV (z, z̃, t) = α dµ dλ .
R R (z̃ − µ)(z − λ)

Therefore, we can recover V̄ from SV using the inversion formula derived in [10] and used in [4] ,
1
V̄ (µ, λ, t) = lim ℜ [SV (λ − i ε, µ + i ε, t) − SV (λ − i ε, µ − i ε, t)] ,
ε→0+ 2π 2 α ρ̃(µ, t) ρ(λ, t)

for any µ in the support of ρ̃(·, t) and λ in the support of ρ(·, t) . Similarly, we have
Z Z √
µλ W̄ (µ, λ, t) ρ̃(µ, t) ρ(λ, t)
SW (z, z̃, t) = α dµ dλ ,
R R (z̃ − µ)(z − λ)

and
1
W̄ (µ, λ, t) = lim √ ℜ [SW (λ − i ε, µ + i ε, t) − SW (λ − i ε, µ − i ε, t)] .
ε→0+ 2π 2 α µλ ρ̃(µ, t) ρ(λ, t)

The case of Ū is a bit trickier, as we need to split SU into four parts. Indeed, one has in the

8
scaling limit
β
Z Z
Ū (µ, λ, t) ρ̃(µ, t) ρ(λ, t) q −α Z Ū (0, λ, t) ρ(λ, t)
SU (z, z̃, t) = α dµ dλ + dλ
R R (z̃ − µ)(z − λ) z̃ R (z − λ)

1
q −1 Z Ū (µ, 0, t) ρ̃(µ, t) ( βq − α) ( 1q − 1)
+α dµ + Ū (0, 0, t) .
z R (z̃ − µ) z z̃

Therefore, we need to use four different inversion formulas to extract Ū in each case:

• For µ in the support of ρ̃(·, t) and λ in the support of ρ(·, t) , we use the same inversion
than for SV and SW ,
1
Ū (µ, λ, t) = lim ℜ [SU (λ − i ε, µ + i ε, t) − SU (λ − i ε, µ − i ε, t)] .
ε→0+ 2π 2 α ρ̃(µ, t) ρ(λ, t)

• For µ = 0 and λ in the support of ρ(·, t) , we use the classical Sokhotski-Plemelj formula
already introduced for the Stieltjes transforms (2.3),
1
Ū (0, λ, t) = lim ℑ [i ε SU (λ − i ε, i ε, t)] .
ε→0+ π ( βq − α) ρ(λ, t)

Taking z̃ = i ε , multipliying it with SU , and sending ε to 0 causes all other integrals to

vanish because the supports of ρ and ρ̃ are included in R∗+ .

• We use the same method for µ in the support of ρ̃(·, t) and λ = 0 ,

1
Ū (µ, 0, t) = lim ℑ [i ε SU (i ε, µ − i ε, t)] .
ε→0+ π α ( 1q − 1) ρ̃(µ, t)

• The last case is simpler, we take z = z̃ = i ε and send ε to 0 which gives

1
Ū (0, 0, t) = lim (i ε)2 SU (i ε, i ε, t) .
ε→0+ ( βq − α) ( 1q − 1)

We are now ready to state our formulas for V̄ , Ū and W̄ for a general initial condition A ,
from which one can always compute SV (·, ·, 0) , SU (·,
·, 0) and SW (·, ·, 0) . Using the notations
yt := 1 − t v(λ, t) − iπt ρ(λ, t) , yt′ := λ yt − q −1 − 1 t , ỹt := 1 − αt ṽ(µ, t) − iαπt ρ̃(µ, t) and
ỹt′ := µ ỹt −(β / q − α) t , along with xA (t) := Sx (yt yt′ , ỹt ỹt′ , 0) and x∗A (t) := Sx (yt yt′ , ỹt∗ (ỹt′ )∗ , 0)
for x ∈ {V , U , W } , we have the following explicit formulas for µ , λ > 0 (more precisely in the
respective supports of ρ̃(·, t) and ρ(·, t)):

  
y ỹ ∗V∗
V̄ (µ, λ, t) = ℜ1 t y ỹ V
t A
− t t A
,

2 2
Z (1−t WA∗ ) −yt yt′ ỹt∗ (ỹt′ )∗ VA∗ UA∗ t2 (1−t WA ) −yt yt′ ỹt ỹt′ VA UA t2








  
yt′ (ỹt′ )∗ UA∗

1 yt′ ỹt′ UA
Ū (µ, λ, t) = µ λ Z ℜ 2 − 2 , (3.2)

 (1−t WA∗ ) −yt yt′ ỹt∗ (ỹt′ )∗ VA∗ UA∗ t2 (1−t WA ) −yt yt′ ỹt ỹt′ VA UA t2




  
∗ 1−t W ∗ +y y ′ ỹ ∗ (ỹ ′ )∗ V ∗ U ∗ t
( A) WA (1−t WA )+yt yt′ ỹt ỹt′ VA UA t

 1 WA t t t t A A
W̄ (µ, λ, t) = √ ℜ − .

 2
µλ Z (1−t WA∗ ) −yt yt′ ỹt∗ (ỹt′ )∗ VA∗ UA∗ t2 (1−t WA )2 −yt yt′ ỹt ỹt′ VA UA t2


9
where Z is the normalisation 2απ 2 ρ̃(µ, t) ρ(λ, t) . For the other cases, we introduce Gt =
limε→0+ G(i ε, t) = − R ρ(λ,t)
R R ρ̃(µ,t) 1 β
λ dλ as well as G̃t = − R µ dµ . Setting c := q −1 and c̃ := q −α
we have

−yt′ UA (t) t
 
1
Ū (0, λ, t) = ℑ ,
π λ ρ(λ, t) (1 − t WA (t))2 + c̃ yt yt′ (1 − αt G̃t ) VA (t) UA (t) t3
 
where xA (t) := Sx yt yt′ , −(1 − αt G̃t ) c̃ t , 0 for x ∈ {V , U , W } ,

−ỹt′ UA (t) t
 
1
Ū (µ, 0, t) = ℑ ,
π α µ ρ̃(µ, t) (1 − t WA (t))2 + c ỹt ỹt′ (1 − t Gt ) VA (t) UA (t) t3

where xA (t) := Sx (−(1 − t Gt ) c t , ỹt ỹt′ , 0) for x ∈ {V , U , W } ,

UA (t) t2
Ū (0, 0, t) = ,
(1 − t WA (t))2 − c c̃ VA (t) UA (t) t4
 
where xA (t) := Sx −(1 − t Gt ) c t , −(1 − αt G̃t ) c̃ t , 0 for x ∈ {V , U , W } .
h i h i
2 2
Hence, we have been able to compute the exact limits of E N ⟨ṽit |vjt ⟩ , E N ⟨ũti |utj ⟩ and
h i
E N ⟨ũti |utj ⟩ ⟨ṽit |vjt ⟩ for eigenvectors in the bulk. These formulas are completely explicit given
the initial condition A .

3.2 The Marchenko-Pastur Case

In this subsection, we show that our formulas simplify when A ≡ 0 . We have already seen in this
case that the distributions ρ and ρ̃ have explicit forms as they are Marchenko-Pastur densities.
We also know their Hilbert transforms (see Appendix F). Furthermore, since all the eigenvalues
are null at t = 0 , we have
α
SV (z, z̃, 0) = ,
z z̃
β
SU (z, z̃, 0) = ,
q z z̃
SW (z, z̃, 0) = 0 .
Therefore, we obtain 
α zt z̃t
SV (z, z̃, t) =
zt zt′ z̃t z̃t′ − αβ t2



 q




zt′ z̃t′

SU (z, z̃, t) = q βz z̃
 zt zt z̃t z̃t′ − αβ
′
q
t2





αβ t

SW (z, z̃, t) = q .


′ ′ αβ 2
zt zt z̃t z̃t − t
q

These forms are explicit and we are able to apply the previous inversion formulas to them, to
obtain simplified forms for our goal functions V̄ , Ū and W̄ .

10
  q 2  q 2   2  2 
β
We get for µ ∈ 1− αq α t , 1 + αβq α t and λ ∈ 1 − √1
q t, 1 + √1
q t ,

(1−α) t µ̄+α (1−β) t λ̄+(1−αβ) (α+ 1q ) t2


V̄ (µ, λ, t) = q

 (1−αβ)2 t2 +q (λ̄−µ̄) (αβ λ̄−µ̄)







(1−β) t µ̄+β (1−α) t λ̄+(1−αβ) (1+ βq ) t2
Ū (µ, λ, t) = q (3.3)

 (1−αβ)2 t2 +q (λ̄−µ̄) (αβ λ̄−µ̄)





 √

W̄ (µ, λ, t) = q (1−αβ) t λ µ
2 2
(1−αβ) t +q (λ̄−µ̄) (αβ λ̄−µ̄)
,
   
where λ̄ := λ − 1 + 1q t and µ̄ := µ − α + βq t . This is the most important result of our
paper. The calculations leading to these simplifications can be found in Appendix F. We note
that these are Cauchy-like functions in λ or µ , as observed in the Wigner setup of [4], as well
as in [2] and [31]. We made our computations in the case M ≥ N and m ≥ n , but these three
expressions are still valid in any other case. Moreover, they are not affected by a specific choice
of bases for the null spaces as they correspond to limiting overlaps between singular vectors
associated with non-zero singular values.
Figure 3.2 shows a comparison of these formulas with simulated rescaled mean squared
overlaps. The fit is excellent.
The other cases for Ū are also simplified into

(1−α) t


 Ū (0, λ, t) = α λ+(1−α) ( 1q −α) t






(1−β) t
Ū (µ, 0, t) = µ+(1−β) ( 1q −α) t
(3.4)







Ū (0, 0, t) = q ,

1−α q

q q
using limε→0+ G(i ε, t) = − (1−q) t and limε→0+ G̃(i ε, t) = − (β−α q) t in the Marchenko-Pastur
setup. The forms (3.4) are specific to our choice of structure for the null spaces made in the
introduction, and to the situation M ≥ N , m ≥ n . Note that numerically there can be some
differences with the overlaps obtained with simulation for certain choices of the parameters q , α
and β , due to the finite matrix size and unexchangeability of the singular vectors.

5
V 1.75 V V
2.0 U U U
W 1.50 W 4 W
1.5 1.25
3
1.00
1.0 0.75 2
0.50 1
0.5
0.25
0.0 0.00 0
0 2 4 6 8 10 12 0 2 4 6 8 10 12 0 2 4 6 8 10 12

Figure 1: Comparison of our formulas for V̄ , Ū and W̄ with numerical simulations of N E [Vij (t)]
(red plain curve for theory and red circles for data), N E [Uij (t)] (blue plain curve for theory
and blue triangles for data) and N E [Wij (t)] (green plain curve for theory and green squares
for data) for M = 300 , q = 0.9 , α = 0.4 , β = 0.8 and t = 3 as a function of λ for a fixed
µ = µ(x, t) . Left: x = 0.9 . Middle: x = 0.5 . Right: x = 0.1 .

11
Appendices
A Correlation
Let 1 ≤ i ≤ m , 1 ≤ l ≤ n , 1 ≤ j ≤ M and 1 ≤ k ≤ N , we have:
M X
N
! m n !
X XX
⟨utj |dBt vkt ⟩ ⟨ũti |dB̃t ṽlt ⟩ = utjr vks
t
dBtrs ũtir ṽls
t
dBtrs
r=1 s=1 r=1 s=1

m X
X n
= utjr ũtir vks
t t
ṽls dt .
r=1 s=1

t = 0 if s > n and for 1 ≤ i ≤ m , ũt = 0 if r > m . Thus we

We recall that for 1 ≤ l ≤ n , ṽls ir
indeed have:
⟨utj |dBt vkt ⟩ ⟨ũti |dB̃t ṽlt ⟩ = ⟨ũti |utj ⟩ ⟨ṽlt |vkt ⟩ dt .

B Burgers Equation
B.1 Deriving the Equation
Applying Itô’s lemma gives:
q
2
N
X λtj dbj (t) M X
N
1
dGN (z, t) = √ + dt
N N (z − λtj )2 N 2 (z − λtj )2
j=1 j=1

N N
1 X λtj + λtk 4 X λtj
+ dt + dt .
N2 (λtj − λtk )(z − λtj )2 N2 (z − λtj )3
j,k=1 j=1
k̸=j

The first and last sum go to 0 in the scaling limit, and the second one converges to − ∂z G(z,t)
q dt .
We need to perform some manipulations to deal with the third sum, that we denote by Σ dt .
We first split it into Σ / 2 + Σ / 2 and invert the indices in the second term. Regrouping the two
sums and applying the identity
1 1 (bj − bi )(2z − bj − bi )
2
− 2
= ,
(z − bj ) (z − bi ) (z − bj )2 (z − bi )2
we get
N
1 X (λtj + λtk ) (2z − λtj − λtk )
Σ=
2N 2 (z − λtj )2 (z − λtk )2
j,k=1
k̸=j

N N
1 X λtj + λtk 1 X λtj + λtk
= +
2N 2 (z − λtj )(z − λtk )2 2N 2 (z − λtj )2 (z − λtk )
j,k=1 j,k=1
k̸=j k̸=j

N
1 X λtj + λtk
= .
N2 (z − λtj )(z − λtk )2
j,k=1
k̸=j

12
We split this forms into two sums:
N N
1 X λtj 1 X λtk
Σ= t t + ,
N2 (z − λj )(z − λk )2 N 2 (z − λj )(z − λtk )2
t
j,k=1 j,k=1
k̸=j k̸=j

where the first sum equals, using λ / (z − λ) = −1 + z / (z − λ) and adding the missing diagonal
terms,
 
N N N
!
z X 1 1 X 1 1 X λtj
−1 +  − ,
N z − λtj N (z − λtk )2 N2 (z − λtj )3
j=1 k=1 j=1

which converges to
(1 − z G(z, t)) ∂z G(z, t) ,
and the second sum equals
 
N N N N
!
1 X 1 1 X 1 z X 1 1 X λtj
  − + − 2 ,
N z − λtj N z − λtk N (z − λtk )2 N (z − λtj )3
j=1 k=1 k=1 j=1

which converges to
G(z, t) (−G(z, t) − z ∂z G(z, t)) .
Finally, regrouping all the terms leads to the announced limiting equation (2.4)
 
1
∂t G(z, t) = 1 − − 2z G(z, t) ∂z G(z, t) − G2 (z, t) .
q

B.2 Solving the Equation

In order to obtain the implicit equation (2.5) satisfied by G , we use the method of characteristics.
We introduce two functions of a new variable s : z(s) and t(s) . We define Ĝ(s) := G(z(s), t(s)) ,
so that the chain rule gives us
dĜ dz dt
= ∂z G(z(s), t(s)) + ∂t G(z(s), t(s))
ds ds ds
   
dz 1 dt dt
= + 1 − − 2z(s) Ĝ ∂z G(z(s), t(s)) − Ĝ2 .
ds q ds ds
Therefore, if we choose the functions z and t such that
(
dt
ds = 1
dz 1
ds = 2z(s) Ĝ(s) + q − 1 ,

then dĜ / ds = −Ĝ2 , meaning Ĝ(s) = Ĝ(0) / (1 + s Ĝ(0)) . This simplifies the differential
equation on z which allows us to obtain
    1  
z(s) = 1 + s Ĝ(0) z(0) 1 + s Ĝ(0) + −1 s .
q
Finally, the solution Ĝ gives Ĝ(0) = Ĝ(s) / (1 − s Ĝ(s)) , i.e. for any s ,
G(z(s), t(0) + s)
G(z(0), t(0)) = .
1 − s G(z(s), t(0) + s)
Evaluating this at s = −t(0) and noticing z(0) and t(0) are free parameters, we obtain the
announced implicit equation (2.5)
G (1 − t G(z, t)) z (1 − t G(z, t)) + 1 − q −1 t , 0




G(z, t) = .
1 + t G ((1 − t G(z, t)) (z (1 − t G(z, t)) + (1 − q −1 ) t) , 0)

13
C Itô Dynamics of the Squared Overlaps
We compute here the Itô dynamics of the different squared overlaps Vij (t) , Uij (t) and Wij (t) .
We detail the calculations for Vij and state the dynamics for the other cases, as the calculations
are similar. For readability, we use the notation [ · ] by

[ailjk ] := ailjk + ailkj + alijk + alikj . (C.1)

In addition, we define oviljk := ⟨ṽit |vjt ⟩ ⟨ṽlt |vkt ⟩ , ouiljk := ⟨ũti |utj ⟩ ⟨ũtl |utk ⟩ and ow t t t t
iljk := ⟨ṽi |vj ⟩ ⟨ũl |uk ⟩ .
Let 1 ≤ i ≤ n and 1 ≤ j ≤ N , we first compute

d ⟨ṽit |vjt ⟩ = ⟨dṽit |vjt ⟩ + ⟨ṽit |dvjt ⟩ + ⟨dṽit |dvjt ⟩

p q
n
1 X µti + µtl 1 X
n µti dw̃il (t) + µtl dw̃li (t)
t t
=− ⟨ṽ |v ⟩ dt + √ ⟨ṽlt |vjt ⟩
2N (µti − µtl )2 i j N l=1 µti − µtl
l=1
l̸=i l̸=i

q q
N
1 X λtj + λtk 1 X
N λtj dwjk (t) + λtk dwkj (t)
t t
− ⟨ṽ |v ⟩ dt + √ ⟨ṽit |vkt ⟩
2N (λtj − λtk )2 i j N k=1 λtj − λtk
k=1
k̸=j k̸=j

n N
1 XX Atiljk
+ ⟨ṽ t |v t ⟩ ,
N (µti − µtl )(λtj − λtk ) l k
l=1 k=1
l̸=i k̸=j

where for any l ̸= i in {1 ; ... ; n} and any k ̸= j in {1 ; ... ; N } ,

q q  q q 
t t t t t
Ailjk := µi dw̃il (t) + µl dw̃li (t) λj dwjk (t) + λk dwkj (t)

hq i
= µti λtj ow
likj dt .

Now, we can compute the dynamics of the squared overlaps:

2
dVij (t) = 2 ⟨ṽit |vjt ⟩ d ⟨ṽit |vjt ⟩ + d ⟨ṽit |vjt ⟩

p q
1
n
X µti + µtl 2
n
X µti dw̃il (t) + µtl dw̃li (t)
=− Vij dt + √ ⟨ṽlt |vjt ⟩ ⟨ṽit |vjt ⟩
N (µti − µtl )2 N µti − µtl
l=1 l=1
l̸=i l̸=i

q q
1
N
X λtj + λtk 2
N
X λtj dwjk (t) + λtk dwkj (t)
− Vij dt + √ ⟨ṽit |vkt ⟩ ⟨ṽit |vjt ⟩
N (λtj − λtk )2 N λtj − λtk
k=1 k=1
k̸=j k̸=j

n N
2 XX Atiljk
+ ⟨ṽ t |v t ⟩ ⟨ṽit |vjt ⟩
N (µti − µtl )(λtj − λtk ) l k
l=1 k=1
l̸=i k̸=j

14
 n N
1 X µti + µtl 1 X λtj + λtk
+ Vlj dt + Vik dt
N (µti − µtl )2 N (λtj − λtk )2
l=1 k=1
l̸=i k̸=j

n N
2 XX Atiljk
+ ⟨ṽ t |v t ⟩ ⟨ṽit |vkt ⟩ ,
N (µti − µtl )(λtj − λtk ) l j
l=1 k=1
l̸=i k̸=j

which can be rewritten as

n N
1 X µti + µtl 1 X λtj + λtk
dVij (t) = (Vlj − Vij ) dt + (Vik − Vij ) dt
N (µti − µtl )2 N (λtj − λtk )2
l=1 k=1
l̸=i k̸=j

 hq t t v i
µi λj oilkj ow

2
n X
X N Vij W̄lk + likj
+ dt
N (µti − µtl )(λtj − λtk )
l=1 k=1
l̸=i k̸=j

p q
2
n
X µti dw̃il (t) + µtl dw̃li (t)
+√ oviljj
N µti − µtl
l=1
l̸=i

q q
2
N
X λtj dwjk (t) + λtk dwkj (t)
+√ oviijk ,
N λtj − λtk
k=1
k̸=j
q
where W̄ij (t) := µti λtj Wij (t) .
Similarly, one finds that
m M
1 X µti + µtl 1 X λtj + λtk
dUij (t) = (Ulj − Uij ) dt + (Uik − Uij ) dt
N (µti − µtl )2 N (λtj − λtk )2
l=1 k=1
l̸=i k̸=j

  hq t t u w
i
2
m X
X M U ij W̄ lk + µ λ o o
i j ilkj iljk
+ dt
N (µi − µl )(λtj − λtk )
t t
l=1 k=1
l̸=i k̸=j

p q
2
m
X µti dw̃li (t) + µtl dw̃il (t)
+√ ouiljj
N µti − µtl
l=1
l̸=i

q q
2
M
X λtj dwkj (t) + λtk dwjk (t)
+√ ouiijk .
N λtj − λtk
k=1
k̸=j

15
 Finally, for Wij , the form is quite heavy as we mix sums with indices ending at four different
bounds: n , m , N and M . We find
q
µti µtl Wlj − µti + µtl Wij


n 2
1 X n−m
dWij (t) = t t dt + Wij dt
N (µi − µl ) 2 2N µti
l=1
l̸=i

q  
1
N
X 2 λtj λtk Wik − λtj + λtk Wij N −M
+ + Wij dt
N (λtj − λtk )2 2N λtj
k=1
k̸=j

q q
1
n X
X N µti λtj [Vij Ulk ] + µtl λtk [Wij Wlk ]
+ dt
N (µti − µtl )(λtj − λtk )
l=1 k=1
l̸=i k̸=j

q h i q h i
n
1 XX
N µti λtk ovlljk ouiijk + µtl λtj ovilkk ouiljj
+ dt
N (µti − µtl )(λtj − λtk )
l=1 k=1
l̸=i k̸=j

q q
1
m
X N
X λtj (Vij Ulk + Vik Ulj ) + 2 λtk oviijk oulljk
+ p dt
N
l=n+1 k=1
µti (λtj − λtk )
k̸=j

p q
n M
1 X X µti (Vij Ulk + Vlj Uik ) + 2 µtl oviljj ouiijk
+ q dt
N (µti − µtl ) λtj
l=1 k=N +1
l̸=i

p q
1
m
X M
X 1
n
X µti dw̃il (t) + µtl dw̃li (t)
+ Vij Ulk dt + √ ow
lijj
µti − µtl
q
N µti λtj i=n+1 j=N +1
N l=1
l̸=i

p q
1
n
X µti dw̃li (t) + µtl dw̃il (t) 1 Xm
+√ ow
iljj +p dw̃li (t) ow
iljj
N µti − µtl N µti l=n+1
l=1
l̸=i

q q
1
N
X λtj dwjk (t) + λtk dwkj (t)
+√ ow
iikj
N λtj − λtk
k=1
k̸=j

q q
1 X
N λtj dwkj (t) + λtk dwjk (t) 1 XM
+√ ow
iijk + √ dwkj (t) ow
iijk .
N k=1 λtj − λtk N λtj k=N +1
k̸=j

16
D System of Partial Differential Equations on the Double Stieltjes Trans-
forms
D.1 First Properties
Here, we define certain tools that will be our main manipulations to derive the system of par-
tial differential equations. We make use of the notations introduced in Appendix C. We first
introduce four symmetrisation properties on sums:
X akp 1 X akp
If akp = apk , = , (S1)
(bk − bp )(z − bk ) 2 (z − bk )(z − bp )
k,p k,p
p̸=k p̸=k
X X
akp + apk = 2 akp , (S2)
k,p k,p
p̸=k p̸=k
X akp + apk X akp
= , (S3)
(bk − bp )(z − bk ) (z − bk )(z − bp )
k,p k,p
p̸=k p̸=k
XX [ailjk ] XX ailjk
t t t t t t = . (S4)
(µi − µl )(z̃ − µi )(λj − λk )(z − λj ) (z̃ − µi )(z̃ − µtl )(z − λtj )(z − λtk )
t
i,l j,k i,l j,k
l̸=i k̸=j l̸=i k̸=j

These properties can be easily proved:

• For (S1), we separate the left sum S into S / 2 + S / 2 and invert the indices in the second
term. Then, we apply the identity
1 1 1
− = . (I)
(bk − bp )(z − bk ) (bk − bp )(z − bp ) (z − bk )(z − bp )

• (S2) is easily obtained by expanding into two sums and inverting the indices in the second
one.

• (S3) is an application of the two previous properties. Indeed, akp + apk is symmetric so we
can use (S1) and obtain
1 X akp + apk
.
2 (z − bk )(z − bp )
k,p
p̸=k

Symmetrisation (S2) then gives the final result.

• Symmetrisation (S4) is an application of (S3) to each double sum separately, i.e. to indices
i and l and then to j and k .

Finally, we prove a reduction property that exploits the specific structure of a certain type
of sum that we will encounter several times in our computation. It shows that despite the fact
sum appears to be of order O(1) given the order of magnitude of the overlaps in the
that this √
bulk (1 / N ), it is in fact going to zero in the scaling limit at least as 1 / N :

cij ⟨ṽit |vkt ⟩ ⟨ṽlt |vjt ⟩ ⟨ṽlt |vkt ⟩

 
1 X X 1
√ =O , (R)
N2 N (z̃ − µti )p1 (z̃ − µtl )p2 (z − λtj )p3 (z − λtk )p4 N
(i,l)∈I˜ (j,k)∈I

for any p1 , p2 , p3 , p4 ≥ 0 and z , z̃ ∈ C \ R and cij = O(1) . With I and I˜ both in {1 ; ... ; N }2
such that the summands are well defined.

17
Proof. We introduce the notations

cij X X ⟨ṽit |vkt ⟩ ⟨ṽlt |vjt ⟩ ⟨ṽlt |vkt ⟩

aij := , bij := ,
t
(z̃ − µi ) 1 (z − λtj )p3
p (z̃ − µtl )p2 (z − λtk )p4
l k
(i,l)∈I˜ (j,k)∈I

that are considered null if the indices (i, j) do not allow the correct definition of the terms. We
have
1 X
Σ= √ aij bij .
N 2 N i,j
Using the Cauchy-Schwarz inequality we get
  
1 X X
|Σ|2 ≤ 5  |aij |2   |bij |2  .
N
i,j i,j

Let us treat both sums separately. First,

X X |cij |2
|aij |2 = = O(N 2 ) ,
|z̃ − µti |2p1 |z − λtj |2p3
i,j i,j

and secondly,

X N X
X N
|bij |2 ≤ |bij |2
i,j i=1 j=1

N X
N
X X X ⟨ṽit |vkt ⟩ ⟨ṽlt |vjt ⟩ ⟨ṽlt |vkt ⟩ ⟨ṽit |vkt ′ ⟩ ⟨ṽlt′ |vjt ⟩ ⟨ṽlt′ |vkt ′ ⟩
≤ .
(z̃ − µtl )p2 (z̃ ∗ − µtl′ )p2 (z − λtk )p4 (z ∗ − λtk′ )p4
i=1 j=1 l,l′ k,k′
(i,l)∈I˜ (j,k)∈I
′
(i,l′ )∈I˜ (j,k )∈I

Since v1t , ... , vN

t is an orthonormal basis of RN , we have

N
X
⟨ṽlt |vjt ⟩ ⟨ṽlt′ |vjt ⟩ = ⟨ṽlt |ṽlt′ ⟩ = δll′ ,
j=1

and similarly
N
X
⟨ṽit |vkt ⟩ ⟨ṽit |vkt ′ ⟩ = δkk′ ,
i=1
so that
X X Vlk (t)
|bij |2 ≤ t = O(N ) .
|z̃ − 2p
µl | 2 |z − λtk |2p4
i,j l,k

Therefore, we get  
2 1
|Σ| = O ,
N2
which means Σ = O(1 / N ) .

Note that this property is also satisfied if we replace the overlaps ⟨ṽ|v⟩ by ⟨ũ|u⟩ , working
with indices in {1 ; ... ; M } . Similarly, we can replace ⟨ṽlt |vjt ⟩ ⟨ṽit |vkt ⟩ with ⟨ũtl |utj ⟩ ⟨ṽit |vkt ⟩ or
⟨ṽlt |vjt ⟩ ⟨ũti |utk ⟩ for example.

18
D.2 Deriving the System
The system of partial differential equations is obtained by applying Itô’s lemma to each of the
three functions. Therefore, we detail how we obtain the equation on SV (the method for SU is
almost identical) and the equation on SW . We work with fixed t and fixed z , z̃ ∈ C \ R , all
three independent of M , N , m , n .

(N )
First Equation Itô’s formula on SV gives
n N n N
(N ) 1 XX dVij (t) 1 XX Vij (t)
dSV = + dµt
N t t
(z̃ − µi )(z − λj ) N (z̃ − µti )2 (z − λtj ) i
i=1 j=1 i=1 j=1

n N n N
1 XX Vij (t) t 1 XX dVij (t)
+ dλj + dµt
N t t
(z̃ − µi )(z − λj ) 2 N (z̃ − µti )2 (z − λtj ) i
i=1 j=1 i=1 j=1

n N n N
1 XX dVij (t) t 1 XX Vij (t)
+ dλ + dµt dλt
N (z̃ − µti )(z − λtj )2 j N (z̃ − µti )2 (z − λtj )2 i j
i=1 j=1 i=1 j=1

n N n N
1 XX Vij (t) t 2

1 XX Vij (t) t 2


+ dµi + dλj .
N (z̃ − µti )3 (z − λtj ) N (z̃ − µti )(z − λtj )3
i=1 j=1 i=1 j=1

Based on the correlations derived in Section 2, we have:

√ t N qλt ow +√λt ow
4 µ P j iikj k iijk v 1
• dVij (t) dµti = N i


λt −λt
oiijk dt = O N2
.
j k
k=1
k̸=j
q
n
√ √
4 λtj µti ow µtl ow
lijj + 1
dVij (t) dλtj oviljj dt = O
P iljj


• = N µti −µtl N2
.
l=1
l̸=i
q q
4 4 1
• dµti dλtj = µti λtj db̃i (t) dbj (t) = µti λtj Wij (t) dt = O


N N N2
.

2
• dµti = N4 µti dt = O 1


N .
 2
• dλtj = N4 λtj dt = O 1


N .

Therefore, using the fact that Vij vanishes as 1 / N in the scaling limit, we can rewrite our
previous Itô formula as
n N n N
(N ) 1 XX dVij (t) 1 XX Vij (t)
dSV = + dµt
N t t
(z̃ − µi )(z − λj ) N (z̃ − µti )2 (z − λtj ) i
i=1 j=1 i=1 j=1

n N
1 XX Vij (t)
+ dλt + o(1) .
N (z̃ − µti )(z − λtj )2 j
i=1 j=1

We denote the sums on the right-hand side respectively by dΣV , dΣµ and dΣλ . We can expand
the first sum dΣV using the dynamics of Vij for Appendix C as
dΣV = (Iµ + Iλ + Iµλ ) dt + dIw̃ + dIw ,

19
where:
n N
µti + µtl (Vlj − Vij )


1 X X
Iµ := 2 ,
N (µti − µtl )2 (z̃ − µti )(z − λtj )
i,l=1 j=1
l̸=i

 
1
n
X N
X λtj + λtk (Vik − Vij )
Iλ := ,
N2 (λtj − λtk )2 (z̃ − µti )(z − λtj )
i=1 j,k=1
k̸=j

h q i
n N V W̄ + µ t λ t ov o w
2 X X ij lk i j ilkj likj
Iµλ := 2 ,
N (µti − µtl )(λtj − λtk )(z̃ − µti )(z − λtj )
i,l=1 j,k=1
l̸=i k̸=j

p q
2
n X
X N µti dw̃il (t) + µtl dw̃li (t)
dIw̃ := √ oviljj ,
N N (µti − µtl )(z̃ − µti )(z − λtj )
i,l=1 j=1
l̸=i

q q
2
n
X N
X λj dwjk (t) + λtk dwkj (t)
t
dIw := √ oviijk .
N N (λtj − λtk )(z̃ − µti )(z − λtj )
i=1 j,k=1
k̸=j

Our goal is to prove the following convergences:

 
• Iµ dt + dΣµ → α − βq − 2αz̃ G̃(z̃, t) ∂z̃ SV dt − α G̃(z̃, t) SV dt ,
 
1
• Iλ dt + dΣλ → 1 − q − 2z G(z, t) ∂z SV dt − G(z, t) SV dt ,

• Iµλ dt → 2 SW SV dt ,

• dIw̃ → 0 and dIw → 0 .

We start
by manipulating Iµ , applying symmetrisation (S3) to indices i and l with ail =
µti + µtl Vlj / (µti − µtl ) , we transform it into
n N
µti + µtl Vlj


1 X X
Iµ = 2 .
N (µti − µtl )(z̃ − µti )(z̃ − µtl )(z − λtj )
i,l=1 j=1
l̸=i

Now, using the dynamics of µti , we have

 
n X
N n
1 X Vij m
 dt + 1
X µti + µtl 
dΣµ = t t dt
N (z̃ − µi )(z − λj )  N N µti − µtl 
i=1 j=1 l=1
l̸=i

n X
N p
2 X Vij µti db̃i (t)
+ √ ,
N N (z̃ − µti )2 (z − λtj )
i=1 j=1

20
which, by inverting the indices i and l in the double sum, can be rewritten as
n N
µti + µtl Vlj


m (N ) 1 X X
dΣµ = − ∂z̃ SV dt − 2 dt + o(1) .
N N (µti − µtl )(z̃ − µtl )2 (z − λtj )
i,l=1 j=1
l̸=i

Using identity (I), we obtain

n N
µti + µtl Vlj


1 X X m (N )
Iµ dt + dΣµ = 2 t 2 t t dt − ∂z̃ SV dt + o(1) ,
N (z̃ − µl ) (z̃ − µi )(z − λj ) N
i,l=1 j=1
l̸=i

where we can add the diagonal terms l = i (that are well defined since we got rid of the µti − µtl
denominators) as their are vanishing in the scaling limit because of the factor 1 / N 2 and of the
order of magnitude of Vlj . We can expand (µti + µtl ) in the sum, the first sum we obtain is
n N n
! n X N

1 X X µti Vlj 1 X µti 1 X V lj
=  
N2 (z̃ − µtl )2 (z̃ − µti )(z − λtj ) N z̃ − µti N (z̃ − µtl )(z − λtj )
i,l=1 j=1 i=1 l=1 j=1

 n n  
(N )
= − + z̃ G̃(z̃, t) −∂z̃ SV
N N
n n 
(N )
= − z̃ G̃(z̃, t) ∂z̃ SV ,
N N
and the second one is
n N n
! n X
N

1 X X µtl Vlj 1 X 1 1 X µ tV
l lj
=  
N2 (z̃ − µtl )2 (z̃ − µti )(z − λtj ) N z̃ − µti N (z̃ − µtl )2 (z − λtj )
i,l=1 j=1 i=1 l=1 j=1

n 
(N ) (N )

= G̃(z̃, t) −SV − z̃ ∂z̃ SV .
N
Since n / N → α and m / N → β / q , we obtain the announced convergence
 
β
Iµ dt + dΣµ → α − − 2αz̃ G̃(z̃, t) ∂z̃ SV − α G̃(z̃, t) SV .
q
The method for the convergence of Iλ dt + dΣλ is identical.
We now derive the limit of Iµλ , applying symmetrisation (S4) to it we obtain
q
2 X n XN V ij W̄ lk + µti λtj ovilkj ow
likj
Iµλ = 2 .
N (z̃ − µi )(z̃ − µl )(z − λj )(z − λtk )
t t t
i,l=1 j,k=1
l̸=i k̸=j

(1) (2)
Expanding the numerator we get two sums Iµλ + Iµλ . Adding the diagonal terms l = i and
k = j to the first one, since they vanish in the scaling limit, gives
 
n X N n X N
!
(1) 1 X V ij 1 X W̄ lk
Iµλ = 2   + o(1)
N (z̃ − µti )(z − λtj ) N (z̃ − µtl )(z − λtk )
i=1 j=1 l=1 k=1

= 2 SW SV + o(1) .

21
Adding to the second sum its diagonal terms that are of order 1 / N , we have
q
n
2 X X
N µti λtj ⟨ṽit |vkt ⟩ ⟨ṽlt |vjt ⟩ ⟨ṽlt |vkt ⟩ ⟨ũti |utj ⟩  
1
(2)
Iµλ = 2 t t t t +O
N (z̃ − µi )(z̃ − µl )(z − λj )(z − λk ) N
i,l=1 j,k=1
√ q (2)
Applying (R) with cij = N µti λtj ⟨ũti |utj ⟩ = O(1) , we get Iµλ → 0 .
Finally, we prove that the Brownian terms dIw and dIw̃ go to zero in the scaling limit. We
detail the method for dIw only. The independence of dw with respect to the other random
variables in the sum indicates that dIw is centered. Furthermore, we can apply symmetrisation
(S3) to the indices j and k which leads to
q
2 X n X N λtj dwjk (t)
dIw = √ ov ,
N N (z̃ − µti )(z − λtj )(z − λtk ) iijk
i=1 j,k=1

and we can write its variance as

 
q
λtj λtj ′ dwjk (t) dwj ′ k′ (t) oviijk ovllj ′ k′
 
n N
h i 4 X X 
E |dIw |2 = 3 E
 
N  t ∗ t t t ∗ t ∗ t
(z̃ − µi )(z̃ − µl )(z − λj )(z − λk )(z − λj ′ )(z − λk′ ) 

i,l=1 j,k,j ′ ,k′ =1 
k̸=j
k′ ̸=j ′

 
n N
4  X X λtj ⟨ṽit |vjt ⟩ ⟨ṽit |vkt ⟩ ⟨ṽlt |vjt ⟩ ⟨ṽlt |vkt ⟩ 
= E 
2
 dt ,
N3  2
i,l=1 j,k=1 (z̃ − µti )(z̃ ∗ − µtl ) z − λtj z − λtk
k̸=j
h i √
which gives E |dIw |2 = O(1 / N 2 ) using (R) with cij = N λtj ⟨ṽit |vjt ⟩ = O(1) . Since the
variances are summable with respect to N , Borel-Cantelli’s lemma indicates that dIw → 0
almost surely.
(N )
We have proved that randomness vanishes almost surely in the equation on SV and leads
to  
∂t SV = g(z, t) ∂z SV + g̃(z̃, t) ∂z̃ SV + 2 SW − G(z, t) − α G̃(z̃, t) SV ,
1 β
with g(z, t) := 1 − q − 2z G(z, t) and g̃(z̃, t) := α − q − 2αz̃ G̃(z̃, t) .

Second Equation The equation on SU is obtained with the same method. The only difference
(N )
comes from the fact that instead of obtaining the term GN (z, t) SU , we get
M
1 X 1 (N )
t SU ,
N z − λj
j=1

which is equal to (recalling that we introduced the notations λtN +1 = ... = λtM = 0 for simplicity)
 
M −N (N )
+ GN (z, t) SU ,
Nz
that converges to
1
!
q −1
+ G(z, t) SU .
z

A similar modification is obtained for the G̃(z̃, t) SU term.

22
Third Equation For the equation on SW , once summed, the Brownian terms almost surely
vanish in the scaling limit using the same argument as for SV . Likewise, the sums of the form
h i
c o v ou
1 X X ij lljk iijk
N2 (z̃ − µti )(µti − µtl )(z − λtj )(λtj − λtk )
i,l=1 j,k=1
l̸=i k̸=j

go to zero (using arguments similar to (R)). Therefore, we focus on the transformation

q of the non
vanishing terms (we recall that in Appendix C we introduced the notation W̄ij = µti λtj Wij ):
q
1
n X
N µti λtj dWij (t)
(N )
X
dSW =
N (z̃ − µti )(z − λtj )
i=1 j=1

n N
!
1 XX W̄ij (t) W̄ij (t)
+ + dµti
N 2µti (z̃ − µti )(z − λtj ) (z̃ − µti )2 (z − λtj )
i=1 j=1

n N
!
1 XX W̄ij (t) W̄ij (t)
+ + dλtj + o(1) .
N 2λtj (z̃ − µti )(z − λtj ) (z̃ − µti )(z − λtj )2
i=1 j=1

We denote by ΣW , Σµ and Σλ the three sums on the right hand side, in their respective order.
From what we said, most of the terms vanish in ΣW so we can write

ΣW = (Iµ + Iλ + IV U + IW ) dt + o(1) ,

where:
n N n P
N
1 P P 2µti W̄lj −(µti +µtl ) W̄ij n−m P W̄ij
• Iµ := N2 (µti −µtl )2 (z̃−µti )(z−λtj )
+ 2N 2 µti (z̃−µti )(z−λtj )
,
i,l=1 j=1 i=1 j=1
l̸=i

n N 2λtj W̄ik −(λtj +λtk ) W̄ij n P

N
1 P P N −M P W̄ij
• Iλ := N2 (λtj −λtk )2 (z̃−µti )(z−λtj )
+ 2N 2 λtj (z̃−µti )(z−λtj )
,
i=1 j,k=1 i=1 j=1
k̸=j
q
n N µti µtl λtj λtk [Wij Wlk ]
1 P P
• IW := N2 (µti −µtl )(z̃−µti )(λtj −λtk )(z−λtj )
,
i,l=1 j,k=1
l̸=i k̸=j

and
n N
1 X X µti λtj [Vij Ulk ]
IV U := 2
N (µti − µtl )(z̃ − µti )(λtj − λtk )(z − λtj )
i,l=1 j,k=1
l̸=i k̸=j

n m N
1 X X X λtj (Vij Ulk + Vik Ulj )
+ 2
N (z̃ − µti )(λtj − λtk )(z − λtj )
i=1 l=n+1 j,k=1
k̸=j

23
 n N M
1 X X X µti (Vij Ulk + Vlj Uik )
+ 2
N (µti − µtl )(z̃ − µti )(z − λtj )
i,l=1 j=1 k=N +1
l̸=i

n m N M
1 X X X X Vij Ulk
+ 2 .
N (z̃ − µti )(z − λtj )
i=1 l=n+1 j=1 k=N +1

We are going to prove the following convergences:

 
• Iµ dt + Σµ → α − βq − 2αz̃ G̃(z̃, t) ∂z̃ SW dt ,
 
1
• Iλ dt + Σλ → 1 − q − 2z G(z, t) ∂z SW dt ,

• IV U → z z̃ SV SU ,
2 .
• IW → SW
We begin with the convergence of Iµ dt + Σµ . We can rewrite the first sum in Iµ as
n N n N
µti W̄lj − W̄ij


1 X X 1 X X µti W̄lj − µtl W̄ij
+ .
N2 (µti − µtl )2 (z̃ − µti )(z − λtj ) N 2 (µti − µtl )2 (z̃ − µti )(z − λtj )
i,l=1 j=1 i,l=1 j=1
l̸=i l̸=i

In the first term, we can replace µti in the numerator by z̃ because the difference between the
two sums is a null sum (the summand is antisymmetric with respect to i and l). Applying
symmetrisation (S3) to both sums, we obtain
n N
1 X X (z̃ + µti ) W̄lj
.
N2 (µti − µtl )(z̃ − µti )(z − λtj )
i,l=1 j=1
l̸=i

Once again, we use identity (I) to transform it into

n N n N
1 X X (z̃ + µti ) W̄lj 1 X X (z̃ + µti ) W̄lj
+ .
N2 (z̃ − µtl )2 (µti − µtl )(z − λtj ) N 2 (z̃ − µtl )2 (z̃ − µti )(z − λtj )
i,l=1 j=1 i,l=1 j=1
l̸=i l̸=i

The second sum converges to  

α − 2αz̃ G̃(z̃, t) ∂z̃ SW ,
and we denote by A the first sum that we will combine with Σµ . We recall that
n
m 1 X µti + µtl
dµti = dt + dt + o(1) ,
N N µti − µtl
l=1
l̸=i

so that
n N n N
m XX W̄ij m XX W̄ij
Σµ = t t t dt + 2 dt
2N 2 µi (z̃ − µi )(z − λj ) N (z̃ − µi )2 (z − λtj )
t
i=1 j=1 i=1 j=1

n N
!
1 X X µti + µtl W̄ij W̄ij
+ 2 + dt + o(1) ,
N µti − µtl 2µi (z̃ − µi )(z − λj ) (z̃ − µi )2 (z − λtj )
t t t t
i,l=1 j=1
l̸=i

24
where the second sum converges to − βq ∂z̃ SW dt and the last sum, if added to A dt (after ex-
changing the indices i and l), equals
n N
n XX W̄ij
− dt .
N 2 2µi (z̃ − µti )(z − λtj )
t
i=1 j=1

Therefore,
n N
m−n X X W̄ij β
Σµ + A dt = 2 t t t dt − ∂z̃ SW dt + o(1)
2N µi (z̃ − µi )(z − λj ) q
i=1 j=1

which cancels out with the second sum in the definition of Iµ dt . Finally, we have proved that
 
β
Iµ dt + Σµ −→ α − − 2αz̃ G̃(z̃, t) ∂z̃ SW dt .
q

The demonstration for the convergence

q of Iλ dt + Σλ is identical.
For IW , we first notice that µi µl λtj λtk [Wij Wlk ] = W̄ij W̄lk . Then, applying symmetri-
t t
 

sation (S4) we get

n N
1 X X W̄ij W̄lk
IW = 2 ,
N (z̃ − µti )(z̃ − µtl )(z − λtj )(z − λtk )
i,l=1 j,k=1
l̸=i k̸=j

which converges to SW 2 .

We now focus on the remaining term IV U . Considering its first sum, one can write µti λtj =
(µti − z̃) λtj + (λtj − z) z̃ + z z̃ to see that we can replace µti λtj by z z̃ in the numerator because the
difference between the two sums are two null sums (antisymmetric with respect to i and l or to
j and k). Therefore, the first sum in IV U equals, after applying symmetrisation (S4),
n N
z z̃ X X Vij Ulk
.
N 2 t
(z̃ − µi )(z̃ − µtl )(z − λtj )(z − λtk )
i,l=1 j,k=1
l̸=i k̸=j

The same type of reasoning can be applied to the other sums composing IV U until we obtain
n N m N
(N ) 1 XX Ulk 1 X X Ulk
IV U = z z̃ SV +
N (z̃ − µtl )(z − λtk ) N z̃(z − λtk )
l=1 k=1 l=n+1 k=1

n M m M
!
1 X X Ulk 1 X X Ulk
+ + + o(1)
N (z̃ − µtl )z N z̃z
l=1 k=N +1 l=n+1 k=N +1

(N ) (N )
= z z̃ SV SU + o(1) .

Thus, IV U converges to z z̃ SV SU .
Finally, we have proven that SW satisfies the announced deterministic differential equation.

25
E Solving the System
Let z , z̃ ∈ C \ R and t ≥ 0 . We introduce a new variable s , as well as functions z(s) , z̃(s)
and t(s) such that z(0) = z , z̃(0) = z̃ and t(0) = t . Moreover, we introduce the notation
SˆV (s) := SV (z(s) , z̃(s) , t(s)) and similarly for our other functions in the equations. Denoting
by c (respectively c̃) the constant 1q − 1 (respectively βq − α), if
′
t (s) = 1

z ′ (s) = 2Ĝ(s) z(s) + c

 ′ ˆ z̃(s) + c̃ ,
z̃ (s) = 2αG̃(s)

then the chain rule gives

 ′ 
ˆ

Sˆ (s) = 2 Sˆ (s) − Ĝ(s) − α G̃(s) SˆV (s)
 V W


′ ˆ
 
SˆU (s) = 2 SˆU (s) − z(s)
c c̃
− Ĝ(s) − z̃(s) − α G̃(s) SˆV (s)

 ˆ ′
 2
SW (s) = SˆW (s) + z(s)z̃(s) SˆV (s) SˆU (s) .

Additionally, under the previous conditions on z(s) , z̃(s) and t(s) we know from equation (2.4)
and its resolution in Appendix B that
• t(s) = t + s ,
  
• z(s) = 1 + s Ĝ(0) z (1 + sĜ(0)) + cs ,

Ĝ(0)
• Ĝ(s) = .
1+s Ĝ(0)

The equation (2.6) on G̃ can give us similarly:

ˆ ˆ
  
• z̃(s) = 1 + αs G̃(0) z̃ (1 + αs G̃(0)) + c̃s ,

ˆ ˆ
G̃(0)
• G̃(s) = ˆ .
1+αs G̃(0)

Therefore, the equations on SˆV and SˆU lead, after integration, to

SˆV (0) Rs
SˆW (u) du
SˆV (s) =   e2 0 ,
ˆ

1 + s Ĝ(0) 1 + αs G̃(0)

z z̃ SˆU (0) Rs
SˆW (u) du
SˆU (s) =   e2 0 .
ˆ

z (1 + sĜ(0)) + cs z̃ (1 + αs G̃(0)) + c̃s

This leaves us with the following differential equation on SˆW :

′ 2 Rs
SˆW (u) du
SˆW (s) = SˆW (s) + z z̃ SˆV (0) SˆU (0) e4 0 .
ˆ
R · are going to solve it explicitly. For readability we introduce the notations f := SW , F :=
We
ˆ ˆ
0 f (u) du and a := z z̃ SV (0) SU (0) . With the change of variable x = F (s) , we get

df
f = f 2 + a e4x ,
dx
so that g := f 2 satisfies
dg
= 2 g + 2a e4x ,
dx

26
which gives,
g(x) = (f 2 (0) − a) e2x + a e4x .
This can be rewritten into an order 1 differential equation on F ,
dF
q
= ± (f 2 (0) − a) e2F + a e4F .
ds
We separate the variables and integrate, which leads to
q p
a + (f 2 (0) − a) e−2F (s) = f 2 (0) ± (f 2 (0) − a) s ,

and finally,
1  p 
log 1 ± 2 f 2 (0) s + (f 2 (0) − a) s2 .
F (s) = −
2
We can now differentiate to obtain
p
∓ f 2 (0) − (f 2 (0) − a) s
f (s) = p .
1 ± 2 f 2 (0) s + (f 2 (0) − a) s2
p
The condition at s = 0 gives f (0) = ∓ f 2 (0) , therefore we end up with

f (0) + (a − f 2 (0)) s
f (s) = .
1 − 2 f (0) s − (a − f 2 (0)) s2

Putting all of this together, we obtain the system

2
  
SˆW (0)+ z z̃ SˆV (0) SˆU (0)−SˆW (0) s
SˆW (s) =

 2
  


 1−2 SˆW (0) s− z z̃ SˆV (0) SˆU (0)−SˆW (0) s2
SˆV (0) 

SˆV (s) = 
ˆ

ˆ ˆ ˆ ˆ 2
 

 ( 1+s Ĝ(0)) 1+αs G̃(0) 1−2 SW (0) s− z z̃ SV (0) SU (0)−SW (0) s
2

ˆ

SˆU (s) = zz̃ SU (0)


  .
ˆ 2
 
(z (1+sĜ(0))+cs) z̃ (1+αs G̃(0))+c̃s 1−2 SˆW (0) s− z z̃ SˆV (0) SˆU (0)−SˆW (0) s2

We denote by D(s) the common denominator

  2
 −1
1 − 2 SˆW (0) s − z z̃ SˆV (0) SˆU (0) − SˆW (0) s2 .

One can solve for D using the previous system of equations, which gives

D(s) = (1 + SˆW (s) s)2 − z(s)z̃(s) SˆV (s) SˆU (s) s2 .

Thus, we can invert the system:

ˆ ˆ SˆV (s) SˆU (s)s



 SˆW (0) = SW (s) (1+SW (s) s)−z(s)z̃(s)
D(s)
ˆ
  
(1+s Ĝ(0)) 1+αs G̃(0) SˆV (s)

ˆ
SV (0) =
 D(s) 
ˆ

(z (1+s Ĝ(0))+cs) z̃ (1+αs G̃(0))+c̃s SˆU (s)


Sˆ (s) =

.
U z z̃ D(s)

Finally, since fˆ(0) = f (z , z̃ , t) , evaluating at s = −t gives the announced result.

27
F Inversion in the Marchenko-Pastur Case
We detail the case of V̄ as the other functions are obtained almost identically. First, we recall
that
lim G(λ ± i ε, t) = v(λ, t) ∓ i π ρ(λ, t)
ε→0+

where r  
(1 + √1 )2 t −λ λ − (1 − √1 )2 t
q q
ρ(λ, t) =
2πλt
and
λ − ( 1q − 1) t
v(λ, t) = .
2λt
We have a similar relation between G̃ and
r
√ q  √ q 
( α + βq )2 t − µ µ − ( α − βq )2 t
ρ̃(µ, t) == ,
2παµt

µ − ( βq − α) t
ṽ(µ, t) = .
2αµt
Therefore if we define SV± := limε→0+ SV (λ − i ε, µ ± i ε, t) , then,
 
α (A − iB) Ã ± iB̃
SV± =    ,
αβ 2
(A − iB) (λ A − ct − i λ B) Ã ± iB̃ µ Ã − c̃t ± i µ B̃ − q t

where:

• A := 1 − t v(λ, t) ,

• B := πt ρ(λ, t) ,

• Ã := 1 − αt ṽ(µ, t) ,

• B̃ := απt ρ̃(µ, t) ,
1
• c := q − 1,
β
• c̃ := q − α.

We can simplify this into

 
α (A − i B) Ã ± i B̃
SV± =   .
αβ 2
(A (λ A − ct) − λ B 2 − i B (2λ A − ct)) Ã (µ Ã − c̃t) − µ B̃ 2 ± i B̃ (2µ Ã − c̃t) − q t

µ+c̃t
This form is very practical since we remark that A = λ+ct
2λ and à = 2µ , therefore 2λ A−ct = λ
and 2µ Ã − c̃t = µ . Furthermore, rewriting B and B̃ leads to

−λ2 + 2 (1 + 1q ) t λ − c2 t2 −µ2 + 2 (α + βq ) t µ − c̃2 t2

B2 = and B̃ 2 = .
4 λ2 4 µ2

28
  
µ̄
Therefore, A (λ A−ct)−λ B 2 = λ̄2 where λ̄ := λ− 1 + 1q t and similarly à (µ Ã−c̃t)−µ B̃ 2 = 2
 
where µ̄ := µ − α + βq t . We end up with
 
α (A − i B) Ã ± i B̃ N±
SV± =    =: .
λ̄ µ̄ αβ 2 D ±
2 − iλB 2 ± i µB̃ − q t

In order to compute V̄ , we need to explicit the real part of SV+ − SV− . We have

N+ D− − N− D+
SV+ − SV− =
D+ D−

∗ D∗
(N+ D− − N− D+ ) D+ −
= ,
|D+ D− |2
t2 λ̄2
so we begin with simplifying the denominator. When needed, we use the fact that λ2 B 2 = q − 4
αβ 2 µ̄2
and µ2 B̃ 2 = q t − 4 .

2  2
µ̄2 α2 β 2
   
2 λ̄ αβ 2 λ̄
|D+ D− | = − iλB + µ2 B̃ 2 − t − iλB µ̄ + 2 t4
2 4 q 2 q

2
λ̄2 α2 β 2
   
αβ 2 αβ 2 λ̄
= − λ2 B 2 − i λ λ̄ B t − t − iλB µ̄ + 2 t4
4 q q 2 q

2
α2 β 2 4 λ̄2 t2 λ̄ µ̄ αβ 2
= t − − + t − i λ λ̄ B + i λ µ̄ B
q2 2 q 2 q

2
α2 β 2 λ̄ (αβ − 1) 2
= 2 t4 (λ̄ − µ̄) + t − i λ (λ̄ − µ̄) B
q 2 q

α2 β 2 2 (αβ − 1)2 4
 
2 λ̄ αβ − 1 2
= 2 t4 2 2
(λ̄ − µ̄) ( + λ B ) + t λ̄ (λ̄ − µ̄) + t
q 4 q q2

α2 β 2 4 t2 (αβ − 1)2 4
 
αβ − 1 2
= t (λ̄ − µ̄)( (λ̄ − µ̄) + t λ̄) + t
q2 q q q2

α2 β 2 6 (αβ − 1)2 2
 
= t (λ̄ − µ̄) (αβ λ̄ − µ̄) + t .
q3 q

Since the denominator is the same for SV , SU and SW , this final form is helpful in all three
computations.
We now focus on the numerator, and more precisely on its real part. The previous compu-
tation gives us
 
∗ ∗ αβ 2 λ̄ αβ − 1 2
D+ D− = t (λ̄ − µ̄) + t + i λ (λ̄ − µ̄) B .
q 2 q

29
Moreover, we have

N+ D− − N− D+ = α (A − i B) (Ã (D− − D+ ) + i B̃ (D− + D+ ))

 
2αβ λ̄
=− t B̃ (A − i B) λ B + i ( + α t) ,
q 2
µ̄ β
after some simplifications using µ Ã − 2 = q t . Also,
 
∗ ∗ λ − ct αβ − 1 2 t αβ − 1 2
D+ D− (A − i B) = t (λ̄ − µ̄) + t A + i (λ̄ − µ̄) B − i t B,
2qλ q q q

using λ̄2 A + λ B 2 = λ−ct λ̄ t

2qλ t and λ A − 2 = q . We can now compute the real part of the entire
numerator, which, after some simplifications, is

2α2 β 2 5
 
1
t B B̃ (1 − α) µ̄ + α(1 − β) λ̄ + (1 − αβ) (α + ) t .
q3 q

We end up with the announced formula,

1
ℜ SV+ − SV−
 
V̄ (µ, λ, t) =
2π 2 α ρ(λ, t) ρ̃(µ, t)

(1 − α) t µ̄ + α (1 − β) t λ̄ + (1 − αβ) (α + 1q ) t2
=q .
(1 − αβ)2 t2 + q (λ̄ − µ̄) (αβ λ̄ − µ̄)

References
[1] Mark Adler, Pierre Van Moerbeke, and Dong Wang. Random matrix minor processes
related to percolation theory. Random Matrices: Theory and Applications, 2(04):1350008,
2013.

[2] Romain Allez, Joël Bun, and Jean-Philippe Bouchaud. The eigenvectors of gaussian matri-
ces with an external source. arXiv preprint arXiv:1412.7108, 2014.

[3] Matej Artac, Matjaz Jogan, and Ales Leonardis. Incremental pca for on-line visual learning
and recognition. In 2002 International Conference on Pattern Recognition, volume 3, pages
781–784. IEEE, 2002.

[4] Elie Attal and Romain Allez. Interlacing eigenvectors of large gaussian matrices. Journal
of Physics A: Mathematical and Theoretical, 2024.

[5] Zhidong Bai and Jack W Silverstein. Spectral analysis of large dimensional random matrices,
volume 20. Springer, 2010.

[6] Jinho Baik, Gérard Ben Arous, and Sandrine Péché. Phase transition of the largest eigen-
value for nonnull complex sample covariance matrices. 2005.

[7] Jean-Philippe Bouchaud and Marc Potters. Financial applications of random matrix theory:
a short review. arXiv preprint arXiv:0910.1205, 2009.

[8] Marie-France Bru. Diffusions of perturbed principal component analysis. Journal of mul-
tivariate analysis, 29(1):127–136, 1989.

30
 [9] Joël Bun, Romain Allez, Jean-Philippe Bouchaud, and Marc Potters. Rotational invari-
ant estimator for general noisy matrices. IEEE Transactions on Information Theory,
62(12):7475–7490, 2016.
[10] Joël Bun, Jean-Philippe Bouchaud, and Marc Potters. Overlaps between eigenvectors of
correlated random matrices. Physical Review E, 98(5):052145, 2018.
[11] T Tony Cai, Tiefeng Jiang, and Xiaoou Li. Asymptotic analysis for extreme eigenvalues of
principal minors of random matrices. The Annals of Applied Probability, 31(6):2953–2990,
2021.
[12] Bart De Ketelaere, Mia Hubert, and Eric Schmitt. Overview of pca-based statistical process-
monitoring methods for time-dependent, high-dimensional data. Journal of Quality Tech-
nology, 47(4):318–335, 2015.
[13] Vrinda Dhingra, Amita Sharma, and Shiv K Gupta. Sectoral portfolio optimization by
judicious selection of financial ratios via pca. Optimization and Engineering, 25(3):1431–
1468, 2024.
[14] AB Dieker and J Warren. On the largest-eigenvalue process for generalized wishart random
matrices. arXiv preprint arXiv:0812.1504, 2008.
[15] Mathias Drton, Hélène Massam, and Ingram Olkin. Moments of minors of wishart matrices.
2008.
[16] T Cabanal Duvillard and A Guionnet. Large deviations upper bounds for the laws of matrix-
valued processes and non-communicative entropies. The Annals of Probability, 29(3):1205–
1261, 2001.
[17] Freeman J Dyson. A brownian-motion model for the eigenvalues of a random matrix.
Journal of Mathematical Physics, 3(6):1191–1198, 1962.
[18] Abel Folch-Fortuny, Francisco Arteaga, and Alberto Ferrer. Pca model building with miss-
ing data: New proposals and a comparative study. Chemometrics and Intelligent Laboratory
Systems, 146:77–88, 2015.
[19] David J Grabiner. Brownian motion in a weyl chamber, non-colliding particles, and random
matrices. In Annales de l’IHP Probabilités et statistiques, volume 35, pages 177–204, 1999.
[20] Peter M Hall, A David Marshall, and Ralph R Martin. Incremental eigenanalysis for
classification. In BMVC, volume 98, pages 286–295. Citeseer, 1998.
[21] Tiefeng Jiang and Yongcheng Qi. Largest eigenvalues of principal minors of deformed
gaussian orthogonal ensembles and wishart matrices. arXiv preprint arXiv:2410.15160,
2024.
[22] Iain M Johnstone. On the distribution of the largest eigenvalue in principal components
analysis. The Annals of statistics, 29(2):295–327, 2001.
[23] Laurent Laloux, Pierre Cizeau, Marc Potters, and Jean-Philippe Bouchaud. Random matrix
theory and financial correlations. International Journal of Theoretical and Applied Finance,
3(03):391–397, 2000.
[24] Olivier Ledoit and Sandrine Péché. Eigenvectors of some large sample covariance matrix
ensembles. Probability Theory and Related Fields, 151(1):233–264, 2011.
[25] Weihua Li, H Henry Yue, Sergio Valle-Cervantes, and S Joe Qin. Recursive pca for adaptive
process monitoring. Journal of process control, 10(5):471–486, 2000.

31
[26] Zeqin Lin and Guangming Pan. Eigenvector overlaps in large sample covariance matrices
and nonlinear shrinkage estimators. arXiv preprint arXiv:2404.18173, 2024.

[27] Satya N Majumdar. Extreme eigenvalues of wishart matrices: application to entangled

bipartite system. arXiv preprint arXiv:1005.4515, 2010.

[28] Vladimir Alexandrovich Marchenko and Leonid Andreevich Pastur. Distribution of eigen-
values for some sets of random matrices. Matematicheskii Sbornik, 114(4):507–536, 1967.

[29] Xavier Mestre. Improved estimation of eigenvalues and eigenvectors of covariance matrices
using their sample estimates. IEEE Transactions on Information Theory, 54(11):5113–5129,
2008.

[30] Philip RC Nelson, Paul A Taylor, and John F MacGregor. Missing data methods in pca
and pls: Score calculations with incomplete observations. Chemometrics and intelligent
laboratory systems, 35(1):45–65, 1996.

[31] Alessandro Pacco and Valentina Ros. Overlaps between eigenvectors of spiked, correlated
random matrices: From matrix principal component analysis to random gaussian land-
scapes. Physical Review E, 108(2):024145, 2023.

[32] Nick Patterson, Alkes L Price, and David Reich. Population structure and eigenanalysis.
PLoS genetics, 2(12):e190, 2006.

[33] Marc Potters and Jean-Philippe Bouchaud. A first course in random matrix theory: for
physicists, engineers and data scientists. Cambridge University Press, 2020.

[34] Terence Tao. Topics in random matrix theory, volume 132. American Mathematical Soc.,
2012.

[35] Antonia M Tulino, Sergio Verdú, et al. Random matrix theory and wireless communications.
Foundations and Trends® in Communications and Information Theory, 1(1):1–182, 2004.

[36] Akihiko Utsugi, Kazusumi Ino, and Masaki Oshikawa. Random matrix theory analysis of
cross correlations in financial markets. Physical Review E—Statistical, Nonlinear, and Soft
Matter Physics, 70(2):026110, 2004.

[37] Jelle Veraart, Els Fieremans, and Dmitry S Novikov. Diffusion mri noise mapping using
random matrix theory. Magnetic resonance in medicine, 76(5):1582–1593, 2016.

[38] Jelle Veraart, Dmitry S Novikov, Daan Christiaens, Benjamin Ades-Aron, Jan Sijbers,
and Els Fieremans. Denoising of diffusion mri using random matrix theory. Neuroimage,
142:394–406, 2016.

[39] Juyang Weng, Yilu Zhang, and Wey-Shiuan Hwang. Candid covariance-free incremental
principal component analysis. IEEE Transactions on Pattern Analysis and Machine Intel-
ligence, 25(8):1034–1040, 2003.

[40] Wei Zhu, Xiaodong Ma, Xiao-Hong Zhu, Kamil Ugurbil, Wei Chen, and Xiaoping Wu.
Denoise functional magnetic resonance imaging with random matrix theory based principal
component analysis. IEEE Transactions on Biomedical Engineering, 69(11):3377–3388,
2022.

32