∎

11institutetext: G. Barbarino 22institutetext: Department of Mathematics and
Systems Analysis, Aalto University, Finland. 22email:
<EMAIL_ADDRESS>33institutetext: V. Noferini 44institutetext:
Department of Mathematics and Systems Analysis, Aalto University, Finland.
44email<EMAIL_ADDRESS>
# The limit empirical spectral distribution of Gaussian monic complex matrix
polynomials

Giovanni Barbarino  Vanni Noferini

###### Abstract

We define the empirical spectral distribution (ESD) of a random matrix
polynomial with invertible leading coefficient, and we study it for complex
$n\times n$ Gaussian monic matrix polynomials of degree $k$. We obtain exact
formulae for the almost sure limit of the ESD in two distinct scenarios: (1)
$n\rightarrow\infty$ with $k$ constant and (2) $k\rightarrow\infty$ with $n$
constant. The main tool for our approach is the replacement principle by Tao,
Vu and Krishnapur. Along the way, we also develop some auxiliary results of
potential independent interest: we slightly extend a result by Bürgisser and
Cucker on the tail bound for the norm of the pseudoinverse of a non-zero mean
matrix, and we obtain several estimates on the singular values of certain
structured random matrices.

###### Keywords:

random matrix polynomial empirical spectral distribution polynomial eigenvalue
problem strong circle law companion matrix

###### MSC:

60B20 15B52 15A15 65F15 65F20

##  1 Introduction

Given $n\times n$ matrices $C_{0},C_{1},\dots,C_{k}\in\mathbb{C}^{n\times n}$,
consider the square matrix polynomial of degree $k$

$P(x)=\sum_{j=0}^{k}C_{j}x^{j};$ (1)

a finite eigenvalue of $P(x)$ is then defined Dopico2018 ; LN20  as a number
$\lambda\in\mathbb{C}$ such that

$\mathrm{rank}_{\mathbb{C}}P(\lambda)<\mathrm{rank}_{\mathbb{C}(x)}P(x).$

The polynomial eigenvalue problem (PEP) is to find all such eigenvalues AT12 ;
DLPVD18 ; Dopico2018 ; GLRBook ; GT17 ; LN20 ; NP15 ; TM01 , possibly (and
depending on the application) together with other objects – such as
eigenspaces, infinite eigenvalues, minimal indices and minimal bases – whose
precise definition is not relevant for this paper. Under the generic
assumption that $\det P(x)\not\equiv 0$, the finite eigenvalues of $P(x)$ are
the roots of its determinant. Polynomial eigenvalue problems are common in
several areas of applied and computational mathematics; their applications
include acoustics, control theory, fluid mechanics, and structural engineering
GLRBook ; GT17 ; TM01 .

Clearly, two very classical mathematical problems arise as special cases of
polynomial eigenvalue problems: finding the roots of a scalar polynomial
corresponds to $n=1$, while finding the eigenvalues of a matrix corresponds to
$k=1$ and $C_{1}=I_{n}$. When randomness enters the game, these two extremes
are well understood. It is known that, when the polynomial coefficients are
i.i.d. normally distributed random variables and in the limit
$k\rightarrow\infty$, then the roots of scalar polynomials are uniformly
distributed on the unit circle. Similarly, classical results in random matrix
theory state that, when the entries of an $n\times n$ matrix are i.i.d.
normally distributed random variables with mean $0$ and variance $n^{-1}$, and
in the limit $n\rightarrow\infty$, then the eigenvalues are uniformly
distributed on the unit disc. Moreover, the phenomenon of universality is well
known: there exist works that, under relatively mild assumptions, extend these
results to several other distributions of coefficients or entries.

To our knowledge, nothing was so far explicitly known about the eigenvalues of
random matrix polynomials, except for the two extremal cases above described.
In this paper, we fill this gap by computing the empirical eigenvalue
distribution of monic ($C_{k}=I)$ square matrix polynomials of size $n$ and
degree $k$, with all but the leading coefficient being i.i.d. complex Gaussian
random matrices, in two different limits: when $n\rightarrow\infty$ with $k$
constant, and when $k\rightarrow\infty$ with $n$ constant. Moreover, our
results can equivalently be interpreted as results on the empirical eigenvalue
distribution of certain structured random matrices: indeed, given a monic
matrix polynomials $P(x)$, a _linearization_ of $P(x)$ is a matrix whose
eigenvalues (as well as their geometric and algebraic multiplicities) coincide
with those of $P(x)$. In the numerical linear algebra literature, numerous
constructions of linearizations are known: see, e.g., DLPVD18 ; GLRBook ; NP15
and the references therein. In particular, the prototype of all linearizations
is the so-called companion matrix, which plays a central role in this paper.

In previous research on matrix polynomials, probability theory was used in the
context of analysing the condition number of PEPs. Namely, in AB19  Armentano
and Beltrán computed the average eigenvalue condition number for Gaussian
random complex matrix polynomials; and in BK19 , Beltrán and Kozhasov extended
the analysis to the case of real Gaussian matrix polynomials. In LN20 , Lotz
and Noferini went beyond the classical idea of condition by imposing a uniform
probability distribution on the sphere for perturbations of a fixed singular
matrix polynomial. However, we are not aware of any previous work where the
exact distribution of the eigenvalues of a random matrix polynomial is
obtained. In addition to being interesting _per se_ , our results can
potentially be valuable to numerical analysts in the context of testing
numerical methods for the solution of the PEP. Indeed, although randomly
generated problems are expected not to be very challenging from the numerical
point of view (by the results in AB19 ; BK19 ), it is common practice to use
them as benchmark for minimal performance requirements; in published research
papers on this subject, tests on random input are in fact often included among
the numerical experiments. The analytic knowledge of the limit eigenvalue
distributions that we obtain in this article can help to predict the behaviour
of randomly generated problems: when scrutinizing a novel algorithm, if the
numerically computed eigenvalues should significantly deviate from the
expectations then this fact can raise legitimate suspicions on the accuracy of
the computations.

The structure of the paper is as follows. In Section 2 we review some
necessary background material on linear algebra, matrix polynomial theory,
probability theory, and random matrix theory. Moreover, we define the
empirical spectral distribution of a random matrix polynomial with invertible
leading coefficient. In Section 3 we obtain our first main result: the almost
sure limit, for $n\rightarrow\infty$, of the empirical spectral distribution
of a random $n\times n$ monic complex Gaussian matrix polynomial of degree
$k$. In Section 4, our second main result is discussed: the almost sure limit,
for $k\rightarrow\infty$, of the empirical spectral distribution of a random
$n\times n$ monic complex Gaussian matrix polynomial of degree $k$. In Section
5 we draw some conclusions and propose new lines of research. To keep the main
part of the paper as easily readable as possible, the proof of some technical
lemmata, needed in Sections 3 and 4, is postponed to Appendix A; however, we
believe that some of those results could have independent interest. In
particular, we slightly improve known results on the tail bounds for
pseudoinverses of random matrices with nonzero mean, and we study the extremal
singular values of certain structured random matrices.

##  2 Mathematical background

###  2.1 Linear algebra

Given an $m\times n$ complex matrix $X$, we denote it singular values by
$\sigma_{1}(X)\geq\dots\geq\sigma_{\min}(X)\geq 0$, having introduced the
shorthand $\sigma_{\min}(X):=\sigma_{\min(m,n)}(X)$. The spectral norm of $X$
is denoted by $\|X\|:=\sigma_{1}(X)$, while the Frobenius norm of $X$ is

$\|X\|_{F}=\left(\sum_{i=1}^{m}\sum_{j=1}^{n}|X_{ij}|^{2}\right)^{1/2}=\sqrt{\mathrm{tr}(X^{*}X)}=\sqrt{\sigma_{1}(X)^{2}+\dots+\sigma_{\min}(X)^{2}}.$

Recall that any $X$ admits a singular value decomposition $U\Sigma V$ with
$U\in\mathbb{C}^{m\times m}$, $V\in\mathbb{C}^{n\times n}$ unitary matrices
and $\Sigma\in\mathbb{R}^{m\times n}$ diagonal real matrix whose diagonal
elements are the singular values $\sigma_{i}(X)\geq 0$. The Moore-Penrose
pseudoinverse of $X$ is the matrix $X^{\dagger}=V^{*}\Sigma^{\dagger}U^{*}$,
where $\Sigma^{\dagger}\in\mathbb{R}^{n\times m}$ is a diagonal real matrix
whose diagonal entries are
$\Sigma_{i,i}^{\dagger}=\Sigma_{i,i}^{-1}=\sigma_{i}(X)^{-1}$ if
$\sigma_{i}(X)>0$ and zero otherwise. Note that, if $X$ has full rank, then
$\|X^{\dagger}\|=1/\sigma_{\min}(X)$. We also use the induced $1$ and $\infty$
matrix norms, defined respectively as

$\|X\|_{1}=\max_{1\leq j\leq
n}\sum_{i=1}^{m}|X_{ij}|,\qquad\|X\|_{\infty}=\max_{1\leq i\leq
m}\sum_{j=1}^{n}|X_{ij}|.$

Since $\mathbb{C}^{m\times n}$ is finite dimensional, the various norms
mentioned above are of course equivalent to each other, and the following
relations will be useful to us:

$\|X\|\leq\|X\|_{F},\qquad\|X\|\leq\sqrt{\|X\|_{\infty}\|X\|_{1}}.$ (2)

An interlacing result for the singular values arises when we consider low rank
perturbation of matrices.

######  Theorem 2.1 (Interlacing Singular Values for Low-Rank Perturbations
THO76 )

Let $A$ and $E$ be $n\times n$ matrices, where $E$ has rank at most $k$. If
$B=A+E$ and the singular values of $A$ and $B$ are, respectively,

$\alpha_{1}\geq\alpha_{2}\geq\dots\geq\alpha_{n},\qquad\beta_{1}\geq\beta_{2}\geq\dots\geq\beta_{n},$

then

$\displaystyle\alpha_{i}\geq\beta_{i+k},$ $\displaystyle i=1,2,\dots,n-k,$
$\displaystyle\beta_{i}\geq\alpha_{i+k},$ $\displaystyle i=1,2,\dots,n-k.$

If the norm of the perturbation, as opposed to its rank, is to be used to
estimate the singular values, then we can appeal to the following result
attributed to Mirsky, which is a corollary of the minimax principle for
singular values.

######  Theorem 2.2 (Perturbation Theorem MIR60 )

Given two ${n\times n}$ matrices $A,B$, with singular values, respectively,

$\alpha_{1}\geq\alpha_{2}\geq\dots\geq\alpha_{n},\qquad\beta_{1}\geq\beta_{2}\geq\dots\geq\beta_{n},$

then

$|\alpha_{i}-\beta_{i}|\leq\|A-B\|,\qquad i=1,2,\dots,n.$

###  2.2 Matrix polynomial theory

Let $P(x)$ be the matrix polynomial defined in (1). We give here a brief
overview of those aspects in the spectral theory of square complex matrix
polynomials that are relevant to this paper. More detailed discussions can be
found, e.g., in AT12 ; Dopico2018 ; GLRBook ; LN20  and the references
therein. As mentioned in the introduction, an element $\lambda\in\mathbb{C}$
is said to be a finite eigenvalue of $P(x)$ if

$\mathrm{rank}_{\mathbb{C}}(P(\lambda))<\mathrm{rank}_{\mathbb{C}(x)}(P(x))=:r,$

where $\mathbb{C}(x)$ is the field of fractions of $\mathbb{C}[x]$, that is,
the field of rational functions with coefficients in $\mathbb{C}$.

If the leading coefficient $C_{k}$ of the matrix polynomial $P(x)$ in (1) is
invertible, then $P(x)$ has $kn$ finite eigenvalues. Under this assumption,
one can define the _companion matrix of $P(x)$_ as (see e.g. AT12 )

$M=\begin{bmatrix}-C_{k}^{-1}C_{k-1}&-C_{k}^{-1}C_{k-2}&\dots&-C_{k}^{-1}C_{1}&-C_{k}^{-1}C_{0}\\\
I_{n}&0&0&\dots&0\\\ 0&I_{n}&0&\dots&0\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\
0&\dots&0&I_{n}&0\end{bmatrix}\in\mathbb{C}^{kn\times kn},$ (3)

where $I_{n}$ and $0$ are, respectively, the $n\times n$ identity and zero
matrices. It is well known that the eigenvalues of $M$, defined in the
classical sense, coincide with the finite eigenvalues of $P(x)$. As a
consequence, under the assumption that $C_{k}$ is invertible, studying the
finite eigenvalues of $P(x)$ is equivalent to studying the eigenvalues of the
structured matrix $M$. Observe that, if $P(x)$ is monic, then $C_{k}=I$ so
that the assumption is automatically satisfied.

We can identify, say via an arbitrary but fixed rearrangement of the real and
imaginary parts of the entries of each coefficient, the (real) vector space of
$n\times n$ complex matrix polynomials of degree up to $k$ with
$\mathbb{R}^{2(k+1)n^{2}}$. In this setting, let
$\mathcal{S}\subset\mathbb{R}^{2(k+1)n^{2}}$ correspond to the subset of
matrix polynomials that are regular and have $kn$ distinct finite eigenvalues.
We conclude this subsection by observing that $\mathcal{S}$ is a nonempty
Zariski open set, and hence, its complement has Lebesgue measure zero: in this
sense, being regular with $kn$ distinct finite eigenvalues is a generic
property of matrix polynomials.

###  2.3 Random Matrix Theory

Often, within our probabilistic arguments it will be crucial to consider
matrices that have some deterministic entries and some other entries
corresponding to (complex) random variables, which in turn can be seen as
pairs of real random variables. We implicitly identify those matrices with a
vector in $\mathbb{R}^{N}$, $N$ being the number of real random variables
involved, and equipping $\mathbb{R}^{N}$ with an appropriate probability
measure. In this context, we will often invoke, without explicit
justification, the well known fact that events that happen in (subsets of)
proper Zariski closed sets of $\mathbb{R}^{N}$ have probability zero: for
example, we may claim that a certain square random matrix is almost surely
invertible. Recalling that any proper Zariski closed set has Lebesgue measure
zero, it follows immediately that the claimed property is true for any
absolutely continuous probability measure (as are all the ones we discuss in
this paper). The verification that, in all the instances where we make such a
claim, the corresponding algebraic set is indeed contained in a proper
algebraic set is a straightforward exercise in linear algebra, and we
therefore omit the details.

####  2.3.1 Empirical spectral distributions

Given a deterministic matrix $A\in M_{m}(\mathbb{C})$ with eigenvalues
$\lambda_{1}(A),\dots,\lambda_{m}(A)$ we say that its empirical spectral
distribution (ESD) is the atomic measure

$\mu_{A}=\frac{1}{m}\sum_{i=1}^{m}\delta_{\lambda_{i}(A)},$

where the eigenvalues are considered with their respective algebraic
multiplicities. A random matrix $A_{n}$ can be seen as a random variable with
values in the appropriate space of matrices, that will usually be
$M_{nk}(\mathbb{C})$, where $n$ and $k$ are fixed parameters. We can extend
the concept of ESDs to random matrices as follows.

###### Definition 2.3

Given a random matrix $A$, its empirical spectral distribution (ESD) is a
random variable with values in the space of probabilities on $\mathbb{C}$,
defined as

$\mu_{A}(\omega):=\mu_{A(\omega)}=\frac{1}{m}\sum_{i=1}^{m}\delta_{\lambda_{i}(A(\omega))}.$

The space of probabilities on $\mathbb{C}$ is a measurable subset of
$\mathcal{M}^{b}(\mathbb{C})$, the space of signed measure of $\mathbb{C}$
with bounded total variation, that is a Hausdorff space when equipped with the
vague (or weak$-^{*}$) convergence of measures.

We will study the spectral distribution for some families of random matrices
$\\{A_{n}\\}_{n}$, and find that in our cases the sequence
$\\{\mu_{A_{n}}\\}_{n}$ always converges almost surely (a.s.) to a constant
random variable, that can be identified with a probability measure
$\mu\in\mathbb{P}(\mathbb{C})$. In this case, we simply write

$\mu_{A_{n}}\xrightarrow{a.s.}\mu.$

The measure $\mu$ will thus be our candidate for the asymptotic spectral
distribution of the family $\\{A_{n}\\}_{n}$.

Finally, having let us consider a random matrix polynomial $P(x;\omega)$ of
shape $n\times n$ and degree $k$, under the assumption that, for all
$\omega\in\Omega$, $P(x;\omega)$ has invertible leading coefficient. This
implies, in particular, that $P(x;\omega)$ has $kn$ finite eigenvalues, that
we denote by $\lambda_{1}(P(x;\omega)),\dots,\lambda_{kn}(P(x;\omega))$.

###### Definition 2.4

Let $P(x;\omega)$ be a random matrix polynomial of size $n$ and degree $k$,
such that its leading coefficient is invertible for all $\omega\in\Omega$. Its
empirical spectral distribution (ESD) is a random variable with values in the
space of probabilities on $\mathbb{C}$, defined as

$\mu_{P}(\omega):=\mu_{P(x;\omega)}=\frac{1}{kn}\sum_{i=1}^{kn}\delta_{\lambda_{i}(P(x;\omega))}.$

It is immediate by Definitions 2.3 and 2.4 that the ESD of a random matrix
polynomial coincides with the ESD of its (random) companion matrix (3).
Indeed, in this paper we will strongly rely on its equivalence.

####  2.3.2 The replacement principle and the circle law

Central to our arguments to derive the empirical spectral distributions is the
so-called replacement principle: a tool in random matrix theory developed by
Tao, Vu and Krishnapur. We recall it below.

######  Theorem 2.5 (Replacement Principle TVK10 )

Let $A_{m},B_{m}$ be two $m\times m$ random matrices. Assume that

  1. 1.

The quantity $\frac{1}{m^{2}}\left(\|A_{m}\|_{F}^{2}+\|B_{m}\|^{2}_{F}\right)$
is bounded a.s.;

  2. 2.

For a.e. $z\in\mathbb{C}$,

$\frac{1}{m}\log\left|\frac{\det(m^{-1/2}A_{m}-zI)}{\det(m^{-1/2}B_{m}-zI)}\right|\xrightarrow{a.s}0.$

Then,
$\mu_{\frac{1}{\sqrt{m}}A_{m}}-\mu_{\frac{1}{\sqrt{m}}B_{m}}\xrightarrow{a.s.}0$.

###### Remark 2.6

The random variable
$\mu_{\frac{1}{\sqrt{m}}A_{m}}-\mu_{\frac{1}{\sqrt{m}}B_{m}}$ takes values in
the space of signed measures on $\mathbb{C}$ with total variation bounded by
$2$.

Thanks to the replacement principle, we will be able to generalize a well-
known result on random Gaussian matrices to the case of monic Gaussian matrix
polynomials.

######  Theorem 2.7 (Strong Circle Law MehtaBook )

Let $A_{m}$ be the ${m\times m}$ random matrix whose entries are iid Gaussian
random variables with mean 0 and variance 1. Then the ESDs of
$\frac{1}{\sqrt{m}}A_{m}$ converges almost surely to the uniform distribution
on the unit disc.

##  3 Empirical spectral distribution for $n\times n$ monic complex Gaussian
matrix polynomials of degree $k$, in the limit $n\rightarrow\infty$

Let $X$ be a complex random variable, normally distributed with mean $0$ and
variance $1$. We consider the $n\times n$ monic matrix polynomial of degree
$k\geq 2$

$P_{n}(x)=I_{n}x^{k}+\sum_{j=0}^{k-1}C_{j}x^{j},$ (4)

where, for $j=0,\dots,k-1$ every coefficient $C_{j}$ is an $n\times n$ random
matrix whose entries are i.i.d. copies of $X$. Note that each $C_{j}$ depends
on $j$ and on $n$, but we omit the dependence on $n$ in the notation. It is
intended moreover that all $C_{j}$ are independent of each other for varying
$j$ and $n$.

The finite eigenvalues of $P_{n}(x)$ coincide with the eigenvalues of its
companion matrix: in particular, substituting $C_{k}=I_{n}$ in (3), we obtain

$M:=\begin{bmatrix}-C_{k-1}&\dots&-C_{1}&-C_{0}\\\ I_{n}&&&\\\ &\ddots&&\\\
&&I_{n}&\end{bmatrix}=:Z+E_{1}C^{T}$ (5)

where $E_{1}^{T}=\begin{bmatrix}I_{n}&0&\dots&0\end{bmatrix}$ and
$C^{T}=-\begin{bmatrix}C_{k-1}&\dots&C_{1}&C_{0}\end{bmatrix}$. Note that the
spectrum of the matrix $E_{1}C^{T}$ consists of the eigenvalues of the random
matrix $-C_{k-1}$, with the addition of the eigenvalue $0$, which appears with
algebraic multiplicity $n(k-1)$. As $C_{k-1}$ is a Gaussian random matrix, the
almost sure limit ESD of $n^{-1/2}C_{k-1}$ follows the circular law (Theorem
2.7), i.e., is distributed with the uniform measure on the unit disc. Hence,
the ESD of $n^{-1/2}E_{1}C^{T}$ converges almost surely, in the limit
$n\rightarrow\infty$, to $\frac{k-1}{k}{\bf 1}_{0}+\frac{1}{k}{\bf 1}_{D}$,
where ${\bf 1}_{0}$, ${\bf 1}_{D}$ denote the uniform probability measures on,
respectively, the set $\\{0\\}$ and the unit disc. Since $n^{-1/2}M$ is a
perturbation of $n^{-1/2}E_{1}C^{T}$, one can expect that the almost sure
limit ESD of $n^{-1/2}M$, and thus the almost sure limit ESD of
$P_{n}(n^{1/2}x)$ coincides with the limit ESD for $n^{-1/2}E_{1}C^{T}$.

This conjecture is also empirically confirmed by the experiments. For example,
in Figure 1, we plotted the complex eigenvalues, multiplied by $n^{-1/2}$, of
$N$ realization of the polynomial $P_{n}(x)$ for different values of the
triple $(k,n,N)$ under the constraint $knN=c$ for some positive integer $c$
(so that the number of the eigenvalues plotted is the same in every image). We
display several subfigures organized as a matrix: the degree of the polynomial
is constant on each row (namely $k=6$ for the first row and $k=4$ for the
second row), while the columns are characterized by different values of $n$,
increasing from left to right. To facilitate the visual comparison with the
above claim, we also superimpose the unit circle on each image.

Figure 1: Scatter plots of the eigenvalues of $P_{n}(x)$ for growing $n$,
multiplied by $\frac{1}{\sqrt{n}}$.

We prove the claim as Theorem 3.1 below. Its proof relies on several technical
lemmata on the behaviour of the singular values of certain matrices: in order
to improve the readability of the paper, these are collected in Appendix A.2
and Appendix A.3. Since, for $k=1$, we recover the well known limit
distribution of the eigenvalues of a Gaussian random matrix, within the proof
we tacitly assume that $k\geq 2$.

###### Theorem 3.1

Let $P_{n}(x)$ be a monic $n\times n$ complex random matrix polynomial of
degree $k$ as in (4), where the entries of each coefficient $C_{j}$ are i.i.d.
complex random variables normally distributed with mean $0$ and variance $1$.
Then, for $n\rightarrow\infty$, the empirical spectral distribution of
$P_{n}(n^{1/2}x)$ converges almost surely to

$\frac{k-1}{k}{\bf 1}_{0}+\frac{1}{k}{\bf 1}_{D},$

where ${\bf 1}_{0}$, ${\bf 1}_{D}$ denote the uniform probability measures on,
respectively, the set $\\{0\\}$ and the unit disc.

###### Proof

The strategy of the proof is to apply the Replacement Principle (Theorem 2.5)
in the special case where $m=kn$, $A_{m}=M$ and $B_{m}=E_{1}C^{T}$, where
$M,E_{1},C$ are the matrices defined in (5) and immediately below. Indeed, by
the observations above, this immediately implies the statement. Thus, we need
to verify that the two assumptions of Theorem 2.5 hold.

  1. 1.

Consider the random variable

$R_{n}=\frac{1}{k^{2}n^{2}}\sum_{i=1}^{2kn^{2}}|X_{i}|^{2},$

where $X_{i}$ are i.i.d. normally distributed complex random variables with
mean $0$ and variance $1$. The $X_{i}$ depend also on $n$, and it is intended
that all $X_{i}$ are i.i.d. for varying $i$ and $n$. Since
$\frac{1}{m^{2}}\left(\|A_{m}\|_{F}^{2}+\|B_{m}\|_{F}^{2}\right)$ has the same
distribution as $R_{n}+\frac{k-1}{k^{2}n}$, it suffices to prove that $R_{n}$
is bounded almost surely. This is tantamount to
$\mathbb{P}(\limsup_{n}R_{n}<\infty)=1$. On the other hand, by the strong law
of large numbers,

$\mathbb{P}\left(\lim_{n}\frac{1}{2kn^{2}}\sum_{i=1}^{2kn^{2}}|X_{i}|^{2}=\operatorname{\mathbb{E}}[|X_{i}|^{2}]=1\right)=1;$

it follows that
$\mathbb{P}\left(\limsup_{n}R_{n}<\infty\right)\geq\mathbb{P}\left(\limsup_{n}R_{n}=\frac{2}{k}\right)=1.$

  2. 2.

Fix a nonzero complex number $w\neq 0$. We need to verify that, for almost
every $w$,

$\frac{1}{kn}\left(\log\left|\det\left(\frac{1}{\sqrt{kn}}E_{1}C^{T}-wI\right)\right|-\log\left|\det\left(\frac{1}{\sqrt{kn}}M-wI\right)\right|\right)\xrightarrow{a.s.}0.$

Defining $z:=w\sqrt{k}$ we readily see that this is equivalent to showing

$\frac{1}{n}\sum_{i=1}^{kn}\left[\log\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)-\log\sigma_{i}\left(\frac{1}{\sqrt{n}}M-zI\right)\right]\xrightarrow{a.s.}0$
(6)

for every $z\neq 0$. Now let $0<\delta<1/2$ and set $f(n):=\lfloor
kn-n^{1-\delta}\rfloor$. Observe that, for any $n$ large enough,
$kn>f(n)>kn-n$. Rather than verifying (6) directly, we will prove a somewhat
stronger statement. Indeed, we claim that the following three facts all hold:

$\frac{1}{n}\sum_{i=f(n)+1}^{kn}\log\sigma_{i}\left(\frac{1}{\sqrt{n}}M-zI\right)\xrightarrow{a.s.}0.$
(7)
$\frac{1}{n}\sum_{i=f(n)+1}^{kn}\log\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)\xrightarrow{a.s.}0.$
(8)
$\frac{1}{n}\sum_{i=1}^{f(n)}\left[\log\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)-\log\sigma_{i}\left(\frac{1}{\sqrt{n}}M-zI\right)\right]\xrightarrow{a.s.}0.$
(9)

It is clear that (7), (8) and (9), together, imply (6). It now remains to
prove each statement separately.

     * •

_Proof of ( 7)._ By Lemma A.6 and Lemma A.8, almost surely, for all $n$
sufficiently large, the following is true:

$\displaystyle\sum_{i=f(n)+1}^{kn}\log\sigma_{i}\left(\frac{1}{\sqrt{n}}M-zI\right)$
$\displaystyle\geq\sum_{i=f(n)+1}^{kn}\log(n^{-a-2})\geq(n^{1-\delta}+1)(-a-2)\log(n),$
$\displaystyle\sum_{i=f(n)+1}^{kn}\log\sigma_{i}\left(\frac{1}{\sqrt{n}}M-zI\right)$
$\displaystyle\leq\sum_{i=f(n)+1}^{kn}\log(d)\leq(n^{1-\delta}+1)\log(d),$

where $a$ and $d$ are the positive constants appearing in Lemma A.6 and Lemma
A.8, and $d$ can be chosen greater than $1$. Hence, dividing by $n$,

$(n^{-\delta}+n^{-1})(-a-2)\log(n)\leq\frac{1}{n}\sum_{i=f(n)+1}^{kn}\log\sigma_{i}\left(\frac{1}{\sqrt{n}}M-zI\right)\leq(n^{-\delta}+n^{-1})\log(d).$
(10)

Thus, (7) follows by the sandwich rule.

     * •

_Proof of ( 8)._ By Lemma A.7 and Lemma A.8, there are positive constants
$\widetilde{a}$ and $d>1$ such that almost surely, for all $n$ sufficiently
large,

$\displaystyle\sum_{i=f(n)+1}^{kn}\log\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)$
$\displaystyle\geq\sum_{i=f(n)+1}^{kn}\log(n^{-\widetilde{a}-2})\geq(n^{1-\delta}+1)(-\widetilde{a}-2)\log(n),$
$\displaystyle\sum_{i=f(n)+1}^{kn}\log\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)$
$\displaystyle\leq\sum_{i=f(n)+1}^{kn}\log(d)\leq(n^{1-\delta}+1)\log(d).$

The latter inequalities imply

$(n^{-\delta}+n^{-1})(-\widetilde{a}-2)\log(n)\leq\frac{1}{n}\sum_{i=f(n)+1}^{kn}\log\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)\leq(n^{-\delta}+n^{-1})\log(d).$
(11)

yielding in turn (8) via the sandwich rule.

     * •

_Proof of ( 9)._ We start by the algebraic manipulation

$\frac{1}{n}\sum_{i=1}^{f(n)}\left[\log\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)-\log\sigma_{i}\left(\frac{1}{\sqrt{n}}M-zI\right)\right]=-\frac{1}{n}\sum_{i=1}^{f(n)}\left[\log\frac{\sigma_{i}\left(\frac{1}{\sqrt{n}}M-zI\right)}{\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)}\right].$

Thanks to Mirsky’s Theorem (Theorem 2.2), we know that, for every $i$,

$\left|\sigma_{i}\left(\frac{1}{\sqrt{n}}M-zI\right)-\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)\right|\leq\frac{1}{\sqrt{n}}\|M-E_{1}C^{T}\|=\frac{1}{\sqrt{n}}$

so, for $i=1,\dots,f(n)$ there exist $d_{i}$ satisfying
$|d_{i}|\leq\frac{1}{\sqrt{n}}$ and such that

$\sigma_{i}\left(\frac{1}{\sqrt{n}}M-zI\right)=\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)+d_{i}.$

Thus,

$\left|\frac{1}{n}\sum_{i=1}^{f(n)}\left[\log\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)-\log\sigma_{i}\left(\frac{1}{\sqrt{n}}M-zI\right)\right]\right|\leq\frac{1}{n}\sum_{i=1}^{f(n)}\left|\log\left(1+\frac{d_{i}}{\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)}\right)\right|$
(12)

Observe now that, using Lemma A.9, we have that, for some positive constants
$t,\varepsilon$, almost surely, for all $n$ sufficiently large and for every
$i\leq f(n)$,

$|x|:=\left|\frac{d_{i}}{\sigma_{i}(n^{-1/2}E_{1}C^{T}-zI)}\right|\leq\left|\frac{n^{-1/2}}{\sigma_{f(n)}(n^{-1/2}E_{1}C^{T}-zI)}\right|\leq
t^{-1}n^{-\varepsilon}.$

For sufficiently large $n$ (i.e. $n>t^{-1/\varepsilon}$), the right hand side
of the latter inequality is bounded above by $1$. Noting that
$|x|<1\Leftrightarrow|\log(1+x)|\leq-\log(1-|x|)$, we obtain the following
upper bound for the right hand side of (12):

$0\leq-\frac{1}{n}\sum_{i=1}^{f(n)}\log\left(1-\frac{|d_{i}|}{\sigma_{i}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)}\right)\leq-\frac{f(n)}{n}\log\left(1-\frac{n^{-1/2}}{\sigma_{f(n)}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)}\right)$

which in turn is bounded above by

$-k\log\left(1-\frac{n^{-1/2}}{\sigma_{f(n)}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)}\right).$

Invoking again Lemma A.9, we have that almost surely

$0\leq\frac{n^{-1/2}}{\sigma_{f(n)}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)}\leq
t^{-1}n^{-\varepsilon}\rightarrow 0$
$\implies-k\log\left(1-\frac{n^{-1/2}}{\sigma_{f(n)}\left(\frac{1}{\sqrt{n}}E_{1}C^{T}-zI\right)}\right)\xrightarrow{a.s.}0,$

and this concludes the proof.

###### Remark 3.2

The relations (7) and (8) still hold if the entries of $C_{i}$ are i.i.d.
copies of any centered random variable with unit variance, using slight
variations of the reported results.

##  4 Empirical spectral distribution for $n\times n$ monic complex Gaussian
matrix polynomials of degree $k$ in the limit $k\rightarrow\infty$

Consider again111The slight notational change with respect to (4) is just to
emphasize that here we will let $k\rightarrow\infty$ rather than
$n\rightarrow\infty$. the monic matrix polynomial

$P_{k}(x)=I_{n}x^{k}+\sum_{j=0}^{k-1}C_{j}x^{j},$ (13)

so that for all $j=0,\dots,k-1$ every coefficient $C_{j}$ is a $n\times n$
random matrix where all the entries are i.i.d. Gaussian complex random
variables with mean zero and variance 1. Note that each $C_{j}$ depends on $j$
and on $k$, but we omit the dependence on $k$ in the notation. It is intended
moreover that all $C_{j}$ are independent of each other for varying $j$ and
$k$.

The finite eigenvalues of $P_{k}(x)$ coincide with those of its companion
matrix $M$ as in (5). However, this time we decompose $M$ as the sum of a
deterministic circulant matrix and a random matrix with rank at most $n$

$M=B+A,\qquad B=\begin{bmatrix}&&&I_{n}\\\ I_{n}&&&\\\ &\ddots&&\\\
&&I_{n}&\end{bmatrix},\qquad
A=\begin{bmatrix}-C_{k-1}&\dots&-C_{1}&-(C_{0}+I_{n})\\\ \phantom{I_{n}}&&&\\\
&\phantom{\ddots}&&\\\ &&\phantom{I_{n}}&\end{bmatrix}=E_{1}\widehat{C}^{T},$
(14)

where $\widehat{C}^{T}=C^{T}-e_{k}^{T}\otimes I_{n}$,
$E_{1}^{T}=\begin{bmatrix}I_{n}&0&\dots&0\end{bmatrix}$ and
$C^{T}=-\begin{bmatrix}C_{k-1}&\dots&C_{1}&C_{0}\end{bmatrix}$. In particular,
$B$ is a circulant matrix DavisBook , with spectrum

$\Lambda(B)=\set{\exp(2\pi\textnormal{i}j/k)}{j=0,1,\dots,k-1}$

where each eigenvalue has multiplicity $n$. It is thus easy to see that the
almost sure limit ESD of $B$ is the uniform (singular) probability measure on
the unit circle $\bf 1_{U}$. The problem can be seen through the lens of the
theory of perturbations for Toeplitz matrices and sequences (see for example
BPZ19 ; VZ21 ), but since the perturbation $E_{1}\hat{C}^{T}$ is a rank $n$
correction to the Teoplitz matrix $B$, this case does not fall in the
classical settings, where it is required that the random perturbation is not
singular with high probability. Since the rank of the perturbation is small
when compared with the growing size $kn$ of the matrices, we may anyway expect
that the ESD of $M$ also converge almost surely to the same distribution $\bf
1_{U}$.

This claim is also empirically supported by the experiments. For example, in
Figure 2 we plot the complex eigenvalues of $N$ realization of the polynomial
$P_{k}(x)$ for different values of the triple $(k,n,N)$. The interpretation of
Figure 2 is the same as Figure 1 after swapping the roles of the matrix sizes
and the degrees: we fix $n$ on each row (to $6$ and $4$ respectively) and we
increase the degree of the polynomial on the columns. The unit circle is drawn
on top of each scatter plot to make easier the comparison with the claim
above; $N$ is always chosen so that the number of eigenvalues plotted, equal
to $knN$, is the same in every image.

Figure 2: Scatter plots of the eigenvalues of $P_{k}(x)$ for growing $k$.

In Theorem 3.1 below, we thus show that the ESD of $M$ also converge almost
surely to $\bf 1_{U}$. Note that the statement includes, as a special case
when $n=1$, the well known limit distribution of random scalar polynomials,
for which we thus provide a novel proof. For the sake of a clearer exposition,
below we focus on the major lines of thought that lead to the proof,
postponing to Appendix A.4 a more detailed analysis of some technicalities
that appear as intermediate steps.

###### Theorem 4.1

Let $P_{k}(x)$ be a monic $n\times n$ complex random matrix polynomial of
degree $k$ as in (4), where the entries of each coefficient $C_{j}$ are i.i.d.
complex random variables normally distributed with mean $0$ and variance $1$.
Then, for $k\rightarrow\infty$, the empirical spectral distribution of
$P_{k}(x)$ converges almost surely to $\bf 1_{U}$, the uniform probability
measure on the unit circumference.

###### Proof

The strategy of the proof follows very closely that of Theorem 3.1: we verify
that the two assumptions of Theorem 2.5 hold in the special case where $m=kn$,
$A_{m}=\sqrt{kn}\ M$ and $B_{m}=\sqrt{kn}\ B$, where $M,B$ are the matrices
defined in (14) and immediately below.

  1. 1.

The first item is treated analogously to Theorem 3.1, and we omit the details.

  2. 2.

Fix a nonzero complex number $z$ such that $|z|\not\in\set{0,1}$. We show
that, for every such $z$,

$\frac{1}{kn}\left(\log\left|\det\left(M-zI\right)\right|-\log\left|\det\left(B-zI\right)\right|\right)\xrightarrow{a.s.}0,$

that is equivalent to

$\frac{1}{k}\sum_{i=1}^{kn}\log(\sigma_{i}\left(M-zI\right))-\log(\sigma_{i}\left(B-zI\right))\xrightarrow{a.s.}0.$
(15)

We claim that the following facts are true:

$\frac{1}{k}\sum_{i=1}^{n}\log(\sigma_{i}\left(M-zI\right))\xrightarrow{a.s.}0,\qquad\frac{1}{k}\sum_{i=kn-n+1}^{kn}\log(\sigma_{i}\left(M-zI\right))\xrightarrow{a.s.}0.$
(16)
$\frac{1}{k}\sum_{i=1}^{n}\log(\sigma_{i}\left(B-zI\right))\xrightarrow{a.s.}0,\qquad\frac{1}{k}\sum_{i=kn-n+1}^{kn}\log(\sigma_{i}\left(B-zI\right))\xrightarrow{a.s.}0.$
(17)
$\frac{1}{k}\sum_{i=n+1}^{kn-n}\log(\sigma_{i}\left(M-zI\right))-\log(\sigma_{i}\left(B-zI\right))\xrightarrow{a.s.}0.$
(18)

It is clear that (16), (17) and (18), together, imply (15). It now remains to
prove each statement separately.

     * •

_Proof of ( 16)._ By Lemma A.10, almost surely, for all $k$ sufficiently
large, the following are true almost surely for some positive constant $r$:

$\sigma_{1}(M-zI)\leq r\sqrt{k}+1+|z|,\qquad\sigma_{n}(M-zI)\geq|1-|z||.$

These facts are enough to conclude that

$\frac{1}{k}\sum_{i=1}^{n}|\log(\sigma_{i}\left(M-zI\right))|\leq\frac{n}{k}\max\left\\{|\log(r\sqrt{k}+1+|z|)|,|\log(|1-|z||)|\right\\}\xrightarrow{a.s.}0$

and thus the first a.s. limit in (16) holds. Moreover, by Lemma A.11, almost
surely, for all $k$ sufficiently large,

$\sigma_{kn}\left(M-zI\right)\geq tk^{-2}$

for some positive constant $t$. Hence, it suffices to estimate

$\frac{1}{k}\sum_{i=kn-n+1}^{kn}|\log(\sigma_{i}\left(M-zI\right))|\leq\frac{n}{k}\max\\{|\log(\sigma_{1}(M-zI))|,|\log(\sigma_{kn}(M-zI))|\\}$
$\leq\frac{n}{k}\max\\{|\log(r\sqrt{k}+1+|z|)|,|\log(tk^{-2})|\\}\xrightarrow{a.s.}0.$

Thus, the second part of (16) also holds.

     * •

_Proof of ( 17)._ Observe that $B-zI$ is a circulant matrix, and hence, in
particular it is normal. Its spectrum is

$\Lambda(B-zI)=\set{\lambda-z}{\lambda\in\Lambda(B)}.$

Since all eigenvalues of $B$ have unitary norm, we can bound the singular
values of $B-zI$ as

$\sigma_{i}(B-zI)=|\lambda_{i}-z|,\qquad|1-|z||\leq|\lambda_{i}-z|\leq 1+|z|.$
(19)

Importantly, these bounds do not depend on $k$. As a consequence,

$\frac{n\log(|1-|z||)}{k}\leq\frac{1}{k}\sum_{i=kn-n+1}^{kn}\log(\sigma_{i}\left(B-zI\right))\leq\frac{n\log(1+|z|)}{k},$
$\frac{n\log(|1-|z||)}{k}\leq\frac{1}{k}\sum_{i=1}^{n}\log(\sigma_{i}\left(M-zI\right))\leq\frac{n\log(1+|z|)}{k},$

and (17) follows by the sandwich rule.

     * •

_Proof of ( 18)._ We start by noting that the statement is implied by

$\frac{1}{k}\sum_{i=n+1}^{kn-n}\left|\log\left(\frac{\sigma_{i}\left(M-zI\right)}{\sigma_{i}\left(B-zI\right)}\right)\right|\xrightarrow{a.s.}0.$

Assume now that $k>2$. Observe that $M-zI$ is a perturbation of rank at most
$n$ of $B-zI$. As a consequence, by Theorem 2.1, we find that

$\sigma_{i+n}(B-zI)\leq\sigma_{i}(M-zI)\leq\sigma_{i-n}(B-zI)$ (20)

for every $n<i\leq nk-n$. Thus,

$\frac{1}{k}\sum_{i=n+1}^{kn-n}\left|\log\left(\frac{\sigma_{i}\left(M-zI\right)}{\sigma_{i}\left(B-zI\right)}\right)\right|\leq\frac{1}{k}\sum_{i=n+1}^{kn-n}\max\left\\{\left|\log\left(\frac{\sigma_{i-n}\left(B-zI\right)}{\sigma_{i}\left(B-zI\right)}\right)\right|,\left|\log\left(\frac{\sigma_{i+n}\left(B-zI\right)}{\sigma_{i}\left(B-zI\right)}\right)\right|\right\\}$

The singular values of $B-zI$ are the moduli of $\lambda_{i}-z$ where
$\lambda_{i}$ are the eigenvalues of $B$. For the rest of this argument, and
for the sake of notational simplicity, let us now drop the dependence on the
argument matrix and simply refer to the $r$th singular value of $B-zI$ as
$\sigma_{r}$. Since all the eigenvalues of $B$ have multiplicity $n$, then
$\sigma_{i-n}=|\lambda_{j}-z|$, $\sigma_{i}=|\lambda_{i}-z|$ and
$\sigma_{i+n}=|\lambda_{s}-z|$, where necessarily $i,j,s$ are pairwise
distinct; specifically, $j$ and $s$ are determined by $z$ coherently with the
decreasing ordering of the singular values. We conclude that
$\sigma_{i}-\sigma_{i+n}$ and $\sigma_{i-n}-\sigma_{i}$ are both bounded above
by

$\min_{j\neq s\neq i\neq
j}\max\\{||\lambda_{i}-z|-|\lambda_{j}-z||,||\lambda_{i}-z|-|\lambda_{s}-z||\\}\leq\min_{j\neq
s\neq i\neq j}\max\\{|\lambda_{i}-\lambda_{j}|,|\lambda_{i}-\lambda_{s}|.\\}$

In particular, as $k>2$, we can choose

$\lambda_{j}=\lambda_{i}\exp(2\pi\textnormal{i}/k),\qquad\lambda_{s}=\lambda_{i}\exp(-2\pi\textnormal{i}/k),$

and hence,

$\min_{j\neq s\neq i\neq
j}\max\\{|\lambda_{i}-\lambda_{j}|,|\lambda_{i}-\lambda_{s}|\\}\leq|1-\exp(2\pi\textnormal{i}/k)|=2\sin(\pi/k).$

Therefore, for all values of $k$ large enough so that
$0<2\sin(\pi/k)<|1-|z||$, we use (19) to obtain

$\left|\log\left(\frac{\sigma_{i+n}}{\sigma_{i}}\right)\right|=-\log\left(1-\frac{\sigma_{i}-\sigma_{i+n}}{\sigma_{i}}\right)\leq-\log\left(1-2\frac{\sin(\pi/k)}{|1-|z||}\right),$
$\left|\log\left(\frac{\sigma_{i-n}}{\sigma_{i}}\right)\right|=\log\left(1+\frac{\sigma_{i-n}-\sigma_{i}}{\sigma_{i}}\right)\leq\log\left(1+2\frac{\sin(\pi/k)}{|1-|z||}\right).$

and, since $0<x<1\implies-\log(1-x)>\log(1+x)$, we conclude that

$\frac{1}{k}\sum_{i=n+1}^{kn-n}\max\left\\{\left|\log\left(\frac{\sigma_{i-n}\left(B-zI\right)}{\sigma_{i}\left(B-zI\right)}\right)\right|,\left|\log\left(\frac{\sigma_{i+n}\left(B-zI\right)}{\sigma_{i}\left(B-zI\right)}\right)\right|\right\\}\leq-\frac{kn-2n}{k}\log\left(1-2\frac{\sin(\pi/k)}{|1-|z||}\right)$

that goes to zero as $k\to\infty$, implying (18).

##  5 Conclusions

We have rigorously obtained the limit of empirical spectral distribution for
monic complex i.i.d. Gaussian matrix polynomials. To our knowledge, and in
spite of the relatively common use of random matrix polynomials in the context
of numerical experiments to test algorithms for the polynomial eigenvalue
problem, the study in the present paper is the first attempt to study
analytically the distribution of eigenvalues of a class of random matrix
polynomials.

We hope that this work may open the path to further future research on
eigenvalues of random matrix polynomials. In particular we believe that it
would be of interest to extend our results by considering, for instance,
different ways to send $k,n\rightarrow\infty$, non-monic polynomials,
coefficients restricted to be real (and/or otherwise structured) and more
general distributions of the entries.

In a forthcoming document, the authors will show further progress about
proving the convergence of non-monic non-Gaussian polynomial empirical
spectral distribution in both cases $n,k\to\infty$.

### Acknowledgement

We acknowledge the computational resources provided by the Aalto Science-IT
project. This version of the article has been accepted for publication, after
peer review, but is not the Version of Record and does not reflect post-
acceptance improvements, or any corrections. The Version of Record is
available online at: https://doi.org/10.1007/s10959-022-01163-3.

## References

  * [1] M. Al-Ammari and F. Tisseur.  Standard triples of structured matrix polynomials.  Linear Algebra Appl., 437(3):817–834, 2012.
  * [2] D. Armentano and C. Beltrán.  The polynomial eigenvalue problem is well conditioned for random inputs.  SIAM J. Matrix Anal. Appl., 40(1):175–193, 2019.
  * [3] A. Basak, E. Paquette, and O. Zeitouni.  Spectrum of random perturbations of toeplitz matrices with finite symbols.  Transactions of the American Mathematical Society, 373:1, 11 2019\.
  * [4] C. Beltrán and K. Kozhasov.  The real polynomial eigenvalue problem is well conditioned on the average.  Found. Comput. Math., 20(2):291–309, 2020.
  * [5] P. Bürgisser and F. Cucker.  Smoothed analysis of Moore-Penrose inversion.  SIAM J. Matrix Anal. Appl., 31(5):2769–2783, 2010.
  * [6] P. Bürgisser and F. Cucker.  Condition. The Geometry of Numerical Algorithms.  Grundlehren der mathematischen Wissenschaften. Springer-Verlag Berlin Heidelberg, 2013.
  * [7] P. J. Davis.  Circulant Matrices.  AMS Chelsea Publishing. Cambridge University Press, 1994.
  * [8] F. Dopico, P. W. Lawrence, J. Pérez, and P. Van Dooren.  Block Kronecker linearizations of matrix polynomials and their backward errors.  Numer. Math., 140:373–426, 2018.
  * [9] F. Dopico and V. Noferini.  Root polynomials and their role in the theory of matrix polynomials.  Linear Algebra Appl., 584:37–78, 2020.
  * [10] I. Gohberg, P. Lancaster, and L. Rodman.  Matrix Polynomials.  SIAM, 2009.  Unabridged republication of the book first published by Academic Press.
  * [11] S. Güttel and F. Tisseur.  The nonlinear eigenvalue problem.  Acta Numer., 26:1–94, 2017.
  * [12] M. Lotz and V. Noferini.  Wilkinson’s bus: Weak condition numbers, with an application to singular polynomial eigenproblems.  Found. Comput. Math., 2020.
  * [13] M. L. Mehta.  Random matrices and the statistical theory of energy levels.  Academic Press, New York, 1967.
  * [14] L. Mirsky.  Symmetry Gauge Functionsa and Unitarily Invariant Norms.  The Quarterly Journal of Mathematics, 11(1):50–59, 01 1960.
  * [15] V. Noferini and F. Poloni.  Duality of matrix pencils, Wong chains and linearizations.  Linear Algebra Appl., 471:730–767, 2015.
  * [16] T. Tao.  Topics in Random Matrix Theory.  Graduate studies in mathematics. American Mathematical Soc.
  * [17] T. Tao and V. Vu.  Random matrices: The circular law.  Commun. Contemp. Math., 10(02):261–307, 2008.
  * [18] T. Tao, V. Vu, and M. Krishnapur.  Random matrices: Universality of ESDs and the circular law.  Ann. Probab., 38(5):2023–2065, 2010.
  * [19] R. Thompson.  Principal submatrices IX: Interlacing inequalities for singular values of submatrices.  Linear Algebra and its Applications, 5(1):1 – 12, 1972.
  * [20] R. Thompson.  The behavior of eigenvalues and singular values under perturbations of restricted rank.  Linear Algebra and its Applications, 13(1):69 – 78, 1976.
  * [21] F. Tisseur and K. Meerbergen.  The quadratic eigenvalue problem.  SIAM Rev., 43(2):235–286, 2001.
  * [22] M. Vogel and O. Zeitouni.  Deterministic equivalence for noisy perturbations.  Proceedings of the American Mathematical Society, ff10.1090/proc/15499ff, 2021.
  * [23] M. A. Woodbury.  Inverting Modified Matrices.  Technical Report Memorandum Report 42, Department of Statistics, Institute for Advanced Study, Princeton University, 1950.

##  Appendix A Appendix: Technical results

In this appendix, we provide the full details on some technical steps that are
necessary for our analysis. For convenience, we have split the appendix into
various subsections, according to the specific nature of the results contained
therein.

###  A.1 Preliminaries and known results

A result we frequently use in our arguments is the following interlacing
property of the singular values of submatrices.

######  Theorem A.1 (Interlacing Singular Values for Submatrices [19])

Given any matrix $A$ and any $p\times q$ submatrix $B$,

$\|A\|\geq\|B\|,\qquad\sigma_{\min(p,q)}(A)\geq\sigma_{\min}(B).$

Moreover, we recall the useful Woodbury identity.

######  Lemma A.2 (Woodbury [23])

Let $A,B,U,V$ be complex matrices satisfying $B=A+UV$ with $A,B$ square. If
$A$ and $I+VA^{-1}U$ are both invertible, then

$B^{-1}=A^{-1}-A^{-1}U(I+VA^{-1}U)^{-1}VA^{-1}.$

When dealing with sequences of random matrices, one can estimate the
distribution of eigenvalues, singular values, or related quantities such as,
for instance, trace, determinant, norms. Here we collect some of the
estimations we use further on. We do not claim that the bounds we mention
below are the best possible ones; yet, they suffice for our purposes.

First, we provide a probabilistic upper bound for the norm of a Gaussian
random matrix.

######  Theorem A.3 ([16])

Suppose that the coefficients of a random matrix $N$ of size ${n\times n}$ are
i.i.d. copies of a normal random variable. Then there exist absolute constants
$C,c>0$ such that

$\mathbb{P}\left(\|N\|>A\sqrt{n}\right)\leq C\exp(-cAn)$

for all $A\geq C$.

A very different kind of estimate, due to Tao and Vu, is needed for the least
singular value of random matrices having nonzero mean.

######  Theorem A.4 ([17])

Let $c$, $d$ be positive constants, and let $X$ be a complex-valued random
variable with non-zero finite variance. Then there are positive constants $a$
and $b$ such that the following holds: if $N_{n}$ is the $n\times n$ random
matrix whose entries are i.i.d. copies of $X$, and $M$ is an $n\times n$
deterministic matrix with spectral norm at most $n^{c}$, then,

$\mathbb{P}\left(\|(M+N_{n})^{-1}\|\geq n^{a}\right)\leq bn^{-d}.$

Note that, in Theorem A.4, we can always choose $d>1$ so that the probability
is summable. Thus, we can use Borel-Cantelli lemma to obtain a lower bound for
the last singular value valid for all sufficiently large $n$.

A widely used distribution in the theory of probability for real random
variables is the Beta distribution on $[0,1]$, that depends on two positive
parameters $\alpha,\beta$ and it is described by its density function

$Y\sim
B(\alpha,\beta)\implies\mathbb{P}\left(Y\leq\lambda\right)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\int_{0}^{\lambda}x^{\alpha-1}(1-x)^{\beta-1}\,{{\rm
d}}x,\quad\forall\,0\leq\lambda\leq 1,$

where $\Gamma(z)=\int_{0}^{\infty}t^{z-1}e^{-t}dt$ is Euler’s gamma function.
It is known that if $z$ is a positive integer, then $\Gamma(z)=(z-1)!$.

Consider now a real random vector $X$ with $N$ components, uniformly
distributed on the unit sphere
$S^{N-1}:=\\{X\in\mathbb{R}^{N}:\|X\|_{2}=1\\}$. It is known [12, Sec. 4.1.1]
that the squared norm of the projection of $X$ onto a $k$-dimensional space is
Beta distributed with parameters $B(k/2,(N-k)/2)$. Note that a complex random
vector with $N$ complex components, uniformly distributed on the respective
complex spherical surface, can be seen as a real vector with $2N$ real
components and uniformly distributed on $S^{2N-1}$. As a consequence, the
squared norm of the projection of such a vector onto a $k$-dimensional complex
space is Beta distributed with parameters $B(k,N-k)$. In particular, if $k=1$,
then

$\mathbb{P}\left(\|P_{1}(X)\|^{2}\leq\lambda\right)=\frac{\Gamma(N)}{\Gamma(1)\Gamma(N-1)}\int_{0}^{\lambda}(1-x)^{N-2}\,{{\rm
d}}x=1-(1-\lambda)^{N-1},\quad\forall\ \lambda\in[0,1].$ (21)

###  A.2 A variation on a result by Bürgisser and Cucker: a formula for the
tail bounds of the norm of the pseudoinverse of a non-zero mean random matrix

Theorem A.5 yields a tail bound on the norm of the Moore-Penrose pseudoinverse
of a random Gaussian complex rectangular matrix _with nonzero mean_. It is a
modification of the results obtained by Bürgisser and Cucker in [5, Sec. 3]
and [6, Ch. 4], with two differences. A minor one is that we work with
complex, as opposed to real, numbers and random variables (as noted already in
[6], this extension is not at all difficult). A more significant one is that
we are interested in the limit case where $\lambda=\frac{n-1}{N}\rightarrow
1$, and we therefore state the result in such a way that it covers that case,
unlike [5, 6] that provide formulae for the regime $\lambda<1$. For these
reasons, as well as for the sake of self-containedness, we provide a full
proof (which still follows very closely the lead of [5, 6]).

###### Theorem A.5

Let $G$ be an $n\times N$ (with $N\geq n$) complex random matrix with i.i.d.
normally distributed entries with mean $0$ and variance $n^{-1}$. Suppose
$R=R_{D}+G$ where $R_{D}\in\mathbb{C}^{n\times N}$ is a deterministic matrix,
and let $R^{\dagger}$ be the Moore-Penrose pseudoinverse of $R$. Then, given
$\tau>0$,

$\mathbb{P}\left(\sigma_{n}(R)\leq\tau\right)\leq\frac{\tau^{2(N-n+1)}}{\sqrt{2\pi}}\frac{(nNe^{2})^{N-n+1}}{(N-n+1)^{2(N-n+1)+1/2}}.$

###### Proof

We know that there exists an unit vector $u\in\mathbb{C}^{n}$ such that

$\|R^{\dagger}u\|=\|R^{\dagger}\|$

and that, for almost every222The set of matrices $R$ can be parametrized by a
random vector in $\mathbb{R}^{2nN}$. The subset of matrices for which this
property fails corresponds to a subset of the proper algebraic set described
by $\mathrm{discriminant}(\det(RR^{*}-xI))=0$.  $R$, $u$ is unique up to
multiplication by a unit of $\mathbb{C}$. If $v\in\mathbb{C}^{n}$ is any unit
vector, then for some $\beta\in\mathbb{C}$ it can be decomposed as

$v=(u^{*}v)u+\beta u^{\perp}$

where $u^{\perp}\in\mathbb{C}^{n}$ is a unit vector orthogonal to $u$. On the
other hand, since $R^{\dagger}u$ is orthogonal to $R^{\dagger}u^{\perp}$, we
have

$\|R^{\dagger}v\|^{2}=|u^{*}v|^{2}\|R^{\dagger}u\|^{2}+\|\beta
R^{\dagger}u^{\perp}\|^{2}\geq|u^{*}v|^{2}\|R^{\dagger}\|^{2}$
$\implies\|R^{\dagger}v\|\geq|u^{*}v|\|R^{\dagger}\|.$

Thus, for any $s\in(0,1)$ and $t>0$ we have

$\mathbb{P}_{R,v}\left(\|R^{\dagger}v\|\geq
t\sqrt{1-s^{2}}\right)\geq\mathbb{P}_{R,v}\left(|u^{*}v|\geq\sqrt{1-s^{2}}\text{
and }\|R^{\dagger}\|\geq t\right)$

or equivalently

$\mathbb{P}_{R,v}\left(\|R^{\dagger}v\|\geq
t\sqrt{1-s^{2}}\right)\geq\mathbb{P}_{R}\left(\|R^{\dagger}\|\geq
t\right)\mathbb{P}_{R,v}\left(|u^{*}v|\geq\sqrt{1-s^{2}}\Bigm{|}\|R^{\dagger}\|\geq
t\right),$

where we choose $v$ uniformly over the unit vectors on $\mathbb{C}^{n}$.
Observe that, by unitary invariance,

$\mathbb{P}_{R,v}\left(|u^{*}v|\geq\sqrt{1-s^{2}}\Bigm{|}\|R^{\dagger}\|\geq
t\right)=\mathbb{P}_{v}\left(|e_{1}^{*}v|\geq\sqrt{1-s^{2}}\right)=\mathbb{P}_{v}\left(|v_{1}|\geq\sqrt{1-s^{2}}\right).$

The vector $v$ can be then seen as a unit real random vector with $2n$
entries, uniformly distributed on the unit sphere $S_{2n-1}$. From (21),

$\mathbb{P}_{v}\left(|v_{1}|\geq\sqrt{1-s^{2}}\right)=(n-1)\int_{1-s^{2}}^{1}(1-x)^{n-2}=(s^{2})^{n-1}.$

Therefore,

$\mathbb{P}_{R,v}\left(\|R^{\dagger}v\|\geq
t\sqrt{1-s^{2}}\right)\geq\mathbb{P}_{R}\left(\|R^{\dagger}\|\geq
t\right)s^{2n-2}$ $\implies\mathbb{P}_{R}\left(\|R^{\dagger}\|\geq
t\right)\leq\mathbb{P}_{R,v}\left(\|R^{\dagger}v\|\geq
t\sqrt{1-s^{2}}\right)s^{2-2n}.$ (22)

Note now that

$\mathbb{P}_{R}\left(\|R^{\dagger}v\|\geq
t\sqrt{1-s^{2}}\right)=\mathbb{P}_{R}\left(\|R^{\dagger}e_{1}\|\geq
t\sqrt{1-s^{2}}\right)$ (23)

since any unitary action on $R$ does not change its property to be the sum of
a Gaussian matrix with mean 0 and variance $n^{-1}I$, plus a deterministic
matrix. Moreover, $w=R^{\dagger}e_{1}$ is the first column of $R^{\dagger}$,
and from $RR^{\dagger}=I$ (which is true with probability $1$ since $R$ is
almost surely full rank), we know that $w^{*}$ is orthogonal to all the rows
of $R$ but $r^{T}=e_{1}^{T}R$. Furthermore,

$r^{T}w=1.$

Let $r_{\perp}$ be the component of $r$ orthogonal to the vector space
$\mathcal{V}$ generated by the other rows of $R$. Since $RR^{*}$ is almost
surely invertible, with probability $1$ we have
$w^{*}=e_{1}^{T}(RR^{*})^{-1}R$, so $w^{*}$ belongs to the row space of $R$.
Since it is also orthogonal to $\mathcal{V}$, we conclude that $w^{*}$ and
$r_{\perp}^{T}$ are parallel, and

$1=r^{T}w=r_{\perp}^{T}w\implies 1=\|r_{\perp}\|\cdot\|w\|.$

Let $\varphi$ be the density of the random matrix $R$, that can be split into
$\varphi=\psi\rho$ where $\psi$ is the density of $r$, and $\rho$ is the
density of the rest of the rows, say, $R_{1}$. We have

$\mathbb{P}_{R}\left(\|R^{\dagger}e_{1}\|\geq p^{-1}\right)=\int_{\|w\|\geq
p^{-1}}\varphi(R)\,{{\rm d}}R=\int_{R_{1}}\rho\int_{\|r_{\perp}\|\leq
p}\psi\,{{\rm d}}r\,{{\rm d}}R_{1}$ (24)

Let us focus on the integral over $r$. Fix $R_{1}$ as a set of $n-1$ linearly
independent (almost surely) vectors, with span $\mathcal{V}$. On the other
hand, $r_{\perp}$ is the projection of $r$ over $\mathcal{V}^{\perp}$. If we
split $r=r_{N}+r_{D}$ where $r_{N}$(the first row of $G$) is a random Gaussian
vector $N(0,n^{-1}I)$ and $r_{D}$ (the first row of $R_{D}$) is a
deterministic vector of bounded norm, then we conclude

$UU^{*}r=r_{\perp}=UU^{*}r_{N}+UU^{*}r_{D}.$

Here, $U$ is a $N\times N-n+1$ matrix such whose columns are an orthonormal
basis of $\mathcal{V}^{\perp}$. We can rewrite $U=QE$ where $Q$ is unitary and
$E^{T}=[I\,\,0]$, so that

$UU^{*}=Q\begin{pmatrix}I&0\\\ 0&0\end{pmatrix}Q^{*}.$

The vector $Q^{*}r_{N}$ is still a random Gaussian vector with mean $0$ and
variance $n^{-1}I$, and hence,

$Q^{*}r_{\perp}=\begin{pmatrix}\widetilde{r}_{N}\\\
0\end{pmatrix}+\begin{pmatrix}\widetilde{r}_{D}\\\ 0\end{pmatrix}$

where $\widetilde{r}_{N}$ is a random Gaussian vector, of length $N-n+1$, with
mean $0$ and variance $n^{-1}$, while $\|\widetilde{r}_{D}\|\leq\|r_{D}\|$.
Setting $\widetilde{r}_{\perp}=\widetilde{r}_{N}+\widetilde{r}_{D}$, then

$Q^{*}r=\begin{pmatrix}\widetilde{r}_{\perp}\\\ *\end{pmatrix}$

so

$\int_{\|r_{\perp}\|\leq p}\psi\,{{\rm
d}}r=\int_{\|\widetilde{r}_{\perp}\|\leq p}\widetilde{\psi}\,{{\rm
d}}\widetilde{r}$ (25)

where $\widetilde{\psi}$ is the distribution of $\widetilde{r}_{\perp}$, and
$\widetilde{r}$ is a generic vector of length $N-n+1$. In other words, we are
restricting to the first $N-n+1$ coordinates of $r$, up to a unitary
transformation. Now, observe that, from (24) and (25),

$\mathbb{P}_{R}\left(\|R^{\dagger}e_{1}\|\geq
p^{-1}\right)=\int_{R_{1}}\mathbb{P}_{v\sim\mathcal{N}(v_{\star},n^{-1}I)}\left(\|v\|\leq
p\right)\,\rho\,{{\rm
d}}R_{1}\leq\int_{R_{1}}\mathbb{P}_{v\sim\mathcal{N}(0,2I)}\left(\|v\|^{2}\leq
2np^{2}\right)\,\rho\,{{\rm d}}R_{1}$ (26)

where $v_{\star}$ is a generic complex vector of dimension $N-n+1$ and $v$ is
a random complex vector of the same rank. Note further that a complex normally
distributed vector $v$ of length $N-n+1$ and variance $2$ can be seen as a
real normally distributed vector $v_{\mathbb{R}}$ of length $2(N-n+1)$ and
variance $1$. Taking this viewpoint, we write

$\mathbb{P}_{v\sim\mathcal{N}(0,2I)}\left(\|v\|\leq\sqrt{2n}p\right)=\mathbb{P}_{v_{\mathbb{R}}\sim\mathcal{N}(0,I)}\left(\|v_{\mathbb{R}}\|\leq\sqrt{2n}p\right)=\frac{1}{(2\pi)^{N-n+1}}\int_{\|X\|\leq\sqrt{2n}p}e^{-\|X\|^{2}/2}\,{{\rm
d}}X$ $\leq\frac{1}{(2\pi)^{N-n+1}}\int_{\|X\|\leq\sqrt{2n}p}1\,{{\rm
d}}X=\frac{(\sqrt{2n}p)^{2(N-n+1)}\pi^{N-n+1}}{(2\pi)^{N-n+1}(N-n+1)!}=\frac{(np^{2})^{N-n+1}}{(N-n+1)!}.$

Recalling Stirling’s bound

$(N-n+1)!\geq\sqrt{2\pi(N-n+1)}\left(\frac{N-n+1}{e}\right)^{N-n+1},$

we get the estimate

$\mathbb{P}_{v\sim
N(0,2I)}\left(\|v\|\leq\sqrt{2n}p\right)\leq\frac{1}{\sqrt{2\pi}}\frac{(np^{2}e)^{N-n+1}}{(N-n+1)^{N-n+3/2}}.$

The latter upper bound does not depend on $R_{1}$, so plugging it into (26) we
obtain

$\mathbb{P}_{R}\left(\|R^{\dagger}e_{1}\|\geq
p^{-1}\right)\leq\int_{R_{1}}\frac{1}{\sqrt{2\pi}}\frac{(np^{2}e)^{N-n+1}}{(N-n+1)^{N-n+3/2}}\,\rho\,{{\rm
d}}R_{1}=\frac{1}{\sqrt{2\pi}}\frac{(np^{2}e)^{N-n+1}}{(N-n+1)^{N-n+3/2}}$
(27)

and, using (22) and (23),

$\mathbb{P}_{R}\left(\|R^{\dagger}\|\geq
t\right)\leq\frac{1}{\sqrt{2\pi}}\frac{\left(\frac{ne}{t^{2}}\right)^{N-n+1}}{(N-n+1)^{N-n+3/2}}\frac{1}{(1-s^{2})^{N-n+1}(s^{2})^{n-1}}$

for every $s\in(0,1)$ and $t>0$. We now sharpen the bound by optimizing in the
parameter $s^{2}$: to this goal, we need to find a maximum over $(0,1)$ of

$q(Y)=(1-Y)^{N-n+1}Y^{n-1}.$

If $n>1$, then a straightforward computation shows that the maximum is
achieved at

$Y=\frac{n-1}{N}=:\lambda\in(0,1)$

and since

$\lambda^{-(n-1)}=\lambda^{-N\lambda}=(\lambda^{-\frac{\lambda}{1-\lambda}})^{N(1-\lambda)}\leq
e^{N(1-\lambda)},$

we have that

$\frac{1}{q(\lambda)}=\frac{1}{(1-\lambda)^{N-n+1}\lambda^{(n-1)}}\leq\frac{e^{N(1-\lambda)}}{(1-\lambda)^{N-n+1}}=\frac{(Ne)^{N-n+1}}{(N-n+1)^{N-n+1}}.$

On the other hand, if $n=1$, then $\sup_{Y}q(Y)=1$ and

$1\leq e^{N}=\frac{(Ne)^{N-n+1}}{(N-n+1)^{N-n+1}}.$

Thus,

$\mathbb{P}_{R}\left(\|R^{\dagger}\|\geq
t\right)\leq\frac{1}{\sqrt{2\pi}t^{2(N-n+1)}}\frac{\left(ne\right)^{N-n+1}}{(N-n+1)^{N-n+3/2}}\frac{(Ne)^{N-n+1}}{(N-n+1)^{N-n+1}},$

and hence, by taking $\tau=t^{-1}$,

$\mathbb{P}\left(\sigma_{n}(R)\leq\tau\right)\leq\frac{\tau^{2(N-n+1)}}{\sqrt{2\pi}}\frac{(nNe^{2})^{N-n+1}}{(N-n+1)^{2(N-n+1)+1/2}}.$

This concludes the proof.

###  A.3 Estimates on the singular values of certain random matrices in
Section 3

In Lemma A.6 below, we obtain (in probability) a lower bound for the smallest
singular value of the matrix $n^{-1/2}M-zI$, where $M$ is defined in in (5).

###### Lemma A.6

Let $M$ be the $kn\times kn$ matrix defined as in (5) and $0\neq
z\in\mathbb{C}$. There exist constants $a,b>0$ such that, for every large
enough $n$,

$\mathbb{P}\left({\sigma_{kn}(n^{-1/2}M-zI)}<n^{-a-2}\right)\leq 2bn^{-2}$

and in particular, with probability 1,

${\sigma_{kn}(n^{-1/2}M-zI)}\geq n^{-a-2}$

for all large enough $n$.

###### Proof

Using the same notation introduced in (5), let us rewrite

$n^{-1/2}M-zI=n^{-1/2}E_{1}C^{T}+n^{-1/2}Z-zI$

and denote $N:=Z-n^{1/2}zI$. recall that the inverse of the least singular
value of an invertible square matrix $X$ is equal to the spectral norm of
$X^{-1}$. By Woodbury Lemma (Lemma A.2), we see that

$(n^{-1/2}E_{1}C^{T}+n^{-1/2}Z-zI)^{-1}=n^{1/2}(N+E_{1}C^{T})^{-1}=n^{1/2}\left[N^{-1}-N^{-1}E_{1}(I+C^{T}N^{-1}E_{1})^{-1}C^{T}N^{-1}\right]$

Here we used that $(I+C^{T}N^{-1}E_{1})$ is invertible with probability 1 and
that $N$ is invertible since $z\neq 0$. As a consequence

$\frac{1}{\sigma_{kn}(n^{-1/2}M-zI)}\leq
n^{1/2}\|N^{-1}\|\left(1+\|(I+C^{T}N^{-1}E_{1})^{-1}\|\|C^{T}\|\|N^{-1}\|\right).$

Observe now that $N=Z-n^{1/2}zI$ is a block Toeplitz matrix lower triangular
matrix. It follows that its inverse is also a block Toeplitz lower triangular
matrix, and it is easily verified that the first block column of $N^{-1}$ is

$N^{-1}=\begin{pmatrix}(\sqrt{n}z)^{-1}I&\\\ -(\sqrt{n}z)^{-2}I&\ddots\\\
\vdots&\ddots\\\ (-1)^{k-1}(\sqrt{n}z)^{-k}I&\ddots\end{pmatrix}.$

Using $\|N^{-1}\|\leq\sqrt{\|N^{-1}\|_{1}\|N^{-1}\|_{\infty}}=\|N^{-1}\|_{1}$
we can bound the norm from above with

$\|N^{-1}\|\leq\|N^{-1}\|_{1}=\sum_{i=1}^{k}(\sqrt{n}|z|)^{-i}=(\sqrt{n}|z|)^{-1}\frac{(\sqrt{n}|z|)^{-k}-1}{(\sqrt{n}|z|)^{-1}-1}=\frac{1-(\sqrt{n}|z|)^{-k}}{\sqrt{n}|z|-1}\leq\frac{2}{\sqrt{n}|z|}$

where we are assuming $n\geq 4/|z|^{2}$. Hence,

$\frac{1}{\sigma_{kn}(n^{-1/2}M-zI)}\leq\frac{2}{|z|}\left(1+\frac{2}{\sqrt{n}|z|}\|(I+C^{T}N^{-1}E_{1})^{-1}\|\|C^{T}\|\right).$

Note that $C^{T}$ is a $n\times kn$ matrix, and can be seen as a submatrix of
a $kn\times kn$ random matrix $\widetilde{C}$ where every entry is an i.i.d
copy of a Gaussian complex random variable $X$. Thanks to an interlacing
theorem for singular values (Theorem A.1), we have that
$\|C^{T}\|\leq\|\widetilde{C}\|$. On the other hand, by Theorem A.3,

$\mathbb{P}(\|C^{T}\|>r\sqrt{n})\leq\mathbb{P}(\|\widetilde{C}\|>r\sqrt{n})\leq
s\exp(-crn)$ (28)

where $c,s,r>0$ are absolute constants. As a consequence, with high
probability (at least $1-s\exp(-crn)$),

$\frac{1}{\sigma_{kn}(n^{-1/2}M-zI)}\leq\frac{2}{|z|}(1+\frac{2r}{|z|}\|(I+C^{T}N^{-1}E_{1})^{-1}\|).$
(29)

Consider now the matrix

$C^{T}N^{-1}E_{1}=\sum_{i=1}^{k}(-1)^{i}(\sqrt{n}z)^{-i}C_{d-i}.$

Clear, each of its entry is a linear combination of i.i.d. Gaussian variables
all having mean $0$, and this is still a normally distributed variable with
mean $0$. Moreover, the variance is

$\sum_{i=1}^{k}(\sqrt{n}|z|)^{-2i}=(\sqrt{n}|z|)^{-2}\frac{(\sqrt{n}|z|)^{-2k}-1}{(\sqrt{n}|z|)^{-2}-1}=\frac{1-(\sqrt{n}|z|)^{-2k}}{(\sqrt{n}|z|)^{2}-1}=:\frac{c(n)^{2}}{n}=\Theta(\frac{1}{n}).$

Hence,

$\|(I+C^{T}N^{-1}E_{1})^{-1}\|^{-1}=\sigma_{n}(I+C^{T}N^{-1}E_{1})=\frac{c(n)}{\sqrt{n}}\sigma_{n}\left(I\frac{\sqrt{n}}{c(n)}+G\right)$

where now $G$ is a matrix where all entries are i.i.d copies of a complex
Gaussian random variable $X$ having mean $0$ and variance $1$. Since for large
values of $n$, $\|I\frac{\sqrt{n}}{c(n)}\|\leq n$, we can apply Theorem A.4
and conclude that there exist positive constants $a,b$ such that

$\mathbb{P}\left(\sigma_{n}\left(I\frac{\sqrt{n}}{c(n)}+G\right)\leq
n^{-a}\right)\leq bn^{-2}$

meaning that, with high probability (at least $1-bn^{-2}$), it holds

$\|(I+C^{T}N^{-1}E_{1})^{-1}\|^{-1}\geq\frac{c(n)}{\sqrt{n}}n^{-a}\geq
n^{-a-1}$ (30)

for sufficiently large $n$. As a consequence, from (29) and (30),

$\frac{1}{\sigma_{kn}(n^{-1/2}M-zI)}\leq\frac{2}{|z|}(1+\frac{2r}{|z|}n^{a+1})$

and thus

${\sigma_{kn}(n^{-1/2}M-zI)}\geq\frac{1}{\frac{2}{|z|}(1+\frac{2r}{|z|}n^{a+1})}\geq
n^{-a-2}$

with probability at least

$1-bn^{-2}-s\exp(-crn)\geq 1-2bn^{-2}$

for any large enough values of $n$. In particular, we can conclude the proof
by invoking the Borel-Cantelli lemma.

Lemma A.7 is the analogue of A.6 when the matrix $E_{1}C^{T}$ is considered.

###### Lemma A.7

Let $E_{1}$ and $C^{T}$ be the matrices defined in (5) and immediately after
it. There exists a constant $\widetilde{a}>0$ such that, for all large enough
$n$,

$\mathbb{P}\left({\sigma_{kn}(n^{-1/2}E_{1}C^{T}-zI)}<n^{-\widetilde{a}-2}\right)\leq
2n^{-2}$

and in particular, with probability $1$,

${\sigma_{kn}(n^{-1/2}E_{1}C^{T}-zI)}\geq n^{-\widetilde{a}-2}$

for all sufficiently large $n$.

###### Proof

The proof is very similar to Lemma A.6, so we only sketch it. In this case,
set $N:=-n^{1/2}zI$ and

$(n^{-1/2}E_{1}C^{T}-zI)^{-1}=n^{1/2}\left[N^{-1}-N^{-1}E_{1}(I+C^{T}N^{-1}E_{1})^{-1}C^{T}N^{-1}\right]$

where $N$ is invertible since $z\neq 0$ and $I+C^{T}N^{-1}E_{1}$ is almost
surely invertible. We have $\|N^{-1}\|=n^{-1/2}/|z|$, and hence, with
probability at least $1-s\exp(-crn)$,

$\frac{1}{\sigma_{kn}(n^{-1/2}E_{1}C^{T}-zI)}\leq\frac{1}{|z|}(1+\frac{r}{|z|}\|(I-\frac{1}{n^{1/2}z}C_{k-1})^{-1}\|).$

Again, we write

$\|(I-\frac{1}{n^{1/2}z}C_{k-1})^{-1}\|^{-1}=\sigma_{n}(I-\frac{1}{n^{1/2}z}C_{k-1})=\frac{1}{n^{1/2}|z|}\sigma_{n}(-n^{1/2}zI+C_{k-1})$

Note that $\|-In^{1/2}z\|\leq n$ for $n$ big enough, so we can apply Theorem
A.4 and find that with high probability (greater than $1-n^{-2}$),

$\sigma_{n}(-n^{1/2}zI+C_{k-1})\geq n^{-\widetilde{a}}$

so that

$\frac{1}{\sigma_{kn}(n^{-1/2}E_{1}C^{T}-zI)}\leq\frac{1}{|z|}(1+rn^{\widetilde{a}+1})$

and

$\sigma_{kn}(n^{-1/2}E_{1}C^{T}-zI)\geq\frac{|z|}{(1+rn^{\widetilde{a}+1})}\geq
n^{-\widetilde{a}-2}.$

By the Borel-Cantelli Lemma, the statement follows.

In Lemma A.8, we control the spectral norms of $M-zI$ and $E_{1}C^{T}-zI$.

###### Lemma A.8

Let $M,E_{1},C^{T}$ be defined as in (5) and immediately after it. There exist
constants $r,s>0$ such that for any $n$ large enough,

$\mathbb{P}\left(\left\|n^{-1/2}E_{1}C^{T}-zI\right\|>r+|z|\right)\leq
s\exp(-r^{2}n)$
$\mathbb{P}\left(\left\|n^{-1/2}M-zI\right\|>r+|z|+1\right)\leq
s\exp(-r^{2}n)$

and in particular, there exists $d>0$ such that, with probability $1$,

${\sigma_{1}(n^{-1/2}E_{1}C^{T}-zI)}\leq
d,\qquad{\sigma_{1}(n^{-1/2}M-zI)}\leq d$

for any sufficiently large $n$.

###### Proof

Note that

$\left\|n^{-1/2}M-zI\right\|\leq\|n^{-1/2}E_{1}C^{T}\|+|z|+n^{-1/2},$
$\left\|n^{-1/2}E_{1}C^{T}-zI\right\|\leq\|n^{-1/2}E_{1}C^{T}\|+|z|,$

and $\|E_{1}C^{T}\|=\|C^{T}\|$, implying (see proof of Lemma A.6)

$\mathbb{P}(\|E_{1}C^{T}\|>r\sqrt{n})\leq s\exp(-crn)$

where $c,s,r>0$ are absolute constants. The statement follows immediately.

Lemma A.9, whose proof relies on Theorem A.5, yields a probabilistic lower
bound on the $f(n)$th singular value of $n^{-1/2}E_{1}C^{T}-zI$.

###### Lemma A.9

Let $0<\delta<1/2$, $E_{1},C^{T}$ be defined as in (5) and immediately after
it (so that $E_{1}C^{T}$ is a $kn\times kn$ matrix with $k>1$), $f(n)=\lfloor
kn-n^{1-\delta}\rfloor$, and $0\neq z\in\mathbb{C}$. Then there exists a
positive constant $t\leq 1$ and a positive constant $\varepsilon>0$ such that
lower bound

$\sigma_{f(n)}\left(n^{-1/2}E_{1}C^{T}-zI\right)\geq tn^{\varepsilon-1/2}$

holds almost surely for all sufficiently large values of $n$.

###### Proof

Note first that, denoting by $\widetilde{T}$ the matrix composed by the first
$f(n)$ rows of $n^{-1/2}E_{1}C^{T}-zI$, then by the interlacing theorem for
singular values (Theorem A.1) we have

$\sigma_{f(n)}(n^{-1/2}E_{1}C^{T}-zI)\geq\sigma_{f(n)}(\widetilde{T})$

so it suffices to study $\widetilde{T}$. To this goal, since $f(n)\geq n$, we
partition

$\widetilde{T}=\begin{pmatrix}H-zI&L&P\\\ &zI&\end{pmatrix}$

where $H$ is a $n\times n$ matrix, $L$ is a $n\times(f(n)-n)$ matrix and $P$
is a $n\times(kn-f(n))$ matrix. Since a permutation of the columns does not
change the singular values, we can equivalently study the matrix

$T=\begin{pmatrix}L&H-zI&P\\\ zI&&\end{pmatrix}=\begin{pmatrix}L&R\\\
zI&\end{pmatrix}$

where $R:=[H-zI\,\,\,P]$. $T$ is an $f(n)\times kn$ matrix, and its least
singular value $\sigma_{\min}(T)$ satisfies

$\sigma_{\min}(T)^{2}=\inf_{\|v\|=1}\|v^{*}T\|^{2}=\inf_{\|v\|=1}\|v_{1}^{*}L+zv_{2}^{*}\|^{2}+\|v_{1}^{*}R\|^{2}\geq\inf_{\|v\|=1}(\|v_{1}^{*}L\|-\|zv_{2}^{*}\|)^{2}+\|v_{1}^{*}R\|^{2},$

where $v^{*}=[v_{1}^{*}\,v_{2}^{*}]$. There are three possibilities. If
$v_{2}$ is zero, then the minimum is attained as
$\sigma_{\min}([H-zI\,\,L\,\,P])^{2}$ and if $v_{1}$ is zero, then the minimum
is simply $|z|^{2}$. If neither is true, we can minimize the expression over
$\|v_{1}\|\neq 0,1$. To this goal, denote $y=\|v_{1}\|^{2}$, so that
$1-y=\|v_{2}\|^{2}$ and define $w_{1}=:v_{1}/\sqrt{y}$,
$w_{2}:=v_{2}/\sqrt{1-y}$. Then

$(\|v_{1}^{*}L\|-\|zv_{2}^{*}\|)^{2}+\|v_{1}^{*}R\|^{2}=(\sqrt{y}\|w_{1}^{*}L\|-\sqrt{1-y}|z|)^{2}+y\|w_{1}^{*}R\|^{2}.$

Define now $\alpha:=\|w_{1}^{*}L\|$, $\beta:=|z|$,
$\gamma^{2}:=\|w_{1}^{*}R\|^{2}$, so that we end up with the problem of
minimizing the function

$g(y)=(\sqrt{y}\alpha-\sqrt{1-y}\beta)^{2}+y\gamma^{2}=\beta^{2}+y(\gamma^{2}+\alpha^{2}-\beta^{2})-2\alpha\beta\sqrt{y-y^{2}}.$

To further simplify the notation, it is convenient to introduce
$a:=\gamma^{2}+\alpha^{2}-\beta^{2}$ and $b:=2\alpha\beta$. This trick yields

$g(y)=\beta^{2}+ya-b\sqrt{y-y^{2}};$
$g^{\prime}(y)=a+b\frac{2y-1}{2\sqrt{y-y^{2}}};$
$g^{\prime\prime}(y)=\frac{b}{4}(y-y^{2})^{-3/2}\geq 0.$

In particular, the computation of the second derivative shows that $g(y)$ is a
convex function, and hence, the roots of its derivative must correspond to
minima. If $b>0$, then the minimum is also unique. Observe that by assumption
$\beta\neq 0$, so $b=0$ implies $w_{1}^{*}L=0$ and $\sigma_{\min}(T)^{2}$ is
surely greater than $|z|^{2}$ or $\min_{\|w_{1}\|=1}\gamma^{2}$. Thus, we
assume instead $b>0$ and find a root of the derivative.

$g^{\prime}(y)=0\implies
b(2y-1)=-2a\sqrt{y-y^{2}}\implies\mathrm{sign}(a)=\mathrm{sign}(1-2y)$
$b^{2}(4y^{2}-4y+1)=4a^{2}(y-y^{2})\implies
4(a^{2}+b^{2})(y^{2}-y+\frac{1}{4})-a^{2}=0\implies(y-\frac{1}{2})^{2}=\frac{a^{2}}{4(a^{2}+b^{2})}$

so the unique root of $g^{\prime}(y)$ is

$y_{\star}=\frac{1}{2}-\frac{a}{2\sqrt{a^{2}+b^{2}}}\in(0,1).$

Moreover,

$g(y_{\star})=\beta^{2}+\frac{a}{2}-\frac{a^{2}}{2\sqrt{a^{2}+b^{2}}}-\frac{b^{2}}{2\sqrt{a^{2}+b^{2}}}=|z|^{2}+\frac{a-\sqrt{a^{2}+b^{2}}}{2}$
$=\frac{\gamma^{2}+\alpha^{2}+|z|^{2}-\sqrt{(\gamma^{2}+\alpha^{2}+|z|^{2})^{2}-4|z|^{2}\gamma^{2}}}{2}\geq\frac{|z|^{2}}{2}\frac{\gamma^{2}}{\gamma^{2}+\alpha^{2}+|z|^{2}}$

To minimize the last expression, we can maximize the denominator by

$\gamma^{2}+\alpha^{2}=\|w_{1}^{*}[L\,R]\|^{2}\leq\|[L\,R]\|^{2}\leq(\frac{1}{\sqrt{n}}\|C^{T}\|+|z|)^{2}$

which is, with high probability (see (28)), bounded by $(r+|z|)^{2}$. We can
thus say that there exists a constant $t>0$ such that

$\sigma_{\min}(T)^{2}\geq\min\\{\min_{\|w_{1}\|=1}t^{2}\gamma^{2},\sigma_{\min}([H-zI\,\,L\,\,P])^{2},|z|^{2}\\}.$

Yet, $\min_{\|w_{1}\|=1}\|w_{1}^{*}R\|=\sigma_{\min}(R)$, and again by the
interlacing theorem for singular values (Theorem A.1)
$\sigma_{\min}([H-zI\,\,L\,\,P])\geq\sigma_{\min}(R)$. We can therefore
conclude that

$\sigma_{f(n)}(n^{-1/2}E_{1}C^{T}-zI)\geq\sigma_{\min}(T)\geq\min\\{t\sigma_{\min}(R),|z|\\}$
(31)

for some absolute constant $0<t\leq 1$ with high probability, and almost
surely for all $n$ sufficiently large.

Almost surely, $R$ is full rank, and thus
$\sigma_{\min}(R)=\|R^{\dagger}\|^{-1}$. In particular, Theorem A.5 yields the
tail bound

$\mathbb{P}_{R}\left(\sigma_{n}(R)\leq\tau\right)\leq\left(2\pi(n^{1-\delta}+1)\right)^{-1/2}\left(\frac{k^{2}n^{2}\tau^{2}e^{2}}{(n^{1-\delta}+1)^{2}}\right)^{\lceil
n^{1-\delta}\rceil+1},$

where we used that $R$ is a $n\times(kn+n-f(n))$ matrix with both dimension
less than $kn$, and

$n^{1-\delta}+1\leq\lceil n^{1-\delta}\rceil+1=(kn+n-f(n))-n+1.$

Now, fix any $\varepsilon$ such that $0<\varepsilon<1/2-\delta$ (as
$0<\delta<1/2$, this is surely possible). Choosing $\tau=n^{\varepsilon-1/2}$,
this implies that for $n$ big enough,

$\mathbb{P}_{R}\left(\sigma_{n}(R)\leq
n^{\varepsilon-1/2}\right)\leq\frac{1}{\sqrt{2\pi}}n^{\delta/2-1/2}\left(\frac{(ke)^{2}n^{1+2\varepsilon}}{n^{2-2\delta}}\right)^{n^{1-\delta}+1}.$

Setting $c:=1-2\delta-2\varepsilon>0$, we conclude that

$\mathbb{P}_{R}\left(\sigma_{n}(R)\leq
n^{\varepsilon-1/2}\right)\leq\frac{(ke)^{2}}{\sqrt{2\pi}}n^{\delta/2-1/2-c}\left((ke)^{-2/c}n\right)^{-cn^{1-\delta}},$

and the right hand side goes exponentially to zero. We conclude that
$\sigma_{\min}(R)$ is at least of the order $n^{\varepsilon-1/2}$ with high
probability, and almost surely for sufficiently large $n$. From (31), we get
that for some absolute constant $0<t\leq 1$ and some $\varepsilon>0$ then
almost surely, for all $n$ sufficiently large,

$\sigma_{f(n)}(n^{-1/2}E_{1}C^{T}-zI)\geq tn^{\varepsilon-1/2}.$

###  A.4 Estimates on the singular values of certain random matrices in
Section 4

In Lemma A.10 below, we obtain probabilistic bounds for the first $n$ singular
values of the matrix $M-zI$, where $M$ is defined in (14).

###### Lemma A.10

Let $M$ be the $kn\times kn$ matrix defined as in (14) and $z\in\mathbb{C}$
with norm different from 1. There exist constants $r,s,c>0$ such that, for
every $k>2$,

$\mathbb{P}\left(\sigma_{1}(M-zI)>r\sqrt{k}+1+|z|\right)\leq
s\exp(-crk),\qquad\sigma_{n}(M-zI)\geq|1-|z||.$

In particular, with probability 1,

$\sigma_{1}(M-zI)\leq r\sqrt{k}+1+|z|,\qquad\sigma_{n}(M-zI)\geq|1-|z||$

for all large enough $k$.

###### Proof

From the proof of Lemma A.6, by exchanging the roles of $k$ and $n$, we know
that

$\mathbb{P}(\|C^{T}\|>r\sqrt{k})\leq s\exp(-crk),$ (32)

where $r,s,c$ are absolute positive constants. As a consequence, with high
probability,

$\sigma_{1}(M-zI)\leq\|C^{T}\|+1+|z|\leq r\sqrt{k}+1+|z|.$

Moreover, from (19) and (20), we know that, if $k>2$,

$\sigma_{n}(M-zI)\geq\sigma_{n+1}(M-zI)\geq\sigma_{2n+1}(B-zI)\geq|1-|z||>0.$

In Lemma A.11 below, we obtain a lower bound (valid with probability $1$) for
the smallest singular value of the matrix $M-zI$, where $M$ is defined in
(14).

###### Lemma A.11

Let $M$ be the $kn\times kn$ matrix defined as in (14) and $z\in\mathbb{C}$
with norm different from 1. There exist a positive constant $t$ such that,
with probability 1,

$\sigma_{kn}\left(M-zI\right)\geq tk^{-2}$

for all sufficiently large values of $k$.

###### Proof

If we use Woodbury Lemma (Lemma A.2) on the splitting
$M-zI=(B-zI)+E_{1}\widehat{C}^{T},$ then333$B-zI$ is invertible, and here we
are assuming $I+\widehat{C}^{T}(B-zI)^{-1}E_{1}$ is also invertible, which is
true with probability $1$.

$(M-zI)^{-1}=(B-zI)^{-1}-(B-zI)^{-1}E_{1}(I+\widehat{C}^{T}(B-zI)^{-1}E_{1})^{-1}\widehat{C}^{T}(B-zI)^{-1}$

and

$\|(M-zI)^{-1}\|\leq\|(B-zI)^{-1}\|\left[1+\|(B-zI)^{-1}\|\|\widehat{C}^{T}\|\|(I+\widehat{C}^{T}(B-zI)^{-1}E_{1})^{-1}\|\right].$
(33)

$B-zI$ is a circulant matrix, so its inverse is still a circulant matrix and
its norm can be estimated with (19) by

$\|(B-zI)^{-1}\|=\sigma_{kn}(B-zI)^{-1}\leq\frac{1}{|1-|z||}.$ (34)

More specifically, $(B-zI)^{-1}$ is a block circulant matrix with first block
column

$(B-zI)^{-1}=\begin{pmatrix}\frac{z^{k-1}}{1-z^{k}}I_{n}&\\\
\frac{z^{k-2}}{1-z^{k}}I_{n}&\ddots\\\ \vdots&\ddots\\\
\frac{1}{1-z^{k}}I_{n}&\ddots\end{pmatrix}.$

Consider now the matrix

$I+\widehat{C}^{T}(B-zI)^{-1}E_{1}=\frac{z^{k}}{z^{k}-1}I-\sum_{i=1}^{k}\frac{z^{k-i}}{1-z^{k}}C_{k-i}.$

It consists of the sum of a constant matrix, and a linear combination of
i.i.d. Gaussian variables all having mean $0$. Such a linear combination is
still a normally distributed variable with mean $0$ and variance

$\sum_{i=1}^{k}\left|\frac{z^{k-i}}{1-z^{k}}\right|^{2}=\frac{1}{|1-z^{k}|^{2}}\sum_{i=0}^{k-1}|z|^{2i}=\frac{1}{1-|z|^{2}}\frac{1-|z|^{2k}}{|1-z^{k}|^{2}}:=c(k)^{2}$

where $c(k)\to|1-|z|^{2}|^{-1/2}$ for $k\to\infty$. We can thus write

$I+\widehat{C}^{T}(B-zI)^{-1}E_{1}=\frac{z^{k}}{z^{k}-1}I+c(k)G$

where $G$ is a Gaussian random $n\times n$ matrix where each entry has mean
zero and unit variance. In Theorem A.5, we proved that for any deterministic
matrix $S$, we have

$\mathbb{P}\left(\sigma_{n}(n^{-1/2}G+S)\leq
t\right)\leq\frac{t^{2}n^{2}e^{2}}{\sqrt{2\pi}},$ (35)

so

$\displaystyle\mathbb{P}\left(\|(I+\widehat{C}^{T}(B-zI)^{-1}E_{1})^{-1}\|\geq
p\right)$
$\displaystyle=\mathbb{P}\left(\sigma_{n}(I+\widehat{C}^{T}(B-zI)^{-1}E_{1})\leq\frac{1}{p}\right)$
$\displaystyle=\mathbb{P}\left(\sigma_{n}\left(\frac{z^{k}}{z^{k}-1}I+c(k)G\right)\leq\frac{1}{p}\right)$
$\displaystyle=\mathbb{P}\left(\sigma_{n}\left(\frac{n^{-1/2}z^{k}}{c(k)(z^{k}-1)}I+n^{-1/2}G\right)\leq(n^{1/2}c(k)p)^{-1}\right)$
$\displaystyle\leq\frac{ne^{2}}{\sqrt{2\pi}c(k)^{2}p^{2}}.$

Since $ne^{2}/\sqrt{2\pi}$ does not depend on $k$, we can choose $p=k$ and
conclude by the Borel-Cantelli Lemma that with probability 1 and for any $k$
sufficiently large,

$\|(I+\widehat{C}^{T}(B-zI)^{-1}E_{1})^{-1}\|<k.$ (36)

Finally, from (32), we have that almost surely, for all $k$ sufficiently
large,

$\|\widehat{C}^{T}\|\leq\|C^{T}\|+1\leq 1+r\sqrt{k}$ (37)

for some absolute constant $r>0$. Gathering all the bounds (34), (36), (37)
and substituting into (33), we find that

$\|(M-zI)^{-1}\|\leq\frac{1}{|1-|z||}\left[1+\frac{1}{|1-|z||}(r\sqrt{k}+1)k\right]\leq
t^{-1}k^{2}\implies\sigma_{kn}(M-zI)\geq tk^{-2}$

for sufficiently large $k$ with probability 1, where $t$ is a positive
constant depending only on $n$ and $z$.