The Lanczos Algorithm Under Few Iterations: Concentration and Location of the Output

TL;DR

This paper analyzes the behavior of the Lanczos algorithm with a randomly chosen initial vector, showing that for few iterations, its output is highly concentrated and nearly deterministic, with implications for spectral approximation and outlier detection.

Contribution

It provides theoretical bounds on the concentration and location of Lanczos outputs after few iterations, including asymptotic results and failure probabilities for outlier detection.

Findings

Lanczos outputs concentrate around medians with high probability after O(log n) iterations.

The spectral density approximation via Jacobi coefficients is justified asymptotically.

Lanczos fails to detect outliers with high probability within O(log n) iterations.

Abstract

We study the Lanczos algorithm where the initial vector is sampled uniformly from $\mathbb{S}^{n-1}$. Let $A$ be an $n \times n$ Hermitian matrix. We show that when run for few iterations, the output of Lanczos on $A$ is almost deterministic. More precisely, we show that for any $ \varepsilon \in (0, 1)$ there exists $c >0$ depending only on $\varepsilon$ and a certain global property of the spectrum of $A$ (in particular, not depending on $n$) such that when Lanczos is run for at most $c \log n$ iterations, the output Jacobi coefficients deviate from their medians by $t$ with probability at most $\exp(-n^\varepsilon t^2)$ for $t<\Vert A \Vert$. We directly obtain a similar result for the Ritz values and vectors. Our techniques also yield asymptotic results: Suppose one runs Lanczos on a sequence of Hermitian matrices $A_n \in M_n(\mathbb{C})$ whose spectral distributions converge in Kolmogorov distance with rate $O(n^{-\varepsilon})$ to a density $\mu$ for some $\varepsilon > 0$. Then we show that for large enough $n$, and for $k=O(\sqrt{\log n})$, the Jacobi coefficients output after $k$ iterations concentrate around those for $\mu$. The asymptotic setting is relevant since Lanczos is often used to approximate the spectral density of an infinite-dimensional operator by way of the Jacobi coefficients; our result provides some theoretical justification for this approach. In a different direction, we show that Lanczos fails with high probability to identify outliers of the spectrum when run for at most $c' \log n$ iterations, where again $c'$ depends only on the same global property of the spectrum of $A$. Classical results imply that the bound $c' \log n$ is tight up to a constant factor.

Equations

$$\mathrm{Kol}(\mu ,\nu ):=\underset{t\in \mathbb{R}}{\sup }\left|\mu ((-\infty ,t])-\nu ((-\infty ,t])\right|.$$

$$⟨f,g⟩={\int }_{\mathbb{R}}f(x)g(x)d\mu (x).$$

$$\int p_k(x)p_l(x)d\mu (x)={\delta }_{kl}.$$

$${\gamma }_k={\left({\int }_{\mathbb{R}}{\pi }_k^2(x)d\mu (x)\right)}^{-\frac{1}{2}}.$$

$${\int }_{\mathbb{R}}{\pi }_k^2(x)d\mu (x)=\underset{q\in {\Gamma }_k}{\min }{\int }_{\mathbb{R}}{q}^2(x)d\mu (x),$$

$$x{p}_k(x)={\beta }_{k-1}{p}_{k-1}(x)+{\alpha }_k{p}_k(x)+{\beta }_k{p}_{k+1}(x),k\ge 1,$$

$$\text{and}x{p}_0(x)={\alpha }_0{p}_0(x)+{\beta }_0{p}_1(x),k=0.$$

$${\gamma }_k={\left(\prod _{i=0}^{k-1}{\beta }_i\right)}^{-1}.$$

$${\beta }_k=\frac{\sqrt{D_{k-1}{D}_{k+1}}}{D_k}\text{and}{\gamma }_k=\sqrt{\frac{D_{k-1}}{D_k}},k\ge 0,$$

$$D_k=\prod _{i=0}^k{\gamma }_i^{-2}.$$

$$p_k(x)=\frac{\det {\tilde{M}}_k(x)}{\sqrt{D_{k-1}{D}_k}}.$$

$${\mu }^u=\sum _{i=1}^n{u}_i^2{\delta }_{{\lambda }_i}.$$

$$\frac{1}{n}\sum _{i=1}^n{\delta }_{{\lambda }_i}.$$

$$t_i^{(k)}=\underset{p\in {\mathcal{P}}_{k-1}^{(i)}}{\min }\underset{j:j\ne i}{\max }|p({\lambda }_j)|,$$

$$\left|\left\{\lambda \in \Lambda :\frac{1}{|T|}\sum _{t\in T}\log |\lambda -t|\ge \log \omega \right\}\right|\ge \delta n.$$

$$2\exp \left\{-\frac{\min \{\delta ,1/50{\}}^2}{32}n\right\}+2\exp \left\{-\frac{1}{64}{\left(\frac{\omega }{4\Vert A\Vert }\right)}^{2i}{\delta }^2{t}^{2}n\right\}.$$

$$k=\min \left\{j,\frac{1}{2\log \frac{M}{\omega }}\left(c\log n+\log \frac{\kappa \delta }{2mg}\right)\right\}$$

$$2\exp \left\{-\frac{\min \{\delta ,1/50{\}}^2}{32}n\right\}+2\exp \left\{-\frac{1}{16}{n}^{1-2c}\right\}$$

$$\Vert J_{k_n}(u)-J_{k_n}(\mu )\Vert {⟶}_{P}0,$$

$$\Vert {\vec{r}}_{k_n}(u)-{\vec{r}}_{k_n}(\mu ){\Vert }_{L^{\infty }({\mathbb{R}}^{k_n})}{⟶}_{P}0,$$

$$\mathbf{P}[|f(u)-\tilde{f}|>\epsilon ]\le \mathbf{P}[u\in {\mathbb{S}}^{n-1}\setminus \Omega ]+2\exp \{-4n{\epsilon }^2/L^2\},$$

$$d_2(f(x),f(y))\le c^{\prime }{d}_1(x,y)\forall x,y\in U.$$

$$d_2(f(w_i),f(w_{i+1}))\le (c+\epsilon )d_1(w_i,w_i+1).$$

$$\Vert {\mathrm{Proj}}_{W_x^{\perp }}(A{v}_j(x))-{\mathrm{Proj}}_{W_y^{\perp }}(A{v}_j(y))\Vert$$

$$\le \Vert {\mathrm{Proj}}_{W_x^{\perp }}(A(v_j(x)-v_j(y)))\Vert +\Vert ({\mathrm{Proj}}_{W_x^{\perp }}-{\mathrm{Proj}}_{W{y}^{\perp }})(A{v}_j(y))\Vert$$

$$=\Vert {\mathrm{Proj}}_{W_x^{\perp }}(A(v_j(x)-v_j(y)))\Vert +\Vert ({\mathrm{Proj}}_{W_x}-{\mathrm{Proj}}_{W_y})(A{v}_j(y))\Vert .$$

$$\Vert {\mathrm{Proj}}_{W_x^{\perp }}(A(v_j(x)-v_j(y)))\Vert \le \Vert A\Vert \Vert v_j(x)-v_j(y)\Vert \le \Vert A\Vert ((4\Vert A\Vert {)}^j{\gamma }_j(\tilde{u})+\epsilon )\Vert x-y\Vert .$$

$$\sum _{i=0}^j(4\Vert A\Vert {)}^i{\gamma }_i(\tilde{u})\le \sum _{i=0}^j{4}^i\Vert A{\Vert }^j{\gamma }_j(\tilde{u})\le \frac{4^{j+1}\Vert A{\Vert }^j{\gamma }_j(\tilde{u})}{3}.$$

$$\Vert {\mathrm{Proj}}_{W_x}(w)-{\mathrm{Proj}}_{W_y}(w)\Vert \le \sum _{i=0}^j\Vert ⟨v_i(x),w⟩v_i(x)-⟨v_i(y),w⟩v_i(y)\Vert$$

$$\Vert ⟨v_i(x),w⟩v_i(x)-⟨v_i(y),w⟩v_i(y)\Vert$$

$$\le \Vert v_i(x)-v_i(y)\Vert \Vert w\Vert +\Vert v_i(x)-v_i(y)\Vert \Vert v_i(y)\Vert \Vert w\Vert$$

$$\le 2(4\Vert A\Vert {)}^i{\gamma }_i(\tilde{u})\Vert x-y\Vert .$$

$$\Vert {\mathrm{Proj}}_{W_x}(w)-{\mathrm{Proj}}_{W_y}(w)\Vert \le \frac{2}{3}\cdot 4^{j+1}\Vert A{\Vert }^j{\gamma }_j(\tilde{u})\Vert x-y\Vert ,$$

Full text

In a different direction, we show that Lanczos fails with high probability to identify outliers of the spectrum when run for at most $c^{\prime}\log n$ iterations, where again $c^{\prime}$ depends only on the same global property of the spectrum of $A$. Classical results imply that the bound $c^{\prime}\log n$ is tight up to a constant factor.

1 Introduction

Eigenvalue problems are ubiquitous in science and engineering. However, most applications require analyzing matrices whose large dimension makes it impractical to exactly compute any important feature of their spectrum. It is for this reason that iterative randomized algorithms have proliferated in numerical linear algebra [Saa11, TBI97].

In this context, iterative randomized algorithms provide an approximation of the spectrum of the matrix in question, where the accuracy of the approximation improves as the number of iterations increases. For any such algorithm, it is natural to ask the following questions:

1. (Q1)

How much does the random output vary? 2. (Q2)

How many iterations are necessary and sufficient to obtain a satisfactory approximation?

The present work, theoretical in nature, addresses the above questions for one of the most widely used algorithms for eigenvalue approximation, namely, the Lanczos algorithm. Throughout the paper we assume exact arithmetic.

1.1 The Lanczos algorithm

Recall that when run for $k$ iterations, the Lanczos algorithm outputs a $k\times k$ tridiagonal matrix called the Jacobi matrix. The nontrivial entries of the the Jacobi matrix are the Jacobi coefficients, which we denote by $\alpha_{k}$ and $\beta_{k}$, and its eigenvalues are the Ritz values, which we denote by $r_{k}$. Oftentimes, the Jacobi coefficients and the Ritz values provide important information about the spectrum of the matrix. In particular, when $k=n$, the Ritz values are exactly the eigenvalues of $A$, and hence the full spectrum is recovered. However, in practice it is usually too expensive to perform $\Theta(n)$ iterations.

Outlying eigenvalues and Ritz values. The success of the Lanczos algorithm resides to some extent in its ability to find the outliers of the spectrum of the matrix $A$ with very few iterations. By outliers, we mean the eigenvalues distant from the region in which the majority of the spectrum accumulates (the bulk). Hence, the algorithm is of particular interest in most applications in science and engineering [Saa11].

Bulk spectrum and Jacobi coefficients. Lanczos-type methods can also be used to approximate the global spectral density of large matrices, also known as density of states; for a survey of techniques see [LSY16]. In applied mathematics, large matrices can arise as discretizations of infinite-dimensional operators such as the Laplacian or as finite-dimensional representations of an infinite-dimensional Hamiltonian. Computing the eigenvalues and Jacobi coefficients of the finite-dimensional operator then yields information about the infinite-dimensional operator and the underlying continuous system. For an example, see [SdJL*+*18], or Section 7 of [VDHVDV01] for numerical experiments and bounds for the Lanczos algorithm applied to an explicit discretized Laplace operator.

In the setting described above, the Jacobi coefficients contain all the information of the spectral density of the infinite-dimensional operator in question and even the first few coefficients are of use. To give an example, in [Hay80] the Haydock method (as it is termed today) was introduced. This method exploits the fact the resolvent of an operator admits a continued fraction expansion where the coefficients are precisely the Jacobi coefficients, and hence knowing these quantities is fundamental to understanding the spectral density of the operator—see 3.2.2 [LSY16] for a summary of the Haydock method.

Using a slightly different perspective, note that from the $k\times k$ Jacobi matrix of an operator one can obtain the $[k-1,k]$ Padé approximation of its resolvent [VA06]. In particular, knowing the $k\times k$ Jacobi matrix is enough to compute the first $2k-1$ moments of the spectral density of the infinite-dimensional operator.

In applications, sophisticated modifications of the Lanczos algorithm are used [GU77, CRS94, LXV*+*16]. Since the goal of the present paper is to introduce proof techniques and theoretical tools that have not been exploited previously, we only deal with the simplest version of the Lanczos algorithm and do not strive to obtain optimal constants in our bounds and theorems when providing answers for questions (1) and (2).

1.2 Question (1): Our contributions

As far as we are aware, there is no previous work addressing this question for the Lanczos algorithm. In this paper we show that there is a $c>0$ such that for $n$ large enough, the output of the Lanczos procedure is almost deterministic when run for at most $c\log n$ iterations. More precisely, in Theorem 2.4 we show that for $k\leq c\log(n)$ and $\varepsilon\in(0,1/2)$, deviations of the order $n^{-\varepsilon}$ of the Jacobi coefficients $\alpha_{k}$ and $\beta_{k}$ computed by Lanczos occur with exponentially small probability. For an illustration, see Figure 1. The strength of our probability bound deteriorates as $k$ grows. The constants in the theorem depend only on an easily computed global property of the spectrum which we call equidistribution.

From the point of view of random matrix theory, the problem treated in the present paper is atypical. In random matrix theory, most of the studied models have a rich probabilistic structure that can be exploited to obtain results about the eigenvalue distribution of the matrix. By contrast, in our case, the Jacobi matrix output by the Lanczos algorithm is a random matrix obtained by running a complicated deterministic dynamic over a minimal source of randomness—a single uniform random unit vector. Hence, in order to obtain results similar to the ones presented in this article, the structure of the algorithm needs to be exploited in an involved way. We use the ubiquitous concentration of measure phenomenon for Lipschitz functions in high dimension, together with a careful control of the variables appearing in the Lanczos algorithm and their Lipschitz constants as functions of the random input. Roughly speaking, the Lipschitz constant is exponential in the number of iterations, which yields concentration in the regime of at most $c\log n$ iterations for sufficiently small $c$. Throughout the analysis we use elementary results in the theory of orthogonal polynomials.

In view of the fact that the output of the Lanczos algorithm is sharply concentrated under few iterations, one may ask which values the output is concentrated around. Toward the end of this introduction we give an overview of our results in this direction.

1.3 Question (2): Previous work

Regarding the problem of detecting outliers of the spectrum via the Lanczos algorithm, theoretical answers to the sufficiency part of Question (2) posed above appeared decades ago. Most of them in essence give an upper bound on the number of iterations required to obtain an accurate approximation when the input is an $n$-dimensional matrix $A$—see [Kan66, Pai71, Saa80]. Roughly speaking, previous literature provides inequalities that state that $k\geq C\log n$ iterations suffice for the Lanczos algorithm to approximate the true extreme eigenvalues of $A$ very well, making the use of $O(\log n)$ iterations common in practice—see [KW94] or [VDHVDV01] for examples of inequalities that give this bound. The constant $C$ in the results mentioned above is determined by features of the spectrum of $A$; typically, these features are the diameter of the spectrum and the gaps between the outliers and the bulk. In recent years, more refined arguments have yielded inequalities in which other features of the spectrum are considered, see [YGL18] for an example or [BSS10] for a survey.

Regarding the necessity part of Question (2), to the best of our knowledge, the only existing negative result regarding detection of outliers is the one given in the recent work [SEAR18]. There, a query complexity bound was proven for any algorithm that is allowed to make queries of matrix-vector products, which in particular applies to the Lanczos algorithm.

1.4 Question (2): Our contributions

In the present paper we study the Lanczos algorithm in the context of approximation of outliers and answer the necessity part of Question (2). That is, we show that if run for at most $k\leq c\log n$ iterations, the Lanczos algorithm fails to approximate outliers with overwhelming probability. Thus, in essence we provide a lower bound on the number of iterations required for accuracy. As in our contribution for Question (1), the constant $c$ depends only on the equidistribution of the spectrum.

To give some rough context, the result in [SEAR18] discussed above shows that if the empirical spectral distribution of a matrix is close to the semicircle distribution plus an outlying “spike,” any algorithm in their class will fail to identify the spike with overwhelming probability, unless given at least $c\log n$ queries. In contrast, our result applies exclusively to the Lanczos algorithm, but shows that outliers are missed for a far more general class of measures than just the semicircle.

In order to analyze asymptotic behavior, we adopt a framework similar to that used in [Kui00] and [Bec00], in which a sequence of Hermitian matrices $A_{n}$ with convergent spectra was considered. These papers studied the behavior of the Lanczos algorithm in the regime of $\Theta(n)$ iterations.

To show that the Lanczos algorithm misses outliers when run for at most $c\log n$ iterations, we use the elementary theory of orthogonal polynomials and standard techniques in high-dimensional probability. Roughly speaking, using a variational principle, we show that for small enough $k$, the roots of the $k$th orthogonal polynomial with respect to a certain random measure are contained in a small blow-up of the convex hull of the bulk of the true spectrum. See Theorem 2.8 or Proposition 2.10 for a precise statement and Figure 2 for an illustration.

1.5 Our result on the locations of the output

One may ask if finer statements about the location of the Jacobi coefficients and Ritz values can be made. Previously, tools from potential theory have been used to answer this question in the regime of $k=\Theta(n)$ iterations [Bec00, Kui00, Kui06]. In the regime of $k$ fixed as $n\to\infty$, a deterministic convergence result for orthogonal polynomials [Gau68, Theorem 4] can be used to show that the Ritz values converge almost surely to the roots of the $k$th orthogonal polynomial of the limiting eigenvalue distribution; see Remark 2.14 for details.

In the present work we use determinantal formulas for orthogonal polynomials and concentration of measure results to locate the Jacobi coefficients and Ritz values in the regime of $k=O(\sqrt{\log n})$ iterations. In particular, we prove that the Ritz values concentrate around the roots of the $k$th orthogonal polynomial for the limiting eigenvalue distribution. See Figure 3 for an illustration. Moreover, also when $k=O(\sqrt{\log n})$, we show that the Jacobi matrix obtained after $k$ iterations is concentrated around the $k$th Jacobi matrix of the limiting measure.

These results may be of particular relevance in applications where an infinite-dimensional operator is discretized with the goal of computing its density. In essence, Theorem 2.12 below states that in this situation the first iterations of the Lanczos algorithm are an accurate approximation of the true Jacobi coefficients of the spectral measure of the infinite-dimensional operator, and hence the procedure is giving valuable information for recovering the limiting measure.

1.6 Organization of the paper

The article is organized as follows. In Section 2, we review the classical background of the Lanczos procedure and orthogonal polynomials and formally state our main theorems. In Section 3, we develop machinery that in Section 4 will be used to prove concentration for the output of the Lanczos algorithm. In Section 5, we prove our complementary results about the location of the Ritz values and Jacobi coefficients. Finally, in Section 6 we discuss further research directions that may be of interest.

2 Preliminaries and statements of theorems

Throughout this paper only elementary facts about orthogonal polynomials are used. For the reader’s convenience in Section 2.1 we include a concise survey of the results that will be used in what follows. Chapter 2 in [Sze39] and Chapters 2 and 3 in [Dei99] are introductory references containing these results.

In order to establish context and notation, in Section 2.2 we describe the Lanczos algorithm and its interpretation in terms of orthogonal polynomials. Some standard references for this matter are Chapter 6 in [TBI97] and Chapter 6 in [Saa11].

In Section 2.3 we introduce the framework in which this paper is developed and formally state the main contributions of our work.

In this paper we use the following notation. We use $\rightarrow_{P}$ to denote convergence in probability. For a sequence of events $E_{n}$, we say $E_{n}$ occurs with overwhelming probability if $\textbf{P}[E_{n}]\geq 1-C\exp\{-n^{c}\}$ for some $c,C>0$. For an $n\times n$ matrix $A$ with eigenvalues $\lambda_{i}$, we say that the empirical spectral distribution of $A$ is the atomic probability measure $\frac{1}{n}\sum_{i=1}^{n}\delta_{\lambda_{i}}$, where $\delta_{x}$ denotes the Dirac mass at $x$. We also let $\|A\|$ denote the spectral norm of $A$ (i.e. the $\ell^{2}\to\ell^{2}$ operator norm). The notation $\|\cdot\|$ applied to a vector will refer to the standard Euclidean norm. We will let $\mathbb{S}^{n-1}$ denote the unit sphere in $\mathbb{R}^{n}$ and denote the uniform probability measure on $\mathbb{S}^{n-1}$ by $\mathrm{Unif}(\mathbb{S}^{n-1})$. Finally, we will let $\mathrm{Kol}(\cdot,\cdot)$ denote the Kolmogorov distance between two measures, namely,

[TABLE]

2.1 Orthogonal polynomials

For now, let $\mu$ be a finite Borel measure on $\mathbb{R}$ and assume that its support, which we denote as $\mathrm{supp}(\mu)$, is compact and has infinitely many points. The set of square integrable functions $L^{2}(\mathbb{R},d\mu)$ becomes a Hilbert space when endowed with the inner product

[TABLE]

The hypothesis that $|\mathrm{supp}(\mu)|=\infty$ implies that the monomials $\{1,x,x^{2},\dots\}$ are linearly independent in $L^{2}(\mathbb{R},d\mu)$. Hence, we can use the Gram-Schmidt procedure to obtain an infinite sequence of polynomials $p_{k}(x)$ with $\mathrm{deg}(p_{k}(x))=k$ and

[TABLE]

The leading coefficient of $p_{k}(x)$ is a quantity of interest in this paper and will be denoted by $\gamma_{k}$. We will denote the monic orthogonal polynomials by $\pi_{k}(x)$. That is, $\pi_{k}(x)=\gamma_{k}^{-1}p_{k}(x)$ and clearly

[TABLE]

Since $\pi_{k}(x)$ is orthogonal to all polynomials with degree less than $k$, the polynomial $x^{k}-\pi_{k}(x)$ is the orthogonal projection of $x^{k}$ onto the span of $\{1,\dots,x^{k-1}\}$. Hence,

[TABLE]

where $\Gamma_{k}$ denotes the space of monic polynomials of degree $k$.

Favard’s theorem ensures that there is a sequence of real numbers $\alpha_{k}$ and a sequence of positive real numbers $\beta_{k}$ such that the following three-term recurrence holds:

[TABLE]

It is clear from the three-term recurrence that the following identity holds:

[TABLE]

These so-called Jacobi coefficients $\alpha_{k}$ and $\beta_{k}$ encode all the information of the measure $\mu$. In fact, since the Stieltjes transform of $\mu$ has a continued fraction expansion in terms of its Jacobi coefficients, knowing the few first elements in these sequences allows one to approximate the measure. See Chapter 4.3 in [Dei99] for an example.

We denote by $J_{k}$ the $k\times k$ Jacobi matrix of $\mu$; that is, $J_{k}$ is the tridiagonal symmetric matrix with $(J_{k})_{ii}=\alpha_{i-1}$ and $(J_{k})_{i+1,i}=(J_{k})_{i,i+1}=\beta_{i-1}$. It is a standard fact that $\pi_{k}(x)=\det(xI-J_{k})$ and that in particular, the roots of $p_{k}(x)$ are exactly the eigenvalues of $J_{k}$, which are real since $J_{k}$ is symmetric.

Another object of importance in this theory is the Hankel matrix of a measure. We will denote $M_{k}$ the $(k+1)\times(k+1)$ Hankel matrix of $\mu$; in other words, if $m_{i}$ denotes the $i$th moment of $\mu$, then $(M_{k})_{ij}=m_{i+j-2}$ for every $1\leq i,j\leq k+1$. From the elementary theory it is known (see [Dei99], Section 3.1) that if we define $D_{k}=\det M_{k}$, then

[TABLE]

where we define $D_{-1}=1$. Note that the second identity in (2.3) implies

[TABLE]

Moreover, if $\tilde{M}_{k}(x)$ denotes the matrix obtained by replacing the last row of $M_{k}$ by the row $(1\>x\>x^{2}\>\cdots\>x^{k})$, we have the following useful identity:

[TABLE]

Note that in the case in which $\mathrm{supp}(\mu)$ has $n$ points, for $n$ a positive integer, the set of monomials $\{1,x,x^{2},\dots\}$ is not linearly independent in $L^{2}(\mathbb{R},d\mu)$. Moreover, the Gram-Schmidt procedure stops after $n$ iterations, and hence it only makes sense to talk about the orthogonal polynomials $p_{k}(x)$ for $k\leq n-1$. However, sometimes it is convenient to define the $n$th monic orthogonal polynomial as the unique monic polynomial of degree $n$ whose roots are the elements of $\mathrm{supp}(\mu)$. In this case, the facts mentioned previously still hold for $k\leq n$.

2.2 The Lanczos algorithm

We understand the Lanczos algorithm as a randomized procedure that takes three inputs: an $n\times n$ Hermitian matrix $A$, a random vector $u$ distributed uniformly in $\mathbb{S}^{n-1}$, and an integer $1\leq k\leq n$. Then, the procedure outputs a $k\times k$ symmetric tridiagonal matrix $J_{k}$ whose diagonal entries will be denoted by $\alpha_{i}$ for $i=0,\dots,k-1$ and whose subdiagonal and superdiagonal entries will be denoted by $\beta_{i}$, for $i=0,\dots,k-2$. The eigenvalues of $J_{k}$ are called the Ritz values and we will usually denote them as $r_{1}\geq\cdots\geq r_{k}$. On the other hand, the eigenvectors of $J_{k}$ give rise (after an orthonormal change of basis determined by the $v_{j}$ below) to the Ritz vectors, that is, the approximations for the eigenvectors of $A$. Algorithm 2.1 below describes how the procedure generates the Jacobi coefficients $\alpha_{i}$ and $\beta_{i}$.

This algorithm has a natural interpretation in terms of orthogonal polynomials. To every $u\in\mathbb{S}^{n-1}$ we can associate a measure supported on the spectrum of $A$ as follows. Let $\lambda_{1}\geq\cdots\geq\lambda_{n}$ be the eigenvalues of $A$ and $u_{1},\dots,u_{n}$ be the coordinates of $u$ when writen in the eigenbasis of $A$. We define the probability measure

[TABLE]

In the language of functional analysis, $\mu^{u}$ is the spectral measure of the operator $A$ induced by the vector state $u$; that is, $\langle f(A)u,u\rangle=\int f(x)\,d\mu^{u}(x)$ for all (say) polynomials $f$. Note that the expectation of the random measure $\mu^{u}$ is just the empirical spectral distribution of $A$, namely,

[TABLE]

It is not hard to see that if $p_{j}(x)$ are the orthogonal polynomials with respect to $\mu^{u}$, then $v_{j}=p_{j}(A)u$. Hence, the coefficients $\alpha_{j}$ and $\beta_{j}$ outputed by the Lanczos algorithm are the Jacobi coefficients of the measure $\mu^{u}$, and the Ritz values after $k$ iterations are the roots of $p_{k}(x)$.

As a last remark, observe that the output of Algorithm 2.1 scales linearly with $A$. Hence, to simplify notation, in some of the proofs below we will start by assuming that $\lVert A\rVert=1$.

2.3 Statement of results

In the remainder of the paper the input matrix, which is assumed to be Hermitian, will be fixed and denoted by $A$. We will use $n$ to denote the dimension of $A$ and usually the number of iterations of the procedure will be denoted by $k$. Note that the Jacobi coefficients $\alpha_{i}$ and $\beta_{i}$ are assigned during the $i$th iteration of the algorithm and are unchanged during future iterations.

Since for our analysis it is necessary to compare outputs of the algorithm resulting from different input vectors $u\in\mathbb{S}^{n-1}$, we will stress this dependence by viewing the respective quantities as a function of $u$ and denoting them by $\alpha_{i}(u),\beta_{i}(u),r_{i}(u),\gamma_{i}(u)$, $p_{k}^{u}(x)$, $v_{i}(u)$, and $J_{k}(u)$. Depending on the context, these quantities will also be thought of as random variables, random polynomials, random vectors, and random matrices, respectively.

Most of our results make a technical assumption, which we have named equidistribution, about the geometry of the spectrum of $A$. We define and motivate this concept below.

Equidistribution. The behaviour of the Lanczos algorithm depends on the spectrum of the input matrix $A$, and what might be true for typical matrices can fail for particular choices of $A$. It is hence challenging, when trying to obtain general theoretical statements, to identify a quantifiable feature of the spectrum that can serve as an assumption but is not too restrictive. This has been done in different ways in previous work, most notably in the seminal work of Saad [Saa80]—see [Saa11, Section 6.6] for a succinct exposition. For example, many of the results in [Saa80] are stated in terms of the quantities

[TABLE]

where $\lambda_{1}\geq\cdots\geq\lambda_{n}$ are the eigenvalues of $A$ and $\mathcal{P}_{k-1}^{(i)}$ denotes the set of all polynomials of degree not exceeding $k-1$ and satisfying $p(\lambda_{i})=1$. In the present work we use a more geometric notion to state our results.

Definition 2.1** (Equidistribution).**

Let $\Lambda$ be any finite set of $n$ real numbers. Let $\delta$ and $\omega$ be positive real numbers and let $j$ be a natural number. We say that $\Lambda$ is $(\delta,\omega,j)$-equidistributed if for any finite set $T$ of at most $j$ real numbers it holds that

[TABLE]

Intuitively, the spectrum is equidistributed if it is not grouped in a small number of tight clusters (see Examples 2.2 and 2.3 below). As we will show in Section 4, the family of well equidistributed point sets includes, but is not limited to, those sets obtained by discretizing an absolutely continuous distribution.

Example 2.2**.**

Let $\Lambda$ be the set of $n$ equally spaced points from $1/n$ to $1$, inclusive. This represents a discretization of the uniform measure $\mu=\mathrm{Unif}([0,1])$. In Section 4.1 we will show that for $j\leq\frac{n}{16}$, the set $\Lambda$ is $(\delta,\omega,j)$-equidistributed for $\delta=1/4$ and $\omega=4e^{-2}$.

Example 2.3**.**

Now consider a set (or multiset) $\Lambda$ of $n>0$ points grouped in $m$ equally spaced small clusters. To make this precise, fix two parameters $\varepsilon,g>0$ and consider $-1=a_{1}\leq b_{1}<a_{2}\leq b_{2}<\cdots<a_{m}\leq b_{m}=1$ such that for every $i=1,\dots,m$ we have $b_{i}-a_{i}=\varepsilon$ and $a_{i+1}-b_{i}=g$. We think of $\varepsilon$ as small with respect to $g$ and of $m$ as small with respect to $n$. If $\Lambda\subset\bigcup_{i=1}^{m}[a_{i},b_{i}]$ with $|\Lambda\cap[a_{i},b_{i}]|\geq\lfloor\frac{n}{m}\rfloor$ for every $i=1,\dots,m$, then $\Lambda$ is $(\frac{m-j}{m},g,j)$-equidistributed and $g\approx 2/m$.

Note that in this case we have good equidistribution parameters unless $j\approx m$. In Section 4 we give a generalization of this assertion in Observation 4.10.

Concentration of the output. We state our main result about concentration of the first few Jacobi coefficients in terms of the equidistribution parameters $\delta,\omega$.

Theorem 2.4** (Concentration of Jacobi coefficients after $i$ iterations).**

Suppose the spectrum of $A$ is $(\delta,\omega,i)$-equidistributed for some $\delta,\omega>0$ and $i\in\mathbb{N}$. Let $\tilde{\alpha}_{i}$ and $\tilde{\beta}_{i}$ denote the medians of the Jacobi coefficients $\alpha_{i}(u)$ and $\beta_{i}(u)$, respectively. Then for all $t>0$, the probabilities $\mathbf{P}[|\alpha_{i}(u)-\tilde{\alpha_{i}}|>t\|A\|]$ and $\mathbf{P}[|\beta_{i}(u)-\tilde{\beta_{i}}|>t\|A\|]$ are both bounded above by

[TABLE]

To clarify the advantages and limitations of the above result we include some remarks.

Remark 2.5**.**

The equidistribution parameters $\delta,\omega$ appearing in the above theorem are typically quite moderate in magnitude and are easy to compute if one can obtain explicit bounds for certain integrals with respect to the spectral distribution of $A$. Note that $\omega\leq\|A\|$ (by taking $T=\{0\}$ in Definition 2.1) and that $\omega$ scales linearly with $A$. As a result, $\omega/\|A\|$ is typically of constant size independent of $n$ in applications. Since $\omega/\|A\|<1$, Theorem 2.4 yields concentration for $i$ at most logarithmic in $n$. Besides Examples 2.2 and 2.3 above, in Section 4.1 we give more examples and a detailed discussion on how to compute these parameters.

Remark 2.6**.**

Ultimately because the Jacobi coefficients $\alpha_{i}(u),\beta_{i}(u)$ turn out to be Lipschitz only on a proper subset of the sphere, we only obtain concentration of the coefficients about their medians, not their means.

Remark 2.7**.**

Using the same techniques one can prove a result analogous to Theorem 2.4 in the case where $A$ is not Hermitian, and even not normal. In the non-Hermitian case, the Lanczos algorithm is called the Arnoldi algorithm and has similar applications to those of the Lanczos algorithm.

In Section 4 we show how Theorem 2.4 can be used to obtain concentration, in the same regime of number of iterations, of the Ritz values and Ritz vectors. See Propositions 4.16 and 4.18 for precise statements.

The rest of the paper focuses on studying the location of the outputs of the Lanczos algorithm.

Undetected outliers of the spectrum. In Section 5.1, we show that if $k$ is a certain fraction of $\log n$, the Ritz values obtained after $k$ iterations are contained in a small blow-up of the convex hull of the bulk of the spectrum of $A$. This complements classical guarantees which show that for some multiple of $\log n$, say, $K$, the Lanczos algorithm approximates with high accuracy the outliers of the spectrum of $A$ when $K$ iterations are performed. Our results are quantitative and use our notion of equidistribution.

Theorem 2.8**.**

Suppose the spectrum of $A$ is $(\delta,\omega,j)$-equidistributed for some $\delta,\omega>0$ and $j\in\mathbb{N}$. Let $M$ be the diameter of the spectrum of $A$. Let $R$ be a real number and let $0<c<1/2$, and suppose there are at most $m\leq\min\{0.02n,2n^{\alpha}\}$ “outliers,” eigenvalues of $A$ lying above $R$, for some $\alpha<1-c$. Let $g=\max_{1\leq i\leq n}\{\lambda_{i}-R\}$ and let $\kappa>0$. Then for up to

[TABLE]

iterations, the probability that the top Ritz value exceeds $R+\kappa$ is at most

[TABLE]

for $n>e^{\frac{1}{1-c-\alpha}}$.

The strength of the above result might be obscured by the appearance of several unintuitive parameters. For the reader’s benefit we include an example below, and to provide a slightly different perspective, we include an asymptotic version of the above result, namely, Proposition 2.10.

Example 2.9**.**

Let $n>0$ and let $A$ be a matrix whose spectrum consists of $n-1$ equally spaced points from $2/n$ to $1$ inclusive, together with an outlier of value 1.1 (compare with Figure 2). In Section 4.1 we will show that for $j\leq n/16$ the spectrum of $A$ is $(1/4,4e^{-2},j)$-equidistributed.

In order to apply Theorem 2.8, we also note that in this case $M=1.08$, $m=1$, and $g=10^{-1}$. Take $\kappa=10^{-4}$. Then, for any $0<c<1/2$, the Ritz values of the Lanczos algorithm on $A$ after $\lfloor\frac{7c}{10}\log n-7/2\rfloor$ iterations will be contained in the interval $[2/n,1+10^{-4}]$ with overwhelming probability.

Proposition 2.10**.**

Let $(A_{n})_{n=1}^{\infty}$ be a sequence of $n\times n$ Hermitian matrices with uniformly bounded norm. Assume their empirical spectral distributions $\mu_{n}$ converge in distribution to a measure $\mu$ with nontrivial absolutely continuous part, and further assume $\text{Kol}(\mu_{n},\mu)=O(1/\log n)$. Suppose there exists $m\in\mathbb{N}$ such that each $A_{n}$ has at most $m$ eigenvalues (“outliers”) greater than $R$, where $R$ denotes the right edge of the support of $\mu$.

Then there exists $c>0$ such that for every $\kappa>0$, the Ritz values of Lanczos applied to $A_{n}$ after $c\log n$ iterations are bounded above by $R+\kappa$ with overwhelming probability for $n$ sufficiently large (depending on how small the gap $\kappa$ is chosen.)

Remark 2.11**.**

Suppose all eigenvalues of $A$ lie in $[-1,1]$ except for one outlier $\lambda_{1}=1+\varepsilon$. Suppose we wish for a Ritz vector $p(A)u=\sum_{i=1}^{n}p(\lambda_{i})u_{i}$ to approximate $u_{1}$, the top eigenvector. Then $p(\lambda_{1})$ must be polynomially larger than $p(\lambda_{2}),\dots,p(\lambda_{n})$, since $|u_{i}|\sim 1/\sqrt{n}$ for all $1\leq i\leq n$ with high probability. If one further imposes that $p(\lambda_{1})$ exceeds $p$ by a polynomial factor on the entire interval $[-1,1]$, then the degree of $p$ (and hence the number of iterations $k$) must be at least $O(\log n/\sqrt{\varepsilon}),$ by the Markov brothers’ inequality and properties of Chebyshev polynomials. However, if the $\lambda_{i}$ are in $j$ tight clusters for $j$ small, one could certainly have $p(\lambda_{2}),\dots,p(\lambda_{n})$ small and $p(\lambda_{1})$ large for some $p$ of degree $j$. Thus, an assumption like equidistribution is natural for Theorem 2.8.

**Location of Jacobi coefficients and Ritz values ** We will now work on the setting of Proposition 2.10, that is, we consider a probability measure $\mu$ and a sequence of matrices $A_{n}$ whose empirical spectral distributions $\mu_{n}$ converge to $\mu$. We give a result about the locations of the Ritz values and Jacobi coefficients when at most $d\sqrt{\log n}$ iterations are performed, with $d$ depending only on $\mu$ and the speed of convergence of the sequence $\mu_{n}$. Essentially, we show that in this regime the Jacobi matrix after $k$ iterations is sharply concentrated around the $k$th Jacobi matrix of the measure $\mu$.

Theorem 2.12** (Location of Jacobi coefficients).**

Let $(A_{n})_{n=1}^{\infty}$ be a sequence of $n\times n$ Hermitian matrices with uniformly bounded operator norm. Assume their empirical spectral distributions $\mu_{n}$ converge in distribution to a measure $\mu$ with nontrivial absolutely continuous part, and further assume $\text{Kol}(\mu_{n},\mu)=O(n^{-c})$ for some $c>0$.

Then there is a constant $d>0$ dependent on $\mu$ and $c$, such that for any sequence of integers $1\leq k_{n}\leq d\sqrt{\log n}$ we have

[TABLE]

where $J_{k_{n}}(u)$ denotes the Jacobi matrix output by the Lanczos algorithm applied to $A_{n}$ under the input $u\sim\mathrm{Unif}(\mathbb{S}^{n-1})$ after $k_{n}$ iterations, and where $J_{k_{n}}(\mu)$ is the $k_{n}$-th Jacobi matrix of the measure $\mu$.

Note that Theorem 2.12 may be of particular relevance in applications where an infinite-dimensional operator is discretized with the goal of computing its density. In essence, Theorem 2.12 states that, in this situation, the first iterations of the Lanczos algorithm are an accurate approximation of the true Jacobi coefficients of the measure $\mu$, and hence the procedure gives valuable information to recover the limiting measure.

From the above proposition, a standard application of the Weyl eigenvalue perturbation inequality yields the following proposition.

Proposition 2.13** (Location of the Ritz values).**

Using the same notation as in Theorem 2.12, let $\vec{r}_{k_{n}}(u)=(r_{1}(u),\dots,r_{k_{n}}(u))$, where $r_{1}(u)\geq\cdots\geq r_{k_{n}}(u)$ are the random Ritz values of the Lanczos algorithm after $k_{n}$ iterations are performed. Then under the assumptions in Theorem 2.12, we have that

[TABLE]

where $\vec{r}_{k_{n}}(\mu)$ is the vector whose entries are the roots of the $k_{n}$-th orthogonal polynomial with respect to $\mu$ in decreasing order.

It remains an open question if similar results can be obtained when $O(\log n)$ iterations are performed. See Section 6 for open questions and further research.

Remark 2.14**.**

For fixed $k$, [Gau68, Theorem 4] states that in the deterministic setting of a weakly convergent sequence of measures $\mu_{n}\to\mu$ with each $\mu_{n}$ supported on $n$ points, the Jacobi coefficients $\alpha_{i},\beta_{i}$ of $\mu_{n}$ for $i\leq k$ converge to those for $\mu$ (and therefore, the same holds for the $i$th orthogonal polynomial and its roots.) To apply this to the Lanczos algorithm, one must use the fact that for independent initial vectors $u_{n}\in\mathbb{S}^{n-1}$, almost surely the measures $\mu_{n}^{u_{n}}$ (defined in 2.6) converge weakly to $\mu$. One therefore obtains the same convergence as in the deterministic setting, almost surely, for the random $\alpha_{i},\beta_{i}$ output by Lanczos. To make the convergence quantitative, one must assume something about the rate of convergence of the $\mu_{n}$ in the hypothesis; this is what is done in Theorem 2.12. Our result holds in the wider range $k=O(\sqrt{\log n})$. We also trade almost-sure convergence for a polynomial rate of convergence holding with overwhelming probability; this rate appears in the proof.

3 Applying the local Lévy lemma

3.1 Strategy

The well known Lévy lemma states, in a quantitative way, that if $f:\mathbb{S}^{n-1}\to\mathbb{R}$ is a Lipschitz function, then $f(u)$ is a random variable concentrated around its median. See Chapter 5.1 in [Ver18] for a detailed discussion. In this direction, the main obstacle for showing concentration of the random variables $\alpha_{i}(u)$ and $\beta_{i}(u)$ is that the functions $\alpha_{i},\beta_{i}:\mathbb{S}^{n-1}\to\mathbb{R}$ are not Lipschitz on the entire sphere. However, we will be able to show that these functions are Lipschitz in a large region of the sphere, which is a common idea in geometric functional analysis. We will use a local version of Lévy’s lemma, which is recorded as Corollary 5.35 in [AS17], and which we restate below with explicit universal constants.

Lemma 3.1** (Local Lévy lemma).**

Let $\Omega\subset\mathbb{S}^{n-1}$ be a subset of measure larger than $3/4$. Let $f:\mathbb{S}^{n-1}\to\mathbb{R}$ be a function such that the restriction of f to $\Omega$ is Lipschitz with constant $L$ (with respect to the geodesic metric on the sphere). Then, for every $\varepsilon>0$,

[TABLE]

where $\tilde{f}$ is the median of $f(u)$ and where $u\sim\mathbb{S}^{n-1}$.

As the function $f$ is allowed to blow up on the small subset $\mathbb{S}^{n-1}\setminus\Omega$, one cannot expect a similar result to hold for concentration around the mean.

In order to identify the correct region of the sphere in which the functions $\alpha_{i}$ and $\beta_{i}$ are Lipschitz, we need a local version of the notion of Lipschitz constant. In what might be a slight departure from standard definitions, we will define local Lipschitz continuity as follows.

Definition 3.2**.**

Let $(X_{1},d_{1})$ and $(X_{2},d_{2})$ be metric spaces. A function $f:X_{1}\to X_{2}$ is said to be locally Lipschitz continuous with constant $c$ at $x_{0}\in X_{1}$ if for every $c^{\prime}>c$ there is a neighborhood $U\subset X_{1}$ of $x_{0}$ such that

[TABLE]

Remark 3.3**.**

For $f$ defined on an open subset of $\mathbb{S}^{n-1}$, we have that $f$ is locally Lipschitz continuous with constant $c$ with respect to the geodesic metric if and only if it is locally Lipschitz continuous with the same constant with respect to the Euclidean (“chordal”) metric.

It is obvious that if a function is locally Lipschitz with constant $c$ on every point of a convex set, then the function is globally Lipschitz on the set with the same constant $c$. However, if the convexity assumption is dropped, a similar conclusion is not guaranteed in general and in order to obtain a global Lipschitz constant the geometry of the set should be analyzed.

Definition 3.4**.**

Let $K>0$ and $(X,d)$ be a metric space. We say that $S_{1}\subset X$ is $K$-connected in $S_{2}$ with $S_{1}\subset S_{2}\subset X$ if for every $x,y\in S_{1}$ there is a rectifiable Jordan arc $\alpha:[0,1]\to S_{2}$ with $\alpha(0)=x$ and $\alpha(1)=y$, such that the length of the trace of $\alpha$ is less than or equal to $Kd(x,y)$.

Now that we have introduced the notion of $K$-connected set we can generalize what we observed for convex sets.

Lemma 3.5**.**

Let $(X_{1},d_{1})$ and $(X_{2},d_{2})$ be metric spaces. Assume that $S_{1}\subset X_{1}$ is $K$-connected in $S_{2}\subset X_{1}$ and let $f:X_{1}\to X_{2}$ satisfy that for every $x_{0}\in S_{2}$, $f$ is locally Lipschitz at $x_{0}$ with constant $c$. Then $f$ is globally Lipschitz on $S_{1}$ with constant $cK$.

Proof.

Fix $x,y\in S_{1}$ and $\varepsilon>0$. We will show that $d_{2}(f(x),f(y))\leq(c+\varepsilon)Kd_{1}(x,y)$. Consider a rectifiable Jordan arc $\alpha:[0,1]\to X_{1}$, such that $\alpha(0)=x$, $\alpha(1)=y$, $\alpha([0,1])\subset S_{2}$ and the length of $\alpha$ is at most $Kd_{1}(x,y)$.

Since the trace of $\alpha$ is contained in $S_{2}$, for every $w\in\alpha([0,1])$ we can take an open ball $U_{w}$ containing $w$ such that $f$ is $(c+\varepsilon)$-Lipschitz on $U_{w}$. Moreover, observe that since $\alpha$ is continuous and injective, for every $w\in\alpha([0,1])$ we can take $U_{w}$ small enough such that $\alpha^{-1}(U_{w})$ is connected and hence an open interval in $[0,1]$.

By compactness of $\alpha([0,1])$ we may take $w_{1},\dots,w_{n}\in\alpha([0,1])$ such that $\{U_{w_{i}}\}_{i=1}^{n}$ is a minimal cover for $\alpha([0,1])$. Now, since each $\alpha^{-1}(U_{w_{i}})$ is connected, and the cover is minimal, we have that $\alpha^{-1}(U_{w_{i}})\cap\alpha^{-1}(U_{w_{i+1}})\neq\emptyset$ for every $1\dots,n-1$.

Furthermore, we will now see that we can modify the sequence of $w_{i}$ such that $w_{i+1}\in U_{w_{i}}$ for every $i=1,\dots,n-1$. Assume that this does not hold and let $i$ be the smallest index for which $w_{i+1}\notin U_{w_{i}}$. Now take some $t\in\alpha^{-1}(U_{w_{i}})\cap\alpha^{-1}(U_{w_{i+1}})$ and define $w^{\prime}=\alpha(t)$. We construct a new sequence $\tilde{w}_{1},\dots,\tilde{w}_{n+1}\in\alpha([0,1])$ by taking $\tilde{w}_{j}=w_{j}$ for $j<i$, $\tilde{w}_{i}=w^{\prime}$, $\tilde{w}_{j+1}=w_{j}$ for $j\geq i$, and $U_{\tilde{w}_{i}}$ to be equal to $U_{w_{i+1}}$. Observe that for the new sequence of points $(\tilde{w}_{i})_{i=1}^{n+1}$ in $\alpha([0,1])$ and sequence of open balls $U_{\tilde{w}_{i}}$ it holds that $\tilde{w}_{j+1}\in U_{\tilde{w}_{j}}$ for all $j\leq i$. By iterating this process we will obtain a finite sequence with the desired property. So, in what follows we can assume without loss of generality that $w_{i+1}\in U_{w_{i}}$ for every $i=1,\dots,n-1$. We then will have

[TABLE]

Using the triangle inequality and the fact that $\sum_{i}d_{1}(w_{i},w_{i+1})$ is bounded by the length of the trace of $\alpha$ the result follows. ∎

In the following section the local Lipschitz constants of the functions $\alpha_{i}(u)$ and $\beta_{i}(u)$ are shown to be related to the orthogonal polynomials of the measure $\mu^{u}$.

3.2 Local Lipschitz constants for Jacobi coefficients

As can be seen from Algorithm 2.1, the dependence of the quantities $\alpha_{i}(u)$, $\beta_{i}(u)$, and $v_{j}(u)$ on $u$ is highly nonlinear, which makes it complicated to show that such quantities are stable under perturbations of the input vector $u$. Here we exploit the fact that during every iteration of the Lanczos algorithm only locally Lipschitz operations are performed. The analysis of the compound effect of iterating the procedure yields a bound on the local Lipschitz constant of the quantities of interests. This bound is exponential in the number of iterations, which is enough to obtain concentration results when $O(\log(n))$ iterations are performed. In what follows, recall that $\gamma_{i}(u)$ denotes the leading coefficient of the $i$th orthonormal polynomial with respect to the measure $\mu^{u}$ defined in (2.6).

Proposition 3.6**.**

Fix $\tilde{u}\in\mathbb{S}^{n-1}$ and let $v_{j}(u)$ be as in Algorithm 2.1. Then, for any $0\leq j\leq n-1$, the functions $v_{j}(u)$ are locally Lipschitz at $\tilde{u}$ with constant $(4\lVert A\rVert)^{j}\gamma_{j}(\tilde{u})$.

Proof.

We proceed by induction. For $j=0$, recall $v_{0}(u)=u$ and $\gamma_{0}(\tilde{u})=1$; the statement follows. Now assume the proposition is true for some $j\geq 0$. For every $x\in\mathbb{S}^{n-1}$ denote $W_{x}=\mathrm{span}\{v_{0}(x)=x,v_{1}(x),\dots,v_{j}(x)\}$ and for any subspace $W\leq\mathbb{R}^{n}$ by $\mathrm{Proj}_{W}$ we mean the orthogonal projection onto $W$.

Take $x,y\in\mathbb{S}^{n-1}$ in a neighborhood $\mathcal{U}$ of $\tilde{u}$ to be determined and note that

[TABLE]

From the induction hypothesis we have that, for any $\varepsilon>0$, we can choose $\mathcal{U}$ small enough so that

[TABLE]

On the other hand, from Algorithm 2.1 it follows that $\beta_{i}(\tilde{u})\leq\lVert A\rVert$ for every $i=0,\dots,n-1$, so in view of (2.2), the $\lVert A\rVert^{i}\gamma_{i}(\tilde{u})$ form an increasing sequence. It then follows that

[TABLE]

For any unit vector $w$, by the triangle inequality, we have that

[TABLE]

and we can bound each term on the right-hand side of (3.3) as follows:

[TABLE]

Hence, adding over $i$ we obtain

[TABLE]

which implies that $\lVert\mathrm{Proj}_{W_{x}}-\mathrm{Proj}_{W_{y}}\rVert\leq\frac{2}{3}\cdot 4^{j+1}\lVert A\rVert^{j}\gamma_{j}(\tilde{u})\lVert x-y\rVert$ and hence

[TABLE]

Putting together inequalities (3.2), (3.2), and (3.4), we get for any $x,y\in\mathcal{U}$ that

[TABLE]

With this we have established that the function $u\mapsto\mathrm{Proj}_{W_{u}^{\perp}}(Av_{j}(u))$ is locally Lipschitz at $\tilde{u}$ with constant $(4\lVert A\rVert)^{j+1}\gamma_{j}(\tilde{u})$. Now consider the function $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ defined by $f(x)=x/\lVert x\rVert$. It is easy to show that for any $x_{0}\neq 0$, $f$ is locally Lipschitz at $x_{0}$ with constant $1/\lVert x_{0}\rVert$. Now recall that by definition $\beta_{j}(\tilde{u})=\lVert\mathrm{Proj}_{W_{\tilde{u}}^{\perp}}(Av_{j}(\tilde{u}))\rVert$. Since the composition of locally Lipschitz functions is locally Lipschitz with the constant being the product of the constants of each of the functions in the composition, we have that the function

[TABLE]

is locally Lipschitz at $\tilde{u}$ with constant $\frac{(4\lVert A\rVert)^{j+1}\gamma_{j}(\tilde{u})}{\beta_{j}(\tilde{u})}=(4\lVert A\rVert)^{j+1}\gamma_{j+1}(\tilde{u})$, where this equality follows from (2.2). ∎

Proposition 3.7**.**

For any $0\leq j\leq n-1$ and any $\tilde{u}\in\mathbb{S}^{n-1}$, the function $\alpha_{j}(u)$ is locally Lipschitz at $\tilde{u}$ with constant $\frac{1}{2}\cdot(4\lVert A\rVert)^{j+1}\gamma_{j}(\tilde{u})$, while $\beta_{j}(u)$ is locally Lipschitz at $\tilde{u}$ with constant $(4\lVert A\rVert)^{j+1}\gamma_{j}(\tilde{u})$.

Proof.

We will use the same notation as in Proposition 3.6. Recall from Algorithm 2.1 that $\alpha_{j}(u)=\langle Av_{j}(u),v_{j}(u)\rangle$. Note that the local Lipschitz constant of the function $u\mapsto Av_{j}(u)$ is obtained by multiplying the local Lipschitz constant of $v_{j}(u)$ by $\lVert A\rVert$ . Then, for any $\varepsilon$ we can pick $\mathcal{U}$ to be a small enough neighborhood of $\tilde{u}$ such that for any $x,y\in\mathcal{U}$ we have

[TABLE]

On the other hand, since $\beta_{j}(u)=\lVert\mathrm{Proj}_{W_{u}^{\perp}}(Av_{j}(u)))\rVert$ and we established in the proof of Proposition 3.6 that this function is locally Lipschitz with constant $(4\lVert A\rVert)^{j+1}\gamma_{j}(\tilde{u})$, the proof is concluded. ∎

Remark 3.8**.**

The local Lipschitz constants presented in the above statements can be improved; the term $4^{j}$ next to $\lVert A\rVert^{j}\gamma_{j}(\tilde{u})$ was chosen for the sake of exposition. Nevertheless, it seems complicated to show that the quantities $v_{j}(u)$ are locally Lipschitz at $\tilde{u}$ with a constant of the form $C_{j}\lVert A\rVert^{j}\gamma_{j}$ and $C_{j}$ subexponential. In any case, the term $\lVert A\rVert^{j}\gamma_{j}$ is typically exponential in $j$, so an improvement on $C_{j}$ would not yield an asymptotic improvement to the final result if the same level of generality is considered. However, as we point out in Section 6, sharpening our constants is of relevance for applications.

3.3 Incompressibility

In Section 4, we will see that our upper bounds for the local Lipschitz constants of the Jacobi coefficients go to infinity if $u$ becomes too close to a sparse vector, roughly speaking. So we only have a good local Lipschitz constant in a certain region of the unit sphere that avoids sparse vectors. In order to upgrade our local Lipschitz constant to a global Lipschitz constant, we must prove

1. (1)

that this region is large enough to apply the local Lévy lemma (Lemma 3.1) and 2. (2)

that this region is $K$-connected for a small enough $K$.

First we give this region a name. Loosely inspired by the compressed sensing literature (see, for example, [Ver09]), we say that a vector $u$ in $\mathbb{S}^{n-1}$ is $(\delta,\varepsilon)$-incompressible if each set of at least $\delta n$ coordinates carries at least $\varepsilon$ of its “$\ell^{2}$ mass.” Otherwise, we say that $u$ is $(\delta,\varepsilon)$-compressible. We denote the set of $(\delta,\varepsilon)$-incompressible vectors in $\mathbb{S}^{n-1}$ by $I_{n}(\delta,\varepsilon)$ and record the formal definition below.

Definition 3.9**.**

[TABLE]

For incompressible $u$ we prove an adequate bound on the local Lipschitz constant in Proposition 4.2. Fortunately, a uniform random unit vector $u$ is incompressible with high probability, as we will now show.

Proposition 3.10**.**

Let $u\in\mathbb{S}^{n-1}$ be a uniform random unit vector, and let $0<\varepsilon<\delta$. Then

[TABLE]

Corollary 3.11**.**

Let $u\in\mathbb{S}^{n-1}$ be a uniform random unit vector, and let $0<\delta\leq 1/50$. Then

[TABLE]

Proof.

Set $\varepsilon=\delta/2$ in Proposition 3.10. Note that $\varepsilon^{2}/8=\delta^{2}/32$ and $2\delta(1+\log 1/\delta)-(1/2)^{2}<-1/32$ for $0<\delta\leq 1/50$. ∎

The proof of the Proposition 3.10 consists of two parts. First, we prove a similar proposition where instead of the $u_{i}$ we have independent Gaussian random variables with the same variance $1/n$. We then use a coupling argument to conclude the desired bound for $u$ drawn uniformly from the unit sphere.

We will need upper and lower tail bounds on the $\chi^{2}$ distribution. One can get good enough bounds using the Chernoff method, but rather than develop these from scratch we will cite the following corollary of Lemma 1 from Section 4.1 of [LM00].

Lemma 3.12**.**

Let $Y$ be distributed as $\chi^{2}(k)$ for a positive integer $k$. Then the following upper and lower tail bounds hold for any $t\geq 0$:

[TABLE]

[TABLE]

Proof of Proposition 3.10.

Let $X_{1},\dots,X_{n}$ denote independent Gaussian random variables each with variance $1/n$, and let $X=(X_{1},\dots,X_{n})$. If we set $u=X/\|X\|$, then $u$ is uniformly distributed on the unit sphere; see e.g. [Mul59].

We seek to upper bound the probability of compressibility $\{u\not\in I_{n}(\delta,\varepsilon)\}$, which is the event that $\sum_{i\in S}u_{i}^{2}<\varepsilon$ for some subset $S$ of coordinates with $|S|\geq\delta n$. This event is contained in the union of the following two events:

$E$, the event that $\sum_{i\in S}X_{i}^{2}\leq 2\varepsilon$ for some $|S|\geq\delta n$, and 2. 2.

$F$, the event that $\sum_{i\in S}X_{i}^{2}\geq\varepsilon+\sum_{i\in S}u_{i}^{2}$ for some $|S|\geq\delta n$.

Indeed, if neither of these events hold, then for all $|S|\geq\delta n$ we have

[TABLE]

so $u$ is incompressible.

To upper bound the probability of $E$, we use the union bound over all sets of size $k=\lceil n\delta\rceil$:

[TABLE]

where in the last step we apply the lower tail bound in Lemma 3.12 with $t$ being the solution to $k-2\sqrt{kt}=2n\varepsilon$. To avoid the bookkeeping of ceiling and floor functions we use the extremely crude inequality $n\delta\leq k\leq 2n\delta$ (valid as long as $\delta n\geq 1$), which will suffice for our purposes:

[TABLE]

We now upper bound the probability of $F$:

[TABLE]

Since $Y=n\|X\|^{2}$ is distributed as $\chi^{2}(n)$, we may apply the upper tail bound in Lemma 3.12 with $t=n\varepsilon^{2}/8$ to obtain

[TABLE]

To conclude, we have $\mathbf{P}[u\not\in I_{n}(\delta,\varepsilon)]\leq\mathbf{P}[E]+\mathbf{P}[F]$, and substituting the bounds we just derived, we obtain the desired inequality. ∎

3.4 K-connectedness of the incompressible region

Having proven that the incompressible region $I_{n}(\delta,\varepsilon)$, where we have a good local Lipschitz constant, is almost the entire sphere, we now turn to proving that the region is $K$-connected for a small enough $K$.

One could try to show that any two points in $I_{n}(\delta,\varepsilon)$ can be connected by a short path contained in $I_{n}(\delta,\varepsilon)$, but for our purposes it is okay to let the path venture out into the larger region $I_{n}(4\delta,\varepsilon/\sqrt{2})$. When upgrading to a global Lipschitz constant, we will have to use the slightly worse upper bound for the local Lipschitz constant in this larger region, but this will still be good enough.

Proposition 3.13**.**

$I_{n}(\delta,\varepsilon)$* is $\sqrt{2/\varepsilon}$-connected in $I_{n}(4\delta,\varepsilon/\sqrt{2}).$*

Proof.

Let $x$ and $y$ be any two endpoints in $I_{n}(\delta,\varepsilon)$. The construction will proceed in two steps. First, we will construct a path from $x$ to $y$ in $\mathbb{R}^{n}$ consisting of $\lceil\delta^{-1}\rceil$ pairwise orthogonal line segments. Then we will project this path radially onto the unit sphere and show that the result indeed lies in $I_{n}(4\delta,\varepsilon/2)$ and has length at most $(2/\sqrt{\varepsilon})\|x-y\|$, which is at most $(2/\sqrt{\varepsilon})d(x,y)$, where $d$ denotes the geodesic distance on $\mathbb{S}^{n-1}$.

Roughly speaking, we will partition the coordinates of $x$ into $1/\delta$ blocks of $\delta n$ coordinates and move the entries of each block linearly from $x$ to $y$ in parallel, one block at a time.

Because basic quantities such as $1/\delta$ and $\delta n$ may not be integers, we will be content to split up $\mathbb{R}^{n}$ as the direct sum $\bigoplus_{i=1}^{m}\mathbb{R}^{n_{i}}$, where $\delta n\leq n_{i}\leq 2\delta n$ for all $i$.111This is possible as long as $n/2\geq\delta n\geq 1$, which will be true in our regime. Note also that this implies $m\geq\frac{1}{\delta}$. Similarly, for any vector $z\in\mathbb{R}^{n}$, we will write $z=\bigoplus_{i=1}^{m}z^{(i)}$, where $z^{(i)}\in\mathbb{R}^{n_{i}}.$

Now we may formally define the path $P_{i}$ to be the line segment

[TABLE]

and define $P$ to be the concatenation of the segments $P_{1},\dots,P_{m}$. The length of $P$ is

[TABLE]

by the Cauchy-Schwarz inequality. Also, $\|P(t)\|\geq\sqrt{\varepsilon/2\delta}$, because

[TABLE]

where we use that $x$ and $y$ are $(\delta,\varepsilon)$-incompressible.

Furthermore, note that $P$ lies inside the closed ball of radius $\sqrt{2}$, because for any $i$ and $t$,

[TABLE]

The path $P$ currently does not lie in the unit sphere, so we project it onto the unit sphere along radii to get our final path $P^{\prime}$. We now show that $P^{\prime}$ indeed lies in $I_{n}(4\delta,\varepsilon/\sqrt{2})$.

At this stage, we will dispense with the direct sum decomposition and use ordinary coordinates $z=(z_{1},\dots,z_{n})$.

Consider any set $S$ of at least $4\delta n$ coordinates, and consider any point $P_{i}(t)$ in our path $P$ (before projection). The $i$th block of coordinates is in motion, and all of the other coordinates are either frozen at their initial value (from $x$) or their final value (from $y$).

The $i$th block consists of at most $2\delta n$ coordinates. Besides these, there are at least $4\delta n-2\delta n=2\delta n$ remaining coordinates in our set $S$. At least $\delta n$ of them are from $x$ or at least $\delta n$ of them are from $y$. By incompressibility of $x$ and $y$, the sum of the squares of these $\delta n$ coordinates is at least $\varepsilon$.

After projecting onto the unit sphere, the sum of the same coordinates is still at least $\varepsilon/\sqrt{2}$, because as we saw, the original path had norm at most $\sqrt{2}$ at every point.

Finally, when projecting onto the unit sphere, the length of the path increases by at most a factor of $1/\sqrt{\varepsilon/2\delta}$, because as we saw earlier, originally each segment lay outside the smaller sphere of radius $\sqrt{\varepsilon/2\delta}$. The verification is an exercise in plane geometry (using the fact that $\tan\theta>\theta$ for $0<\theta<\pi/2$) and also follows from the arc length formula $d\textbf{s}=\sqrt{r^{2}+(dr/d\theta)^{2}}\,d\theta\geq r\,d\theta$.

Thus, finally, we have shown that the path $P^{\prime}$ is contained in $I_{n}(4\delta,\varepsilon/\sqrt{2})$ and has length at most

[TABLE]

∎

4 Concentration of the output

We now analyze the local Lipschitz constant for the entries $\alpha_{i}$ and $\beta_{i}$ of the Jacobi matrix. To simplify notation, in what follows we assume that $\lVert A\rVert=1$ by rescaling $A$. Recall that this will also rescale the Ritz values and Jacobi coefficients by a factor $1/\lVert A\rVert$.

By Corollary 3.7, the function $\alpha_{i}(u)$ has local Lipschitz constant $2\cdot 4^{i}\gamma_{i}(u)$, and $\beta_{i}(u)$ has local Lipschitz constant $4^{i+1}\gamma_{i}(u)$. Thus we are naturally led to the question of finding upper bounds for $\gamma_{k}(u)$. Recall that $\gamma_{k}(u)$ is defined as the leading coefficient of the $k$th orthogonal polynomial with respect to the measure $\mu^{u}=\sum_{i=1}^{n}u_{i}^{2}\delta_{\lambda_{i}}$ and that $\pi_{k}^{u}$ is the monic orthogonal polynomial with respect to the same measure.

The Equations (2.1) and (2.6) imply

[TABLE]

We seek to upper bound $\gamma_{k}(u)$ in terms of $u$, so we need to lower bound the quantity

[TABLE]

where $r_{1}(u),\dots,r_{k}(u)$ are the roots of $\pi_{k}^{u}(z)$, i.e. the Ritz values.

Now, if it happens to be the case that the $n$ eigenvalues $\lambda_{i}$ are all clustered very close to the $k$ Ritz values $r_{j}$, then we won’t get a good lower bound. However, if $k\ll n$ and if the $\lambda_{i}$ are reasonably spread out, we expect to get a good lower bound for most $i$. To make this precise, we are led to the notion of equidistribution, which was stated in Section 2.3 and which we restate below.

Definition 4.1**.**

(Restatement of Definition 2.1). Let $\Lambda$ be any finite set of $n$ real numbers. Let $\delta$ and $\omega$ be positive real numbers and let $j$ be a natural number. We say that $\Lambda$ is $(\delta,\omega,j)$-equidistributed if for any finite set $T$ of at most $j$ real numbers,

[TABLE]

We will show in Section 4.1 that a wide range of spectra are equidistributed.

Now we apply the definition. Returning to our effort to upper bound $\gamma_{j}(u)$, we see that if we assume the spectrum of $A$ is $(\delta,\omega,k)$-equidistributed, then

[TABLE]

where $S$ is some subset of $\{1,\dots,n\}$ of size at least $\delta n$. However, for an arbitrary unit vector $u$ and an arbitrary subset $S$, we have no lower bound on the sum $\sum_{i\in S}u_{i}^{2}$—it could even be zero. This leads to our definition of incompressibility in Section 3, which is satisfied by $u$ with high probability.

Indeed, if we assume that the unit vector $u$ is $(\delta,\varepsilon)$-incompressible, then the right hand side expression above is greater than $\varepsilon\omega^{2k}$. Putting together the last few equations, we have $\gamma_{k}(u)\leq(\varepsilon\omega^{2k})^{-1/2}$. We summarize the result in the following proposition.

Proposition 4.2**.**

Suppose the spectrum of $A$ is $(\delta,\omega,k)$-equidistributed and suppose that $u$ is $(\delta,\varepsilon)$-incompressible for some $\delta,\omega,\varepsilon>0$ and $k\in\mathbb{N}$. Then

[TABLE]

4.1 Equidistribution

In this section we establish sufficient conditions for equidistribution that apply to a wide range of spectra. First, we present an immediate generalization of the notion of equidistribution which applies to measures $\mu$ instead of finite sets $\Lambda$. The definitions coincide for finite sets if one identifies $\Lambda$ with the uniform probability distribution on $\Lambda$.

Definition 4.3** (Equidistribution for measures).**

Let $\mu$ be a probability measure on $\mathbb{R}$. Let $\delta,\omega>0$ and $j$ be a natural number. We say that $\mu$ is $(\delta,\omega,j)$-equidistributed if for any finite set $T$ of at most $j$ real numbers,

[TABLE]

If a measure is $(\delta,\omega,j)$-equidistributed for every $j\in\mathbb{N}$, we will just say that it is $(\delta,\omega)$-equidistributed.

For absolutely continuous measures, we have the following general equidistribution result.

Proposition 4.4** (Absolutely continuous measures are equidistributed).**

Let $\nu$ be a compactly supported probability measure on $\mathbb{R}$ with a nontrivial absolutely continuous part. Then there exist constants $\delta,\omega>0$ such that $\nu$ is $(\delta,\omega)$-equidistributed.

Proof.

By the assumption, we may write $\nu$ = $\nu_{1}+\nu_{2}$, where $\nu_{1}$ is absolutely continuous with respect to Lebesgue measure. By cutting off the portion where the density of $\nu_{1}$ is greater than some large $M>0$ and assigning that mass to $\nu_{2}$ instead, we may assume without loss of generality that the density function of $\nu_{1}$ is bounded.

We now utilize a Markov inequality type argument. Let $T$ be any set of $j$ real numbers. Define the logarithmic potential

[TABLE]

Since $\nu_{1}$ has a bounded density function, $\log|x-t|$ is integrable against $\nu_{1}$ for all $t$, so the integral $\int_{-\infty}^{\infty}V_{t}(x)\,d\nu_{1}(x)$ is finite for each $t\in T$. Averaging over all $t\in T$, we find that

[TABLE]

for some constant $a<\infty$. Then

[TABLE]

Relating this back to the definition of equidistribution, we have

[TABLE]

Hence we may take $\delta=\frac{1}{2}\nu_{1}(\mathbb{R})$ and $\omega=e^{-2a}$. ∎

Given our framework, it will be useful to have a statement relating the equidistribution of an absolutely continuous measure to a discretization of that measure. If the two measures are close in Kolmogorov distance, then we can prove such a statement.

Proposition 4.5**.**

Let $\mu$ and $\nu$ be probability measures. If $\mu$ is $(\delta,\omega,j)$-equidistributed for some $\delta,\omega>0$ and $j\in\mathbb{N}$, then $\nu$ is $(\delta-\varepsilon,\omega,j)$-equidistributed, where $\varepsilon=4j\mathrm{Kol}(\mu,\nu)$.

Proof.

Let $T$ be any set of at most $j$ real numbers. Since $p(x)=\prod_{t\in T}|x-t|$ is the absolute value of a polynomial of degree $j$, each of its level sets is a union of at most $2j$ intervals. Hence,

[TABLE]

∎

Thus, to prove equidistribution for an atomic measure, it suffices to prove equidistribution for a nearby absolutely continuous measure.

The above propositions immediately yield a useful corollary for analyzing the Lanczos procedure in the regime of $O(\log n)$ iterations.

Corollary 4.6**.**

Let $\mu$ be a compactly supported probability measure with nontrivial absolutely continuous part. Let $\{\mu_{n}\}$ be a sequence of probability measures such that $\mathrm{Kol}(\mu_{n},\mu)\leq\frac{C}{\log n}$ for some $C>0$. Then for all $n$, for all $j\leq\frac{1}{2C}\log n$ we have that $\mu_{n}$ is $(\delta,\omega,j)$-equidistributed for some $\delta,\omega>0$.

Remark 4.7**.**

If $\mu$ is $(\delta,\omega,j)$-equidistributed and $\nu$ is the pushforward of $\mu$ under the affine map $x\mapsto ax+b$, then $\nu$ is $(\delta,a\omega,j)$-equidistributed.

We now compute the equidistribution for a few example measures, following the proof of Proposition 4.4.

Example 4.8**.**

Let $\mu$ denote the uniform measure on $[0,1]$. Then

[TABLE]

Thus, $\mu$ is $(1/2,4e^{-2})$-equidistributed.

Example 4.9**.**

Let $\nu$ denote the semicircle law $d\nu=\frac{1}{2\pi}\sqrt{(4-x^{2})_{+}}\,dx$. Then

[TABLE]

Thus, $\nu$ is $(1/2,e^{-1})$-equidistributed.

With the above the claims made in the examples of Section 2.3 are now trivial.

Proof of Example 2.2 and Example 2.9.

It is enough to put together Proposition 4.5 and Example 4.8. ∎

Note that for a given set of points that does not resemble a discretization of an absolutely continuous distribution, it will still be likely that the equidistribution parameters are well behaved (relative to their scale) provided that the points are somewhat spread out. On the other hand, if the points are clustered in a few small clusters the analysis becomes trivial.

Observation 4.10**.**

Let $\Lambda$ be a set (or multiset) of $n$ points. Let $a_{1}\leq b_{1}<a_{2}\leq b_{2}<\cdots<a_{m}\leq b_{m}$ be such that $\Lambda\subset\bigcup_{i=1}^{m}[a_{i},b_{i}]$. Define $n_{i}=|\Lambda\cap[a_{i},b_{i}]|$ and let $g$ the minimal gap between clusters, namely, $g=\min_{1\leq i\leq m-1}a_{i+1}-b_{i}$. Then $\Lambda$ is $(\frac{k_{j}}{n},\frac{g}{2},j)$-distributed, where $k_{j}=\min_{S}\sum_{i\in S^{c}}n_{i}$ and $S$ runs over all subsets of $\{1,\dots,m\}$ of size $j$.

Proof.

The proof follows directly from the definition of equidistribution. ∎

Remark 4.11**.**

A particular case of Observation 4.10 is when $n_{i}\geq\lfloor\frac{n}{m}\rfloor$ and $g=a_{i+1}-b_{i}$ for every $i=1,\dots,m$, which yields Example 2.3 above. More generally, if each $n_{i}$ is roughly $n/m$, then $k_{j}$ will be roughly $m-j$, and hence the $\delta$ parameter for the equidistribution of $\Lambda$ will only degrade when $j\approx m$. In other words, Theorem 2.4 is still strong for matrices whose spectrum consists of small clusters if the number of such clusters exceeds the number of iterations of the Lanczos procedure. On the other hand, if the number of iterations exceeds the number of clusters it is not hard to show that the Lanczos procedure will output (with overwhelming probability) at least one Ritz value per cluster.

4.2 Jacobi coefficients

We now have the necessary tools to prove concentration for the entries of the Jacobi matrix.

Proposition 4.12** (Jacobi coefficients are globally Lipschitz).**

Suppose the spectrum of $A$ is $(4\delta,\omega,i)$-equidistributed for some $\delta,\omega>0$ and $i\in\mathbb{N}$. Then for any $0<\varepsilon<\delta$, functions $\alpha_{i}(u)$ and $\beta_{i}(u)$ are globally Lipschitz on $I_{n}(\delta,\varepsilon)$ with constant $L_{i,\varepsilon}\leq\frac{4^{i+2}\lVert A\rVert^{i+1}}{\omega^{i}\varepsilon}$.

Proof.

Proposition 3.7 says that $\alpha_{i}(u)$ and $\beta_{i}(u)$ both have local Lipschitz constant at most $4^{i+1}\|A\|^{i+1}\gamma_{i}(u)$ for all $u\in\mathbb{S}^{n-1}$. Proposition 4.2 says that because the spectrum of $A$ is $(4\delta,\omega,i)$-equidistributed, $\gamma_{i}(u)\leq\frac{1}{\omega^{i}\sqrt{\varepsilon/\sqrt{2}}}$ for all $u\in I_{n}(4\delta,\varepsilon/\sqrt{2})$. Combining these, we have that $\alpha_{i}(u)$ and $\beta_{i}(u)$ are locally Lipschitz with constant

[TABLE]

for all $u\in I_{n}(4\delta,\varepsilon/\sqrt{2})$. Proposition 3.13 says that $I_{n}(\delta,\varepsilon)$ is $\sqrt{2/\varepsilon}$-connected in the larger set $I_{n}(4\delta,\varepsilon/\sqrt{2})$, so Lemma 3.5 implies that $\alpha_{i}(u)$ and $\beta_{i}(u)$ are globally Lipschitz on $I_{n}(\delta,\varepsilon)$ with constant

[TABLE]

∎

We now have the tools to prove our first main theorem, which quantifies the concentration of the Jacobi coefficients around their medians.

Theorem 4.13** (Restatement of Theorem 2.4).**

Suppose the spectrum of $A$ is $(\delta,\omega,i)$-equidistributed for some $\delta,\omega>0$ and $i\in\mathbb{N}$. Let $\tilde{\alpha}_{i}$ and $\tilde{\beta}_{i}$ denote the medians of the Jacobi coefficients $\alpha_{i}(u)$ and $\beta_{i}(u)$, respectively. Then for all $t>0$, the quantities $\mathbf{P}[|\alpha_{i}(u)-\tilde{\alpha_{i}}|>t\|A\|]$ and $\mathbf{P}[|\beta_{i}(u)-\tilde{\beta_{i}}|>t\|A\|]]$ are both bounded above by

[TABLE]

Proof.

The local Lévy lemma (Lemma 3.1) yields that $\mathbf{P}[|\alpha_{i}(u)-\tilde{\alpha_{i}}|>t\|A\|]$ and $\mathbf{P}[|\beta_{i}(u)-\tilde{\beta_{i}}|>t\|A\|]$ are both at most

[TABLE]

where $L_{i,\varepsilon}$ is the global Lipschitz constant on $I_{n}(\delta,\varepsilon)$ obtained in Proposition 4.12. Note that if $\delta>1/50$, then $A$ is still $(1/50,\omega,i)$-equidistributed, so we may set $\varepsilon=\delta/7$ and apply Corollary 3.11 to bound $\mathbf{P}[u\not\in I_{n}(\delta,\varepsilon)]$. We obtain the upper bound

[TABLE]

[TABLE]

as desired. ∎

Combining the previous theorem with Corollary 4.6 we get convergence in probability of the Jacobi matrices in the regime $k=O(\log n)$.

Proposition 4.14**.**

Let the spectra $\mu_{n}$ of $A_{n}$ converge to the spectrum $\mu$ of $A$ in Kolmogorov distance with rate $O(1/\log n)$. Suppose $\mu$ has a nontrivial absolutely continuous part. Then there exists $c_{2}>0$ and a sequence $k_{n}\geq c_{2}\log n$ such that the Jacobi matrices $J_{k_{n}}$ output by the Lanczos algorithm after $k_{n}$ iterations converge to entrywise in probability to deterministic constants.

Proof.

By Corollary 4.6, we have that $\mu_{n}$ is $(\delta,\omega,k)$-equidistributed for all $k\leq c_{1}\log n$. Picking $c_{2}<c_{1}$ and applying Theorem 2.4, for $i\leq c_{2}\log n$ this yields the bound

[TABLE]

so as long as $2c_{2}\log(\omega/4)+1>0$, we have convergence in probability of the Jacobi coefficients as $n\to\infty$. But this is certainly true for small enough $c_{1}$. The $\beta_{i}$ have the same bound as the $\alpha_{i}$, so we are done. ∎

As mentioned in the introduction, convergence for fixed $k$ to the infinite Jacobi matrix $J$ of $\mu$ for deterministic $\mu_{n}$ (with no hypothesis on the rate of convergence of $\mu_{n}$) is proven in [Gau68, Theorem 4]. In Proposition 4.14 we leave it open to prove that the limit is actually $J$ (see Section 6), but if we reduce the number of iterations from $k=O(\log n)$ to $k=O(\sqrt{\log n})$, we can indeed prove that the limit is $J$. This is the content of Theorem 2.12, proven in Section 5.

4.3 Ritz values

Theorem 2.4 yields concentration of the entries of the random matrix $J_{k}(u)$. In general, controlling the entries of a random matrix does not yield control over its random eigenvalues, but, since $J_{k}(u)$ is Hermitian we know that its spectrum is stable with respect to small perturbations of the entries. More precisely, we will use the well known Weyl’s inequality—see [HJ12, Theorem 4.3.1] for a reference.

Lemma 4.15** (Weyl).**

For every matrix $X$, let $\lambda_{1}(X)\geq\dots\geq\lambda_{n}(X)$ denote the eigenvalues of $X$. If $A$ and $B$ are $n\times n$ Hermitian matrices, then for all $1\leq i\leq n$ we have

[TABLE]

Following the notation in Theorem 2.4, let $\tilde{J}_{k}$ be the $k\times k$ Jacobi matrix with entries $\tilde{\alpha}_{i}$ and $\tilde{\beta}_{i}$, and denote the eigenvalues of $\tilde{J}_{k}$ by $\tilde{r}_{1}\geq\dots\geq\tilde{r}_{k}$.

Proposition 4.16** (Concentration of the Ritz values).**

Assume that the spectrum of $A$ is $(\delta,\omega,k)$-equidistributed for some $\delta,\omega>0$ and $k\in\mathbb{N}$. With the notation described above, let $\vec{r}=(\tilde{r}_{1},\dots,\tilde{r}_{k})$ and let $\vec{r}(u)=(r_{1}(u),$ $\dots,r_{k}(u))$ be the vector of Ritz values after $k$ iterations. Then the probability $\mathbf{P}[\lVert\vec{r}(u)-\vec{r}\rVert_{\infty}\geq t\|A\|]$ is bounded above by

[TABLE]

Proof.

Since $\tilde{J}_{k}$ and $J_{k}(u)$ are tridiagonal matrices, we may split $J_{k}-\tilde{J}_{k}$ into the sum of three matrices consisting of the diagonal, the subdiagonal, and the superdiagonal and then use the triangle inequality to obtain

[TABLE]

Hence, we deduce that

[TABLE]

where the first inequality follows from Lemma 4.15 and the second inequality from (4.2). Now observe that the event $\{\max_{0\leq i\leq k-1}\{|\alpha_{i}(u)-\tilde{\alpha}_{i}|\}+2\max_{0\leq i\leq k-2}\{|\beta_{i}(u)-\tilde{\beta}_{i}|\}\geq t\}$ is contained in the event

[TABLE]

which in turn is contained in the event

[TABLE]

Using a union bound and applying Theorem 2.4, we obtain the desired result. ∎

4.4 Ritz Vectors

Here we will use the same notation as in Section 4.3. Let $\tilde{w}_{i}$ be the eigenvector of $\tilde{J}_{k}$ corresponding to $\tilde{r}_{i}$ and let $w_{i}(u)$ be the eigenvector of $J_{k}(u)$ corresponding to $r_{i}(u)$. We will use the fact that $J_{k}(u)$ concentrates around $\tilde{J}_{k}$, together with the Davis-Kahan theorem [DK69] to establish the concentration of the vectors $w_{i}(u)$.

Theorem 4.17** (Davis-Kahan).**

Here we use the notation of Lemma 4.15. Fix $i\in\{1,\dots,n\}$ and assume that $\lambda_{i}(A)$ has multiplicity 1. Define

[TABLE]

and let $\theta\in[0,\pi/2]$ denote the angle between the $i$-th eigenvectors of $A$ and $A+B$. Then

[TABLE]

Under the assumption that $\tilde{r}_{i}(u)$ is not close to the other Ritz values, we get the following result.

Proposition 4.18** (Concentration of the Ritz vectors).**

Assume that the spectrum of $A$ is $(\delta,\omega,k)$-equidistributed for some $\delta,\omega>0$ and $k\in\mathbb{N}$ and fix some $i\in\mathbb{N}$ with $1\leq i\leq k$. With the notation described above, let $\theta\in[0,\pi/2]$ be the angle between $w_{i}(u)$ and $\tilde{w}_{i}$ and let $\varepsilon=\min_{j:j\neq i}|\tilde{r}_{i}-\tilde{r}_{j}|$. Then for any $0\leq c<1/2$, the probability $\mathbf{P}\left[\sin\theta\geq 2\|A\|/\varepsilon n^{c}\right]$ is bounded above by

[TABLE]

Note. The same result holds for the Ritz vectors, since these are obtained by applying an isometry to the $w_{i}(u)$.

Proof.

From Theorem 4.17 we have that

[TABLE]

and hence

[TABLE]

where the latter inequality was established in the proof of Proposition 4.16. Using the bound obtained in the aforementioned proof and substituting $t=\frac{2}{\varepsilon n^{c}}$ we obtain the desired result. ∎

5 Proofs of Proposition 2.10 and Theorem 2.12

5.1 Proof of Proposition 2.10

We now prove our theorem about the Lanczos algorithm missing outliers in the spectrum.

Proof of Proposition 2.10.

By Proposition 4.5, we have that $\mu_{n}$ is $(\delta,\omega,j)$-equidistributed for some $\delta,\omega>0$ and all $j<c\log n$. Suppose $u\in I_{n}(\delta,\varepsilon)$, which happens with overwhelming probability by Proposition 3.10. Then by Proposition 4.2, we have an upper bound on the leading coefficient of the $j$th orthogonal polynomial: $\gamma_{j}(u)\leq\frac{1}{\omega^{j}\sqrt{\varepsilon}}.$ Equivalently, this is a lower bound on the $L^{2}$ norm of the $j$th monic orthogonal polynomial: $\|\pi_{j}^{u}\|_{L^{2}(\mu^{u})}\geq\omega^{j}\sqrt{\varepsilon}.$ As mentioned in the preliminaries in Section 2, it is a classical fact that the monic orthogonal polynomial of any given degree has minimal $L^{2}$ norm over all monic polynomials of that degree. Thus, we in fact have

[TABLE]

for all monic polynomials $q$ of degree $j$, with equality when $q(x)$ is the $k$th orthogonal polynomial $p_{k}^{u}(x)$.

For all unit vectors $u$, let $\rho(u)$ denote the top Ritz value, i.e. the maximum root of $p_{k}^{u}(x)$. We wish to show that $\rho(u)<R+\kappa$ with high probability.

Take $p_{k}^{u}(x)$ and replace its top root by $t$ to form the monic polynomial $P_{t}$. By the first-order condition for the variational characterization of $p_{k}^{u}$ mentioned above, to show $\rho(u)\leq R+\kappa$ it suffices to show that $\|P_{t}\|_{L^{2}(\mu^{u})}$ is strictly increasing in $t$ for $t>R+\kappa$. We have

[TABLE]

where we let $r_{2},\dots,r_{k}$ denote the roots of $p_{k}^{u}(x)$ besides the maximum root $\rho(u)$, and we omit the argument $u$ for brevity. We calculate the derivative

[TABLE]

We wish to show that this quantity is positive whenever $t\geq R+\kappa$. We have assumed that there are only $m$ outliers, so assume $\lambda_{i}\leq R$ for all $i>m$. Then $t-\lambda_{i}\geq\kappa$ for every $m<i\leq n$.

Thus,

[TABLE]

where in the last step we used the inequality (5.1) on the degree $k-1$ polynomial $p_{k}^{u}(x)/(x-\rho(u))$. Simplifying, we have

[TABLE]

By uniform boundedness of the spectra, there exists $M$ large such that $\lambda_{i}-r_{j}\leq M$ for all $1\leq i\leq m$. Let $g$ be the maximum of the outlier gaps $\lambda_{i}-R$ over all $1\leq i\leq m$. Recall that $t\geq R+\kappa$, so $\lambda_{i}+\kappa-t\leq\lambda_{i}-R\leq g$ for all $1\leq i\leq m$. Finally, we have with overwhelming probability $\sum_{i=1}^{m}u_{i}^{2}<n^{-c}$ for any positive $c<1/2$; we will defer the proof to Lemma 5.3 below. Putting this all together, we have

[TABLE]

This quantity is strictly positive when

[TABLE]

Rearranging, we get

[TABLE]

for $n$ large. Note that $\omega<M$, because $\omega$ is a lower bound on geometric means of distances that are all less than $M$. In conclusion, with high probability, $\frac{d}{dt}\|P_{t}\|_{L^{2}(\mu^{u})}^{2}>0$ for all $t>R+\kappa$ when

[TABLE]

For $n$ large, we may absorb the constants $m,g,\kappa,\varepsilon,\omega$ (which do not depend on $n$) into a single constant $c^{\prime}>0$, and we get the desired $k\leq c^{\prime}\log n$.

∎

Remark 5.1**.**

We have focused on the right side of the spectrum for ease of exposition. Similar results hold for outliers on both sides.

Remark 5.2**.**

There are several parameters that can be tuned in the above proof. For example, one could envision a situation in which $\kappa$ converges to zero as $n\to\infty$, at the expense of some other parameter.

Lemma 5.3**.**

Let $0<c<1/2$ and suppose $m\leq n^{\alpha}$, where $\alpha<1-c$. Then $\sum_{i=1}^{m}u_{i}^{2}<n^{-c}$ with overwhelming probability. To be precise,

[TABLE]

Proof.

We proceed just as in the proof of Proposition 3.10. Define $X_{i}$ as in that proof. Then

[TABLE]

Using Lemma 3.12, we solve for the parameter $\sqrt{t}=\frac{-2\sqrt{m}+\sqrt{2}n^{\frac{1}{2}-\frac{c}{2}}}{4}$ (which requires $\alpha<1-c$) and then we get

[TABLE]

which is an overwhelmingly small probability because $\frac{1}{2}-\frac{c}{2}+\frac{\alpha}{2}<1-c$ when $\alpha<1-c$.

Now following the same coupling argument in the proof of Proposition 3.10 and using Lemma 3.12 again, we get

[TABLE]

∎

Proof of Theorem 2.8.

From the proof of Proposition 2.10, setting $\varepsilon=\delta/2$ we have that the Ritz values are contained in the desired interval for

[TABLE]

as long as $k\leq j$, $u\in I_{n}(\delta,\delta/2)$ and $\sum_{i=1}^{m}u_{i}^{2}>n^{-c}$. Applying Corollary 3.11, the probability that $u$ violates either condition is at most

[TABLE]

where in the last step, we apply Lemma 5.3 and note that for $n\geq e^{\frac{1}{1-c-\alpha}}$ we have $4\sqrt{2}n^{\frac{1-c+\alpha}{2}}\leq n^{1-c}$.

∎

5.2 Proof of Theorem 2.12

For $C>0$ let $\mathcal{P}_{C}$ denote the space of Borel probability measures supported on $[-C,C]$. In order to prove Theorem 2.12 we will show that the Jacobi coefficients of a measure are locally Lipschitz quantities on the space $\mathcal{P}_{C}$ equipped with the Kolmogorov metric. Note that in Section 3 similar results were obtained in the case in which the space of measures in consideration is restricted to atomic measures supported on $n$ fixed points, namely, the eigenvalues of $A_{n}$. Since $\mathcal{P}_{C}$ is a much larger and complicated space we are not able to obtain results as strong as in Proposition 3.7. It remains an open question if a better rate can be achieved at this level of generality; see the concluding remarks for some natural directions to pursue.

We will use the following well known result which, for convenience of the reader, we restate as it appears in Lemma 1.1 in [God17].

Lemma 5.4**.**

Let $A$ and $B$ be two $k\times k$ matrices. Then $\det(A+B)$ is equal to the sum of the determinants of the $2^{k}$ matrices obtained by replacing each subset of the columns of $A$ by the corresponding subset of the columns of $B$.

Proof.

The result follows directly from the fact that the determinant is multilinear in the columns of the matrix. ∎

Lemma 5.5**.**

Let $A$ and $B$ be two $k\times k$ matrices. For $1\leq i\leq k$, let $A^{(i)}$ and $B^{(i)}$ be the $i$th columns of $A$ and $B$, respectively. Let $C,\varepsilon>0$ and assume that

[TABLE]

Then

[TABLE]

Proof.

By the assumption in (5.3) we can write $B=A+E$, where $E$ is a matrix with columns of norm less than or equal to $\varepsilon$. Then, using Lemma 5.4, the inequalities in (5.3), and the fact that the determinant of a matrix is bounded by the product of the Euclidean norms of its columns, we obtain

[TABLE]

where the last inequality follows from the mean value theorem. ∎

We now argue that the moments of a measure are Lipschitz quantities in $\mathcal{P}_{C}$, where the constant is exponential in the order of the moment. With this end fix a Borel measure $\mu$ on $\mathbb{R}$ and denote

[TABLE]

A standard application of Fubini’s theorem yields that if $\mu$ is a finite positive Borel measure supported in $[0,\infty)$, then

[TABLE]

This identity is enough to obtain the following bound.

Lemma 5.6**.**

Let $\mu,\nu\in\mathcal{P}_{C}$ and $k>0$, then $|m_{k}(\mu)-m_{k}(\nu)|\leq 2C^{k}\mathrm{Kol}(\mu,\nu)$.

Proof.

Start by decomposing $\mu$ into $\mu_{+}$ and $\mu_{-}$ as follows:

[TABLE]

Hence $\mu(A)=\mu_{+}(A)+\mu_{-}(-A)$. Define $\nu_{+}$ and $\nu_{-}$ analogously. Note that these new measures are supported on $[0,\infty)$.

Observe that $m_{k}(\mu)=m_{k}(\mu_{+})+(-1)^{k}m_{k}(\mu_{-})$ and that the analogous formula holds for $m_{k}(\nu)$. Hence

[TABLE]

Now, for $t\geq 0$ define $F_{\mu_{+}}(t)=\mu_{+}(t,\infty)$ and $F_{\nu_{+}}(t)=\nu_{+}(t,\infty)$. By definition of Kolmogorov distance we have that

[TABLE]

On the other hand, by (5.4) we have that

[TABLE]

In the exact same way we can bound $|m_{k}(\mu_{-})-m_{k}(\nu_{-})|$ to conclude the proof. ∎

Given $\mu\in\mathcal{P}_{C}$ we denote the $(k+1)\times(k+1)$ Hankel matrix of $\mu$ by $M_{k}(\mu)$ and define $D_{k}(\mu)=\det M_{k}(\mu)$. We will denote the Jacobi coefficients of $\mu$ by $\alpha_{i}^{\mu}$ and $\beta_{i}^{\mu}$. For the proof of the following results, many of the facts stated in Section 2.1 will be used.

Proposition 5.7**.**

Let $\mu,\nu\in\mathcal{P}_{C}$ and let $s_{k}>0$ be constants satisfying

[TABLE]

for $j=1,\dots,k$. Then

[TABLE]

for some $g>0$ dependent of $\mu$ and $\nu$ but independent of $k$.

Proof.

To shorten notation let $x_{j}=D_{j}(\mu)$ and $y_{j}=D_{j}(\nu)$. Without loss of generality $C>1$. A direct application of Lemma 5.6 yields a rough bound between the distance in the Euclidean norm of the corresponding columns of the matrices $M_{j}(\mu)$ and $M_{j}(\nu)$. Namely, the columns are at distance less than $\sqrt{j+1}C^{2j-1}\mathrm{Kol}(\mu,\nu)$. The same reasoning yields that the norm of any column in $M_{j}(\mu)$ or $M_{j}(\nu)$ is bounded by $\sqrt{j+1}C^{2j-1}$. Hence, using Lemma 5.5 we get

[TABLE]

for some $g>0$ independent of $k$.

In what follows we will bound two other terms whose logarithm is also $O(k^{2})$. The implied constants depend only on $\mu$ and $\nu$, so we can modify $g$ to be big enough for the following inequalities to hold as well. By the first expression in (2.3) we have that

[TABLE]

To bound the first term on the right-hand side of the above inequality we see that

[TABLE]

which yields

[TABLE]

On the other hand,

[TABLE]

The result then follows from combining the previous inequalities (5.2), (5.6), and (5.7). ∎

Remark 5.8**.**

The constants $s_{k}$ have already been studied with sophisticated techniques for some families of measures; see [Sze77] for an example. However, using results only from Section 4 it will be easy to show that for measures with an absolutely continuous part we have $|\log(s_{k})|=O(k^{2})$, where the implied constant depends only on $\mu$, which is enough for the proof of Theorem 2.12.

In a similar fashion we can show that the coefficients of $p_{k}^{\mu}(x)$ are locally Lipschitz.

Proposition 5.9**.**

Fix a positive integer $k$. Let $\mu,\nu$ and $s_{k}$ be as in Proposition 5.7. Denote the coefficients of $x^{i}$ in $p_{k}^{\mu}(x)$ and $p_{k}^{\nu}(x)$ by $a_{i}^{\mu}$ and $a_{i}^{\nu}$ respectively. Then

[TABLE]

for some $g>0$ dependent on $\mu$ and $\nu$ but independent of $k$.

Proof.

For $1\leq i\leq k$ let $M_{k}^{(i)}(\mu)$ be the matrix obtained by removing the $k$th row and $i$th column of $M_{k}(\mu)$ and let $d_{i}(\mu)=\det(M_{k}^{(i)}(\mu))$. From identity (2.5) we have

[TABLE]

Using the same notation as in the proof of Proposition 5.7 we have that

[TABLE]

As before $\frac{1}{\sqrt{x_{k-1}x_{k}}}\leq\frac{1}{s_{k}}$, while $|d_{i}(\mu)-d_{i}(\nu)|\leq 2\mathrm{Kol}(\mu,\nu)\exp\{gk^{2}\}$ for some $g>0$ dependent on $\mu$ and $\nu$ only. To bound the second term on the right-hand side of the above inequality note that $d_{i}(\nu)\leq\exp\{gk^{2}\}$ and that

[TABLE]

where the last inequality is a consequence of (5.6). The result follows. ∎

Corollary 5.10**.**

Let $\mu,\nu,s_{k}$ be as in Proposition 5.7. Then

[TABLE]

Proof.

Recall that

[TABLE]

As mentioned above, the quantities $a_{i}^{\mu},a_{i}^{\nu}$, and $m_{i}(\mu),n_{i}(\nu)$ are of size $O(\exp\{gk^{2}\})$. Putting this together with Proposition 5.9 and Lemma 5.6 we get that

[TABLE]

By adding over $i,j$ and modifying $g$ the result follows. ∎

In order to prove Theorem 2.12 and Proposition 2.13 we need one final lemma, which states that with overwhelming probability, the random measure $\mu_{n}^{u}$ is close in Kolmogorov distance to $\mu_{n}$.

Lemma 5.11**.**

For $n$ large enough we have that

[TABLE]

Proof.

We must show that

[TABLE]

for all $1\leq k\leq n$ with probability at least $1-\exp\{-n^{1/4}/8\}$.

Fix $1\leq k\leq n$. As in Section 3.3 start by considering $X_{1},\dots,X_{k}$ independent centered Gaussian random variables of variance $\frac{1}{n}$ and let $Z_{k}=\sum_{i=1}^{k}X_{i}^{2}$. Then by Lemma 3.12 we have that

[TABLE]

where $t_{1}$ and $t_{2}$ are the solutions to

[TABLE]

respectively. Since $k\leq n$ it is clear from (5.8) that $\min\{t_{1},t_{2}\}\geq\frac{n^{\frac{1}{4}}}{4}$. This implies that

[TABLE]

Now, letting $k$ run from 1 to $n$, a union bound yields that

[TABLE]

where the last equality holds for $n$ large enough. Now, as in the proof of Proposition (3.10) we can show by a standard coupling argument that if we take $u_{i}=X_{i}/\sqrt{Z_{n}}$, we will have that

[TABLE]

and the result follows. ∎

Proof of Theorem 2.12.

From Lemma 5.11, for $n$ large enough, we have that $\mathrm{Kol}(\mu^{u},\mu_{n})\leq n^{-\frac{1}{4}}$ with overwhelming probability. By the assumption $\mathrm{Kol}(\mu_{n},\mu)=n^{-c}$ we then have that $\mathrm{Kol}(\mu^{u},\mu)\leq n^{-c^{\prime}}$ also with overwhelming probability for $c^{\prime}=\min\{1/4,c\}$. Hence, under the event $\{\mathrm{Kol}(\mu^{u},\mu)\leq n^{-c^{\prime}}\}$ we can apply Proposition 5.7 and Corollary 5.10 and use the fact that the Jacobi matrices are tridiagonal to obtain that

[TABLE]

Since $\mu$ has an absolutely continuous part we know from Proposition 4.4 and Corollary 4.6 that $|\log(\gamma_{k}^{\mu})|=O(k)$. Hence, from (2.4) we get $|\log s_{k}|=O(k^{2})$, which makes it clear that there exists $d>0$ and a sequence $k_{n}\leq d\sqrt{\log n}$ satisfying the theorem statement. ∎

Proof of Proposition 2.13.

As mentioned in Section 2, this proposition is a direct consequence of Theorem 2.12 and Lemma 4.15. ∎

Remark 5.12**.**

Observe that the above proofs repeatedly use the fact that moments are Lipschitz quantities on $\mathcal{P}_{C}$ and that the Jacobi coefficients are an explicit function of the moments. However, going from moments to Jacobi coefficients is an expensive process that we pay for by getting a rate of $O(\sqrt{\log n})$ instead of $\Theta(\log n)$. At first glance, it may seem that the results in Section 3.2 may be used in a similar fashion to obtain a better rate; however, even if we have strong concentration results for the Jacobi coefficients of the random measures $\mu_{n}^{u}$, it is a difficult task to control the location of the medians (or means) of $\alpha_{j}(u)$ and $\beta_{j}(u)$ and hence it is hard to show that these quantities converge at a good enough rate to the Jacobi coefficients of $\mu$.

6 Concluding remarks

Several directions can be pursued to expand the results presented throughout this paper. Currently, we have only analyzed the Lanczos algorithm in its prototypical form, but have not analyzed the more sophisticated variants that are used in practice. Obtaining similar concentration results and negative results for these modifications, and more generally for Krylov subspace methods, would be of great interest.

The Lanczos algorithm is used in practice for non-Hermitian matrices and even nonnormal matrices, despite these cases being far less understood. In this incarnation, the algorithm is referred to as the Arnoldi algorithm. Extending the results of this paper to the Arnoldi algorithm is a natural direction to pursue. As mentioned in Remark 2.7, it is easy to extend Theorem 2.4 to the non-Hermitian setting, but no longer so easy to prove concentration of the Ritz values or to say anything about their location.

The concentration guarantees of the output of the Lanczos algorithm in the present work hold only for up to $C\log(n)$ iterations, where $C$ is a function of the equidistribution parameters of the spectrum of the input matrix. However, we do not know what is the optimal function of $n$ for which a result of this sort holds. Proving a general concentration statement that holds for $\Omega(\log(n))$ iterations would be interesting and require essentially different ideas. A particular case, which might have a simpler solution but is of great interest in applications, is to show concentration of the eigenvectors that correspond to outlying eigenvalues. A possible approach to the latter problem is to take classical guarantees for approximation of outlying eigenvectors, such as [Saa80], and convert them into probabilistic statements like the one given in Theorem 2.4.

In the same direction, another interesting task is to extend Theorem 2.12 to hold in the regime $k=\Omega(\log n)$ instead of the current setting of $O(\sqrt{\log n})$. Some key difficulties are discussed in Remark 5.12.

Another direction is to translate the concentration of the Jacobi coefficients (Theorem 2.4) into some quantitative statement about the quality of the approximate spectral density obtained from the first few Jacobi coefficients. This would require analyzing the conditioning of the Haydock method or Padé approximation mentioned in the introduction. The Jacobi coefficients are also used for estimating matrix functionals via quadrature; see [Gau04] for a comprehensive account. It would be interesting to understand the implications of concentration in this setting as well.

Finally, we pose a question that is of independent mathematical interest, but which would also allow us to remove the square root on the $\log{n}$ in the statement of Theorem 2.12 and Proposition 2.13.

Question 6.1**.**

For $C>0$ let $\mathcal{P}_{C}$ be the set of probability measures with support contained in the interval $[-C,C]$. Is there a natural metric on $\mathcal{P}_{C}$ inducing a topology for which the set of atomic measures is a dense subset of $\mathcal{P}_{C}$, and the Jacobi coefficients

[TABLE]

have local Lipschitz constant at most exponential in $j$?

Acknowledgments

We would like to thank Nikhil Srivastava for many helpful discussions, for introducing us to this topic, and for providing valuable comments on a draft of this paper. We would also like to thank Stanisław Szarek, Lin Lin, and Max Simchowitz for helpful discussions. Finally, we would like to thank the anonymous referees for comments on an earlier version of this paper which helped improve the presentation.