# Right large deviation principle for the top eigenvalue of the sum or product
of invariant random matrices

Pierre Mergny<EMAIL_ADDRESS>LPTMS, CNRS, Univ. Paris-Sud, Université
Paris-Saclay, 91405 Orsay, France  Chair of Econophysics $\&$ Complex Systems,
Ecole Polytechnique, 91128 Palaiseau Cedex, France  Marc Potters  Capital Fund
Management, 23 rue de l’Université, 75007 Paris, France

###### Abstract

In this note we study the right large deviation of the top eigenvalue (or
singular value) of the sum or product of two random matrices $\mathbf{A}$ and
$\mathbf{B}$ as their dimensions goes to infinity. The matrices $\mathbf{A}$
and $\mathbf{B}$ are each assumed to be taken from an invariant (or bi-
invariant) ensemble with a confining potential with a possible _wall_ beyond
which no eigenvalues/singular values are allowed. The introduction of this
wall puts different models in a very generic framework. In particular, the
case where the wall is exactly at the right edge of the limiting spectral
density is equivalent, from the point of view of large deviations, to
considering a fixed diagonal matrices, as studied previously in Ref. [1]. We
show that that the tilting method introduced in Ref. [1] can be extended to
our general setting and is equivalent to the study of a spherical spin glass
model specific to the operation - sum of symmetric matrices / product of
symmetric matrices / sum of rectangular matrices - we are considering.

###### Contents

  1. 1 Introduction
  2. 2 Reminder on right large deviation for one random matrix: The pulled Coulomb gas approach
    1. 2.1 Classical Rotationally Invariant Ensemble
    2. 2.2 Invariant Ensemble with a wall
    3. 2.3 Coulomb gas approach
      1. 2.3.1 The pulled Coulomb gas ($x>\mathrm{a}_{+}$)
      2. 2.3.2 A word on the pushed Coulomb gas ($x<\mathrm{a}_{+}$)
    4. 2.4 The case of one rectangular random matrices
      1. 2.4.1 Singular value decomposition
      2. 2.4.2 Bi-invariant ensemble
      3. 2.4.3 Bi-invariant ensemble with a wall and fixed diagonal rectangular matrices
      4. 2.4.4 Right large deviation of the top singular value
  3. 3 The tilting method
    1. 3.1 Notation
    2. 3.2 Idea of the tilting method and general expression for the rate function
    3. 3.3 Spherical Integrals as tilting functions
  4. 4 Large deviation for the sum of symmetric matrices
    1. 4.1 The sum of invariant matrices: free convolution and the R-transform
    2. 4.2 The additive spherical integral and the SSK model
    3. 4.3 Asymptotic behavior of the annealed and quenched free energies of the SSK model
    4. 4.4 Retrieving the right large deviation for one random matrix
    5. 4.5 Optimal inverse temperature for the sum
    6. 4.6 Expression for the rate function
  5. 5 Large deviation for the product of symmetric matrices
    1. 5.1 Product of invariant matrices: the multiplicative free convolution and the S-transform
    2. 5.2 The multiplicative spherical integral and the LSSK model
    3. 5.3 Asymptotic behavior of the quenched and annealed free energies of the LSSK model
    4. 5.4 Optimal temperature for the product
    5. 5.5 Expression for the rate function
  6. 6 Large deviation for the top singular value of sum of rectangular matrices
    1. 6.1 Rectangular free convolution
    2. 6.2 Rectangular spherical integral and the BSSK model
    3. 6.3 Asymptotic behavior of the annealed and quenched free energy of the BSSK model
    4. 6.4 Optimal temperature for the rectangular case
    5. 6.5 Expression for the rate function
  7. 7 Conclusion
  8. A Properties of transforms of RMT and Free Probability
    1. A.1 The Stieltjes transform and its inverse
    2. A.2 monotonic behavior of the linearizing transforms
      1. A.2.1 for the R-transform
      2. A.2.2 for the S-transform
      3. A.2.3 for the rectangular C-transform
    3. A.3 Inequalities at the edge of the free convolutions
  9. B Asymptotic of the quenched and annealed free energy
    1. B.1 For the SSK model (additive spherical integral)
      1. B.1.1 Complex integral representation and quenched free energy
      2. B.1.2 Computation of the annealed free energy
    2. B.2 For the LSSK model (multiplicative spherical integral)
    3. B.3 For the BSSK model (rectangular spherical integral)
      1. B.3.1 complex integral representation
      2. B.3.2 saddle-point computations
  10. C Large deviation for rank-one deformation
    1. C.1 Additive rank-one deformation
      1. C.1.1 The tilting method and its failure for left large deviation.
      2. C.1.2 Large deviation beyond the threshold and the Sherman-Morrison formula
      3. C.1.3 Expression for the rate function
    2. C.2 Multiplicative case
    3. C.3 Rectangular case
  11. D Rank-one plus rank-one

##  1 Introduction

Since the pioneer work of Wishart [2] and Wigner [3] in the 60s, Random Matrix
Theory (RMT) has found applications in many domains of research, ranging from
the theory of disordered system [4, 5] to telecommunications [6] and finance
[7] and even more recently to statistical learning theory [8, 9], to cite few
examples. In particular the top eigenvalue/singular value of a random matrix
plays a fundamental role in many fields: for example in statistics, Principal
Component Analysis (PCA) [10] is often used to reduce the dimensions of a raw
matrix based on the values of first top eigenvalues; in the study of the
stability of a complex system [11, 12, 13, 14, 15] the top eigenvalue of the
(opposite of the) stability matrix of a randomly linear system indicates
whether the system is locally stable or not; in the theory of disordered
system, the law of the top eigenvalue is related to the so-called complexity
of the system, see Refs. [16, 17, 18].

Therefore, a natural question in RMT and related fields is the following:
given a random symmetric (respectively rectangular) matrix $\mathbf{A}$, as
the dimensions of this matrix grows, what is the typical value of the top
eigenvalue (resp. singular value) and what is the probability to find it at a
position, say $x$, far from this typical value? When this probability is
exponentially small, one says that the top eigenvalue satisfies a _large
deviation principle_ and the goal is to calculate both the _speed of
convergence_ (what is the power in $N$ in the exponential decay?) and the
_rate function_ (what is the leading prefactor and how does it depend on the
position $x$?). The first natural case is to consider, $\mathbf{A}$ as taken
from an _invariant ensemble_ , a family containing both the famous _Gaussian
orthogonal ensemble_ (GOE) and the _Wishart_ ensemble of RMT. Based on a
direct analogy between the ensemble of eigenvalues and a Coulomb gas of
particles restricted on the real line, the large deviation can be computed
explicitly, see Refs. [19, 20, 21, 22, 23, 24] and [25] for a review. In
particular, both the speed of convergence and the rate functions depend on
whether $x$ is above or below the typical value. The former is known as the
right large deviation principle and the latter is known as the left large
deviation principle. Several studies have gone beyond the invariant ensemble
case by looking for example at generalized Wigner matrices, see Refs. [26, 27,
28, 29] or looking at small-rank deformation of an invariant ensemble, see
Refs. [30, 31, 32]. Recently, there has been an interest in the so-called
full-rank deformation:

  1. 1.

In Ref. [1], the authors studied the right large deviation for the sum
$\mathbf{C}=\tilde{\mathbf{A}}+\mathbf{O}\tilde{\mathbf{B}}\mathbf{O}^{\mathsf{T}}$,
where $\mathbf{O}$ is a random (uniform) orthogonal matrix, and
$\tilde{\mathbf{A}}$ and $\tilde{\mathbf{B}}$ are two ’fixed diagonal’
matrices, based on a tilting method with the additive spherical integral.

  2. 2.

In Ref. [33], based on similar ideas, the author studied the (right) large
deviation of the top eigenvalue of the sum
$\mathbf{C}=\tilde{\mathbf{A}}+\mathbf{B}$ where $\mathbf{B}$ is a (slight
generalization of a) GOE matrix.

  3. 3.

In Ref. [34] and based again on similar methods, the case of the large
deviation of the top eigenvalue of the product
$\mathbf{C}=\sqrt{\tilde{\mathbf{A}}}\mathbf{B}\sqrt{\tilde{\mathbf{A}}}$,
where $\tilde{\mathbf{A}}$ is a diagonal matrix and $\mathbf{B}$ is a Wishart
matrix, was studied.

In this paper, we obtain _explicitly_ the right large deviation principle for
the top eigenvalue (or singular value if the matrices are rectangular) of the
sum or the product of two _arbitrary_ matrices; by which we mean that each
matrix can be either taken from an invariant ensemble or can be a (randomly
rotated) diagonal matrix. In particular, our results allows one to recover the
three specific cases considered previously in Refs. [1, 33, 34] and to obtain
new results for cases that have not been previously considered. To obtain our
result, we first introduce an _invariant ensemble with a wall_ , a natural
generalization of the classical invariant ensemble. From the point of view of
large deviations, we show that for a proper choice of the position of the wall
(namely the wall is exactly at the edge of the spectrum), we recover the cases
of deterministic matrices. Second, we extend the tilting method introduced in
Ref. [1]. We argue that when one is considering respectively the product of
two symmetric matrices or the sum of two rectangular matrices, one should
replace the additive spherical integral with respectively the multiplicative
spherical integral and the rectangular spherical integral. In each case, we
give a natural interpretation of those spherical integrals as the partition
function of a disordered system. Based on ideas develop in Ref. [35], we give
the precise asymptotic behavior for the _annealed free energies_ of any
invariant random matrix, going beyond the case of GOE and Wishart matrices
which can tackled by direct Gaussian integration. Combining this result with
the asymptotic behavior of the _quenched free energy_ of those three
(additive, multiplicative, rectangular) spherical integrals derived in Refs.
[36, 37, 38, 39, 40] allows us to get the rate function in each case.

The rest of the paper is organized as follows: in Sec. 2, we recall the main
results concerning large deviations of the top eigenvalue of _one_ random
matrix from a (classical) invariant ensemble. This allows us to introduce the
main concepts and notations used in this paper. In particular, we define the
so-called invariant ensemble with a wall and also described the case of bi-
invariant rectangular random matrices. In Sec. 3, we introduce the main tool
to compute the rate function: the tilting method with spherical integrals. Our
description of the tilting method has been made in a general framework in
order to describe the idea of the computation for each of the three cases (sum
of symmetric matrices, product of symmetric matrices, sum of rectangular
matrices) simultaneously. We then go into detail for each case separately. We
consider the case of the sum of two symmetric matrices in Sec. 4, the product
of two symmetric matrices in Sec. 5 and the sum of rectangular matrices in
Sec. 6. In each Section, we give explicitly the expression for the rate
function together with concrete examples. This rate function admits up to
three different regimes. Based on the similarity with the rate function of the
simpler model of a rank-one deformation of one invariant random matrix,
described in App. C and the rate function of the toy model of the sum of two
rank-one matrix, described in App. D, we give a natural interpretation for
each regime. In App. A, we recall classical properties of RMT and _free
probability_ that are used in the main text and App. B contains derivation of
the quenched and annealed free energies.

##  2 Reminder on right large deviation for one random matrix: The pulled
Coulomb gas approach

Before considering the addition or multiplication of two random matrices,
let’s first briefly recall the simpler case of _one_ random matrix in an
invariant ensemble. We first start with _symmetric_ random matrices.

###  2.1 Classical Rotationally Invariant Ensemble

For an analytic confining potential $V(.)$, we say that a matrix is drawn from
a (classical rotationally) invariant ensemble if the probability to observe
the matrix $\mathbf{A}\equiv\mathbf{A}_{N}$ in a region $R$ in the space of
$(N\times N)$ symmetric matrices is defined as:

$\displaystyle\mathbb{P}_{V}\left[\mathbf{A}\in R\right]$
$\displaystyle=\frac{1}{Z_{N,V}}\int_{R}\mathrm{e}^{-\frac{N}{2}{\rm
Tr}\,V(\mathbf{A})}\mathrm{d}\mathbf{A}\,,$ (1)

where $\mathrm{d}\mathbf{A}$ is the Lebesgue measure over the space of
$(N\times N)$ symmetric matrices and $Z_{N,V}$ is a constant ensuring that
this probability is normalized to one. When considering positive semi-definite
matrices, we will implicitly assume the potential to be defined on
$\mathbb{R}_{+}$. For any orthogonal matrix $\mathbf{O}$ we have
$V(\mathbf{O}\mathbf{A}\mathbf{O}^{\mathsf{T}})=\mathbf{O}V(\mathbf{A})\mathbf{O}^{\mathsf{T}}$
which together with the cyclical property of the trace gives:

$\displaystyle\mathbb{P}_{V}\left[\mathbf{O}\mathbf{A}\mathbf{O}^{\mathsf{T}}\in
R\right]$ $\displaystyle=\mathbb{P}_{V}\left[\mathbf{A}\in R\right]$
$\displaystyle\quad\left(\forall\mathbf{O}\in\mathsf{O}(N)\right)\,,$ (2)

where $\mathsf{O}(N)$ is the group of $(N\times N)$ orthogonal matrix; hence
the name (rotationally) invariant ensemble.

Example (GOE matrices): If one considers the elements $A_{ij}$ of the matrix
$\mathbf{A}$ to be independent (up to the symmetry) Gaussian random variables
with mean zero and variance $\sigma^{2}/N$ for the off-diagonal elements and
$2\sigma^{2}/N$ for the diagonal elements, then this ensemble corresponds to
the famous _Gaussian Orthogonal Ensemble_ (GOE) for which the potential is
equal to:

$\displaystyle V(\lambda)$ $\displaystyle=\frac{\lambda^{2}}{2\sigma^{2}}\,.$
(3)

Example (Wishart matrices): Consider a $(N\times M)$ rectangular matrix
$\mathbf{X}$ with $M\geq N$, where the elements of $\mathbf{X}$ are Gaussian
independent random variables with mean zero and variance one, from which we
construct the $(N\times N)$ square matrix
$\mathbf{A}:=\frac{1}{M}\mathbf{X}\mathbf{X}^{\mathsf{T}}$. The matrix
$\mathbf{A}$ is a (Gaussian White) _Wishart_ matrix. It is rotationally
invariant with potential:

$\displaystyle V_{N}(\lambda)$
$\displaystyle=\frac{M}{N}\lambda+\left(1-\frac{M}{N}+\frac{1}{N}\right)\log\lambda\to
V(x)=\frac{\lambda}{q}+\left(1-\frac{1}{q}\right)\log\lambda\,,$ (4)

where the asymptotic behavior in Eq. (4) corresponds to the double scaling
limit $N\to\infty$ and $M\to\infty$ with $\frac{N}{M}\to q\in(0,1)$. The
invariant ensemble with the potential $V(.)$ in Eq. (4) is sometimes known as
the _Laguerre Orthogonal Ensemble_.

From there, it is a standard result of RMT that $\mathbf{A}$ admits the
following spectral decomposition
$\mathbf{A}=\mathbf{V}\mathrm{Diag}(\lambda_{1}(\mathbf{A}),\dots,\lambda_{N}(\mathbf{A}))\mathbf{V}^{\mathsf{T}}$,
where the matrix of eigenvectors $\mathbf{V}$ is taken uniformly over
$\mathsf{O}(N)$. To ease notation when it is needed, we write
$\lambda_{i}\equiv\lambda_{i}(\mathbf{A})$ for the eigenvalues of the matrix
$\mathbf{A}$. The joint density of the eigenvalues is given by:

$\displaystyle\mathcal{P}_{N}(\lambda_{1},\dots,\lambda_{N})$
$\displaystyle=\frac{1}{Z_{N}}\mathrm{exp}\left[-\frac{N}{2}\sum_{i=1}^{N}V(\lambda_{i})+\frac{1}{2}\sum_{ij|i\neq
j}\log|\lambda_{i}-\lambda_{j}|\right]\,.$ (5)

The term $\sum_{ij|i\neq j}\log|\lambda_{i}-\lambda_{j}|$ is the Jacobian of
the change of variable $\mathbf{A}\to(\mathbf{V},\\{\lambda_{i}\\})$ and is
exactly the pairwise repulsive interaction of the $2d$-Coulomb gas. In the
large $N$ limit, the empirical (random) spectral density converges to a non-
random smooth density $\mu_{A}(.)$:

$\displaystyle\mu_{\mathbf{A}}(\lambda)$
$\displaystyle:=\frac{1}{N}\sum_{i=1}^{N}\delta(\lambda-\lambda_{i}(\mathbf{A}))\underset{N\to\infty}{\to}\mu_{A}(\lambda)\,.$
(6)

where $\mu_{A}(.)$ is the solution of the Tricomi problem:

$\displaystyle\mathrm{P.V.}\int\frac{\mu_{A}(\lambda^{\prime})}{\lambda-\lambda^{\prime}}\mathrm{d}\lambda^{\prime}$
$\displaystyle=\frac{V^{\prime}(\lambda)}{2}\,,$ (7)

where $\mathrm{P.V.}$ stands for Principal Value. Example (GOE matrices and
the semi-circle distribution):  For GOE matrices with a quadratic potential
given by Eq. (3), the solution of this Tricomi equation is given by the famous
_semi-circle distribution_ :

$\displaystyle\mu_{\mathrm{sc}}(\lambda)$
$\displaystyle=\frac{\sqrt{4\sigma^{2}-\lambda^{2}}}{2\pi\sigma^{2}}\mathbb{I}_{[-2\sigma,2\sigma]}\,.$
(8)

where $\mathbb{I}_{[a,b]}$ is the indicator function, it is equal to $1$ if
$\lambda\in[a,b]$ and zero otherwise. Example (Wishart matrices and the
Marčenko-Pastur distribution): For Wishart matrices with potential given by
Eq. (4), one obtains the _Marčenko-Pastur distribution_ :

$\displaystyle\mu_{\mathrm{MP}_{q}}(\lambda)$
$\displaystyle=\frac{\sqrt{(\mathrm{a}_{+}-\lambda)(\lambda-\mathrm{a}_{-})}}{2\pi
q\lambda}\mathbb{I}_{[\mathrm{a}_{-},\mathrm{a}_{+}]}\,.$ (9)

where the edges are given by $\mathrm{a}_{\pm}=(1\pm\sqrt{q})^{2}$.

###  2.2 Invariant Ensemble with a wall

Importantly, Eq. (5) for the joint density of eigenvalues still makes sense if
we introduce a _wall_ at a position $w_{A}\geq\mathrm{a}_{+}$ beyond which the
potential is infinite. We say that a matrix is taken from an invariant
ensemble with a wall at $w_{A}$, which we denote by
$\mathbf{A}\sim\mathbb{P}_{V,w_{A}}$, if
$\mathbf{A}=\mathbf{O}\mathrm{Diag}\left(a_{1},\dots,a_{N}\right)\mathbf{O}^{\mathsf{T}}$
with $\mathbf{O}$ uniform over $\mathsf{O}(N)$ and the $\\{\lambda_{i}\\}$
follow the joint law of Eq. (5) with $V(.)$ a confining potential such that
$V(x>w_{A})=\infty$. By construction, we still have the property (2) for
random matrices taken from this ensemble. It is important to notice that the
introduction of this wall does not change the limiting equilibrium density
$\mu_{A}(.)$ since $w_{A}\geq\mathrm{a}_{+}$ and the solution of the Tricomi
problem of Eq. (7) only depends on the values of the potential between the two
edges $\mathrm{a}_{\pm}$. The introduction of these invariant ensembles with a
wall might seem odd at first, but as we will see later on, this construction
allows the study of the sum of two matrices where one (or both111When
considering two diagonal matrices, we will be considering the sum
$\mathbf{C}=\tilde{\mathbf{A}}+\mathbf{O}\tilde{\mathbf{B}}\mathbf{O}^{\mathsf{T}}$
with $\mathbf{O}$ uniform over $\mathsf{O}(N)$, such that $\tilde{\mathbf{A}}$
and $\mathbf{O}\tilde{\mathbf{B}}\mathbf{O}^{\mathsf{T}}$ are asymptotically
free and the spectrum of the eigenvalues of the sum (respectively of the
product) is given asymptotically by the _free convolution_ described in Sec.
4.1 (resp. the _multiplicative free convolution_ described in Sec. 5.1). )
matrix is a fixed diagonal matrix:

$\displaystyle\tilde{\mathbf{A}}=\begin{pmatrix}\tilde{a}_{1}&0&\dots&0\\\
0&\ddots&&\vdots\\\ \vdots&&\ddots&0\\\
0&\dots&0&\tilde{a}_{N}\end{pmatrix}\,.$ (10)

The $\\{\tilde{a}_{i}\\}$ are ’frozen’ sequences of numbers such that as we
increase the size of the matrix, their empirical distribution converges to the
same $\mu_{A}(.)$:

$\displaystyle\mu_{\tilde{\mathbf{A}}}(\lambda)$
$\displaystyle:=\frac{1}{N}\sum_{i=1}^{N}\delta(\lambda-\tilde{a}_{i})\underset{N\to\infty}{\to}\mu_{A}(\lambda)\,.$
(11)

and importantly the minimum and the maximum of the $\\{\tilde{a}_{i}\\}$
converge to the edges $\mathrm{a}_{-}$ and $\mathrm{a}_{+}$ of $\mu_{A}(.)$
(that is, in the large $N$ limit, no outlier survives). To be more concrete, a
typical example for the choice of the $\\{\tilde{a}_{i}\\}$ is to take them as
the quantiles of the distribution $\mu_{A}(.)$, sometimes called the
”classical positions” of the particles. For a given $N$ and for each $i$
ranging from $1$ to $N$, this means that $\tilde{a}_{i}$ is solution of the
integral equation:

$\displaystyle\int_{\mathrm{a}_{-}}^{\tilde{a}_{i}}\mu_{A}(\lambda)\mathrm{d}\lambda$
$\displaystyle=\frac{i}{N+1}\,,$ (12)

and by construction their empirical distribution converge to $\mu_{A}(.)$.
Another natural choice is to draw independently each $\tilde{a}_{i}$ from the
distribution $\mu_{A}(.)$. Now in general, diagonal matrices of the form of
Eq. (10) are very different from random matrices taken from an invariant
ensemble since their eigenvalues are fixed and don’t have a true repulsion as
in Eq. (5). However, when one consider the problem of large deviation of the
top eigenvalue, we will see that they behave as an invariant ensemble with a
wall exactly at the edge ($w_{A}=\mathrm{a}_{+}$). This will be made more
precise later on, but one can already see that if there is a wall exactly at
the edge, then at finite $N$ the top eigenvalue cannot fluctuate outside the
support of the bulk density $\mu_{A}(.)$ and hence it is in a sense fixed.
This construction is superfluous when considering just one matrix (in this
case the question of large deviation for the top eigenvalue of a deterministic
matrix as in Eq. (10) is trivial), but turns out to be very convenient when
considering the sum or product of two matrices.

Remark (wall at infinity and classical invariant ensemble): One may note that
we recover the case of classical invariant ensemble of Sec. 2.1 by sending
$w_{A}\to\infty$.

Figure 1: Potential and limiting density for a random matrix taken from an
invariant ensemble with a wall. Beyond the wall, the potential is infinite.
The limiting density is the same as if there were no wall $(w_{A}\to\infty)$
since $w_{A}\geq\mathrm{a}_{+}$. The black dots represent a ’typical’
configurations of the eigenvalues at finite $N$.

###  2.3 Coulomb gas approach

Coming back to random symmetric matrices from an invariant ensemble, if we
order the eigenvalues in decreasing order,
$\lambda_{1}(\mathbf{A})\geq\dots\geq\lambda_{N}(\mathbf{A})$ then we have:

$\displaystyle\lambda_{1}(\mathbf{A})$
$\displaystyle\underset{N\to\infty}{\to}\mathrm{a}_{+}\,.$ (13)

This result is valid in the formal $N\to\infty$ and a natural question is to
estimate the probability at large but finite $N$,
$\mathbb{P}\left[\lambda_{1}(\mathbf{A})\simeq x\right]$, of having the top
eigenvalue $\lambda_{1}(\mathbf{A})$ at a position $x$ different from its
typical value $\mathrm{a}_{+}$. If the density $\mu_{A}(.)$ is ’non-critical’,
by which we mean that it behaves as a square-root near the edge
$\mathrm{a}_{+}$,

$\displaystyle\mu_{A}(x)$
$\displaystyle\underset{x\nearrow\mathrm{a}_{+}}{\sim}\frac{\gamma_{0}^{3/2}}{\pi}\sqrt{\mathrm{a}_{+}-x}\quad\,,$
(14)

then the small deviations of $\lambda_{1}(\mathbf{A})$ around the limiting
value $\mathrm{a}_{+}$ are given by the Tracy-Widom law for fluctuations of
order $O(N^{-\frac{2}{3}})$, see Refs. [41, 42]. For a value of $x$ far from
the edge $\mathrm{a}_{+}$, one is outside the scope of the Tracy-Widom regime
describing typical fluctuations and one needs to estimate a very rare event
dictated by a _large deviation principle_. This can be summarized (see Ref.
[25]) by the following set of equations:

$\displaystyle\mathbb{P}\left[\lambda_{1}(\mathbf{A})\simeq x\right]$
$\displaystyle\approx\left\\{\begin{array}[]{lll}\mathrm{exp}\left[-N^{2}\Psi_{-}(x)+o(N^{2})\right]&\mbox{for
}x<\mathrm{a}_{+}\mbox{ and }|x-\mathrm{a}_{+}|\sim O(1)\,,\\\ \\\
\gamma_{0}N^{2/3}(\mathcal{F}^{(1)})^{\prime}\left(\gamma_{0}N^{2/3}(x-\mathrm{a}_{+})\right)&\mbox{for
}|x-\mathrm{a}_{+}|\sim O(N^{-\frac{2}{3}})\,,\\\ \\\
\mathrm{exp}\left[-N\Psi(x)+o(N)\right]&\mbox{for }x>\mathrm{a}_{+}\mbox{ and
}|x-\mathrm{a}_{+}|\sim O(1)\,.\end{array}\right.$ (20)

The function $\mathcal{F}^{(1)}$ in Eq. (20) is the $\beta=1$ Tracy-Widom
function. The scaling (or speed of convergence) of the large deviation
principle is different if $x$ is either above or below the edge
$\mathrm{a}_{+}$, and can be naturally interpreted thanks to the $2d$-Coulomb
gas picture.

####  2.3.1 The pulled Coulomb gas ($x>\mathrm{a}_{+}$)

Let’s consider the case $x>\mathrm{a}_{+}$ which is the main subject of this
note. Integrating the joint density in Eq. (5), one has that computing the
(logarithm of the) probability $\mathbb{P}[\lambda_{1}(\mathbf{A})\simeq x]$
is equivalent to computing the difference of energy between the configuration
of a $2d$-Coulomb gas where the top particle is _pulled_ at the position $x$
and the configuration of the unperturbed $2d$-Coulomb gas. For the perturbed
gas, we are just moving one particle away from the bulk, and thus we expect
that this perturbation does not change the equilibrium density $\mu_{A}$ of
the $N-1$ other particles inside the bulk. This induces the scaling in Eq.
(20) for $x>\mathrm{a}_{+}$, and the right rate function is given by:

$\displaystyle\Psi(x)$
$\displaystyle=\frac{1}{2}\left[V(x)-V(\mathrm{a}_{+})-2\int\log(x-\lambda)\mu_{A}(\lambda)\mathrm{d}\lambda+2\int\log(\mathrm{a}_{+}-\lambda)\mu_{A}(y)\mathrm{d}\lambda\right]\,,$
(21)

which can be written in integral form as:

$\displaystyle\Psi(x)$
$\displaystyle=\int_{\mathrm{a}_{+}}^{x}\left(\frac{V^{\prime}(t)}{2}-g_{A}(t)\right)\mathrm{d}t\,,$
(22)

where $g_{A}(.)$ is the _Stieltjes transform_ of $\mu_{A}$:

$\displaystyle g_{A}(z)$
$\displaystyle:=\int_{\mathrm{a}_{-}}^{\mathrm{a}_{+}}\frac{\mu_{A}(\lambda)}{z-\lambda}\mathrm{d}\lambda\,.$
(23)

It will be convenient to introduce the _second branch of the Stieltjes
transform_ , (see App. A.1), which, in the case of invariant ensemble, is
given for $x>\mathrm{a}_{+}$ by:

$\displaystyle\bar{g}_{A}(x)$ $\displaystyle:=V^{\prime}(x)-g_{A}(x)\,,$ (24)

so that from Eq. (22) we can interpret the rate function as (half) the area
between the two branches of the Stieltjes transform up to the position $x$:

$\displaystyle\Psi(x)$
$\displaystyle=\frac{1}{2}\int_{\mathrm{a}_{+}}^{x}\left(\bar{g}_{A}(t)-g_{A}(t)\right)\mathrm{d}t\,.$
(25)

Note that for $g_{A}(\mathrm{a}_{+})=\infty$, the rate function is finite.
Remark (Recovering the potential): From the eigenvalue density
$\mu_{A}(\lambda)$ and the rate function $\Psi(x)$, one can recover the
potential $V(x)$ (up to an arbitrary constant) over the whole range of
possible eigenvalues. Indeed, from the density one can compute the Stieltjes
transform using Eq. (23) and rewrite Eqs. (7) and (22) as222The potential for
$x<\mathrm{a}_{-}$ can be recovered similarly using the rate function for
$-A$.

$\displaystyle\frac{V^{\prime}(x)}{2}=\left\\{\begin{array}[]{lll}\mathfrak{Re}\,g_{A}(x-\mathrm{i}0^{+})&\mbox{for
}\mathrm{a}_{-}\leq x\leq\mathrm{a}_{+}\\\ \\\
\Psi^{\prime}(x)+g_{A}(x)&\mbox{for }x\geq\mathrm{a}_{+}\,.\\\
\end{array}\right.$ (29)

For random matrices that are not necessarily drawn from an invariant ensemble,
we can use this formula to define an effective potential, i.e. the potential
of the invariant ensemble that has the same density and rate functions. In
particular, for a fixed diagonal matrix with eigenvalue density $\mu_{A}(x)$
the rate function is infinite, and the effective potential is infinite beyond
$w=\mathrm{a}_{+}$ exactly as described in Sec. 2.2.

(a)  (b)

Figure 2: On the left, a representation of a ’typical’ configuration of a
pushed-to-the-origin Coulomb gas in a harmonic potential. To have the top
eigenvalue at the position $x=0$, one needs to also push all the other
eigenvalues to the right and in the large $N$ limit, this creates a different
limiting density (with black solid line) compared to the unperturbed semi-
circle distribution (with dotted line). On the right, a ’typical’
configuration corresponding to a Coulomb gas in a harmonic potential pulled at
the position $x=2.5>\mathrm{a}_{+}=2$. Only the top eigenvalue pops out of the
limiting distribution.

Example (Rate function for GOE matrices): for a GOE matrix, whose limiting
spectrum is the semi-circle distribution of Eq. (8), the Stieltjes transform
is given by:

$\displaystyle g_{\mathrm{sc}}(z)$
$\displaystyle=\frac{z-\sqrt{z^{2}-4\sigma^{2}}}{2\sigma^{2}}$
$\displaystyle(\mbox{for }z\geq 2\sigma)\,,$ (30)

and the second branch of the Stieltjes is given by:

$\displaystyle\bar{g}_{\mathrm{sc}}(z)$
$\displaystyle=\frac{z+\sqrt{z^{2}-4\sigma^{2}}}{2\sigma^{2}}$
$\displaystyle(\mbox{for }z\geq 2\sigma)\,,$ (31)

and therefore integrating according to Eq. (22), the right large deviation of
the top eigenvalue is given by the rate function:

$\displaystyle\Psi_{\mathsf{GOE}}(x)$
$\displaystyle=\frac{x\sqrt{x^{2}-4\sigma^{2}}}{4\sigma^{2}}+\log\left(\frac{2\sigma}{\sqrt{x^{2}-4\sigma^{2}}+x}\right)\,.$
(32)

The two branches of the Stieltjes transform and the rate function are given in
Fig. 3 (Left). Example (Rate function for Wishart matrices): For a (Gaussian
White) Wishart matrix, whose limiting spectrum is the Marčenko-Pastur
distribution of Eq. (9), the Stieltjes transform is given by:

$\displaystyle g_{\mathrm{MP}_{q}}(z)$
$\displaystyle=\frac{z-(1-q)-\sqrt{z-\mathrm{a}_{-}}\sqrt{z-\mathrm{a}_{+}}}{2qz}$
$\displaystyle(\mbox{for }z\geq\mathrm{a}_{+})\,,$ (33)

with $\mathrm{a}_{\pm}=(1\pm\sqrt{q})^{2}$, and the second branch is given by:

$\displaystyle\bar{g}_{\mathrm{MP}_{q}}(z)$
$\displaystyle=\frac{z-(1-q)+\sqrt{z-\mathrm{a}_{-}}\sqrt{z-\mathrm{a}_{+}}}{2qz}$
$\displaystyle(\mbox{for }z\geq\mathrm{a}_{+})\,.$ (34)

The two branches of the Stieltjes transform are illustrated in Fig. 3 (Right).
The right rate function is given by:

$\displaystyle\Psi_{\mathsf{Wish}}(x)$
$\displaystyle=\int_{\mathrm{a}_{+}}^{x}\frac{\sqrt{(t-\mathrm{a}_{+})(t-\mathrm{a}_{-})}}{2qt}\mathrm{d}t\,.$
(35)

and is also represented in Fig. 3 (Right). Note that in integral in Eq. (35)
can be computed analytically, but the result is not very enlightening. For
$q=1$ the rate function simplifies considerably, and we have:

$\displaystyle\Psi_{\mathsf{Wish}}(x)$
$\displaystyle=\frac{\sqrt{x(x-4)}}{2}+\log\left(\frac{x-2-\sqrt{x(x-4)}}{2}\right)$
$\displaystyle(\mbox{for }q=1)\,.$ (36)

In this case, one may notice the following identity:

$\displaystyle\Psi_{\mathsf{Wish}}(x^{2})$
$\displaystyle=2\Psi_{\mathsf{GOE}}(x)$ $\displaystyle(\mbox{for }q=1)\,,$
(37)

which will be discussed in more detail in Sec. 2.4.2 concerning bi-invariant
rectangular matrices.

(a)  (b)

Figure 3: On the left, the rate function (in red) for the largest eigenvalue
of a GOE matrix with $\sigma=1$, whose expression is given by Eq. (32). On the
right, the rate function (in red) for the largest eigenvalue of a Wishart
matrix with $q=1$, given by Eq. (36). In each case, the rate function is
infinite for values below the edge of the limiting distribution and is
otherwise given as half the area between the curve of the second branch of the
Stieltjes transform (in blue) and the curve of the Stieltjes transform (in
cyan), see Eq. (25).

Remark (Behavior near the edge and the Tracy-Widom ’3/2’ scaling): If one is
looking at a non-critical density satisfying the condition of Eq. (14) near
the edge, then both the Stieltjes and its second branch behaves near the top
edge as:

$\displaystyle\begin{array}[]{ll}g_{A}(\mathrm{a}_{+}+\epsilon)&=g_{A}(\mathrm{a}_{+})-\gamma_{0}^{\frac{3}{2}}\sqrt{\mathrm{a}_{+}+\epsilon}+o(\sqrt{\epsilon})\,,\\\
\bar{g}_{A}(\mathrm{a}_{+}+\epsilon)&=g_{A}(\mathrm{a}_{+})+\gamma_{0}^{\frac{3}{2}}\sqrt{\mathrm{a}_{+}+\epsilon}+o(\sqrt{\epsilon})\,.\end{array}$
(40)

so approximating the integral of Eq. (25) by Euler method, at first order one
has for the rate function:

$\displaystyle\Psi(\mathrm{a}_{+}+\epsilon)$
$\displaystyle=\frac{2\gamma^{3/2}}{3}\epsilon^{3/2}+o\left(\epsilon^{\frac{3}{2}}\right)\,,$
(41)

and hence the probability behaves as

$\displaystyle\mathbb{P}\left[\lambda_{1}(\mathbf{C})\simeq
x\right]\approx\mathrm{exp}\left[-\frac{2}{3}u^{3/2}+o(N)\right]\quad\quad\text{where
}\,u=\gamma N^{2/3}(x-\mathrm{c}_{+})\,.$ (42)

The scaling of the reduced variable $u$ and the asymptotic behavior matches
the large argument behavior of the Tracy-Widom regime which describes the
probability of finding an eigenvalue near the edge, see Ref. [25].

####  2.3.2 A word on the pushed Coulomb gas ($x<\mathrm{a}_{+}$)

For $x<\mathrm{a}_{+}$, one can still use the Coulomb gas analogy, but now the
perturbed $2d$-Coulomb gas is compressed such that its top particle is at the
position $x$. Unlike the case $x>\mathrm{a}_{+}$, the equilibrium measure in
the bulk is modified, since one needs to _push_ a large fraction of the
particles to satisfy this constraint. This explains the different scaling in
Eq. (20) in this scenario, and we refer to Ref. [25] for a description of the
rate function in this case. This left large deviation will not be discussed in
the rest of this note.

###  2.4 The case of one rectangular random matrices

In this paragraph, we describe the case of rectangular random matrices. The
reader only interested in symmetric random matrices may skip this section and
jump directly to Sec. 3.

####  2.4.1 Singular value decomposition

Let $\mathbf{A}\equiv\mathbf{A}_{N,M}$ be a $(N\times M)$ real matrix. Without
loss of generality, we consider $M\geq N$. We will be interested in the double
scaling limit where $N\to\infty$ and $M\to\infty$ but the ratio stays
finite:333In the rest of this section and when we are considering rectangular
matrices, all limits are assumed to be taken under this setting.

$\displaystyle\frac{N}{M}$ $\displaystyle\to q\in(0,1)\,.$ (43)

The singular value decomposition (SVD) of $\mathbf{A}$ is given by:

$\displaystyle\mathbf{A}$
$\displaystyle=\mathbf{U}_{1}\mathbf{\Sigma}\mathbf{V}_{1}^{\mathsf{T}}\,,$
(44)

where $\mathbf{U}_{1}$ (resp. $\mathbf{V}_{1}$) is a $(N\times N)$ (resp.
$(M\times M)$) orthogonal matrix, $\mathbf{\Sigma}$ is a diagonal rectangular
matrix:

$\displaystyle\mathbf{\Sigma}$
$\displaystyle=\begin{pmatrix}s_{1}(\mathbf{A})&0&0\\\
0&\ddots&0&\mathbf{0}_{M-N}\\\ 0&0&s_{N}(\mathbf{A})\end{pmatrix}\,,$ (45)

and the $s_{i}(\mathbf{A})$ are the _singular values_ of $\mathbf{A}$. They
are related to the eigenvalues of the matrix
$\mathbf{A}\mathbf{A}^{\mathsf{T}}$ by the relation:

$\displaystyle s_{i}(\mathbf{A})$
$\displaystyle=\sqrt{\lambda_{i}(\mathbf{A}\mathbf{A}^{\mathsf{T}})}\,.$ (46)

####  2.4.2 Bi-invariant ensemble

We say that a rectangular matrix is taken from a bi-invariant ensemble if its
law can be written as:

$\displaystyle\mathbb{P}_{V}\left[\mathbf{A}\in R\right]$
$\displaystyle=\frac{1}{Z_{N,V}}\int_{R}\mathrm{e}^{-\frac{N}{2}{\rm
Tr}\,V(\mathbf{A}\mathbf{A}^{\mathsf{T}})}\mathrm{d}\mathbf{A}\,,$ (47)

for an analytic potential $V(.)$, since it satisfies the property:

$\displaystyle\mathbb{P}_{V}\left[\mathbf{U}\mathbf{A}\mathbf{V}^{\mathsf{T}}\in
R\right]$ $\displaystyle=\mathbb{P}_{V}\left[\mathbf{A}\in R\right]\,,$ (48)

for any $\mathbf{U}\in\mathsf{O}(N)$ and any $\mathbf{V}\in\mathsf{O}(M)$. The
change of variable $\mathbf{A}\to\left(\mathbf{U},\\{s_{i}\equiv
s_{i}\left(\mathbf{A}\right)\\},\mathbf{V}\right)$ introduces a Jacobian [43]
which is equal to:

$\displaystyle\mathrm{Jac}\left[\mathbf{A}\to\left(\mathbf{U},\\{s_{i}\\},\mathbf{V}\right)\right]$
$\displaystyle\propto\mathrm{e}^{\frac{1}{2}\sum_{i,j|i\neq
j}\log|s_{i}^{2}-s_{j}^{2}|}\prod_{i=1}^{N}s_{i}^{M-N}\,.$ (49)

In this case, the matrix $\mathbf{A}$ admits the SVD of Eq. (44) where
$\mathbf{U}_{1}$ (resp. $\mathbf{V}_{1}$) is taken uniformly over
$\mathsf{O}(N)$ (resp. over $\mathsf{O}(M)$) and the singular values $s_{i}$
admit the following joint density:

$\displaystyle\mathcal{P}_{N}(s_{1},\dots,s_{N})$
$\displaystyle=\frac{1}{Z_{N}}\mathrm{exp}\left[-\frac{N}{2}\sum_{i=1}^{N}\left(V(s_{i}^{2})+\left(1-\frac{M}{N}\right)\log(s_{i}^{2})\right)+\frac{1}{2}\sum_{i,j|i\neq
j}\log|s_{i}^{2}-s_{j}^{2}|\right]\,.$ (50)

If we introduce the modified potential

$\displaystyle\tilde{V}_{q}(x)$
$\displaystyle:=V(x)+\left(1-\frac{1}{q}\right)\log x\,,$ (51)

with $q=N/M$, this can be written as:

$\displaystyle\mathcal{P}_{N}(s_{1},\dots,s_{N})$
$\displaystyle=\frac{1}{Z_{N}}\mathrm{exp}\left[-\frac{N}{2}\sum_{i=1}^{N}\tilde{V}_{q}(s_{i}^{2})+\frac{1}{2}\sum_{i,j|i\neq
j}\log|s_{i}^{2}-s_{j}^{2}|\right]\,.$ (52)

Note that if we do the change of variable
$\\{s_{i}(\mathbf{A})\\}\to\\{\lambda_{i}\left(\mathbf{A}\mathbf{A}^{\mathsf{T}}\right)\\}$
given by Eq. (48) in the joint law of Eq. (52), we have that the matrix
$\mathbf{A}\mathbf{A}^{\mathsf{T}}$ is taken from an invariant ensemble with
the modified potential $\tilde{V}_{q}(.)$ (plus a vanishing term coming from
the change of variable).444Taking into account the Jacobian amounts to do the
change $\tilde{V}_{q}(x)\to\tilde{V}_{q}(x)+\frac{1}{N}\log x$. In particular,
this vanishing term prevents the eigenvalues of the matrix
$\mathbf{A}\mathbf{A}^{\mathsf{T}}$ to be negative for $q=1$ where the
coefficient in front of the logarithm in Eq. (51) is null in this case. For
$q=1$, one has to consider $\tilde{V}_{1}(.)$ as being infinite for negative
values. As a consequence, in the large double scaling limit of Eq. (43), the
empirical singular value distribution converges to a smooth limit:

$\displaystyle\rho_{\mathbf{A}}(s)$
$\displaystyle:=\frac{1}{N}\sum_{i=1}^{N}\delta(s-s_{i}(\mathbf{A}))\to\rho_{A}(s)\,,$

where the limiting singular value distribution (LSVD) $\rho_{A}$ is given as:

$\displaystyle\rho_{A}(s)$ $\displaystyle=2s\mu_{AA^{\mathsf{T}}}(s^{2})\,,$
(53)

where $\mu_{AA^{\mathsf{T}}}$ is solution of the Tricomi problem of Eq. (7)
with the potential $V(.)$ replaced by $\tilde{V}_{q}(.)$ of Eq. (51):

$\displaystyle\mathrm{P.V.}\int\frac{\mu_{AA^{\mathsf{T}}}(\lambda^{\prime})}{\lambda-\lambda^{\prime}}\mathrm{d}\lambda^{\prime}$
$\displaystyle=\frac{\tilde{V}^{\prime}_{q}(\lambda)}{2}\,.$ (54)

The edges $\mathrm{a}_{\pm}$ of the LSVD $\rho_{A}$ are the square root of the
edges of the distribution $\mu_{AA^{\mathsf{T}}}$. Equivalently, the LSVD
$\rho_{A}$ is the solution of the following equation:

$\displaystyle\mathrm{P.V.}\int\frac{\rho_{A}(s^{\prime})}{s^{2}-{s^{\prime}}^{2}}\mathrm{d}s^{\prime}$
$\displaystyle=\frac{\tilde{V}^{\prime}_{q}(s^{2})}{2}\,,$ (55)

which can be directly seen from the joint law of Eq (52).

Example (Gaussian rectangular random matrices and Ginibre matrices): Let’s
consider a $(N\times M)$ matrix $\mathbf{A}$ with Gaussian i.i.d entries with
mean zero and variance $1/M$. This corresponds to $V(x)=\frac{x}{q}$. The
matrix $\frac{1}{\sigma^{2}}\mathbf{A}\mathbf{A}^{\mathsf{T}}$ is a Gaussian
white Wishart matrix with density given by Eq. (9). Since The LSVD of the
matrix $\mathbf{A}$ is related to the spectrum of Marčenko-Pastur distribution
by
$\rho_{A}(x)=2\,\frac{x}{\sigma^{2}}\,\mu_{\mathrm{MP}_{q}}(\frac{x^{2}}{\sigma^{2}})$,
one gets:

$\displaystyle\rho_{A}(s)$
$\displaystyle=\frac{\sqrt{4q\sigma^{4}-(s^{2}-\sigma^{2}(1+q))^{2}}}{\pi\sigma^{2}qs}\mathbb{I}_{[\sigma(1-\sqrt{q}),\sigma(1+\sqrt{q})]}\,.$
(56)

In particular the case $q=1$ corresponds to _Ginibre random matrices_ and Eq.
(56) becomes the so-called _quarter-circle distribution:_

$\displaystyle\rho_{A}(s)$
$\displaystyle=\frac{\sqrt{4\sigma^{2}-s^{2}}}{\pi\sigma^{2}}\mathbb{I}_{[0,2\sigma]}\,.$
(57)

####  2.4.3 Bi-invariant ensemble with a wall and fixed diagonal rectangular
matrices

Similarly to the symmetric case, we say that a rectangular random matrix is
taken from a bi-invariant ensemble _with a wall_ at the position $w_{A}$ if
its SVD is given by Eq. (44) with $\mathbf{U}$ (resp. $\mathbf{V}$) taken
uniformly over $\mathsf{O}(N)$ (resp. over $\mathsf{O}(M)$) and the singular
values follow the law of Eq. (52) with now $\tilde{V}(x>w_{A}^{2})=\infty$,
such that no singular value can be higher than the position $w_{A}$. The
introduction of this wall will allow us to study the sum of rectangular
matrices where one (or both555When considering two diagonal rectangular
matrices, we will be considering the sum
$\mathbf{C}=\tilde{\mathbf{A}}+\mathbf{U}\tilde{\mathbf{B}}\mathbf{V}^{\mathsf{T}}$
with $\mathbf{U}$ and $\mathbf{V}$ uniform over $\mathsf{O}(N)$ and
$\mathsf{O}(M)$ respectively, such that $\tilde{\mathbf{A}}$ and
$\mathbf{U}\tilde{\mathbf{B}}\mathbf{V}^{\mathsf{T}}$ are asymptotically bi-
free, and the spectrum of the singular values of the sum is given
asymptotically by the _rectangular free convolution_ described in Sec. 6.1. )
of them is a _fixed diagonal rectangular matrix_ of the form:

$\displaystyle\tilde{\mathbf{A}}$
$\displaystyle=\begin{pmatrix}\tilde{a}_{1}&0&0\\\
0&\ddots&0&\mathbf{0}_{M-N}\\\ 0&0&\tilde{a}_{N}\end{pmatrix}\,,$ (58)

where the $\\{\tilde{a}_{i}\\}$ are frozen sequence of number such that at
large $N$, the empirical distribution of the $\\{\tilde{a}_{i}\\}$ converges
to the same $\rho_{A}$:

$\displaystyle\rho_{\tilde{\mathbf{A}}}(s)$
$\displaystyle:=\frac{1}{N}\sum_{i=1}^{N}\delta(s-\tilde{a}_{i})\to\rho_{A}(s)\,.$

####  2.4.4 Right large deviation of the top singular value

As in the symmetric case, the goal is to estimate for large $N$, the
probability of finding the top singular value at a position $x$ above its
typical value given by the edge $\mathrm{a}_{+}$ of $\rho_{A}$:

$\displaystyle\mathbb{P}\left[s_{1}\left(\mathbf{A}\right)\simeq x\right]$
$\displaystyle=\mathrm{exp}\left[-N\Phi(x)\right]\,.$ (59)

From the relation (48), we have:

$\displaystyle\mathbb{P}\left[s_{1}\left(\mathbf{A}\right)\simeq x\right]$
$\displaystyle=\mathbb{P}\left[\lambda_{1}(\mathbf{A}\mathbf{A}^{\mathsf{T}})\simeq
x^{2}\right]\,,$ (60)

and hence

$\displaystyle\Phi(x)$ $\displaystyle=\Psi_{AA^{\mathsf{T}}}(x^{2})\,,$ (61)

where $\Psi_{AA^{\mathsf{T}}}(.)$ is the rate function associated with the
invariant matrix $\mathbf{A}\mathbf{A}^{\mathsf{T}}$ in the potential
$\tilde{V}_{q}(.)$. Namely, using Eq. (25), we have:

$\displaystyle\Phi(x)$
$\displaystyle=\frac{1}{2}\int_{\mathrm{a}_{+}^{2}}^{x^{2}}\left(\bar{g}_{AA^{\mathsf{T}}}(t)-g_{AA^{\mathsf{T}}}(t)\right)\mathrm{d}t\,,$
(62)

where $g_{AA^{\mathsf{T}}}$ and $\bar{g}_{AA^{\mathsf{T}}}$ are respectively
the Stieltjes and second branch of the Stieltjes of $\mu_{AA^{\mathsf{T}}}$.
Equivalently, we can write the rate function $\Phi$ as:

$\displaystyle\Phi(x)$
$\displaystyle=\int_{\mathrm{a}_{+}}^{x}t\left(\tilde{V}^{\prime}_{q}(t^{2})-\frac{1}{2}\int_{\mathrm{a}_{-}}^{\mathrm{a}_{+}}\frac{\rho_{A}(s)}{t^{2}-s^{2}}\mathrm{d}s\right)\mathrm{d}t\,.$
(63)

Remark (square matrix and symmetrized density): In the case where $q=1$,
corresponding to the case where $\mathbf{A}$ is an (asymptotic) square matrix
but not necessarily symmetric, there exist a nice relation with the symmetric
case of Sec. 2.1. Let’s consider a _symmetric_ matrix
$\hat{\mathbf{A}}\sim\mathbb{P}_{\hat{V}}(.)$, where its potential
$\hat{V}(.)$ is related to the potential $V(.)$ of Eq. (47) of the square (but
non-symmetric) matrix $\mathbf{A}$ by:

$\displaystyle\hat{V}(\lambda)$ $\displaystyle:=\frac{V(\lambda^{2})}{2}\,.$
(64)

In the large $N$ limit, the limiting density $\mu_{\hat{A}}(.)$ of
$\hat{\mathbf{A}}$ satisfies Eq. (7) with $V(.)$ replaced by $\hat{V}(.)$,
that is using Eq. (64):

$\displaystyle\mathrm{P.V.}\int\frac{\mu_{\hat{A}}(\lambda^{\prime})}{\lambda-\lambda^{\prime}}\mathrm{d}\lambda^{\prime}$
$\displaystyle=\lambda V^{\prime}(\lambda^{2})\,.$ (65)

On the other hand, for $q=1$ since $V_{1}(\lambda^{2})=V(\lambda^{2})$ (see
Eq. (51)), Eq. (55) reads:

$\displaystyle\mathrm{P.V.}\int\frac{\rho_{A}(\lambda^{\prime})}{\lambda^{2}-{\lambda^{\prime}}^{2}}\mathrm{d}\lambda^{\prime}$
$\displaystyle=\frac{V^{\prime}(\lambda^{2})}{2}\,,$ (66)

so using the identity:

$\displaystyle\frac{1}{\lambda^{2}-{\lambda^{\prime}}^{2}}$
$\displaystyle=\frac{1}{2\lambda}\left(\frac{1}{\lambda-{\lambda^{\prime}}}+\frac{1}{\lambda+{\lambda^{\prime}}}\right)\,,$
(67)

we have:

$\displaystyle\mathrm{P.V.}\int\frac{1}{\lambda-\lambda^{\prime}}\left(\frac{\rho_{A}(\lambda^{\prime})+\rho_{A}(-\lambda^{\prime})}{2}\right)\mathrm{d}\lambda^{\prime}$
$\displaystyle=\lambda V^{\prime}(\lambda^{2})\,,$ (68)

and hence by comparing Eq. (65) and (68), we have that the distribution
$\mu_{\hat{A}}(.)$ is the _symmetrized_ distribution of $\rho_{A}(.)$:

$\displaystyle\mu_{\hat{A}}(.)$
$\displaystyle=\frac{\rho_{A}(.)+\rho_{A}(-.)}{2}\,.$ (69)

Furthermore, using again the identity Eq. (67) in Eq. (63) for $q=1$ with the
definition of $\hat{V}(.)$ given by Eq. (64), we can write the rate function
as:

$\displaystyle\Phi(x)$
$\displaystyle=\int_{\mathrm{a}_{+}}^{x}\left(\hat{V}^{\prime}(t)-2g_{\hat{A}}(t)\right)\mathrm{d}t=2\Psi_{\hat{A}}(x)\,.$
(70)

where $g_{\hat{A}}$ is the Stieltjes transform of $\mu_{\hat{A}}$, and
$\Psi_{\hat{A}}$ is the rate function associated to the large deviation of the
largest eigenvalue of $\hat{\mathbf{A}}$. In other words, for a bi-invariant
square matrix, the rate function associated to the largest singular value is
_twice_ the one associated to the largest eigenvalue of the invariant
symmetric matrix whose limiting eigenvalue density is equal to the symmetrized
distribution of $\rho_{A}$. This can be heuristically guessed by remarking
that for $q=1$, we can express Eq. (52) as:

$\displaystyle\mathcal{P}_{N}(s_{1},\dots,s_{N})$
$\displaystyle=\frac{1}{Z_{N}}\mathrm{exp}\left[-\frac{N}{2}\sum_{i=1}^{N}\left[\hat{V}(s_{i})+\hat{V}(-s_{i})\right]+\frac{1}{2}\sum_{i,j|i\neq
j}\log|s_{i}-s_{j}|+\frac{1}{2}\log|s_{i}+s_{j}|\right]\,,$ (71)

where the potential $\hat{V}(x)$ is symmetric by construction. The joint law
of the $N$ (positive) singular values can be interpreted as the law of $2N$
eigenvalues following the usual Eq. (5) where the first $N$ variables are
constraint to be positive and each of the last $N$ is constraint to equal
minus its positive counterpart. In the large $N$ limit, these constraints are
irrelevant: the two problems have the same density (which does not depend on
$N$) and differ by a factor of two for the rate function (which as an explicit
$N$ factor).

Example (Rate function for Ginibre matrices): If $\mathbf{A}$ is a Ginibre
matrix, the LSVD is the quarter-circle law of Eq. (57) and its symmetrized
density is the semi-circle law of Eq. (8). As a consequence, the rate function
$\Phi_{\mathsf{Gin}}$ of the largest singular value of a Ginibre matrix is
given by:

$\displaystyle\Phi_{\mathsf{Gin}}(x)=2\Psi_{\mathsf{GOE}}(x)$
$\displaystyle=\frac{x\sqrt{x^{2}-4\sigma^{2}}}{2\sigma^{2}}+2\log\left(\frac{2\sigma}{\sqrt{x^{2}-4\sigma^{2}}+x}\right)\,.$
(72)

For $\sigma=1$, $\mathbf{A}\mathbf{A}^{\mathsf{T}}$ is a Wishart with shape
parameter $q=1$, so Eq. (72) and Eq. (61) give the relation (37).

Figure 4:  The Rate function (in red solid line) for the largest singular
value of a Ginibre random matrix with $\sigma=1$. This rate function is twice
the rate function of the largest eigenvalue of a GOE matrix with also
$\sigma=1$, see Eq. (72).

##  3 The tilting method

###  3.1 Notation

In the previous section, we have reviewed the results concerning the right
large deviation of the top eigenvalue (or singular value) of _one_ random
matrix. The goal of this section is to describe the general framework to
tackle the case of the large deviation of the largest eigenvalue (resp.
singular value) of a matrix $\mathbf{C}$ given as the (free) sum or product of
two symmetric (resp. rectangular) random matrices $\mathbf{A}$ and
$\mathbf{B}$. In the large $N$ limit, the limiting density $\mu_{C}$ (resp.
$\rho_{C}$) of eigenvalues (resp. singular values) of the matrix $\mathbf{C}$
is given by _free probability theory_ and depends precisely on the operation
(sum of symmetric matrices, product of symmetric matrices or sum of
rectangular matrices666As we will see, the case of the product of two
rectangular matrices can be deduced from the symmetric case.) one is
considering. This will be discussed in each dedicated section, see Secs. 4,
Sec. 5 and Sec. 6 respectively, but we argue that the strategy to get the
right large deviation is the same. To put everything in the same framework we
denote by

$\displaystyle\zeta_{i}\left(\mathbf{C}\right):=\left\\{\begin{array}[]{ll}\lambda_{i}(\mathbf{C})&\mbox{if
}\mathbf{C}\mbox{ is the sum/product of symmetric matrices}\,,\\\ \\\
s_{i}(\mathbf{C})&\mbox{if }\mathbf{C}\mbox{ is the sum of rectangular
matrices}\,.\end{array}\right.$ (76)

In any case, if we denote by $\mathrm{c}_{+}$ the edge of the spectrum of
$\mu_{C}$ (resp. $\rho_{C}$), we have:

$\displaystyle\zeta_{1}\left(\mathbf{C}\right)$
$\displaystyle\underset{N\to\infty}{\to}\mathrm{c}_{+}\,,$ (77)

and similar to the one-random matrix model, the natural question is to
estimate the probability
$\mathbb{P}\left[\zeta_{1}\left(\mathbf{C}\right)\simeq x\right]$ at large but
finite $N$. Since one recovers the original setting by taking the limit
$\mathbf{B}\to\mathbf{0}$ (the null matrix) in the additive case or
$\mathbf{B}\to\mathbf{I}$ (the identity matrix) in the multiplicative case,
one should expect again to have the same different scaling as before for
$x>\mathrm{c}_{+}$. As before, to tackle both the symmetric and rectangular
cases at the same time, let’s denote by:

$\displaystyle\Pi_{C}(x):=\left\\{\begin{array}[]{ll}\Psi_{C}(x)&\mbox{if
}\mathbf{C}\mbox{ is the sum/product of symmetric matrices}\,,\\\ \\\
\Phi_{C}(x)&\mbox{if }\mathbf{C}\mbox{ is the sum of rectangular
matrices}\,,\end{array}\right.$ (81)

such that $\Pi_{C}$ is the rate function we want to compute:

$\displaystyle\mathbb{P}\left[\zeta_{1}(\mathbf{C})\simeq x\right]$
$\displaystyle\approx\mathrm{exp}\left[-N\Pi_{C}(x)+o(N)\right]\quad\quad\text{for}\,x>\mathrm{c}_{+}\,,$
(82)

from the knowledge of the laws of the matrices $\mathbf{A}$ and $\mathbf{B}$.

###  3.2 Idea of the tilting method and general expression for the rate
function

The starting point for the derivation of the large deviation of _one_ random
matrix is the joint law density of the eigenvalues/singular values
$\zeta_{i}(\mathbf{C})$, from which we can use the Coulomb gas approach. For
the sum or the product of matrices, one does not have a simple expression for
the joint density $\mathcal{P}_{N}\left(\zeta_{1},\dots,\zeta_{N}\right)$ of
the eigenvalues/singular values, except in some specific cases, see Ref. [44].
Instead of directly looking at the matrix $\mathbf{C}$, the idea introduced in
Ref. [1] in the context of RMT, is to look at a _weighted_ realization of this
matrix. If we denote by $\mathcal{P}(.)$ the probability density of the random
matrix $\mathbf{C}$ in the space of symmetric/rectangular matrix, let’s
consider _another_ random matrix $\mathbf{C}^{\prime}$ whose probability
density is given by:

$\displaystyle\mathcal{P}^{(\theta)}(\mathbf{C}^{\prime})$
$\displaystyle:=\frac{\mathcal{Z}_{\mathbf{C}^{\prime}}(\theta)}{\mathbb{E}\left[\mathcal{Z}_{\mathbf{C}}(\theta)\right]}\mathcal{P}(\mathbf{C}^{\prime})\,.$
(83)

where
$\mathcal{Z}_{\mathbf{C}^{\prime}}(\theta)\equiv\mathcal{Z}(\zeta_{1}(\mathbf{C}^{\prime}),\dots,\zeta_{N}(\mathbf{C}^{\prime}),\theta)$
is a function whose precise expression will be given later on. For now, let us
note that this function depends on $\mathbf{C}^{\prime}$ only through its
eigenvalues/singular values in addition to the free parameter $\theta$. The
expectation in the denominator of the RHS of Eq. (83) is an average over
$\mathbf{C}\sim\mathcal{P}(.)$ (and hence an average over both $\mathbf{A}$
and $\mathbf{B}$):
$\mathbb{E}\left[\mathcal{Z}_{\mathbf{C}}(\theta)\right]:=\int\mathcal{Z}_{\mathbf{C}}(\theta)\,\mathcal{P}(\mathbf{C})\mathrm{d}\mathbf{C}$
such that $\mathcal{P}^{(\theta)}(.)$ is well normalized. Because
$\mathcal{Z}_{\mathbf{C}^{\prime}}(\theta)$ only depends on the
$\zeta_{i}(\mathbf{C}^{\prime})$ we can relate the joint density
$\mathcal{P}_{N}^{(\theta)}(\zeta_{1},\dots,\zeta_{N})$ of the
eigenvalues/singular values of $\mathbf{C}^{\prime}$ to the (unknown) joint
density $\mathcal{P}_{N}\left(\zeta_{1},\dots,\zeta_{N}\right)$ of
eigenvalues/singular values of $\mathbf{C}$:

$\displaystyle\mathcal{P}_{N}^{(\theta)}(\zeta_{1},\dots,\zeta_{N})$
$\displaystyle=\frac{\mathcal{Z}_{\mathbf{C}^{\prime}}(\theta)}{\mathbb{E}\left[\mathcal{Z}_{\mathbf{C}}(\theta)\right]}\mathcal{P}_{N}\left(\zeta_{1},\dots,\zeta_{N}\right)\,.$
(84)

by integrating Eq. (84) we have:

$\displaystyle\mathrm{Prob.}\left[\zeta_{1}(\mathbf{C}^{\prime})\simeq
x\right]$
$\displaystyle=\frac{1}{\mathbb{E}\left[\mathcal{Z}_{\mathbf{C}^{\prime}}(\theta)\right]}\int\delta(\zeta_{1}-x)\mathcal{Z}_{\mathbf{C}^{\prime}}(\theta)\mathcal{P}_{N}\left(\zeta_{1},\dots,\zeta_{N}\right)\mathrm{d}\zeta_{1}\dots\mathrm{d}\zeta_{N}\,.$
(85)

Now if we _choose_ $\mathcal{Z}_{\mathbf{C}^{\prime}}(\theta)$ such that for
large $N$ this function explicitly depend on the position of the largest
eigenvalue/singular value $\zeta_{1}\equiv\zeta_{1}(\mathbf{C}^{\prime})$ but
is self-averaging with respect to the other $\zeta_{2},\dots,\zeta_{N}$ and
thus independent of them, then the integral over these variables is nothing
else than the probability of Eq. (82) we want to estimate. At large $N$ the
density of eigenvalues (or singular values) of $\mathbf{C}^{\prime}$ is the
same as $\mathbf{C}$, but the position of top eigenvalue/singular value might
different. In other words, we have:

$\displaystyle\mathrm{Prob.}\left[\zeta_{1}(\mathbf{C}^{\prime})\simeq
x\right]$
$\displaystyle=\mathrm{e}^{N\left[J_{C}(x,\theta)-F_{C}(\theta)-\Pi_{C}(x)\right]+o(N)}\,.$
(86)

where we have introduced the _quenched free energy_ as the large $N$ limit

$\displaystyle J_{C}(x,\theta)$
$\displaystyle:\approx\frac{1}{N}\log\mathcal{Z}_{\mathbf{C}|\\{\zeta_{1}=x\\}}(\theta)\,,$
(87)

where $\\{\zeta_{1}=x\\}$ indicates that the limit is taken with the
constraint that the top eigenvalue (or singular value) is fixed at $x$; and
similarly the _annealed free energy_ as:

$\displaystyle F_{C}(\theta)$
$\displaystyle:\approx\frac{1}{N}\log\mathbb{E}\mathcal{Z}_{\mathbf{C}}(\theta)\,.$
(88)

Given a value of the free parameter $\theta$, what is the typical value
$x^{*}(\theta)$ of the largest eigenvalue (or singular value) of the matrix
$\mathbf{C}^{\prime}$ as $N$ goes to infinity? By Jensen inequality, the
function in the bracket of the RHS of Eq. (85) is always non-positive
(otherwise the probability would be higher than one) and so the typical value
$x^{*}(\theta)$ corresponds to the case where this function is exactly zero,
namely $x^{*}(\theta)$ is given as the solution of:

$\displaystyle\Pi_{C}(x^{*}(\theta))-J_{C}(x^{*}(\theta),\theta)-F_{C}(\theta)$
$\displaystyle=0\,.$ (89)

Now the idea of the tilting method is to look at Eq. (89) the other way: fix
any $x>\mathrm{c}_{+}$, can we vary the free parameter $\theta$ until it
reaches a value $\theta^{*}\equiv\theta^{*}(x)$ such that now the event
$\\{\zeta_{1}(\mathbf{C}^{\prime})\simeq x\\}$ becomes typical? If so, we can
relate the unknown rate function $\Pi_{C}$ to the annealed and quenched free
energies:

$\displaystyle\Pi_{C}(x)$
$\displaystyle=J_{C}(x,\theta^{*}(x))-F_{C}(\theta^{*}(x))\,.$ (90)

Since the LHS of (85) is a probability measure, it is always bounded by one,
so this optimal $\theta^{*}$ \- if it exists - is given by:

$\displaystyle\theta^{*}(x)$
$\displaystyle=\underset{\theta>0}{\mathrm{argsup}}\\{I_{x}(\theta):=J_{C}(x,\theta)-F_{C}(\theta)\\}\,.$
(91)

In particular, if the supremum is a maximum, one has that $\theta^{*}$ is
solution of:

$\displaystyle
I^{\prime}_{x}(\theta^{*}):=\partial_{\theta}J_{C}(x,\theta)|_{\theta^{*}}-F^{\prime}_{C}(\theta^{*})$
$\displaystyle=0\,.$ (92)

Next since we have Eq. (77), by construction $\Pi_{C}(\mathrm{c}_{+})=0$. One
can therefore take the derivative of Eq. (90) with respect to $x$ to get the
following integral representation for the rate function:

$\displaystyle\Pi_{C}(x)$
$\displaystyle=\int_{\mathrm{c}_{+}}^{x}\partial_{t}J_{C}(t,\theta^{*}(t))+\theta^{*\prime}(t)(\partial_{\theta}J_{C}(t,\theta)|_{\theta^{*}}-F^{\prime}(\theta^{*}))\mathrm{d}t\,,$
(93)

but by Eq. (92) the second term is null, so that we have the following simple
formula:

$\displaystyle\Pi_{C}(x)$
$\displaystyle=\int_{\mathrm{c}_{+}}^{x}\partial_{t}J_{C}(t,\theta^{*}(t))\mathrm{d}t$
$\displaystyle\mbox{with $\theta^{*}$ solution of Eq.\ \eqref{eq:Theta_opt.2}
}\,.$ (94)

###  3.3 Spherical Integrals as tilting functions

To summarize, we need to find a tilting function $\mathcal{Z}_{\mathbf{C}}$
which in the large $N$ limit only depends on the limiting distribution of the
matrix $\mathbf{C}$ and the position of its top eigenvalue/singular value and
such that:

  1. 1.

we can compute the partial derivatives of the quenched free energy Eq. (87),

  2. 2.

we can compute the derivative of the annealed free energy of Eq. (88),

  3. 3.

we can show that for each $x>\mathrm{c}_{+}$, there is one optimal temperature
$\theta^{*}(x)$ solution of Eq. (92) and we compute it,

then we get the rate function thanks to Eq. (94). Based on this necessary
properties, we argue that a natural candidate for the tilting function is
given by the _spherical integral_ of the operation we are considering. The
reason is twofold:

  * •

On the one hand, the quenched free energy associated to each spherical
integral has been computed before in the literature, see Refs. [36, 37, 38,
39, 40], and is known to satisfy a transition depending on the parameter
$\theta$, between a phase where it does not depend explicitly on the position
of the top eigenvalue/singular value and a phase where it does. The
asymptotics of the quenched free energy is summarized in Sec. 4.3, in Sec. 5.3
and in Sec. 6.3 for respectively the additive spherical integral, the
multiplicative spherical integral, and the rectangular spherical integral.

  * •

On the other hand, the spherical functions satisfy (by construction, as we
will see later on) the following _decomposition property:_

$\displaystyle\mathbb{E}\left[\mathcal{Z}_{\mathbf{C}}(\theta)\right]$
$\displaystyle=\mathbb{E}_{\mathbf{A}}\left[\mathcal{Z}_{\mathbf{A}}(\theta)\right]\,\mathbb{E}_{\mathbf{B}}\left[\mathcal{Z}_{\mathbf{B}}(\theta)\right]\,.$
(95)

Applying the logarithm function and dividing by $N$ Eq. (95), this means that
the annealed free energy is given by:

$\displaystyle F_{C}(\theta)$
$\displaystyle=F_{A}(w_{A},\theta)+F_{B}(w_{B},\theta)\,,$ (96)

where for $\mathbf{A}\sim\mathbb{P}_{V_{A},w_{A}}$ (and similarly for
$\mathbf{B}\sim\mathbb{P}_{V_{B},w_{B}}$), we have defined

$\displaystyle F_{A}(w_{A},\theta)$
$\displaystyle:\approx\frac{1}{N}\log\mathbb{E}_{\mathbf{A}}\mathcal{Z}_{\mathbf{A}}(\theta)\,,$
(97)

and put explicitly the dependence in the position of the wall $w_{A}$ (resp.
$w_{B}$). Since we have a precise description for the law of each matrix
$\mathbf{A}$ and $\mathbf{B}$ separately but not for $\mathbf{C}$, this allows
us to compute the annealed free energy and the results are also given in Sec.
4.3, in Sec. 5.3 and in Sec. 6.3.

As a consequence, once we have an expression for the (derivatives of) the
annealed and quenched free energies, we only need to prove that the optimal
$\theta^{*}(x)$ is well defined for each (attainable) $x$ and compute it. This
is done for each case in Sec. 4.5; Sec 5.4 and Sec. 6.4 . Injecting this
expression in Eq. (94), we can then get an expression for the rate function
and the results are given in in Sec. 4.6, in Sec. 5.5 and in Sec. 6.5.

Remark (Spherical function as partition function of spherical spin glass
model): We have used the notations of statistical physics to denote the
quantities of Eq. (87) and (88) as respectively the quenched and annealed free
energies. The reason is that the spherical integral
$\mathcal{Z}_{\mathbf{C}}(\theta)$ can be seen as the partition function over
configuration $\bm{\sigma}$ living on the sphere, with disorder $\mathbf{C}$
and inverse temperature $\theta$:

$\displaystyle\mathcal{Z}_{\mathbf{C}}(\theta)$
$\displaystyle:=\langle\mathrm{e}^{N\theta\mathcal{H}(\bm{\sigma})}\rangle\,.$
(98)

The precise description of the spin glass model will depend on the nature of
the operation we are considering and be given explicitly in Sec. 4.2 for the
case of the sum, in 5.2 for the case of the product, and in Sec. 6.2 for the
rectangular case. Let us mention that the rate function is given as the
highest difference between the quenched and the annealed free energy for a
value of the top eigenvalue/singular value being fixed. Since for such systems
it is known that both the quenched and annealed free energies are equal (this
is the so-called ’paramagnetic phase’) below a certain threshold of the
inverse temperature, the optimal inverse temperature $\theta^{*}$ in Eq. (91)
is necessarily attained after this threshold (in the so-called ’spin-glass
phase’).

Remark (the case of diagonal matrices and wall at the edge): If one of the
matrix, say $\mathbf{A}$, is replaced by a diagonal (resp. rectangular
diagonal) matrix, $\tilde{\mathbf{A}}$ as described in Sec. 2.2 (resp. as
described in Sec. 2.4.3), then the decomposition property (95) is now given
by:

$\displaystyle\mathbb{E}\left[\mathcal{Z}_{\mathbf{C}}(\theta)\right]$
$\displaystyle=\mathcal{Z}_{\tilde{\mathbf{A}}}(\theta)\,\mathbb{E}_{\mathbf{B}}\left[\mathcal{Z}_{\mathbf{B}}(\theta)\right]\,.$
(99)

Since there is no constraint on the top eigenvalue/singular value, and the
limiting density of $\tilde{\mathbf{A}}$ is the same as the one of
$\mathbf{A}$, we have for the free energy

$\displaystyle F_{C}(\theta)$
$\displaystyle=J_{A}(\mathrm{a}_{+},\theta)+F_{B}(w_{B},\theta)\,.$ (100)

Now we argue that in each case, we have the following property:

$\displaystyle\partial_{\theta}F_{A}(w_{A}=\mathrm{a}_{+},\theta)$
$\displaystyle=\partial_{\theta}J_{A}(\mathrm{a}_{+},\theta)\,,$ (101)

such that the derivative of the annealed free energy of Eq. (96) and Eq. (100)
are equal. Since no other quantity is modified, from the point of view of
large deviation, it is equivalent to replacing the diagonal matrix
$\tilde{\mathbf{A}}$ by a random matrix
$\mathbf{A}\sim\mathbb{P}_{V,w_{A}=\mathrm{a}_{+}}$ from an invariant ensemble
with the same limiting density $\mu_{A}(.)$, but with a wall at the edge. The
same argument holds if one considers the free sum of two diagonal matrices.

##  4 Large deviation for the sum of symmetric matrices

In this section, we consider the case where the matrix $\mathbf{C}$ is given
as

$\displaystyle\mathbf{C}$ $\displaystyle=\mathbf{A}+\mathbf{B}\,,$ (102)

where $\mathbf{A}\sim\mathbb{P}_{V_{A},w_{A}}$ and
$\mathbf{B}\sim\mathbb{P}_{V_{B},w_{B}}$ are two _symmetric_ matrices, each
taken from an invariant ensemble with a wall as defined in Sec. 2.2. We aim at
computing the rate function777We recall that since we are studying the
symmetric case, one needs to replace the notation
$\zeta_{i}\left(\mathbf{C}\right)$ and $\Pi_{C}(x)$ in Sec. 3 by
$\lambda_{i}\left(\mathbf{C}\right)$ and $\Psi_{C}(x)$ respectively.
$\Psi_{C}(x)$:

$\displaystyle\mathbb{P}\left[\lambda_{1}(\mathbf{C})\simeq x\right]$
$\displaystyle=\mathrm{exp}\left[-N\,\Psi_{C}(x)+o(N)\right]\quad\quad\text{for}\,x>\mathrm{c}_{+}\,,$
(103)

where $\mathrm{c}_{+}$ is the edge of the limiting spectrum $\mu_{C}$ of
$\mathbf{C}$, and $\mu_{C}$ is described by _free probability_ , see the next
paragraph.

###  4.1 The sum of invariant matrices: free convolution and the R-transform

The limiting spectral distribution $\mu_{C}(.)$ is a smooth density with right
edge $\mathrm{c}_{+}$ and is given as the unique probability distribution
solution of:

$\displaystyle\mathcal{R}_{C}(y)$
$\displaystyle=\mathcal{R}_{A}(y)+\mathcal{R}_{B}(y)\,,$ (104)

for all $y$ in the complex plane close enough to the origin, where
$\mathcal{R}_{A}$ (resp. $\mathcal{R}_{B},\mathcal{R}_{C}$) is the R-transform
of the limiting spectral distribution $\mu_{A}$ (resp. $\mu_{B},\mu_{C}$):

$\displaystyle\mathcal{R}_{A}(y)$
$\displaystyle:=g_{A}^{\langle-1\rangle}(y)-\frac{1}{y}\,,$ (105)

where $g_{A}^{\langle-1\rangle}(.)$ denotes the (functional) inverse of the
Stieltjes transform of Eq. (23). This inverse, and hence also the R-transform,
is only defined between $(0,g_{A}(\mathrm{a}_{+}))$. However, for invariant
ensemble, one can extend analytically this function for values beyond
$g_{A}(\mathrm{a}_{+})$. In this case, this corresponds to inverting the
second branch of the Stieltjes transform given by Eq. (24) and not the
Stieltjes branch itself, see App. A.1. In the following, the R-transform has
to be understood with this analytical continuation procedure. Eq. (104) shows
that the R-transform linearizes the sum, and hence it is the random matrix
analogous of the cumulant generating function of classical probability. The
limiting spectral distribution is therefore called the free (additive)
convolution of $\mu_{A}$ and $\mu_{B}$ and is usually denoted by
$\mu_{C}:=\mu_{A}\boxplus\mu_{B}$ in the RMT literature. We refer the reader
to Ref. [45, 46, 47] for more detail on the free convolution. We will be using
the two following properties of the additive free convolution:

  * •

If we denote by

$\displaystyle\mathrm{r}_{A}$
$\displaystyle:=\lim_{w_{A}\to\infty}\bar{g}_{A}(w_{A})\,,$ (106)

then for $\theta\in(0,\mathrm{r}_{A})$, $\mathcal{R}_{A}(.)$ is an increasing
function of $\theta$, see App. A.2.1.

  * •

The Stieltjes transform at the edge of the spectrum of the matrix $\mathbf{C}$
satisfies the following inequality:

$\displaystyle g_{C}(\mathrm{c}_{+})\leq$
$\displaystyle\min\left(g_{A}(\mathrm{a}_{+}),g_{B}(\mathrm{b}_{+})\right)\,,$
(107)

see App. A.3.

Example (R-transform of GOE matrices): For a GOE matrix inverting Eq. (30),
one has that the inverse of the Stieljes is given by:

$\displaystyle g_{\mathrm{sc}}^{\langle-1\rangle}(y)$
$\displaystyle=\sigma^{2}y+\frac{1}{y}\,,$ (108)

and so the R-transform is given by:

$\displaystyle R_{\mathrm{sc}}(y)=\sigma^{2}y\,.$ (109)

Example (R-transform of Wishart matrices): For a Wishart matrix inverting Eq.
(33), one has its R-transform is given by:

$\displaystyle R_{\mathrm{MP}_{q}}(y)=\frac{1}{1-qy}\,.$ (110)

###  4.2 The additive spherical integral and the SSK model

Based on what we have discussed in Sec. 3.3, we argue that a natural candidate
for the tilt function of Eq. (83) is given by the _additive spherical
function_ $\mathcal{Z}_{\mathbf{C}}(\theta)$, where for _any_ symmetric matrix
$\mathbf{M}$ and $\theta$ positive, it is defined by:

$\displaystyle\mathcal{Z}_{\mathbf{M}}(\theta)$
$\displaystyle:=\int_{\mathbb{S}^{N-1}}\mathrm{exp}\left[\frac{N\theta}{2}\bm{\sigma}^{\mathsf{T}}\mathbf{M}\bm{\sigma}\right]\,\mathrm{d}\omega(\bm{\sigma})\,,$
(111)

where $\mathbb{S}^{N-1}:=\\{\bm{\sigma}\in\mathbb{R}^{N}\text{ s.t
}\|\bm{\sigma}\|=1\\}$ is the hypersphere of radius one and
$\mathrm{d}\omega(\bm{\sigma})$ is the compact notation for the uniform
measure over $\mathbb{S}^{N-1}$ normalized to one. It is clear that this
function only depends on the eigenvalues of $\mathbf{M}$, because we can
always absorb the matrix of eigenvectors in the variable of integration
$\bm{\sigma}$ by a change of variable. Equivalently, the additive spherical
integral can be written as

$\displaystyle\mathcal{Z}_{\mathbf{M}}(\theta)$
$\displaystyle=\int_{\mathsf{O}(N)}\,\mathrm{e}^{\frac{N}{2}{\rm
Tr}\,\mathbf{O}\mathbf{M}\mathbf{O}^{\mathsf{T}}(\theta\mathbf{e}\mathbf{e}^{\mathsf{T}})}\,\mathrm{d}\mathbf{O}\,,$
(112)

where $\mathrm{d}\mathbf{O}$ is the normalized uniform Haar measure over
$\mathsf{O}(N)$. In RMT, $\mathcal{Z}_{\mathbf{M}}(\theta)$ is also known as
the rank-one Harish-Chandra-Itzykson-Zuber (HCIZ) integral, see Ref. [48, 49].
The asymptotic behavior of this spherical integral is related to the
$R$-transform, as we will see later. Note that if one replaces the rank one
matrix $\theta\mathbf{e}\mathbf{e}^{\mathsf{T}}$ in Eq. (112) by a full rank
matrix, then the asymptotic is governed by a different complex variational
principle that was first derived by Matytsin, see Refs. [50, 51, 52] and it is
an open problem to understand the crossover between the two different regimes.
By Haar invariance, for _any_ symmetric matrices $\mathbf{M}_{1}$ and
$\mathbf{M}_{2}$, we have the following property:

$\displaystyle\int_{\mathsf{O}(N)}\mathcal{Z}_{\mathbf{M}_{1}+\mathbf{O}\mathbf{M}_{2}\mathbf{O}^{\mathsf{T}}}(\theta)\mathrm{d}\mathbf{O}$
$\displaystyle=\mathcal{Z}_{\mathbf{M}_{1}}(\theta)\,\mathcal{Z}_{\mathbf{M}_{2}}(\theta)$
$\displaystyle\,.$ (113)

If we now apply this relation for $\mathbf{M}_{1}=\mathbf{A}$ and
$\mathbf{M}_{2}=\mathbf{B}$, since by Eq. (2) we have
$\mathbf{B}\overset{\text{in
law}}{=}\mathbf{O}\mathbf{B}\mathbf{O}^{\mathsf{T}}$, and $\mathbf{A}$ and
$\mathbf{B}$ are independent, we have after taking the average over both
$\mathbf{A}$ and $\mathbf{B}$, the desired property (95).

As we have argued in Sec. 3.3, we can interpret the additive spherical as the
partition function of spherical model. Indeed, we can always write Eq. (111)
as

$\displaystyle\mathcal{Z}_{\mathbf{C}}(\theta)$
$\displaystyle:=\langle\mathrm{e}^{\frac{N}{2}\theta\mathcal{H}^{\mathrm{SSK}}(\bm{\sigma})}\rangle\,,$
(114)

where to follow the standard notation in statistical physics, we have denoted
by
$\langle.\rangle\equiv\int_{\mathbb{S}^{N-1}}.\,\mathrm{d}\omega(\bm{\sigma})$
the uniform average over the spins living in the sphere of radius one and the
_Hamiltonian_ $\mathcal{H}^{\mathrm{SSK}}(.)$ is given by the quadratic form:

$\displaystyle\mathcal{H}^{\mathrm{SSK}}(\bm{\sigma})$
$\displaystyle:=\sum_{i}^{N}\lambda_{i}\left(\mathbf{M}\right)\sigma_{i}^{2}\,.$
(115)

This model is known (see Refs. [53, 54, 55, 56]) as the ($p=2$) _Spherical
Sherrington-Kirkpatrick_ (SSK in short) model888The SSK is usually introduced
with a different convention by absorbing the $N$ in the spin variable:
$\tilde{\bm{\sigma}}:=\sqrt{N}\bm{\sigma}$ so that the spins lives on a sphere
with radius $\sqrt{N}$, which does not change the spherical integral.. The
matrix $\mathbf{M}$ is the _disordered pairwise interaction_ and the parameter
$\theta$ is the _inverse temperature_ of the model. The SSK model has been
studied in detail in the literature, and in particular the case where
$\mathbf{M}$ is a GOE matrix has received a lot of attention. In this paper,
we are interested in the case where $\mathbf{M}=\mathbf{C}$ given by Eq.
(102). Qualitatively, for high temperature it is known that the system is in a
paramagnetic phase and all eigenvalue of the matrix $\mathbf{C}$ contributes
roughly equally to the partition function while for low temperature the
situation is drastically different and the system is in a spin glass phase
where the partition function is dominated by rare configurations which put
more weight on the top eigenvalue $\lambda_{1}(\mathbf{C})$.

###  4.3 Asymptotic behavior of the annealed and quenched free energies of the
SSK model

As explained in Sec. 3, in order to get the rate function of the largest
eigenvalue of $\mathbf{C}$, we need to compute the derivatives of the quenched
and annealed free energies of Eq. (87) and Eq. (88) with the partition
function given by Eq. (111). For a matrix $\mathbf{C}$ with limiting density
$\mu_{C}$ conditioned to have its largest eigenvalue $\lambda_{1}(\mathbf{C})$
fixed at the position $x$, the partial derivatives of the quenched free energy
are known [36] to be given by:

$\displaystyle\partial_{\theta}J_{C}(x,\theta)=\frac{1}{2}\left\\{\begin{array}[]{ll}\mathcal{R}_{C}(\theta)&\mbox{for
}\theta\leq g_{C}(x)\,,\\\ \\\ x-\frac{1}{\theta}&\mbox{for }\theta\geq
g_{C}(x)\,,\end{array}\right.$ (119)

and by:

$\displaystyle\partial_{x}J_{C}(x,\theta)=\frac{1}{2}\left\\{\begin{array}[]{ll}0&\mbox{for
}\theta\leq g_{C}(x)\,,\\\ \\\ \theta-g_{C}(x)&\mbox{for }\theta\geq
g_{C}(x)\,,\end{array}\right.$ (123)

Thanks to Eq. (96), the computation of the (derivative of the) annealed free
energy $F_{C}(.)$ reduced to the computation of the annealed free energies
$F_{A}(.)$ and $F_{B}(.)$ given by Eq. (97). Those annealed free energies can
be computed based on (an extension of) recent ideas developed in Ref. [35] in
the case of GOE. The proof is left in the App. B. For
$\mathbf{A}\sim\mathbb{P}_{V_{A},w_{A}}$ (and similarly for
$\mathbf{B}\sim\mathbb{P}_{V_{B},w_{B}}$), the annealed free energy is given
by

$\displaystyle\partial_{\theta}F_{A}(w_{A},\theta)=\frac{1}{2}\left\\{\begin{array}[]{ll}\mathcal{R}_{A}(\theta)&\mbox{for
}\theta\leq\bar{g}_{A}(w_{A})\,,\\\ \\\ w_{A}-\frac{1}{\theta}&\mbox{for
}\theta\geq\bar{g}_{A}(w_{A})\,,\end{array}\right.$ (127)

where $\bar{g}_{A}(x)$ is the second branch of the Stieltjes transform from
Eq.(24). Let us mention two important remarks:

Remark (wall at the edge and diagonal matrices): If we choose the position of
the wall to be exactly at the edge: $w_{A}=\mathrm{a}_{+}$, since we have the
relation

$\displaystyle\bar{g}_{A}(\mathrm{a}_{+})$
$\displaystyle=g_{A}(\mathrm{a}_{+})\,,$ (128)

Eq. (127) reads in this case:

$\displaystyle\partial_{\theta}F_{A}(w_{A}=\mathrm{a}_{+},\theta)=\frac{1}{2}\left\\{\begin{array}[]{ll}\mathcal{R}_{A}(\theta)&\mbox{for
}\theta\leq g_{A}(\mathrm{a}_{+})\,,\\\ \\\
\mathrm{a}_{+}-\frac{1}{\theta}&\mbox{for }\theta\geq
g_{A}(\mathrm{a}_{+})\,.\end{array}\right.$ (132)

Comparing Eq. (132) and Eq. (119) (with the index $C$ replaced by the index
$A$ and for $x=\mathrm{a}_{+}$), we see that we have indeed the relation (101)
such that from the point of large deviation we can consider a fixed diagonal
matrix as an invariant matrix with a wall at the edge of its distribution.

Remark (wall at infinity and classical invariant ensemble): Since classical
ensembles are obtained by taking the limit $w_{A}\to\infty$, we have for the
corresponding annealed free energy:

$\displaystyle\partial_{\theta}F_{A}(w_{A}\to\infty,\theta)=\frac{1}{2}\left\\{\begin{array}[]{ll}\mathcal{R}_{A}(\theta)&\mbox{for
}\theta\leq\mathrm{r}_{A}\,,\\\ \\\ \infty&\mbox{for
}\theta\geq\mathrm{r}_{A}\,,\end{array}\right.$ (136)

where we recall that $\mathrm{r}_{A}$ is defined as the limit of the second
branch of Stieltjes transform, see Eq. (106). One may note that the second
line of Eq. (136) is removed if $\mathrm{r}_{A}=\infty$ (which is for example
is the case for a GOE matrix, see Eq. (31) for $z\to\infty$) but is present
otherwise (which is for example is the case for a Wishart matrix for which
$\mathrm{r}_{A}=\frac{1}{q}$, see Eq. (34) for $z\to\infty$).

###  4.4 Retrieving the right large deviation for one random matrix

The next step in order to get the rate function is to show that there exist an
optimal temperature $\theta^{*}(x)$ and compute it. Let’s consider the case of
_one_ random matrix $\mathbf{A}\sim\mathbb{P}_{V}(.)$ in a classical invariant
ensemble, and let’s retrieve the expression of Eq. (25) with the tilting
method as a warm-up exercise. For any $x>\mathrm{a}_{+}$, the optimal inverse
temperature $\theta^{*}$ is given as the supremum of Eq. (91) with the
notation $J_{C}(.)$ and $F_{C}(.)$ replaced by $J_{A}(.)$ and
$F_{A}(w_{A}\to\infty,.)$.

One may notice by integrating Eqs. (119) and (136) with respect to $\theta$,
that for $\theta$ between $0$ and $g_{A}(x)$, we have:

$\displaystyle J_{A}(x,\theta)$
$\displaystyle=F_{A}(w_{A}\to\infty,\theta)=\frac{1}{2}\int_{0}^{\frac{\theta}{2}}\mathcal{R}_{A}(\theta^{\prime})\mathrm{d}\theta^{\prime}$
$\displaystyle(\mbox{for }\theta\in(0,g_{A}(x))\,.$ (137)

This corresponds to the high temperature regime (or _paramagnetic phase_) of
the system where both the annealed and quenched free energy are equal.
Necessarily, the optimal temperature $\theta^{*}(x)$, if there is one, cannot
be in this region since from Eq. (91), we want precisely the difference
between the two free energies to be as high as possible. We can therefore
restrict the range of possible optimal temperature to be in the _spin glass
phase_ , $\theta>g_{A}(x)$:

$\displaystyle\theta^{*}(x)$
$\displaystyle=\underset{\theta>g_{A}(x)}{\mathrm{argsup}}\left\\{I_{x}(\theta):=J_{A}(x,\theta)-F_{A}(w_{A}\to\infty,\theta)\right\\}\,.$
(138)

Let’s compute the derivative with respect to $\theta$ of this function
$I_{x}(\theta)$. According to Eq. (119) and Eq. (136), and the definition of
the $R$-transform given by Eq. (105), it is simply given by:

$\displaystyle I_{x}^{\prime}(\theta)$
$\displaystyle=\frac{1}{2}\left(x-g_{A}^{\langle-1\rangle}(\theta)\right)$
$\displaystyle(\mbox{for }\theta>g_{A}(x))\,.$ (139)

Here the function
$g_{A}^{\langle-1\rangle}(\theta)=\mathcal{R}_{A}(\theta)+1/\theta$ contains
the inverse of both branch of the Stieltjes transform. It is decreasing until
it reaches the value $\theta=g_{A}(\mathrm{a}+)$ and then it is increasing
until it reaches the (possible infinite) value,
$\theta=\mathrm{r}_{A}=\bar{g}_{A}(w_{A}\to\infty)$ where it goes to infinity.
Conversely, the function $I_{x}^{\prime}(\theta)$ of Eq. (139), seen as
function of $\theta$ for $x$ fixed, starts at zero at $\theta=g_{A}(x)$ and
then is increasing until it reaches the point $\theta=g_{A}(\mathrm{a}+)$ and
then decreasing again, and goes to $-\infty$ as $\theta\to\mathrm{r}_{A}$.
Thus, as one varies $\theta$ starting at $g_{A}(x)$, this function is positive
and then negative and only crosses the real axis once. As a consequence, the
supremum in Eq. (138) is a maximum and this maximum is unique. This maximum
$\theta^{*}$ is given at the unique point where the function
$I_{x}^{\prime}(\theta)$ of Eq. (139) crosses the real axis in the region
$\theta>g_{A}(\mathrm{a}_{+})$. In other words, finding $\theta^{*}$ amounts
to solve the equation:

$\displaystyle x$ $\displaystyle=g_{A}^{\langle-1\rangle}(\theta^{*}(x))$
$\displaystyle\mbox{for }\theta^{*}(x)>g_{A}(\mathrm{a}_{+})\,,$ (140)

which is nothing else than the definition of the second branch of the
Stieltjes transform, that is we have:

$\displaystyle\theta^{*}(x)$ $\displaystyle=\bar{g}_{A}(x)\,.$ (141)

If we now use the integral representation Eq. (94) of the rate function, with
the expression of Eq. (123) for the partial derivative of the quenched free
energy, together with Eq. (141) for the expression of $\theta^{*}$, we recover
Eq. (25) as expected. A plot of the function $I_{x}^{\prime}$, for
$\mathbf{A}$ a GOE matrix, is given in Fig. 5.

Figure 5: Derivative of the difference of the free energy in the case of one
GOE random matrix with $\sigma=1$ and $x=3$, given as the argument of Eq.
(138). For $\theta\leq g_{\mathrm{sc}}(x)$, this function (in black) is null
since the two free energies are equals, see Eq. (137), and this corresponds to
the paramagnetic phase. For $\theta\geq g_{\mathrm{sc}}(x)$, this function (in
brown) is increasing and then decreasing with a maximum at
$g_{\mathrm{sc}}(\mathrm{a}_{+}=2\sigma)=1$, and this corresponds to the spin
glass phase. The optimal inverse temperature (in blue) corresponds to the
value where this function crosses the real axis in the spin glass phase.

###  4.5 Optimal inverse temperature for the sum

Let’s now consider the general case given by Eq. (102). Without loss of
generality, we can consider999For general $\mathbf{A}$ and $\mathbf{B}$,
having $w_{A}\leq w_{B}$ does not imply
$\bar{g}_{A}(w_{A})\leq\bar{g}_{B}(w_{B})$ nor its converse. One can even come
up with examples where $\mathbf{A}$ has a no wall ($w_{A}\to\infty$) while
$\mathbf{B}$ has a finite wall $w_{B}$ but still
$\bar{g}_{A}(w_{A})\leq\bar{g}_{B}(w_{B})$.

$\displaystyle\bar{g}_{A}(w_{A})$ $\displaystyle\leq\bar{g}_{B}(w_{B})\,.$
(142)

We further assume the non-trivial condition:

$\displaystyle g_{C}(\mathrm{c}_{+})$ $\displaystyle<\infty\,,$ (143)

as for $g_{C}(\mathrm{c}_{+})=\infty$ (which necessarily implies
$g_{A}(\mathrm{a}_{+})=g_{B}(\mathrm{b}_{+})=\infty$ by Eq. (107)), the right
large deviation is infinite for any $x>\mathrm{c}_{+}$. As in the previous
section, we first want to show that the supremum in Eq. (91) is attained at a
unique point, where $F_{C}$ is given by the sum of Eq. (96) and the annealed
free energies $F_{A}(.)$ and $F_{B}(.)$ are given by Eq. (127). Since the
function $I_{x}(\theta)$ is given as the sum of three piece-wise functions,
let’s first note that we have the following set of inequalities:

$\displaystyle g_{C}(x)\leq g_{C}(\mathrm{c}_{+})\leq g_{A}(\mathrm{a}_{+})$
$\displaystyle\leq\bar{g}_{A}(w_{A})\leq\bar{g}_{B}(w_{B})\,.$ (144)

The first inequality is due to the fact that the Stieltjes is decreasing for
$x>\mathrm{c}_{+}$. The second inequality is the property (107) of the free
convolution. The third is due to the second branch of the Stieltjes transform
being monotonically increasing. The fourth is the previously mentioned
convention of Eq. (142). Using the asymptotics of Eqs. (119) (127) for the
quenched and annealed free energies, together with the linearizing property
(104) of the R-transform, one has the following behavior for the difference
between the derivative of the annealed and quenched free energy:

$\displaystyle
I_{x}^{\prime}(\theta)=\frac{1}{2}\left\\{\begin{array}[]{llll}0&\mbox{for
}\theta\leq g_{C}(x)\,,\\\ \\\ x-g_{C}^{\langle-1\rangle}(\theta)&\mbox{for
}g_{C}(x)\leq\theta\leq\bar{g}_{A}(w_{A})\,,\\\ \\\
x-w_{A}-\mathcal{R}_{B}(\theta)&\mbox{for
}\bar{g}_{A}(w_{A})\leq\theta\leq\bar{g}_{B}(w_{B})\,,\\\ \\\
x-w_{A}-w_{B}+\frac{1}{\theta}&\mbox{for
}\theta\geq\bar{g}_{B}(w_{B})\,.\end{array}\right.$ (152)

where if $\bar{g}_{B}(w_{B})=\infty$ one has to remove the last line and
similarly if $\bar{g}_{A}(w_{A})=\infty$, one has to remove the last two
lines. This function is represented in Fig. 6 for different values of $x$.
Let’s look at each interval separately.

  1. 1.

For $\theta<g_{C}(x)$, we are in the paramagnetic phase where both the
annealed and the free energy are equal. Since for each $x>\mathrm{c}_{+}$, we
want again the difference between the two to be as high as possible, the
optimal inverse temperature is not in this region of the phase space.

  2. 2.

For $g_{C}(x)\leq\theta\leq\bar{g}_{A}(w_{A})$, as we have seen in the simple
case of one invariant matrix, the function $\theta\mapsto
x-g_{C}^{\langle-1\rangle}(\theta)$ is increasing until it reaches
$g_{C}(\mathrm{c}_{+})$ and then decreasing.

  3. 3.

For $\bar{g}_{A}(w_{A})\leq\theta\leq\bar{g}_{B}(w_{B})$, since the
R-transform is increasing, the function $\theta\mapsto
x-w_{A}-\mathcal{R}_{B}(\theta)$ is decreasing.

  4. 4.

For $\bar{g}_{A}(w_{A})\leq\theta\leq\bar{g}_{B}(w_{B})$, the function
$\theta\mapsto x-w_{A}-w_{B}+\frac{1}{\theta}$ is decreasing.

One can easily check that the function of Eq. (152) is continuous at each
point where its behavior changes. At $\bar{g}_{A}(w_{A})$ it is equal to:

$\displaystyle I_{x}^{\prime}(\bar{g}_{A}(w_{A}))$
$\displaystyle=\frac{1}{2}\left(x-g_{C}^{\langle-1\rangle}(\bar{g}_{A}(w_{A}))\right)\,,$
(153)

and at $\bar{g}_{B}(w_{B})$, it is equal to:

$\displaystyle I_{x}^{\prime}(\bar{g}_{B}(w_{B}))$
$\displaystyle=x-w_{A}-w_{B}+\frac{1}{\bar{g}_{B}(w_{B})}\,.$ (154)

To summarize, in the spin glass phase $\theta\geq g_{C}(x)$ the function of
Eq. (152) is continuously increasing until $\theta=g_{C}(\mathrm{c}_{+})$ and
then it is continuously decreasing. For $x>w_{A}+w_{B}$, it is easy to check
that this function never crosses the real axis for values of
$\theta>g_{C}(x)$, as a consequence,

$\displaystyle\theta^{*}(x)$ $\displaystyle=\infty$ $\displaystyle\mbox{for
}x>w_{A}+w_{B}\,.$ (155)

This is expected because for the sum of two matrices we have the classical
inequality:

$\displaystyle\lambda_{1}(\mathbf{C})$
$\displaystyle\leq\lambda_{1}(\mathbf{A})+\lambda_{1}(\mathbf{B})\,,$ (156)

and since by definition of the walls, $\lambda_{1}(\mathbf{A})\leq w_{A}$ and
$\lambda_{1}(\mathbf{B})\leq w_{B}$, the top eigenvalue of $\mathbf{C}$ cannot
exceed $w_{A}+w_{B}$. So we find that the rate function is infinite for
$x>w_{A}+w_{B}$. Otherwise, for values of $x<w_{A}+w_{B}$, this function
always crosses the real axis once in this region. The correct equation for
$\theta^{*}(x)$ \- the point where the function $I^{\prime}_{x}(.)$ touches
the real axis, see Eq. (92) - depends on if the value of this function at
$\bar{g}_{A}(w_{A})$ or $\bar{g}_{B}(w_{B})$ is above or below zero, and hence
on the value of $x$. There exist three possible cases, separated by two
critical points, $x_{c_{1}}$ and $x_{c_{2}}$ defined respectively as the
solution of the RHS of Eq. (153) and the RHS of Eq. (154)) being equal to
zero, that is:

$\displaystyle x_{c_{1}}$
$\displaystyle:=g_{C}^{\langle-1\rangle}(\bar{g}_{A}(w_{A}))=w_{A}+\mathcal{R}_{B}(\bar{g}_{A}(w_{A}))\,,$
(157)

and

$\displaystyle x_{c_{2}}$
$\displaystyle:=w_{A}+w_{B}-\frac{1}{\bar{g}_{B}(w_{B})}\,.$ (158)

Note that $x_{c_{1}}\leq x_{c_{2}}$ as we have postulated Eq. (142). We have:

  1. 1.

for $\mathrm{c}_{+}<x<x_{c_{1}}$, the optimal inverse temperature is attained
in the region $g_{C}(x)\leq\theta\leq\bar{g}_{A}(w_{A})$ and so replacing in
Eq. (92) the expression of the difference of the free energies by the top line
of Eq. (152), it is solution of the same equation (140) as the one in the
simple one invariant random matrix case (with $g_{A}^{\langle-1\rangle}$
replaced by $g_{C}^{\langle-1\rangle}$) and thus we have:

$\displaystyle\theta^{*}(x)$ $\displaystyle=\bar{g}_{C}(x)\,,$ (159)

where $\bar{g}_{C}(.)$ is defined as the inverse of
$g_{C}^{\langle-1\rangle}(.)$ for values beyond $g_{C}(\mathrm{c}_{+})$.

  2. 2.

For $x_{c_{1}}<x<x_{c_{2}}$, the optimal inverse temperature is attained in
the region $\bar{g}_{A}(w_{A})\leq\theta\leq\bar{g}_{B}(w_{B})$ and so from
the expression of the second line of the RHS of Eq. (152), it is solution of:

$\displaystyle\mathcal{R}_{B}\left(\theta^{*}(x)\right)$
$\displaystyle=x-w_{A}\,,$ (160)

since the R-transform is continuously increasing, it has an inverse which we
denote by $\mathcal{R}_{B}^{\langle-1\rangle}$ so that the optimal temperature
is given by:

$\displaystyle\theta^{*}(x)$
$\displaystyle=\mathcal{R}_{B}^{\langle-1\rangle}\left(x-w_{A}\right)\,.$
(161)

  3. 3.

For $x_{c_{2}}<x<w_{A}+w_{B}$, $\theta^{*}(x)$ is attained in the region
$\theta\geq\bar{g}_{B}(w_{B})$ and so solving the third line of Eq. (152)
being equal to zero, we have:

$\displaystyle\theta^{*}(x)$ $\displaystyle=\frac{1}{w_{A}+w_{B}-x}\,.$ (162)

One can check that the piecewise function $\theta^{*}(x)$ is actually
continuously increasing.

(a)  (b)

(c)

Figure 6: Representation of the function of Eq. (152) for $\mathbf{A}$ and
$\mathbf{B}$ two matrices from a GOE ensemble with a wall at their edge, with
$\sigma_{A}=1$ and $\sigma_{B}=9/10$ and for different values of $x$. Each
color represents the correct expression of this piecewise continuous function
in a given interval. In the upper left, $x=2.75$ and the optimal inverse
temperature is attained in the first region where the brown curve crosses the
real axis. In the upper right, $x=2.85$ and the optimal inverse temperature is
attained in the second region where the yellow curve crosses the real axis. In
the center, $x=2.95$ and the optimal inverse temperature is attained in the
third region where the green curve crosses the real axis.

###  4.6 Expression for the rate function

Now that we have the expression for the optimal temperature, we can get the
expression for the (right) rate function $\Psi_{C}(.)$ thanks to Eq. (94).
This gives:

$\displaystyle\Psi_{C}(x)=\left\\{\begin{array}[]{lll}\displaystyle\frac{1}{2}\int_{\mathrm{c}_{+}}^{x}\left(\bar{g}_{C}(t)-g_{C}(t)\right)\mathrm{d}t&\mbox{for
}\mathrm{c}_{+}\leq x\leq x_{c_{1}}\,,\\\ \\\ \displaystyle
K_{1}+\frac{1}{2}\int_{x_{c_{1}}}^{x}\left(\mathcal{R}_{B}^{\langle-1\rangle}\left(t-w_{A}\right)-g_{C}(t)\right)\mathrm{d}t&\mbox{for
}x_{c_{1}}\leq x\leq x_{c_{2}}\,,\\\ \\\ \displaystyle
K_{2}+\frac{1}{2}\log\left(\frac{1}{w_{A}+w_{B}-x}\right)-\frac{1}{2}\int_{x_{c_{2}}}^{x}g_{C}(t)\mathrm{d}t&\mbox{for
}x_{c_{2}}\leq x\leq w_{A}+w_{B}\,.\end{array}\right.$ (168)

and is infinite for values of $x$ outside $(\mathrm{c}_{+},w_{A}+w_{B})$. The
constant $K_{1}$ (resp. $K_{2}$) is given when $x_{c_{1}}<\infty$ (resp. when
$x_{c_{2}}<\infty$) by:

$\displaystyle K_{1}$
$\displaystyle:=\frac{1}{2}\int_{\mathrm{c}_{+}}^{x_{c_{1}}}\left(\bar{g}_{C}(t)-g_{C}(t)\right)\mathrm{d}t\,,$
(169)

and by:

$\displaystyle K_{2}$
$\displaystyle:=K_{1}+\frac{1}{2}\log\left(\frac{1}{\bar{g}_{B}(w_{B})}\right)+\frac{1}{2}\int_{x_{c_{1}}}^{x_{c_{2}}}\left(\mathcal{R}_{B}^{\langle-1\rangle}\left(t-w_{A}\right)-g_{C}(t)\right)\mathrm{d}t\,,$
(170)

such that the rate function is continuous (and actually at least $C^{2}$) at
each critical point. Remark (Effective potential): The effective potential Eq.
(29) can be computed from the rate function and gives that $V(x)/2$ is equal
to the same expression as Eq. (168) with a plus sign instead of a minus sign
in front of $g_{C}(t)$ (including in the constants Eqs. (169) and (169)). The
potential is only defined up to an arbitrary constant chosen here such that
$V(\mathrm{c}_{+})=0$. An example of effective potential in plotted in Fig.7.

Remark (Tracy-Widom ’3/2’-scaling near the edge): Since the first regime
matches the ones of the classical case of one random matrix given by Eq. (25),
we retrieve in particular the Tracy-Widom ’3/2’-scaling of Eq. (42) near the
edge for the free convolution of non-degenerate densities, as recently
investigated in Ref. [57].

Remark (Number of critical points): In general, if one is considering the sum
of two random matrices taken from an invariant ensemble with a wall, then the
rate function has two possible critical points. However, if
$\bar{g}_{A}(w_{A})=\bar{g}_{B}(w_{B})$, then from Eqs. (157) (158), one can
see that two critical points merge, and we have at most one critical point. In
particular, this happens when one is considering the free sum
$\tilde{\mathbf{A}}+\mathbf{O}\tilde{\mathbf{A}}\mathbf{O}^{\mathsf{T}}$,
where $\tilde{\mathbf{A}}$ is a fixed diagonal matrix. Since $x_{c_{1}}$
diverges if $w_{A}\to\infty$ while $x_{c_{2}}$ diverges if either
$w_{A}\to\infty$ or $w_{B}\to\infty$, the free sum of a two fixed diagonal and
a fluctuating matrix (no wall) has at most one critical point while the sum of
two fluctuating matrices never has a critical point, in this case the rate
function is the same as that of an invariant ensemble Eq. (25).

Remark (Interpretation of the three regimes): For $c_{+}\leq x\leq x_{c_{1}}$
the rate function is the same as the one from an invariant ensemble. In this
regime, the eigenvalues (including rare large ones) of the matrix $\mathbf{C}$
behave exactly as an invariant ensemble with the potential (and its analytical
continuation outside the segment $[c_{-},c_{+}]$) $V_{C}(.)$ compatible with
the limiting density $\mu_{C}(.)$, that is, $V_{C}(.)$ is given by Eq. (7). In
particular, walls (if any) do not modify the rate function in this regime. For
$x_{c_{1}}\leq x\leq x_{c_{2}}$, the wall $w_{A}$ starts to matter, the
derivative of the rate function is now the same as if the matrix $\mathbf{A}$
were replaced by a rank-1 matrix with eigenvalue $w_{A}$ but with the still
correct $g_{C}(x)$ (see App. C). Finally, for $x_{c_{1}}\leq x\leq
w_{A}+w_{B}$, both walls matter and the derivative of the rate function is now
the same as for the sum of two rank-one matrices with eigenvalues $w_{A}$ and
$w_{B}$, again up to the correct $g_{C}(x)$ (see App. D). In particular, very
close to the maximal value $w_{A}+w_{B}$, we have:

$\displaystyle\mathbb{P}\left[\lambda_{1}\left(\mathbf{C}\right)\simeq(w_{A}+w_{B})(1-\epsilon)\right]$
$\displaystyle\approx\epsilon^{N/2}$ $\displaystyle\mbox{ for }\epsilon\ll
1\,,$ (171)

which is the asymptotic probability of two random vectors of having a squared
overlap of order $1-\epsilon$ in dimension $N$.

Example (Free sum of two diagonal semi-circle matrices): In this paragraph, we
will compute explicitly the rate function $\Psi_{C}(x)$ for a matrix
$\mathbf{C}=\tilde{\mathbf{A}}+\mathbf{O}\tilde{\mathbf{B}}\mathbf{O}^{\mathsf{T}}$,
where $\tilde{\mathbf{A}}$ and $\tilde{\mathbf{B}}$ are two fixed diagonal
matrices with the semi-circle distribution of Eq. (8) as their limiting
spectrum. Without loss of generality, let’s consider that their variance are
given by $\sigma_{A}=1$ and by $\sigma_{B}^{2}\equiv\sigma^{2}\leq 1$
respectively. The computation is equivalent to the sum of two invariant random
matrices in quadratic potentials $V_{i}(x)=x^{2}/(2\sigma_{i}^{2})$ and a wall
at $w_{i}=2\sigma_{i}$ for $i=A,B$ respectively. The free convolution of two
semi-circle distributions is again a semi-circle distribution, with variance
the sum of the variance. In other words, the limiting law and Stieltjes
transform of the matrix $\mathbf{C}$ is given respectively by Eq. (8) and Eq.
(30) with $\sigma_{C}=\sqrt{1+\sigma^{2}}$. Since we have $g_{A}(2)\leq
g_{B}(2\sigma)$ and the R-transform of the semi-circle distribution is given
by Eq. (109), the optimal inverse temperature is given by:

$\displaystyle\theta^{*}(x)=\left\\{\begin{array}[]{lll}\frac{x+\sqrt{x^{2}-4(1+\sigma^{2})}}{2(1+\sigma^{2})}&\mbox{for
}2\sqrt{1+\sigma^{2}}\leq x\leq 2+\sigma^{2}\,,\\\ \\\
\frac{x-2}{\sigma^{2}}&\mbox{for }2+\sigma^{2}\leq x\leq 2+\sigma\,,\\\ \\\
\frac{1}{2+2\sigma-x}&\mbox{for }2+\sigma\leq x\leq
2+2\sigma\,.\end{array}\right.$ (177)

and so the rate function is given for $x\in[2\sqrt{1+\sigma^{2}},2+2\sigma]$
by

$\displaystyle\Psi_{\mathrm{sc}+\mathrm{sc}}(x)=\left\\{\begin{array}[]{lllll}\displaystyle\frac{x\sqrt{x^{2}-4(1+\sigma^{2})}}{4(1+\sigma^{2})}+\log\left(\frac{2\sqrt{1+\sigma^{2}}}{\sqrt{x^{2}-4(1+\sigma^{2})}+x}\right)&\mbox{for
}2\sqrt{1+\sigma^{2}}\leq x\leq 2+\sigma^{2}\,,\\\ \\\
\displaystyle\frac{(x-2)^{2}}{4\sigma^{2}(1+\sigma)^{2}}+\frac{4\sigma^{2}+x^{2}-8x+12+x\sqrt{x^{2}-4(1+\sigma^{2})}}{8(1+\sigma^{2})}\\\
\displaystyle+\frac{1}{4}\log\left(\frac{x^{2}-2(1+\sigma^{2})-x\sqrt{x^{2}-4(1+\sigma^{2})}}{2(1+\sigma^{2})^{2}}\right)&\mbox{for
}2+\sigma^{2}\leq x\leq 2+\sigma\,,\\\ \\\
\displaystyle\frac{1}{4}\log\left(\frac{(2+\sigma)(2+\sigma-\sqrt{\sigma(4-3\sigma)})-2(1+\sigma^{2})}{2}\right)\\\
\displaystyle+\frac{6(1+\sigma^{2})-x^{2}+x\sqrt{x^{2}-4(1+\sigma^{2})}}{8(1+\sigma^{2})}\\\
\displaystyle+\frac{1}{2}\log\left(\frac{\sigma\left(2+\sigma+\sqrt{\sigma(4-3\sigma)}\right)}{(1+\sigma^{2})(2(1+\sigma)-x)\left(x+\sqrt{x^{2}-4(1+\sigma^{2})}\right)}\right)&\mbox{for
}2+\sigma\leq x\leq 2+2\sigma\,.\end{array}\right.$ (186)

and is infinite otherwise. This function has been plotted in Fig. 7 for
$\sigma=9/10$.

(a)  (b)

Figure 7: One the left, the rate function of the largest eigenvalue of the sum
of $\mathbf{A}$ and $\mathbf{B}$ from a GOE ensemble with a wall at their
edge, with $\sigma_{A}=1$ and $\sigma_{B}=\sigma=9/10$, see Eq. (186). On the
right, the effective potential, Eq. (29), for the same problem. Note that for
$x<x_{c_{1}}$, the effective potential is given by
$V(x)=x^{2}/(2(1+\sigma^{2}))$ and that nothing special happens at
$x=c_{+}=2\sqrt{1+\sigma^{2}}$. The potential is non-analytic at the two
points $x_{c_{1}}=2+\sigma^{2}$ and $x_{c_{2}}=2+\sigma$ (the two vertical
dashed lines) and diverges at $w_{A}+w_{b}=2+2\sigma$ after which it is
infinite (shaded area). Although it is impossible to see from this graph, it
deviates from the quadratic potential (dotted green line) for $x>x_{c_{1}}$.

##  5 Large deviation for the product of symmetric matrices

In this section we consider the case where the matrix $\mathbf{C}$ is given as
the symmetric product:101010One may note that since we are interested in the
eigenvalues, we can equivalently consider the product $\mathbf{A}\mathbf{B}$
since this matrix is similar to the matrix $\mathbf{C}$ and hence has the same
eigenvalues. However, unlike $\mathbf{C}$, the matrix $\mathbf{A}\mathbf{B}$
is a priori not symmetric.

$\displaystyle\mathbf{C}$
$\displaystyle=\sqrt{\mathbf{A}}\mathbf{B}\sqrt{\mathbf{A}}\,,$ (187)

where $\mathbf{A}\sim\mathbb{P}_{V_{A},w_{A}}$ and
$\mathbf{B}\sim\mathbb{P}_{V_{B},w_{B}}$ are two positive semi-definite random
matrices. In the large $N$ limit, the limiting density $\mu_{C}$ of
$\mathbf{C}$ is described by the multiplicative free convolution of the next
section and our goal is to compute the rate function $\Psi_{C}(x)$:111111We
recall once again that since we are studying the case of symmetric matrices,
one needs to replace the notation $\zeta_{i}\left(\mathbf{C}\right)$ and
$\Pi_{C}(x)$ by respectively $\zeta_{i}\left(\mathbf{C}\right)$ and
$\Psi_{C}(x)$ in Sec. 3.

$\displaystyle\mathbb{P}\left[\lambda_{1}(\mathbf{C})\simeq x\right]$
$\displaystyle\approx\mathrm{exp}\left[-N\,\Psi_{C}(x)+o(N)\right]\quad\quad\text{for}\,x>\mathrm{c}_{+}\,,$
(188)

describing the right large deviation of the top eigenvalue of $\mathbf{C}$ far
from its typical value given by the edge $\mathrm{c}_{+}$ of $\mu_{C}$. Let us
make the following important remark concerning the case of the product of
rectangular matrices. Remark (product of rectangular random matrices):  We
argue that the case of the product of two rectangular matrices reduced to the
symmetric cases. Indeed, if we have two rectangular matrices $\mathbf{A}$ and
$\mathbf{B}$ of size $(N\times M)$ and $(M\times K)$ respectively, then the
non-zero singular values of the product are given by:

$\displaystyle s_{i}\left(\mathbf{A}\mathbf{B}\right)$
$\displaystyle=\sqrt{\lambda_{i}\left(\mathbf{A}\mathbf{B}\mathbf{B}^{\mathsf{T}}\mathbf{A}^{\mathsf{T}}\right)}=\sqrt{\lambda_{i}\left(\sqrt{\mathbf{A}^{\mathsf{T}}\mathbf{A}}\mathbf{B}\mathbf{B}^{\mathsf{T}}\sqrt{\mathbf{A}^{\mathsf{T}}\mathbf{A}}\right)}\,.$
(189)

Since $\mathbf{A}^{\mathsf{T}}\mathbf{A}$ and
$\mathbf{B}\mathbf{B}^{\mathsf{T}}$ are two $(M\times M)$ symmetric matrices
we recover the case of Eq. (187).

###  5.1 Product of invariant matrices: the multiplicative free convolution
and the S-transform

The limiting density $\mu_{C}(.)$ is given as the unique probability
distribution solution of:

$\displaystyle\tilde{\mathcal{S}}_{C}(y)$
$\displaystyle=\tilde{\mathcal{S}}_{A}(y)\,\tilde{\mathcal{S}}_{B}(y)\,,$
(190)

for all $z$ in the complex plane close enough to the origin, where
$\tilde{\mathcal{S}}_{A}$ (resp.
$\tilde{\mathcal{S}}_{B},\tilde{\mathcal{S}}_{C}$) is the (modified) _S-
transform_ 121212The standard S-transform $\mathcal{S}$ is usually defined in
the literature as the reciprocal of our S-transform:
$\tilde{\mathcal{S}}(z)=\frac{1}{\mathcal{S}(z)}$. As we will see, the
convention $\tilde{\mathcal{S}}(.)$ will appear more naturally.  of the
limiting spectral distribution $\mu_{A}$ (resp. $\mu_{B},\mu_{C}$):

$\displaystyle\tilde{\mathcal{S}}_{A}(y)$
$\displaystyle:=\frac{y}{y+1}t_{A}^{\langle-1\rangle}(y)\,,$ (191)

where $t_{A}$ (resp. $t_{B},t_{C}$), is the _T-transform_ of the limiting
distribution $\mu_{A}$ (resp. $\mu_{B},\mu_{C}$), which is defined for all
$z\in\mathbb{C}\setminus[\mathrm{c}_{-},\mathrm{c}_{+}]$ by:

$\displaystyle t_{A}(z)$
$\displaystyle:=zg_{A}(z)-1=\int_{\mathrm{a}_{-}}^{\mathrm{a}_{+}}\frac{\lambda\,\mu_{A}(\lambda)}{z-\lambda}\mathrm{d}\lambda\,.$
(192)

Similarly, we define the second branch of the T-transform as:

$\displaystyle\bar{t}_{A}(z)$ $\displaystyle:=z\bar{g}_{A}(z)-1\,.$ (193)

The function $\bar{t}_{A}(.)$ is the inverse of $t_{A}^{\langle-1\rangle}(.)$
continued to values higher than $t_{A}(\mathrm{a}_{+})$. As a consequence of
Eq. (190), the logarithm of the S-transform linearizes the multiplication and
hence it is the random matrix analogous of the (log of the) Mellin transform
of classical probability. The limiting spectral distribution is therefore
called the _free multiplicative convolution_ of $\mu_{A}$ and $\mu_{B}$ and is
usually denoted by $\mu_{C}:=\mu_{A}\boxtimes\mu_{B}$ in the RMT literature.
We will be using the following properties of the multiplicative free
convolution:

  * •

If we denote by

$\displaystyle\mathrm{s}_{A}$
$\displaystyle:=\lim_{w_{A}\to\infty}\bar{t}_{A}(w_{A})\,,$ (194)

then for $\theta\in(0,\mathrm{s}_{A})$, $\tilde{\mathcal{S}}_{A}(.)$ is an
increasing function of $\theta$, see App. A.2.2.

  * •

The T-transform of the matrix $\mathbf{C}$ satisfies the following inequality
at the edge:

$\displaystyle t_{C}(\mathrm{c}_{+})$
$\displaystyle\leq\min\left(t_{A}(\mathrm{a}_{+}),t_{B}(\mathrm{b}_{+})\right)\,,$
(195)

see App. A.3.

Example (S-transform of a Wishart matrix): For a Wishart matrix, from its
Stieltjes transform given by Eq. (33), one can easily get the T-transform
defined by Eq. (192) and invert it to get the S-transform of Eq. (191):

$\displaystyle\tilde{\mathcal{S}}_{\mathrm{MP}_{q}}(y)$
$\displaystyle=1+qy\,.$ (196)

###  5.2 The multiplicative spherical integral and the LSSK model

For the multiplicative case, we argue that a good candidate for the tilt is
the _multiplicative spherical integral_ $\mathcal{Z}_{\mathbf{C}}(\theta)$ ,
where for _any_ positive semi-definite matrix $\mathbf{M}$ and $\theta\geq 0$,
this spherical integral is given by:

$\displaystyle\mathcal{Z}_{\mathbf{M}}(\theta)$
$\displaystyle:=\int_{\mathbb{S}^{N-1}}\mathrm{d}\omega(\bm{\sigma})\,(\bm{\sigma}^{\mathsf{T}}\mathbf{M}\bm{\sigma})^{\frac{N}{2}\theta}\,.$
(197)

where we recall that $\mathbb{S}^{N-1}$ is the hyper-sphere of radius one and
$\mathrm{d}\omega(.)$ is the unit uniform measure over $\mathbb{S}^{N-1}$. By
Haar property, for any two positive semi-definite matrices $\mathbf{M}_{1}$
and $\mathbf{M}_{2}$, this function satisfies:

$\displaystyle\int_{\mathsf{O}(N)}\mathcal{Z}_{\sqrt{\mathbf{M}_{1}}\mathbf{O}\mathbf{M}_{2}\mathbf{O}^{\mathsf{T}}\sqrt{\mathbf{M}_{1}}}(\theta)\mathrm{d}\mathbf{O}$
$\displaystyle=\mathcal{Z}_{\mathbf{M}_{1}}(\theta)\,\mathcal{Z}_{\mathbf{M}_{2}}(\theta)\,,$
(198)

which gives for $\mathbf{M}_{1}=\mathbf{A}$ and $\mathbf{M}_{2}=\mathbf{B}$
taken from an invariant ensemble and after averaging, the desired
decomposition property of Eq. (95).

We can again interpret $\mathcal{Z}_{\mathbf{M}}(\theta)$ as the partition
function of a spherical model:

$\displaystyle\mathcal{Z}_{\mathbf{M}}(\theta)$
$\displaystyle:=\langle\mathrm{e}^{\frac{N}{2}\theta\mathcal{H}^{\mathrm{LSSK}}(\bm{\sigma})}\rangle\,,$
(199)

with the Hamiltonian:

$\displaystyle\mathcal{H}^{\mathrm{LSSK}}(\bm{\sigma})$
$\displaystyle:=\log\left(\sum_{i=1}^{N}\lambda_{i}(\mathbf{M})\sigma_{i}^{2}\right)\,.$
(200)

Since $\mathbf{M}$ is the symmetric product of definite positive matrices, its
eigenvalues are positive, so this Hamiltonian is well defined. In this paper,
we are interested in the case $\mathbf{M}=\mathbf{C}$ given by Eq. (187). Due
to the logarithmic term, we denote this model as the _Logarithmic Spherical
Sherrington-Kirkpatrick_ (LSSK for short) model. The behavior of
$\mathcal{Z}_{\mathbf{C}}(\theta)$ for $N$ large has been recently
investigated in Refs. [38, 39] and are given in the following section. Due to
its similarity with the original SSK model, one should expect to have a
similar behavior, with a paramagnetic phase at high temperature and a spin
glass phase at low temperature.

###  5.3 Asymptotic behavior of the quenched and annealed free energies of the
LSSK model

We summarize here the asymptotic behavior of (the derivatives of) the quenched
and annealed free energies.

For $\mathbf{C}$ conditioned to have its largest eigenvalue fixed at the
position $x$, the partial derivatives of the quenched free energy of Eq. (87)
with $\mathcal{Z}_{\mathbf{C}}(\theta)$ given by Eq. (197) satisfy a phase
transition. In the high temperature regime, it has been shown in Refs. [38,
39] that the derivative with respect to the parameter $\theta$ of the quenched
free energy is for small enough $\theta$ the logarithm of the S-transform. One
can also get the behavior for low temperature where there is a saturation, and
we have:

$\displaystyle\partial_{\theta}J_{C}(x,\theta)=\frac{1}{2}\left\\{\begin{array}[]{ll}\log\tilde{\mathcal{S}}_{C}(\theta)&\mbox{for
}\theta\leq t_{C}(x)\,,\\\ \\\
\log\left(\frac{x\theta}{\theta+1}\right)&\mbox{for }\theta\geq
t_{C}(x)\,.\end{array}\right.$ (204)

Similarly, the partial derivative with respect to $x$ is given by:

$\displaystyle\partial_{x}J_{C}(x,\theta)=\frac{1}{2}\left\\{\begin{array}[]{ll}0&\mbox{for
}\theta\leq t_{C}(x)\,,\\\ \\\ \frac{\theta+1}{x}-g_{C}(x)&\mbox{for
}\theta\geq t_{C}(x)\,.\end{array}\right.$ (208)

For the derivative of the annealed free energy $F_{C}(.)$, one only needs the
of $F_{A}(w_{A},.)$ and $F_{B}(w_{B},.)$ separately since we have the
decomposition of Eq. (96), which are given by (see App. B.2 for a derivation):

$\displaystyle\partial_{\theta}F_{A}(w_{A},\theta)=\frac{1}{2}\left\\{\begin{array}[]{ll}\log\tilde{\mathcal{S}}_{A}(\theta)&\mbox{for
}\theta\leq\bar{t}_{A}(w_{A})\,,\\\ \\\
\log\left(\frac{w_{A}\theta}{\theta+1}\right)&\mbox{for
}\theta\geq\bar{t}_{A}(w_{A})\,,\end{array}\right.$ (212)

and similarly for $F_{B}(w_{B},.)$. Note that if the wall is at the edge
$w_{A}=\mathrm{a}_{+}$, since
$t_{A}(\mathrm{a}_{+})=\bar{t}_{A}(\mathrm{a}_{+})$, we have Eq. (101) and
hence form the large deviation of the top eigenvalue, we can indeed consider
fixed diagonal matrices as invariant matrices with a wall at the edge.
Conversely, classical invariant ensembles correspond to the limit
$w_{A}\to\infty$ from which we see that the annealed free energy is equal to
the logarithm of the S-transform for $\theta<\mathrm{s}_{A}$, with
$\mathrm{s}_{A}$ given by the limit of (194), and is otherwise infinite.

###  5.4 Optimal temperature for the product

Without any loss of generality, let’s assume

$\displaystyle\bar{t}_{A}(w_{A})$ $\displaystyle\leq\bar{t}_{B}(w_{B})\,,$
(213)

and

$\displaystyle\bar{t}_{C}(\mathrm{c}_{+})$ $\displaystyle\leq\infty\,.$ (214)

Our goal is to show that the supremum in Eq. (91) is attained at a unique
point by looking at the derivative $I_{x}^{\prime}(.)$ with respect to
$\theta$. Paying attention to the bounds in Eqs. (204) and (212), one has the
following behavior:

$\displaystyle
I_{x}^{\prime}(\theta)=\frac{1}{2}\left\\{\begin{array}[]{llll}0&\mbox{for
}\theta\leq t_{C}(x)\,,\\\ \\\
\log\dfrac{x}{t_{C}^{\langle-1\rangle}(\theta)}&\mbox{for
}t_{C}(x)\leq\theta\leq\bar{t}_{A}(w_{A})\,,\\\ \\\
\log\dfrac{x}{w_{A}\tilde{\mathcal{S}}_{B}(\theta)}&\mbox{for
}\bar{t}_{A}(w_{A})\leq\theta\leq\bar{t}_{B}(w_{B})\,,\\\ \\\
\log\dfrac{x(\theta+1)}{w_{A}w_{B}\theta}&\mbox{for
}\theta\geq\bar{t}_{B}(w_{B})\,.\end{array}\right.$ (222)

Based on similar monotonous argument as in the additive case of Sec. 4.5 , one
can show that for $\theta\geq t_{C}(x)$ this function is continuously
increasing until it reaches the point $t_{C}(\mathrm{c}_{+})$ and then it is
continuously decreasing. For values of $x>w_{A}\,w_{B}$, it never crosses the
real axis, and we have $\theta^{*}(x)=\infty$. Otherwise, it crosses the real
axis exactly one time and the equation determining $\theta^{*}(x)$ depends on
the position of $x$ with respect to the two critical points $x_{c_{1}}$ and
$x_{c_{2}}$ defined by

$\displaystyle x_{c_{1}}$
$\displaystyle:=t_{C}^{\langle-1\rangle}(\bar{t}_{A}(w_{A}))=w_{A}\,\tilde{\mathcal{S}}_{B}\left(\bar{t}_{A}(w_{A})\right)\,,$
(223)

and by:

$\displaystyle x_{c_{2}}$
$\displaystyle:=w_{A}w_{B}\,\frac{\bar{t}_{B}(w_{B})}{\bar{t}_{B}(w_{B})+1}\,.$
(224)

  1. 1.

for $\mathrm{c}_{+}<x<x_{c_{1}}$, $\theta^{*}$ is attained in
$(t_{C}(\mathrm{c}_{+}),\bar{t}_{A}(w_{A}))$ and so setting the second line of
the RHS of Eq. (222) being equals to zero, gives:

$\displaystyle\theta^{*}(x)$ $\displaystyle=\bar{t}_{C}(x)\,;$ (225)

  2. 2.

for $x_{c_{1}}<x<x_{c_{2}}$, the optimal inverse temperature is attained in
the region $\bar{t}_{A}(w_{A})\leq\theta\leq\bar{t}_{B}(w_{B})$ and
$\theta^{*}(x)$ is solution of the third line of the RHS of Eq.(222) being
equal to zero, that is:

$\displaystyle\theta^{*}(x)$
$\displaystyle=\tilde{\mathcal{S}}_{B}^{\langle-1\rangle}\left(\frac{x}{w_{A}}\right)\,;$
(226)

  3. 3.

for $x_{c_{2}}<x<w_{A}w_{B}$, $\theta^{*}(x)$ is attained in the region
$\theta\geq\bar{t}_{B}(w_{B})$ and so from Eq. (222) it is given by:

$\displaystyle\theta^{*}(x)$ $\displaystyle=\frac{x}{w_{A}w_{B}-x}\,.$ (227)

###  5.5 Expression for the rate function

Using the expressions of the previous section for the optimal temperature
together with the expression of the partial derivative of the quenched free
energy of Eq. (208) in Eq. (94), we have that the rate function is given by:

$\displaystyle\Psi_{C}(x)=\left\\{\begin{array}[]{lll}\displaystyle\frac{1}{2}\int_{\mathrm{c}_{+}}^{x}\left(\bar{g}_{C}(t)-g_{C}(t)\mathrm{d}t\right)&\mbox{for
}\mathrm{c}_{+}\leq x\leq x_{c_{1}}\,,\\\ \\\
K_{1}+\displaystyle\frac{1}{2}\int_{x_{c_{1}}}^{x}\left(\frac{\tilde{\mathcal{S}}_{B}^{\langle-1\rangle}\left(\frac{t}{w_{A}}\right)+1}{t}-g_{C}(t)\right)\mathrm{d}t&\mbox{for
}x_{c_{1}}\leq x\leq x_{c_{2}}\,,\\\ \\\
K_{2}+\displaystyle\frac{1}{2}\log\left(\dfrac{x}{w_{A}w_{B}-x}\right)-\frac{1}{2}\int_{x_{c_{2}}}^{x}g_{C}(t)\mathrm{d}t&\mbox{for
}x_{c_{2}}\leq x\leq w_{A}w_{B}\,.\end{array}\right.$ (233)

and is infinite for other values of $x$. The constant $K_{1}$ is again given
by Eq. (169) with now $x_{c_{1}}$ given Eq. (223) (if $x_{c_{1}}<\infty)$ and
$K_{2}$ is defined by:

$\displaystyle K_{2}$
$\displaystyle:=K_{1}+\frac{1}{2}\log\left(\frac{1}{\bar{t}_{B}(w_{B})}\right)+\frac{1}{2}\int_{x_{c_{1}}}^{x_{c_{2}}}\left(\frac{\tilde{\mathcal{S}}_{B}^{\langle-1\rangle}\left(\frac{t}{w_{A}}\right)+1}{t}-g_{C}(t)\right)\mathrm{t}s$
(234)

if $x_{c_{2}}<\infty$.

Example (Rate Function for Generalized Wishart): Let’s consider the case in
Eq. (187) where $\mathbf{B}$ is a (White) Wishart with shape ratio $q$ and
$\mathbf{A}$ is a fixed diagonal (positive semi-definite) matrix. In this
case, we have $w_{A}=\mathrm{a}_{+}$ and $w_{B}=\bar{t}_{B}(w_{B})=\infty$. As
a consequence, there is just one critical point and using the expression (191)
for the S-transform of the Wishart matrix, it is given by:

$\displaystyle x_{c_{1}}$
$\displaystyle=\mathrm{a}_{+}(1+q\,t_{A}(\mathrm{a}_{+}))\,.$ (235)

and using Eq. (196) for the S-transform, we have for the rate function:

$\displaystyle\Psi_{C}(x)=\left\\{\begin{array}[]{lll}\displaystyle\frac{1}{2}\int_{\mathrm{c}_{+}}^{x}\left(\bar{g}_{C}(t)-g_{C}(t)\right)\mathrm{d}t&\mbox{for
}\mathrm{c}_{+}\leq x\leq x_{c_{1}}\,,\\\ \\\ \displaystyle
K_{1}+\frac{1}{2}\left(-\frac{1}{q}-t_{A}(\mathrm{a}_{+})+\frac{x}{q\mathrm{a}_{+}}+\left(1-\frac{1}{q}\right)\log\left(\frac{x}{x_{c_{1}}}\right)-\int_{x_{c_{1}}}^{x}g_{C}(t)\mathrm{d}t\right)&\mbox{for
}x\geq x_{c_{1}}\,,\\\ \end{array}\right.$ (239)

which is up to a change in the notation, the results obtained in Ref. [34].

##  6 Large deviation for the top singular value of sum of rectangular
matrices

In this section, we consider the case where the matrix $\mathbf{C}$ is given
as

$\displaystyle\mathbf{C}$ $\displaystyle=\mathbf{A}+\mathbf{B}\,,$ (240)

where $\mathbf{A}\sim\mathbb{P}_{V_{A},w_{A}}$ and
$\mathbf{B}\sim\mathbb{P}_{V_{B},w_{B}}$ are two _rectangular_ matrices, each
taken from a bi-invariant ensemble with a wall as defined in Sec. 2.4.3. We
aim at computing the rate function $\Phi_{C}(x)$:131313We recall that since we
are studying the rectangular case, one needs to replace the notation
$\zeta_{i}\left(\mathbf{C}\right)$ and $\Pi_{C}(x)$ in Sec. 3 by
$s_{i}\left(\mathbf{C}\right)$ and $\Phi_{C}(x)$ respectively.

$\displaystyle\mathbb{P}\left[s_{1}(\mathbf{C})=x\right]$
$\displaystyle\approx\mathrm{exp}\left[-N\,\Phi_{C}(x)+o(N)\right]\quad\quad\text{for}\,x>\mathrm{c}_{+}\,,$
(241)

where $\mathrm{c}_{+}$ is the edge of the limiting density of singular values
$\rho_{C}$ of $\mathbf{C}$, described by the _rectangular free convolution_ ,
see the next paragraph. Note that unlike the case of the product of
rectangular matrices, this case does not boil down to consider the sum of
symmetric matrices, since:

$\displaystyle s_{i}\left(\mathbf{A}+\mathbf{B}\right)$
$\displaystyle=\sqrt{\lambda_{i}\left(\mathbf{A}\mathbf{A}^{\mathsf{T}}+\mathbf{B}\mathbf{B}^{\mathsf{T}}+\mathbf{A}\mathbf{B}^{\mathsf{T}}+\mathbf{B}\mathbf{A}^{\mathsf{T}}\right)}\neq\sqrt{\lambda_{i}\left(\mathbf{A}\mathbf{A}^{\mathsf{T}}+\mathbf{B}\mathbf{B}^{\mathsf{T}}\right)}\,.$
(242)

For specific values of $q$ (namely $q=0$, corresponding to _long matrices_ and
$q=1$, corresponding to square matrices), it is known - as we will see - that
this rectangular free convolution is related to the additive free convolution
of Sec. 4.1.

###  6.1 Rectangular free convolution

The LSVD $\rho_{C}$ of $\mathbf{C}$ of Eq. (240) is given by the _rectangular
free convolution (with shape ratio $q$)_ [58, 59], denoted by
$\rho_{A}\boxplus_{q}\rho_{B}$ in the RMT literature, which is – similarly to
the additive and multiplicative free convolution – defined via a linearizing
transform called the (modified)141414The standard convention for the
_rectangular R-transform_ , $\mathcal{R}^{(q)}_{A}$, is related to our
rectangular C-transform $\tilde{\mathcal{C}}^{(q)}_{A}$, via
$\tilde{\mathcal{C}}^{(q)}_{A}(t)=:\mathcal{R}^{(q)}_{A}(t^{2})/t$. As we will
see, the convention $\tilde{\mathcal{C}}^{(q)}_{A}$ will appear more
naturally. _rectangular C-transform (with shape ratio $q$)_:

$\displaystyle\tilde{\mathcal{C}}^{(q)}_{C}(y)$
$\displaystyle=\tilde{\mathcal{C}}^{(q)}_{A}(y)+\tilde{\mathcal{C}}^{(q)}_{B}(y)\,.$
(243)

This transform is given by the formula:

$\displaystyle\tilde{\mathcal{C}}^{(q)}_{A}(y)$
$\displaystyle:=\frac{U\left(y\,d_{A}^{\langle-1\rangle}(y)\right)}{y}\,,$
(244)

where $U(.)$ is defined by:

$\displaystyle U(y)$
$\displaystyle:=\frac{-1-q+\sqrt{(1-q)^{2}+4qy^{2}}}{2q}\,.$ (245)

Note that $U(.)$ is an increasing function of $y$ whose inverse transform is
given by the simple formula:

$\displaystyle U^{\langle-1\rangle}(z)$ $\displaystyle=\sqrt{(1+z)(1+qz)}\,.$
(246)

The function $d_{A}^{\langle-1\rangle}(.)$ in Eq. (244) is the inverse
functional of the _D-transform_ , which is the rectangular counterpart of the
Stieltjes transform (resp. T-transform) in the additive (resp. multiplicative)
case, defined by:

$\displaystyle d_{A}(z)$
$\displaystyle:=\sqrt{\left(\int\frac{z}{z^{2}-s^{2}}\rho_{A}(s)\mathrm{d}s\right)\,\left(q\int\frac{z}{z^{2}-s^{2}}\rho_{A}(s)\mathrm{d}s+\frac{1-q}{z}\right)}\,.$
(247)

Note that if we denote by $g_{AA^{\mathsf{T}}}(.)$ the Stieltjes transform of
the measure $\mu_{AA^{\mathsf{T}}}(.)=\frac{\rho_{A}(\sqrt{.})}{2\sqrt{.}}$,
the limiting Stieltjes of the matrix $\mathbf{A}\mathbf{A}^{\mathsf{T}}$, then
we have:

$\displaystyle d_{A}(z)$
$\displaystyle=\sqrt{qz^{2}\left(g_{AA^{\mathsf{T}}}(z^{2})\right)^{2}+(1-q)g_{AA^{\mathsf{T}}}(z^{2})}\,.$
(248)

Remark (long matrices ($q\to 0$) and additive free convolution): In the limit
$q\to 0$, corresponding to the case of ($N\times M)$ rectangular _long_
matrices with $1\ll N\ll M$, we have for the function $U$ and the D-transform,

$\displaystyle U(y)$ $\displaystyle\underset{q\to 0}{\to}y^{2}-1\,,$ (249)

and

$\displaystyle d_{A}(z)$ $\displaystyle\underset{q\to
0}{\to}\sqrt{g_{AA^{\mathsf{T}}}(z^{2})}\,,$ (250)

and so the inverse of the D-transform is given by:

$\displaystyle d_{A}^{\langle-1\rangle}(y)$ $\displaystyle\underset{q\to
0}{\to}\sqrt{g_{AA^{\mathsf{T}}}^{\langle-1\rangle}(y^{2})}\,.$ (251)

As a consequence, the rectangular C-transform is related to the R-transform by

$\displaystyle\tilde{\mathcal{C}}^{(0)}_{A}(y)$
$\displaystyle=y\,\mathcal{R}_{AA^{\mathsf{T}}}\left(y^{2}\right)\,,$ (252)

and so by the linearizing property of the C-transform and the R-transform we
have:

$\displaystyle\tilde{\mathcal{C}}^{(0)}_{C}(y)$
$\displaystyle=y\,\mathcal{R}_{AA^{\mathsf{T}}+BB^{\mathsf{T}}}\left(y^{2}\right)\,.$
(253)

In other words for long matrices, if one is looking at the LSVD of the sum,
one can replace the symbol ’$\neq$’ in Eq. (242) by an equality! Remark
(square matrices ($q=1$) and symmetrized density): For $q=1$, corresponding to
(asymptotic) square matrices, the function $U$ is simply given by:

$\displaystyle U(y)$ $\displaystyle\underset{q\to 1}{\to}y-1\,,$ (254)

and the D-transform of Eq. (247) considerably simplifies (for
$z>\mathrm{a}_{+}$) into:

$\displaystyle d_{A}(z)$
$\displaystyle=\int\frac{z}{z^{2}-s^{2}}\rho_{A}(s)\mathrm{d}s=\int\frac{\mu_{\hat{A}}(\lambda)}{z-\lambda}\mathrm{d}\lambda=:g_{\hat{A}}(z)$
$\displaystyle(\mbox{for }q=1)\,.$ (255)

where we recall that $\mu_{\hat{A}}$ is the symmetrized density of $\rho_{A}$,
given by Eq. (69). The C-transform of Eq. (244) reads in this case:

$\displaystyle\tilde{\mathcal{C}}^{(1)}_{A}(y)$
$\displaystyle=\mathcal{R}_{\hat{A}}(y)\,,$ (256)

and by linearizing property, we have that the LSVD of the matrix $\mathbf{C}$
is given as the unique probability measure on $\mathbb{R}_{+}$ such that:

$\displaystyle\tilde{\mathcal{C}}^{(1)}_{C}(y)$
$\displaystyle=\mathcal{R}_{\hat{A}+\hat{B}}(y)\,.$ (257)

In other words, the singular values of the sum of two (bi-free) square
matrices is given asymptotically by the additive free convolution of Sec. 4.1
of their respective symmetrized singular value densities. Example (Gaussian
rectangular random matrices): Let’s consider the case of Gaussian rectangular
matrices with LSVD given Eq. (56). Using Eq. (248) with the expression of Eq.
(196) for the Stieltjes transform of the Marčenko-Pastur distribution, one
gets the following expression for the D-transform (for $z>0$):

$\displaystyle d_{A}(z)$
$\displaystyle=\frac{1}{\sigma^{2}}\sqrt{\frac{z^{2}-(1+q)\sigma^{2}-\sqrt{z^{4}-2(1+q)\sigma^{2}z^{2}+(1-q)^{2}\sigma^{4}}}{2q}}\,,$
(258)

whose inverse is given by

$\displaystyle d_{A}^{\langle-1\rangle}(y)$
$\displaystyle=\frac{\sqrt{(1+\sigma^{2}y^{2})(1+q\sigma^{2}y^{2})}}{y}\,.$
(259)

The argument inside the square-root function is nothing else than the inverse
$U^{\langle-1\rangle}$ evaluated at $\sigma^{2}y^{2}$, see Eq. (246). Using
Eq. (244), the rectangular C-transform of the Gaussian rectangular matrix is
given simply by:

$\displaystyle\tilde{\mathcal{C}}^{(q)}_{A}(y)$ $\displaystyle=\sigma^{2}y\,.$
(260)

In the following, we will be using the two following properties of the
rectangular free convolution:

  * •

we argue that the function $\tilde{\mathcal{C}}^{(q)}_{A}(y)$ is a continuous
increasing function see App. A.2

  * •

The D-transform at the edge of $\rho_{C}$ satisfies the following inequality:

$\displaystyle d_{C}(\mathrm{c}_{+})$
$\displaystyle\leq\min\left(d_{A}(\mathrm{a}_{+}),d_{B}(\mathrm{a}_{+})\right)\,$
(261)

see App. A.3.

###  6.2 Rectangular spherical integral and the BSSK model

Similar to the cases of the sum of symmetric matrices, we choose the tilt
function for the problem of the sum of rectangular matrices to be given by the
_rectangular spherical integral_ $\mathcal{Z}_{\mathbf{C}}(\theta)$ defined
for _any_ rectangular matrix $\mathbf{M}$ and $\theta\geq 0$ by:

$\displaystyle\mathcal{Z}_{\mathbf{M}}(\theta)$
$\displaystyle:=\int_{\mathbb{S}^{N-1}}\mathrm{d}w(\bm{\sigma}_{1})\int_{\mathbb{S}^{M-1}}\mathrm{d}w(\bm{\sigma}_{2})\,\mathrm{exp}\left[\sqrt{NM}\theta\bm{\sigma}^{\mathsf{T}}_{1}\mathbf{M}\bm{\sigma}_{2}\right]\,.$
(262)

For any rectangular matrices $\mathbf{M}_{1}$ and $\mathbf{M}_{2}$, this
function satisfies the following property:

$\displaystyle\int_{\mathsf{O}(N)}\int_{\mathsf{O}(M)}\mathcal{Z}_{\mathbf{M}_{1}+\mathbf{U}\mathbf{M}_{2}\mathbf{V}^{\mathsf{T}}}(\theta)\mathrm{d}\mathbf{U}\mathrm{d}\mathbf{V}$
$\displaystyle=\mathcal{Z}_{\mathbf{M}_{1}}(\theta)\,\mathcal{Z}_{\mathbf{M}_{2}}(\theta)\,,$
(263)

which evaluated for $\mathbf{M}_{1}=\mathbf{A}$ and
$\mathbf{M}_{2}=\mathbf{B}$ given as rectangular bi-invariant random matrices,
gives after integration over the laws of $\mathbf{A}$ and $\mathbf{B}$, the
decomposition property of Eq. (95). This rectangular spherical can again be
understood as the partition function of a spherical model with inverse
temperature $\theta$. Indeed, if we denote by

$\displaystyle\bm{\sigma}$
$\displaystyle:=\frac{1}{\sqrt{2}}\left[\bm{\sigma}_{1},\bm{\sigma}_{2}\right]\in\mathbb{S}^{N+M-1}\,,$
(264)

$\mathcal{Z}_{\mathbf{C}}(\theta)$ can be written as:

$\displaystyle\mathcal{Z}_{\mathbf{M}}(\theta)$
$\displaystyle:=\langle\mathrm{e}^{\sqrt{NM}\theta\mathcal{H}^{\mathrm{BSSK}}(\bm{\sigma})}\rangle\,,$
(265)

with the Hamiltonian

$\displaystyle\mathcal{H}^{\mathrm{BSSK}}(\theta)$
$\displaystyle:=\bm{\sigma}^{\mathsf{T}}\begin{pmatrix}\mathbf{0}&\mathbf{M}\\\
\mathbf{M}^{\mathsf{T}}&\mathbf{0}\end{pmatrix}\bm{\sigma}\,,$ (266)

and by bi-invariance, this can be also written as:

$\displaystyle\mathcal{H}^{\mathrm{BSSK}}(\theta)$
$\displaystyle=\sum_{i=1}^{N}s_{i}(\mathbf{M})\sigma_{i}\sigma_{N+i}.$ (267)

This model is known [60, 61, 62] as the _Bipartite Spherical Sherrington-
Kirkpatrick_ (BSSK in short) spin model, due to the graph structure of the
interaction matrix: each coordinate of one family vector interacts only with
members of the other family.

###  6.3 Asymptotic behavior of the annealed and quenched free energy of the
BSSK model

To have the rate function, we first need to compute the derivatives of the
quenched and annealed free energy of the BSSK model of Eq. (87) and Eq. (88)
respectively with the tilt function $\mathcal{Z}_{\mathbf{C}}$ given by Eq.
(262). For $\mathbf{C}$ conditioned to have its largest singular value at the
position $x$, the quenched annealed free energy is known to be related to the
rectangular C-transform for small enough value of the inverse temperature
$\theta$. For higher value, it depends also on the position $x$, see the
derivation in the Appendix, and we have in full generality:

$\displaystyle\partial_{\theta}J_{C}(x,\theta)=\left\\{\begin{array}[]{ll}\tilde{\mathcal{C}}^{(q)}_{C}(\theta)&\mbox{for
}\theta\leq d_{C}(x)\,,\\\ \\\ \frac{U(\theta x)}{\theta}&\mbox{for
}\theta\geq d_{C}(x)\,.\end{array}\right.$ (271)

Similarly, the partial derivative with respect to $x$ is given by:

$\displaystyle\partial_{x}J_{C}(x,\theta)=\left\\{\begin{array}[]{ll}0&\mbox{for
}\theta\leq d_{C}(x)\,,\\\ \\\
\frac{\sqrt{(1-q)^{2}+4q\theta^{2}x^{2}}-\sqrt{(1-q)^{2}+4qd_{A}(x)^{2}x^{2}}}{2qx}&\mbox{for
}\theta\geq t_{C}(x)\,.\end{array}\right.$ (275)

The derivative of the annealed free energy $F_{C}$ is the sum of the
derivative of the annealed free energy $F_{A}$ and $F_{B}$ which are given by:

$\displaystyle\partial_{\theta}F_{A}(w_{A},\theta)=\left\\{\begin{array}[]{ll}\tilde{\mathcal{C}}_{A}^{(q)}(\theta)&\mbox{for
}\theta\leq\bar{d}_{A}(w_{A})\,,\\\ \\\ \frac{U(\theta
w_{A})}{\theta}&\mbox{for }\theta\geq\bar{d}_{A}(w_{A})\,,\end{array}\right.$
(279)

and similarly for $F_{B}$. If the wall is at the edge $w_{A}=\mathrm{a}_{+}$,
since $d_{A}(\mathrm{a}_{+})=\bar{d}_{A}(\mathrm{a}_{+})$, we have Eq. (101)
and hence form the large deviation of the top singular value, we can consider
fixed rectangular diagonal matrices as bi-invariant matrices with a wall at
the edge.

###  6.4 Optimal temperature for the rectangular case

Without loss of generality, we assume:

$\displaystyle\bar{d}_{A}(w_{A})$ $\displaystyle\leq\bar{d}_{B}(w_{B})\,,$
(280)

and the non-trivial condition:

$\displaystyle d_{C}(\mathrm{c}_{+})$ $\displaystyle<\infty\,.$ (281)

In this case, the derivative with respect to $\theta$ of the function
$I_{x}(.)$ in the supremum of Eq. (91) is given by:

$\displaystyle I_{x}^{\prime}(\theta)=\left\\{\begin{array}[]{llll}0&\mbox{for
}\theta\leq d_{C}(x)\,,\\\ \\\ \frac{U(\theta
x)}{\theta}-\tilde{\mathcal{C}}^{(q)}_{C}(\theta)&\mbox{for
}d_{C}(x)\leq\theta\leq\bar{d}_{A}(w_{A})\,,\\\ \\\ \frac{U(\theta
x)}{\theta}-\frac{U(\theta
w_{A})}{\theta}-\tilde{\mathcal{C}}^{(q)}_{B}(\theta)&\mbox{for
}\bar{d}_{A}(w_{A})\leq\theta\leq\bar{d}_{B}(w_{B})\,,\\\ \\\ \frac{U(\theta
x)}{\theta}-\frac{U(\theta w_{A})}{\theta}-\frac{U(\theta
w_{B})}{\theta}&\mbox{for }\theta\geq\bar{d}_{B}(w_{B})\,.\end{array}\right.$
(289)

By property of the rectangular free convolution, for $\theta>d_{C}(x)$, this
function is first increasing with $\theta$ until it reaches the value
$\theta=d_{C}(\mathrm{c}_{+})$ and then it is decreasing with $\theta$. One
can check that it crosses the real axis if $x$ is in the interval
$[\mathrm{c}_{+},w_{A}+w_{B}]$, and in this case, the position of the optimal
inverse temperature $\theta^{*}$ depends on two critical points $x_{c_{1}}$
and $x_{c_{2}}$. The first one is given by

$\displaystyle x_{c_{1}}$
$\displaystyle:=d_{C}^{\langle-1\rangle}(\bar{d}_{A}(w_{A}))\,.$ (290)

Unfortunately, unlike the sum and the product of symmetric matrices, the
expression for $x_{c_{1}}$ in terms of the rectangular C-transform is quite
involved. if we introduce the function:

$\displaystyle f_{q}(z)$
$\displaystyle:=\frac{1}{2}\sqrt{(1-q)^{2}+4qz^{2}}\,,$ (291)

to ease the notation, then we have:

$\displaystyle x_{c_{1}}$
$\displaystyle=\displaystyle\sqrt{w_{A}^{2}+\tilde{\mathcal{C}}_{B}^{(q)}(\bar{d}_{A}(w_{A}))\left(q\tilde{\mathcal{C}}_{B}^{(q)}(\bar{d}_{A}(w_{A}))+\dfrac{f_{q}(w_{A}\bar{d}_{A}(w_{A}))}{\bar{d}_{A}(w_{A})}\right)}\,.$
(292)

The other critical point $x_{c_{2}}$ is given by:

$\displaystyle x_{c_{2}}$
$\displaystyle:=\frac{1}{\bar{d}_{B}(w_{B})}U^{\langle-1\rangle}\left(U(w_{A}\bar{d}_{B}(w_{B}))+U(w_{B}\bar{d}_{B}(w_{B}))\right)\,,$
(293)

where we recall that $U^{\langle-1\rangle}$ is given by Eq. (246). Eq. (293)
can be written in semi-explicit form with the function $f_{q}$ of Eq. (291):

$\displaystyle x_{c_{2}}$
$\displaystyle=\frac{\sqrt{\left(f_{q}\left(w_{A}\bar{d}_{B}(w_{B})\right)+f_{q}\left(w_{B}\bar{d}_{B}(w_{B})\right)-1\right)\left(f_{q}\left(w_{A}\bar{d}_{B}(w_{B})\right)+f_{q}\left(w_{B}\bar{d}_{B}(w_{B})\right)-q\right)}}{\sqrt{q}\bar{d}_{B}(w_{B})}\,.$
(294)

  1. 1.

for $\mathrm{c}_{+}\leq x\leq x_{c_{1}}$, $\theta^{*}$ is attained in
$(d_{C}(\mathrm{c}_{+}),\bar{d}_{A}(w_{A}))$ and hence $\theta^{*}$ is
solution of the second line of the RHS of Eq. (289) being equals to zero.
Using Eq. (244) to express the C-transform in terms of the function $U$, one
gets:

$\displaystyle U(\theta^{*}x)$
$\displaystyle=U(d^{\langle-1\rangle}_{C}(\theta^{*})\theta^{*})\,;$ (295)

and so by applying the (monotonous) function $U^{\langle-1\rangle}(.)$ to this
equation and dividing by $\theta^{*}$, one gets to solve the equation:

$\displaystyle d^{\langle-1\rangle}_{C}(\theta^{*})$ $\displaystyle=x$
$\displaystyle\left(\mbox{for
}\theta^{*}\in(d_{C}(\mathrm{c}_{+}),\bar{d}_{A}(w_{A}))\right)\,.$ (296)

The solution is given by the second branch of the D-transform:

$\displaystyle\theta^{*}(x)$ $\displaystyle=\bar{d}_{C}(x)\,.$ (297)

  2. 2.

for $x_{c_{1}}\leq x\leq x_{c_{2}}$, the optimal inverse temperature is
attained in the region $\bar{d}_{A}(w_{A})\leq\theta\leq\bar{d}_{B}(w_{B})$.
$\theta^{*}\equiv\theta^{*}(x)$ is solution of the third line of the RHS of
Eq.(289) being equal to zero so that it satisfies:

$\displaystyle U(x\theta^{*})$
$\displaystyle:=U(w_{A}\theta^{*})+\tilde{\mathcal{C}}^{(q)}_{B}(\theta^{*})\theta^{*}\,;$
(298)

If one applies the function $U^{\langle-1\rangle}$ on each side, one gets
after simplification an analytical expression for the function
$\theta^{*}\mapsto x(\theta^{*})$:

$\displaystyle x(\theta^{*})$
$\displaystyle=\sqrt{w_{A}^{2}+\tilde{\mathcal{C}}_{B}^{(q)}(\theta^{*})\left(q\tilde{\mathcal{C}}_{B}^{(q)}(\theta^{*})+\frac{f_{q}\left(w_{A}\theta^{*}\right)}{\theta^{*}}\right)}\,,$
(299)

which is by definition the inverse of the function $\theta^{*}(x)$. Now
unfortunately for general values of the parameter $q$, we do not have a simple
analytical formula for the optimal inverse temperature function of the
position $x$ and so we simply denote by $F_{1}$ the solution of Eq. (299) with
unknown $\theta^{*}$ for $x$ between $x_{c_{1}}$ and $x_{c_{2}}$.

  3. 3.

for $x_{c_{2}}\leq x\leq w_{A}+w_{B}$, $\theta^{*}(x)$ is attained in the