Random Perturbations of Matrix Polynomials

A sum of a large-dimensional random matrix polynomial and a fixed low-rank matrix polynomial is considered. The main assumption is that the resolvent of the random polynomial converges to some deterministic limit. A formula for the limit of the resolvent of the sum is derived and the eigenvalues are localised. Three instances are considered: a low-rank matrix perturbed by the Wigner matrix, a product $HX$ of a fixed diagonal matrix $H$ and the Wigner matrix $X$ and a special matrix polynomial. The results are illustrated with various examples and numerical simulations.


Introduction
Motivation. Since the seminal works of Wigner [63] and Marchenko and Pastur [40] spectral theory of random matrices has gathered a huge interest. In particular, studying the limit laws of eigenvalues was considered many times in the literature ( [1,10,11,16,19,20,21,28,54,58,60]). One of the recent techniques in this field is to investigate the limit (in probability) of the resolvent (X N − zI N ) −1 . This was done already for Hermitian matrices, e.g. when X N ∈ C N ×N is a generalised Wigner matrix, see [13,27,31,34,35]. In particular, the local isotropic semicircle law, stated in [13], says that for a suitably chosen family of compact set S N in the upper half-plane sup converges in probability, with a rate O(N − ω 2 ), to zero, see Example 7 for details. Here m W (z) denotes the Stjelties transform of the Wigner semicircle law. As the sets S N approach the real line, the local isotropic semicircle law becomes a tool to study the distribution of the eigenvalues.
Our aim is to investigate the limit of the resolvent for some classes of nonsymmetric matrices and matrix polynomials. Let us recall that already several studies have addressed the canonical forms of nonrandom structured matrices and matrix polynomials [30] and their change under a low-rank perturbation, see e.g. [4,22,24,51,36,42,43,44,45]. However, the theory of random matrix polynomials is yet uncharted.
The additional motivation for the current research lies in noncommutative probability. Recall that deforming a random matrix X N one obtains a deformation of the moment expansion of its limiting resolvent. This was already studied in [52] for X N = H N W N , where W N is a symmetric Wigner matrix and H N = diag(c 1 , 1, . . . , 1), c 1 ∈ R. See also [17,18] for other works on moment deformations.
The results. Let us recall first the basic notions. For a matrix polynomial X(z) = k j=0 X (j) z j , with X (j) ∈ C N ×N , j = 1, . . . n a point λ ∈ C is called an eigenvalue if X(λ)x = 0 for some nonzero x ∈ C N . A polynomial is called regular if det X(z) is a nonzero function. In such case, the matrix X(λ) (λ ∈ C) is invertible if and only if λ is not an eigenvalue. This allows us to define the resolvent of a regular matrix polynomial as which is a matrix valued rational function with poles in the eigenvalues. We will consider eigenvalues and resolvent only for regular polynomials with the leading coefficient being invertible matrix. Hence, we will not investigate the eigenvalue infinity. Let us turn now to the main results of the paper, a further review of necessary linear algebra and probability notions is contained in Section 1.
In Section 2 we will consider a general setting of random matrix polynomials where the degree k of the polynomial is fixed and does not depend on N and the matrices X (i) N are either deterministic or random. The leading assumption is that the polynomias X N (z) are regular and X N (z) is invertible on a set S N ⊆ C, S N ⊆ S N +1 for N = 1, 2, . . . and the resolvent X N (z) −1 converges, pointwise in z, in probability to M (z) on the union of S N (see Definition 5 for details). In such setting we will investigate how these objects behave under a low rank perturbation X N (z) + A N (z). Our first main result, Theorem 11, states precisely how the sets S N , the limit M N (z) and the convergence rate are deformed in this general situation. Further, in Theorem 14, we locate and count the eigenvalues of X N (z) + A N (z), appearing in the union of S N , after such deformation.
Further sections are devoted to the study of concrete ensembles. And so, Section 3 contains a result on low rank non-Hermitian perturbations of Wigner and random sample covariance matrices. We obtain the limit of the resolvent (A N + X N − zI N ) −1 and show the limit points and convergence rate of the non-real points of the spectrum in Theorem 16. Shortly these can be formulated as follows.
Let X N be a generalised Wigner or sample covariance matrix and let A N be deterministic, fixed low rank perturbation, e.g., A N := C ⊕ 0 N −n,N −n ∈ C N ×N . Then the resolvent of A N + X N converges in the maximum norm in probability to M N (z) = m (z)I N − m 2 (z)((C −1 + m (z)I n ) −1 ⊕ 0 N −n,N −n ) .
Furthermore, if z 0 ∈ C + is such that ξ = − 1 m(z0) is an eigenvalue of C with the algebraic multiplicity k ξ and the size of the largest Jordan block equal to p ξ , then the k ξ eigenvalues λ N 1 , . . . , λ N k ξ of X N + A N closest to z 0 are simple and converge to z 0 in probability, with the rate O(N − 1 2p ξ ). Matrices of type HX with H = H * invertible are well known in linear algebra, see e.g [30,42]. Section 4 discusses products H N X N , where H N is a deterministic diagonal matrix with (H N − I N ) being of fixed low rank and X N is a Wigner or a random sample covariance matrix. In Theorem 22 we provide the limit of the resolvent and limit and convergence rate of nonreal eigenvalues for H M X N . It is important to notice that already in this matrix problem it is necessary to apply the main results to nontrivial matrix polynomials of degree one (linear pencils). Namely, we set Section 5 contains a study of matrix polynomials of the form where p, q are polynomials, X N is either a Wigner or a random sample covariance matrix, u N is some deterministic vector and C n , D n ∈ C n×n are diagonal deterministic matrices. This choice is motivated by the fact that matrix polynomials of this form appear in numerical methods for partial differential equations, see [6]. Again, we localise the spectrum of the given above polynomial by means of Theorems 11 and 14 and show difficulties appearing in a particular example.
Relation to existing results. The main novelty of the current paper lies in the formula for the limit of the resolvent and in considering nonlinear eigenvalue problems. So far the limit laws for the resolvent were considered only for matrices and were of isotropic type, i.e., m(z)I N with a scalar analytic function m(z). Our construction leads to limit laws of different type and also for polynomials of degree greater than one. Note that although eigenvalue problems for matrices are special cases of eigenvalue problems for polynomials, many of the methods suitable for matrices, like, e.g., analysis of tr X n , cannot be adapted in the polynomial case. Our general theorems form Section 2 develop a method which works like a 'black box': knowing a limit of the resolvent of a polynomial X N (z) one is able to compute the limit of the resolvent of X N (z) + A N (z). Furthermore, by detecting the sets on which X N (z) + A N (z) is invertible, we localise the eigenvalues.
As it was already said, our technique is applied also to matrices, i.e., to polynomials X N − zI N . Although the low-rank perturbations of Wigner matrices were considered in many papers, see e.g. [7,8,9,10,11,13,16,20,58,60], the authors usually concentrate on Hermitian or symmetric perturbations. The exceptions are the papers [33,58], see the former for the literature on physical motivations. In the latter paper Rochet considered the possibly non Hermitian finite-rank perturbations A of Wigner matrices W , proving results on the limit and convergence rate of non-real eigenvalues. More precisely, the 'Furthermore' part of Theorem 16 (see also the simplified version above) is, generally speaking, a repetition of Theorems 2.3 and 2.10 from [58]. Note that the paper [58] was a continuation of [11], where the authors considered a low-rank perturbation of a random matrix with a distribution invariant under the left and right actions of the unitary group. Analogous convergence rates for outliers O(N −1/(2p ξ ) ) were obtained therein.
In the current paper we show how the resolvent tools can be used to find the convergence rate of the eigenvalues of matrices converging to the non-real limits, repeating the aforementioned result on eigenvalues from [58]. However, in addition to [58], we provide a formula and convergence rate for the limit of the resolvent after perturbation. We also estimate the rate of the convergence (to zero) of the imaginary part non-real eigenvalues which are not outliers, see Example 20. These two aspects were not studied in [58].
The results on the products H N W N also refine the existing ones from [52,64] by showing the limit of the resolvent and considering a much wider class of H N .
The last section on polynomials contains original, up to our knowledge, results on nonlinear eigenvalue problems with random coefficients.
The outcome. There are two main outcomes of the present paper: • extension of the knowledge of limit laws for the resolvents by providing new limit laws for the resolvents of polynomials of type where A N are B N are low rank and non-symmetric matrix and p(z) and q(z) are scalar polynomials, we stress that these limit laws are no longer isotropic; • analysis of limits in N of spectra of polynomials of the above type, with a special emphasis on investigating the convergence rates. In the future, employing results for non-symmetric matrices or structured matrix polynomials would be most desirable, see e.g. [56] for applications in neural networks. However, the limit laws for the resolvent have not been discovered yet, see [15] for a review. Nonetheless, the general scheme we propose in Section 2 is perfectly suited for studying those as well.
1. Preliminaries 1.1. Linear algebra. First let us introduce various norms on spaces of matrices. If b is a vector then by b p we denote the p -norm of b. If A ∈ C k×l then We abbreviate A p,p to A p . Recall that Recall also the following formulas, valid for A = [a ij ] ∈ C k×l , Further, A max denotes the maximum of the absolute values of all entries of A, and clearly By I N we denote the identity matrix of size N . For matrices P ∈ C N ×n and Q ∈ C n×N we define Let us denote the maximal number of nonzero entries in each row of Q by r(Q) and the maximal number of nonzero entries in each column of P by c(P ).

Proposition 1.
For Q ∈ C n×N and P ∈ C N ×n the following inequalities hold we obtain, using formula (4), the following Similarly, The last claim results from the inequalities and the relation (5).
The following elementary result on matrices will be of frequent use. Let A, B ∈ C n×n and let A be nonsingular. Then A + B is nonsingular if and only if I n + BA −1 is nonsingular, and in such case Let · denote any matrix norm. Then Furthermore, In many places of this article we will use the well known Woodbury matrix identity. Let us recall that for invertible matrices X ∈ C N ×N , C ∈ C k×k , and matrices P ∈ C N ×k , Q ∈ C k×N , the matrix X + P CQ is invertible if and only if L := C −1 + QX −1 P is invertible. In such case (11) (X + P CQ) −1 = X −1 − X −1 P L −1 QX −1 .

Probability theory.
In the whole paper we will work with one probability space, which is hidden in the background in the usual manner. By P and E we denote the probability and expectation, respectively. We will use the symbol 'const' to denote a universal constant, independent from N .

Remark 4.
Recall that in the literature (cf. e.g. [13] Definition 2.1) the symbol ≺ denotes the uniform stochastic domination, namely for all ε > 0 and γ > 0 we have (13) sup . It is clear that simultaneous stochastic domination implies uniform stochastic domination and if the variables are Lipschitz continuous then the converse implication also holds, see, e.g., Remark 2.6 of [13] or Corollary 3.19 of [26]. See also Lemma 3.2 of [13] for other properties of stochastic uniform domination. To avoid assuming Lipschitz continuity, we will speak only about simultaneous stochastic domination.
Let us now introduce one of the main objects of our study: a limit law for the resolvent, defined here for random matrix polynomials.
be a random matrix polynomial, i.e. the matrices X +∞] be sequences of deterministic functions such that for any z ∈ N S N the sequence (Ψ N (z)) N converges to zero. We say that the resolvent X N (z) −1 has the limit law M N (z) on sets S N with the rate Ψ N (z) if the eigenvalues of X N (z) are with high probability outside the set S N and Remark 6. Note that the requirement that the eigenvalues are with high probability outside S N should be read formally as: the (parameter-free) event {the eigenvalues are outside S N } holds with probability ≥ 1 − N −γ for N large enough and all γ > 0. This is equivalent to saying that the polynomial X N (z) + A N (z) is invertible with high probability simultaneously in z ∈ S N .
We present main examples, which are the motivation for the above definition.
N be an N × N Hermitian matrix whose entries W ij are independent complex-valued random variables for i ≤ j, such that where const(p) denotes a constant depending on p only. The function is the Stieltjes transform of Wigner semicircle distribution. It was shown in [13] (see also [31]) that for each ω ∈ (0, 1), the resolvent (W N − zI N ) −1 has a limit law Indeed, this can easily be deduced from Theorem 2.12, remark after Theorem 2.15 (see also Remark 2.6) and Lemma 3.2(i) of [13]. The authors call this the isotropic local limit law because of the form M N (z) = m W (z)I N . In the next section we will provide polynomials with the resolvent having limit law of a different type.
Another example of a resolvent having a limit law is given by the same polynomial W N − zI N but now with S N = T, where T is some compact set in the upper halfplane. Observe that in this setting we again have M N (z) = m W (z)I N with the same rate Ψ W N (z), but the estimate (16) can be improved to (17) sup In what follows we will need both constructions presented in this example.
Our second example is the isotropic local Marchenko-Pastur limit law.
whose entries Y ij are independent complex-valued random variables such that Let also As in the previous example, this is an example of an isotropic limit law and can be deduced from the results in [13]: Theorem 2.4, Remark 2.6 and Lemma 3.2(i). Furthermore, one has that As in Example 7 we change the setting by putting S N = T, where T is some compact set in the upper half-plane, which leads to the estimate Further examples of local limit laws in the literature (which are, in particular, limit laws for the resolvent) concern matrices of type ( where w is a complex parameter, applied in the Hermitisation technique, see [14, Theorem 6.1].
Remark 9. It is worth mentioning, that having a limit law for the resolvent of the family of polynomials N it is relatively easy to derive a limit law for the resolvent of the polynomials where X N is a generalised Wigner matrix, is an example of a first order polynomial having a limit law for the resolvent, with a nontrivial leading coefficient.
We conclude this section with an example of a random matrix without a resolvent limit law, to show the difference between the resolvent limit law and stochastic convergence of eigenvalues.
Example 10. Let X N ∈ C N ×N be a diagonal matrix, with elements on the diagonal being i.i.d. standard normal variables. Although the empirical measures of the eigenvalues of X N converge weakly in probability to the normal distribution, the resolvent (X N − zI N ) −1 does not converge in any reasonable sense.

Main results
2.1. The resolvent. In this subsection we will show how a low-dimensional perturbation deforms a resolvent limit law. Recall that r(B), c(B) denote, respectively, the maximal number of nonzero entries in each row and column of a matrix B.
Theorem 11. Let (n N ) N be a nondecreasing sequence and let be a random matrix polynomial. We assume that (a1) X N (z) −1 has a limit law M N (z) on a family of sets S N with rate Ψ N (z), see Definition 5, Then, for any β ∈ (0, α), the eigenvalues of the random polynomial X N (z) + A N (z) are with high probability outside the set under the additional assumption that Ψ N (z) converges to zero for z ∈ N S N .
Proof. We set β = α/2, the proof for arbitrary β < α requires only few technical adjustments. Fix arbitrary γ > 0. Due to (a1) and the definition of stochastic simultaneous domination (Definition holds with probability ≥ 1 − N −γ , for N ≥ N 0 (α, γ) sufficiently large. Note that, if Θ occurs, one has that, for all z ∈ S N , (24), Note that by the assumption (a3), if Θ occurs then Consequently, the above inequality holds with probability ≥ 1 − N −γ , for sufficiently large N ≥ N 0 (α, γ). By the Woodbury matrix equality, the matrix This, together with (7) implies that on the event Θ the matrix X N (z)+P N C N (z)Q N is invertible for sufficiently large N ≥ N 0 (α, γ). As γ was arbitrary we see that the eigenvalues of X N (z) + A N (z) are outside S N with high probability. Now we prove the convergence of (X N (z) + A N (z)) −1 . Let Confronting with (26) and dropping the z-dependence (z ∈ S N ) and the Ndependence we obtain the difference to estimate We will estimate the maximum norm of each summand in the right hand side of the above equation. For this aim we state some preliminary inequalities. Recall that by Proposition 1, assumptions on P N and Q N and (5) one has The stochastic domination below in this proof is simultaneous in z ∈ S N . One has We can now derive the announced estimation of summands of E (2) (z). In the following estimations we again drop the z-dependence and the N -dependence, the stochastic domination below in this proof is simultaneous in z ∈ S N . And so we have (30), (28), (2), (30), (23), (29), , by (23), (23).
Due to the fact that 2.2. The spectrum. In the current subsection the dimension n N will be constant and denoted by n. First let us prove a technical lemma.
Lemma 12. Let the matrices A, B ∈ C k×k . Then Proof. Observe that with The next step in the analysis of the spectra of matrices X N + A N is the following theorem, for its formulation let us introduce a usual technical definition.
Theorem 14. Let n be fixed and let be a random matrix polynomial. We assume that (a1) X N (z) −1 has a limit law M N (z) on a family of sets S N with the rate Ψ N (z), see Definition 5, (a2.1) C(z) is invertible for z ∈ N S N , (a3.1) the following estimate holds with some α > 0, (a4) P N 2 , Q N 2 , c(P N ), r(Q N ) ≤ const, (a5.1) the matrix-valued function z → Q N M N (z)P N is analytic on the interior of N S N and does not depend on N . Let also Assume that the function det K(z) has a zero of order k > 0 at a point z 0 lying in the interior of N S N and let λ N 1 , . . . , λ N k , . . . be the zeros of det L N (z) written down with multiplicities in the order given by their distance to z 0 . Then the first k of them converge to z 0 in the following sense while the k+1-st, and in consequence all following ones, do not converge to z 0 , more precisely, for any β > 0 the set of random variables |λ N k+1 −z 0 | is not simultaneously stochastically dominated by N −β .
Proof. Fix ε, γ > 0, with ε < α. We show that there are exactly k zeros λ large enough and any β ≤ −ε+α k . This will prove both statements. Indeed, setting β = −ε+α k shows that condition (12) in the definition of stochastic simultaneous domination is satisfied for any ε < α, and hence in an obvious way for any ε > 0. Setting β to be arbitrary small shows that |λ N k+1 − z 0 | is not simultaneously stochastically dominated by N −β .
Let us fix an open bounded set T such that z 0 ∈ T and the closure of T is contained in the interior of some S N (N ≥ 1). We may assume without loss of generality that X N (z) is invertible on T. Note that due to (a2') and (a5) the function K(z) is continuous on the closure of T, hence, Moreover, due to Proposition 1 one has Hence, by (a1) and (a4) the probability of the event Hence, the probability of the event is higher than 1−N −γ , for N ≥ N 2 (ε, γ) large enough. As ε 3 −α < 0 and K(z) max is bounded on T the probability of the event Consequently, for sufficiently large N ≥ N 4 (ε, γ) and any β ≤ −ε+α k . Combining (34) and (35) we get with probability higher than 1 − N −γ for N ≥ N 5 (ε, γ) large enough. However, (36) implies, via the Rouché theorem, that det K(z) and det L N (z) have the same number of zeros in B(z 0 , N −β ). Hence, there are exactly k zeros λ 1 , . . . , λ k of det L N (z) in B(z 0 , N −β ) with probability higher than 1 − N −γ for N ≥ N 5 (ε, γ) large enough.
Remark 15. Let us compare Theorems 11 and 14. First note that the latter one has slightly stronger assumptions, which are however necessary for defining the limit point z 0 , also the sequence n N is required there to be constant. Comparing the claims let us note that both theorems state more or less the same fact: the eigenvalues converge to some limit points with a certain convergence rate. If n = 1 the claim of Theorem 14 is slightly than the one of Theorem 11. Namely, the function K(z) is a scalar function in this situation, and the condition |K(z)| ≥ N −α , which constitutes the set S N , locates with high probability the eigenvalues in the (approximate) discs around the points z 0 and with radius equal to N −β , for β < α, while Theorem 14 already states that β = α.
However, already for n = 2 the estimates given by Theorem 14 are weaker, we will see this more clearly in Section 3. The main reason for stating Theorem 14, despite it is giving a weaker estimate, is that it allows us in some situations to count the eigenvalues, while Theorem 11 does not even guarantee that inside each connected component of the complement of S N there is any eigenvalue of X N +A N . Therefore, in what follows, we will use Theorem 11 to get the optimal convergence rate and Theorem 14 to get the number of eigenvalues which converge to z 0 .

Random perturbations of matrices
In this section we will consider the situation where A N (z) = A N is a matrix and X N (z) = X N − zI N . In this subsection, like in Theorem 11, neither n nor C depend on N .
Theorem 16. Let n be fixed and let C ∈ C n×n , A N := P N CQ N ∈ C N ×N , be deterministic matrices, where P N ∈ C N ×n , Q N ∈ C n×N , N = 1, 2, . . . Let X N ∈ C N ×N be a random matrix. We assume that (a1.2) X N is either a Wigner matrix from Example 7 or a random sample covariance matrix from Example 8, so that the resolvent of X N − zI N has a limit law m (z)I N on the family of sets S N,ω with the rate Ψ N (z), where ∈ {W, MP}, respectively, and let T be a compact set that does not intersect the real line.
Then the eigenvalues of X N + A N are with high probability outside the set (37) S N,ω := z ∈ S N,ω : min where ρ < τ < 1 2 , p ξ denotes the size of the largest block corresponding to ξ and σ(D) is the set of eigenvalues of D. The resolvent of the polynomial A N +X N −zI N has on S N,ω and T N the following limit law respectively. Furthermore, if z 0 ∈ C \ R is such that ξ = − 1 m (z0) is an eigenvalue of D with algebraic multiplicity k ξ and the size of the largest Jordan block is equal to p ξ , then the k ξ eigenvalues λ N 1 , . . . , λ N k ξ of X N + A N closest to z 0 are simple and converge to z 0 in the following sense (39) |λ provided that the independent random variables constituting the matrix X N , i.e. W ij in Example 7 or, respectively, Y ij in Example 8, have continuous distributions.
Proof. I: Limit law for the resolvent First we prove that the resolvent of X N + A N contains with high probability the set S N,ω . We fix τ < 1 2 and let τ ∈ (τ, 1/2). Let us note that the assumptions (a1)-(a4) of Theorem 11 are satisfied with α = ω/2 and β = τ ω < α. With K(z) as in Theorem 11 and D = SJ D S −1 , where S is invertible and J D is in Jordan normal form, one has Note that (I n + m (z)J D ) −1 is a block-diagonal matrix with blocks corresponding to possibly different eigenvalues ξ of D, of possibly different sizes r, of the form and the non-indicated entries are zeros. Hence, As m (z) is bounded on N S N,ω we can apply estimates (16) and (21) and Theorem 11 with α = ω 2 and get that the eigenvalues of X N + A N are with high probability outside the set which in turn contains the following sets z ∈ S N,ω : max ξ∈σ(D) where δ is an appropriate constant and N > N 0 (τ, τ , δ, ω) is sufficiently large. The formula for the limit law M N (z) follows straightforwardly, let us discuss now the convergence rate. First consider the case S N = S N,ω and apply Theorem 11. For this aim we need to check the additional assumption that Ψ N (z) = N ω 2 Ψ N (z) converges to zero for each z. This is, however, satisfied due to Ψ N (z) = O(N −1/2 ) and ω < 1. Now consider the case S N = T, we apply Theorem 11 with α = τ and β = ρ. Repeating the arguments from the previous case (with (17), (22) used instead of (16) and (21)) we see that the eigenvalues of X N + A N are with high probability outside T N . The formula for the limit law M N (z) follows straightforwardly, to see the convergence rates observing that the additional assumption on convergence of Ψ N (z) is satisfied due to τ < 1/2.
II: Eigenvalues ('Furthermore' part). First note that due to the fact that m (z 0 ) = 0 the function det K(z) = det C −1 det(I n + m (z)D) has a zero of order k ξ at z 0 . By Theorem 14, det L N (z) has exactly k ξ zeros, counting with multiplicities, λ N 1 , . . . , λ N k ξ , that converge to z 0 as (40) |λ N l − z 0 | ≺ N −α , l = 1, . . . , k ξ with some α > 0. Each of these zeros is an eigenvalue of X N + A N converging to z 0 . We show now that α = 1 2p ξ . Let us fix a compact set T, not intersecting the real line, and such that z 0 is in the interior of T. As d(1+m (z)ξ) dz (z 0 ) = 0 one gets immediately from (40) and the form of the set T N in (38) that α ≥ τ /p ξ with arbitrary τ < 1 2 . Hence, (40) holds with α = 1 2p ξ as well. Let us see that the zeros of L(z) are almost surely simple. Note that det L(z) = det C −1 det(X N −zI+A) det(X N −zI) is a rational function and if it has a double zero then det(X N − zI + A) has a double zero. However, it is a well-known fact that the eigenvalues of the random matrix with the continuously distributed entries are almost surely simple, see, e.g., [61, Exercise 2.6.1]. Thus the zeros λ N 1 , . . . , λ N k ξ are almost surely mutually different.
Remark 17. The Theorem above may be easily generalised to the situation where the resolvent of X N − zI N has a limit law of the form µ(z)I N with µ(z) being a Stieltjes transform of a probability measure. The only exception is the counting of the eigenvalues in (39), namely it is not true that the eigenvalues converging to z 0 need to be simple and that their number has to be precisely k ξ . The proof of this generalisation follows exactly the same lines except the last two paragraphs.  Figure 1. Note that the log-log plot support the conjecture that the exponents in estimates of the convergence rates cannot be in practice improved. Namely, in all four cases the slope of the corresponding group of points are approximately −1/2 in cases (1), (2) and (4) and −1/6 case (3). Furthermore, it is visible that while the exponent in the rate of convergence in the cases (1), (2) and (4) is  Example 20. In this example we will compare the convergence of the eigenvalues of X N +A N to the real axis in three different situations. Here X N is again a Wigner matrix. First let us take the matrix A (1) N from Example 18. In this situation there is one eigenvalue λ N 1 of X N + A N converging to z 0 = 63 i /8, cf. Example 18, and the set S W N,ω is for large N a rectangle with an (approximately) small disc around z 0 removed, similarly as for C (4) in Figure 2. The other eigenvalues converge to the real line with the rate N −1 , i.e. ∆(N ) ≺ N −1 where By definition of S W N,ω one can see that ∆(N ) ≺ N −1+ω with the arbitrary parameter ω ∈ (0, 1). However, by the definition of stochastic domination, it means that ∆(N ) ≺ N −1 , see Remark 3. The numbers ∆(N ) are plotted in Figure 4. One can observe that the plot bends in the direction of the the line given by N −1 , which is still in accordance with the definition of stochastic domination.
The second situation to consider is A  However, z 0 = 0 can be seen as a solution, if we define m W (0) as lim y↓0 m W (y i) = 1i. Hence, the set S W N is a rectangle with an (approximately) half-disc around z 0 = 0 removed, see Figure 3. The half disc has radius of order N − 1 2 , hence (42) ∆(N ) = max {| Im λ| : λ ∈ σ(X N + A N )} converges to zero with the rate N − 1 2 , which can be seen in Figure 4. The last situation to consider is A Here, according to Corollary 19 the sets S W N and S W N coincide. Hence, all the eigenvalues converge to the real axis with the rate N −1 , the same comments concerning ω as in the C (4) case above apply. The plot of ∆(N ), defined as in (42), can be seen in Figure 4.
As in Example 18, the plot suggest that the exponent in the convergence rate cannot be improved in the discussed examples.
The next corollary will concern the class of port Hamiltonian matrices, i.e. matrices of the form A − Z, where A = −A * and Z is positive definite. This class has recently gathered some interest [46,47] due to its role in mathematical modelling. Clearly, the spectrum of A − Z lies in the closed left half-plane. We will consider below the case where A = C ⊕ 0 N −n,N −n is a nonrandom matrix with n fixed and Z = Y * Y is the random sample covariance matrix. For the sake of simplicity we will take a square random sample covariance matrix (N = M ).
Proof. Consider the matrix Z N − A N . As i C is Hermitian, the matrix D = −C from Theorem 16 does not have any Jordan chains longer than one. Note that −z j is the solution of 1 + i t j m 1 (z) = 0. The claim follows now directly from Theorem 16.

Random matrix pencils and H-selfadjoint random matrices
In this section we will employ the setting of random matrix pencils. The theory is aimed on localisation of the spectrum of the products of matrices H N X N . Although the linear pencil appear only in the proof of the main result (Theorem 22), its role here is crucial. In what follows H N is a nonrandom diagonal matrix and X N is a generalised Wigner or random sample covariance matrix. To prove a limit law for the resolvent we need n N to be a slowly increasing sequence, while to count the number of eigenvalues converging to their limits we need n N to be constant. Note that unlike in the case of perturbations X N + A N considered in Theorem 16 we do not need to localise the spectrum near the real line, as the spectrum of X N is symmetric with respect to the real line and contains at most n N points in the upper half plane, see e.g. [30]. The following theorem explains the behavior of all non-real eigenvalues of H N X N . It covers the results on locating the nonreal eigenvalues of H N X N of [52] and [64], where the case n N = 1 was considered. In addition, the convergence rate and formula for the resolvent are obtained.
Theorem 22. Let n N ≤ log N be a sequence of nonrandom natural numbers and let where (c j ) j is a negative sequence such that the sequence (c −1 j ) j is bounded. We also assume that (a1.3) X N is either a Wigner matrix from Example 7 or a random sample covariance matrix from Example 8, so that the resolvent of X N − zI N has a limit law m (z)I N on the compact set T with the rate Ψ N (z), where ∈ {W, MP}, respectively. Then, with high probability, the eigenvalues of the matrix H N X N are outside the set where β < α < 1 2 . Furthermore, the resolvent of H N X N − zI N has a limit law on T N (43) M N (z) = diag(g 1 (z), . . . , g n N (z)) ⊕ m (z)I N −n N , z ∈ T N , with g j (z) = m (z) 1 (c j − 1)zm (z) + c j and the rate n N N α Ψ N (z).
If, additionally, n N = n is constant and if k j denotes the number of repetitions of c j in the sequence c 1 , . . . , c n and if, using the notation of Example 8, , j = 1, . . . , n then, for any j = 1, 2, . . . , n, there are exactly k j eigenvalues λ N 1 , . . . , λ N kj of H N X N converging to z j in probability as Proof. First note that for z ∈ T we have (44) Note that, by (17) and (22), the polynomials X N (z) = X N − zI N , A N (z) = zP N C N Q N satisfy the assumptions of Theorem 11 with any α < 1 2 . Hence, the resolvent of H N (X N − zI N + zP N C N Q N ) has a limit law . . , g n N (z)) ⊕ m (z)I N −n N on the set T N with β < α < 1 2 . The rate of this convergence is n N N α Ψ N (z), due to the fact that H −1 N is a diagonal matrix with bounded entries. Finally, note that n N N α Ψ x N (z) converges to zero pointwise in z. Now let us prove the statement concerning the eigenvalues. By (44) the eigenvalues of H N X N are the eigenvalues of the linear pencil X N − zI N + zP N C N Q N . Define K(z) as in Theorems 11 and 14: As the matrix C is diagonal we have that det K(z) = 0 is equivalent to (45) c j c j − 1 for some j ∈ {1, . . . n}. Consequently, the points z j , j = 1, . . . , n are precisely the zeros of det K(z). If X N is a Wigner matrix, then using the well known equality it is easy to show that the equation (45) has for each j = 1, . . . n two complex solutions If X N is a random sample covariance matrix then direct computations give the formula for z j . By Theorem 14, for each z j , there are k j eigenvalues of X N converging to z j as (47) |λ with some α > 0. By the first part of the theorem, with high probability, the eigenvalues of H N X N are outside the set (47) can be chosen as arbitrary β < 1 2 . Hence by the definition of stochastic domination (see Definition 2), equation (47) holds with α = 1 2 as well.
Remark 23. Let us note two facts about the formula for z j in Theorem 22 above. In the Wigner matrix case the point z j is in the upper half-plane and its complex conjugate is also a limit point of k j eigenvalues of H N X N due to the symmetry of spectrum of H N X N with respect to the real line.
Further, in the random sample covariance matrix case it holds that z j < 0. Further, if the underlying matrix is a square matrix, then the formula simplifies to Moreover, one eigenvalue of X N converges in probability to z 1 = √ 2 2 i and two eigenvalues converge in probability to z 2 = 2 √ 3 3 i.
We conclude the section with a different type of a random pencil.
Remark 25. We recall the method of detecting damped oscillations in a noisy signal via Padé approximations of the Z-transform of the signal, proposed in [5]. The method has found several practical applications [12,32,50,55], its numerical analysis can be found in [5]. Here finding the limit law (if it exists) for the resolvent of the pencil zU 0 − U 1 , where , and s = s 0 , . . . s 2n−1 is a white noise, would substantially contribute to the signal analysis, via Theorem 11 above, see [2,3] for details. The spectral properties of Hankel matrices were studied e.g. in [19]. The investigation of the pencil zU 0 − U 1 is left for subsequent papers.

Analysis of some random quadratic matrix polynomials
In the present section we will consider the spectra of matrix polynomials of the form (48) X N − p(z)I N + q(z)u N u * N , where X N is either a Wigner or a random sample covariance matrix and u N is some deterministic vector and p(z) and q(z) are some (scalar-valued) polynomials and polynomials of type where C n , D n ∈ C n×n are diagonal deterministic matrices. This we see as a step forward in a systematic study of polynomials with random entries. The problem of localising the spectrum of (48) is a clear extension of the usual perturbation problem for random matrices, see e.g. [11,16,20,28,33,42,43,49,54,60]. Indeed, the latter problem can be seen as studying the pencil X N − zI N + u N u * N . In Example 27 below we will demonstrate an essential difference between these two problems. The general case, i.e. following the setting of Theorem 11, is as yet unclear and requires developing new methods for estimating the expression K(z) −1 2 appearing therein. Note that matrix polynomials of type A − p(z)I N + q(z)uu * (A ∈ C n×n , u ∈ C n ) appear in many practical problems connected with modelling, cf. [6].
be a random matrix polynomial, and let neither p(z) nor q(z) depend on N . We assume that (a1.2) X N is either a Wigner matrix from Example 7 or a random sample covariance matrix from Example 8, so that the resolvent of X N has a local limit law m (z)I N on the family of sets S N,ω with the rate Ψ N (z), where ∈ {W, MP}, respectively, and let T be a compact set that does not intersect the real line. (a2.3) p(z) and q(z) are fixed nonzero polynomials, (a4.3) u N is a deterministic vector of norm one, having at most n nonzero entries, where n is fixed and independent from N .
Then the eigenvalues of X N (z) + A N (z) are with high probability outside the set and (49) where β < α < 1 2 . The resolvent of the polynomial X N (z) + A N (z) has on S N,ω and T N the following limit law with the rates N ω 2 Ψ N (z) and N α Ψ N (z), respectively.
Furthermore, for each solution z 0 with p(z 0 ) ∈ C \ R of the equation m (p(z)) + q(z) −1 = 0 there exists eigenvalues λ N j , j = 1, . . . , k, where k ≥ 1, of X N (z) converging to z 0 as |z 0 − λ N j | ≺ N −1/2 , j = 1, . . . , k. Proof. For the proof we note that the assumptions of Theorem 11 are satisfied, X N (z) has a limit law m (p(z))I N and that and so the first part of the claim follows directly.
The statement concerning the convergence of eigenvalues follows from Theorem 14 and from the form of the set T N , cf. the proofs of Theorems 16 and 22.
Example 27. We present an example showing how the techniques above may be useful in localising eigenvalues of second order matrix polynomials. Consider a particular matrix polynomial (48): where W N is the real Wigner matrix. Note that a similar matrix can be found in [6] as Problem acoustic_wave_1d (after substitution z = λ N ), but with the Wigner matrix replaced by  The matrix polynomial in (50) is real and symmetric, hence its spectrum is symmetric with respect to the real line. Let us note that the spectrum is a priori not localised on the real and imaginary axis, for small N we can have eigenvalues with relatively large both real and imaginary parts. However, these eigenvalues will converge to the real and imaginary axis as N grows. To see this, first note that the Figure 5. The spectrum of the quadratic random polynomial (50) with N = 500 (red crosses) and the set (51) (blue). setsS W N,ω ⊂ p −1 (S W N,ω ) are contained in the first and third quadrant of the plane. However, due to the symmetry of the spectrum, we may extend the set to (51) z ∈ C : z ∈S W N,ω orz ∈S W N,ω , which lies in all four quadrants of the complex plane. The spectrum of P (λ) is located with high probability in the complement of the set (51). A sample set (51) is plotted in Figure 5. Note that the complement of (51) contains both the real and the imaginary axis with some neighbourhood, which is not clearly seen in the picture. Hence, for large N , the spectrum is either (approximately) on the real or imaginary axis. Furthermore, the real spectrum concentrates on p −1 ([−2, 2]) = − 1 π , 1 π and there is also one real eigenvalue (near z = 0.5), which converges to a real solution of (52) m W (4π 2 z 2 − 2) + 1 2π 2 z 2 + 2π i z − 1 = 0 (we extend the function m W (z) onto the real line as m W (x) = lim y↓0 m(x + i y), similarly for the matrix A (5) in Example 20). Due to its role in applications, computation of spectra of quadratic polynomials is currently an important task, see e.g. [6,38,41,53,62]. The usual procedure is a linearisation (see [37,39]), which, however, has some limitations [29]. Above we obtained a family of matrix polynomials (i.e. one coefficient is a random matrix) for which we are able to control the real part of the non-real eigenvalues. This property can be useful e.g. in testing particular numerical algorithms, by plotting maximum of | Re λ|, over all non-real eigenvalues λ. If the algorithm works properly then this quantity should converge to zero as approximately N −1 . A preliminary picture in matlab ( Figure 6) of maximum of | Re λ|, over all non-real eigenvalues λ, does not reveal any numerical anomalies for N ≤ 10 3 and the upper bound of order N −1 is visible. We present yet another possible application of the main results. Again, investigating this particular polynomial is motivated by examples from [6].  . . , d n ∈ C are fixed. Let X N (z) = z 2 I + zX N ∈ C N ×N [z] be a random matrix polynomial. We assume that (a1.4) X N is either a Wigner matrix from Example 7 or a random sample covariance matrix from Example 8, so that the resolvent of z 2 I N + zX N has a limit law z −1 m (z)I N on the family of sets S N,ω with the rate Ψ N (z), where ∈ {W, MP}, respectively, and let T be a compact set that does not intersect the real line. (a2.4) z 2 c i + d i = 0 for z ∈ S N,ω or z ∈ T, respectively, i = 1, . . . , n, e.g. c i , d i > 0, for i = 1, . . . , n.
Then the eigenvalues of X N (z) + A N (z) are with high probability outside the set S N,ω = z ∈ S N,ω : min where β < α < 1 2 . The resolvent of the polynomial X N (z) + A N (z) has on S N,ω and T N the following limit law , i = 1, . . . n, with the rates N ω 2 |z| −2 Ψ N (z) and N α Ψ N (z), respectively.
Furthermore, for each solution z 0 with z 0 ∈ C \ R of the equation (z 2 c i + d i ) −1 + z −1 m (z) = 0 there exists eigenvalues λ N j , j = 1, . . . , k, where k ≥ 1, of X N (z) converging to z 0 as |z 0 − λ N j | ≺ N −1/2 , j = 1, . . . , k. Proof. First note that the resolvent of z 2 I N +zX N has indeed the limit law z −1 m (z)I N on the same sets and with the same convergence rate as X N − zI N . Now the proof becomes another application of Theorem 11 with . . .