Spectral density estimation for normal matrices
May 29, 2026
Abstract
The spectral density estimation problem asks for an algorithm that, given an \(n\times n\) matrix \(A\), outputs a probability measure that is a good approximation to the uniform distribution on the eigenvalues of \(A\), called the spectral density of \(A\). This paper considers the setting where \(A\) is a large normal matrix that is accessible only through matrix-vector product queries. We provide an algorithm that makes just \(m\) matrix-vector queries to \(A\) and returns, with high probability, a measure within earth mover’s distance \(O(1/m+\log m/{\sqrt n})\) of the true spectral density of \(A\). We provide a complementary lower bound that any algorithm producing an \(\varepsilon\)-approximation to the true spectral density for large matrices must make \(\Omega(1/\varepsilon)\) matrix-vector queries. The lower bound holds even for the more restricted case of real symmetric input matrices. In combination with our upper bound, it shows that spectral density estimation is essentially no harder for complex normal matrices than for real symmetric matrices.
1 Introduction↩︎
Given a matrix \(A \in \mathbb{C}^{n \times n}\), its spectral density \(\mu_A\) is the uniform distribution on its discrete set of \(n\) complex eigenvalues. Spectral density estimation (SDE) is the problem of producing, given access to \(A\), a measure \(\widehat\mu_A\) that approximates \(\mu_A\) in some metric; we will use the earth mover’s distance (also called the Wasserstein distance, Kantorovich-Rubinstein distance, and optimal transport distance).
Spectral density estimation arises as a natural problem in a variety of contexts: in physics and quantum chemistry (where the spectral density is often referred to as the “density of states”) [1], [2], in the study of the Hessians of neural nets during optimization [3], and as a subroutine in large-scale eigenvalue computation [4]. In many of these contexts, the matrix \(A\) is constrained to be real symmetric or Hermitian, and there are several standard references for spectral density estimation in these cases, e.g. [2], [5], [6]. In this paper, we expand the scope of spectral density estimation to normal matrices. Normal matrices which are not Hermitian arise in several settings. The most important class is unitary matrices. Other examples include circulant matrices, normal random matrix models in physics [7], [8], normally regular digraphs [9], and directed Cayley graphs where the set of generators is contained in its own normalizer. Moreover, several matrix operations preserve normality: applying matrix functions, taking Kronecker products, adding commuting normal matrices, and taking the compound matrix (matrix of minors), to name a few.
The spectral density can in principle be computed to a high degree of precision by diagonalizing \(A\), at cost \(O(n^\omega)\) (or \(O(n^3)\) in practice) [10]. For modern problems where \(n\) is enormous or \(A\) is given implicitly (e.g., as a black box), this is prohibitive. The earth mover’s relaxation—only requiring that \(\mu_A\) and \(\widehat\mu_A\) are close in earth mover’s distance—opens the door to dramatically faster algorithms.
In particular, our algorithm is a matrix-vector query (matvec) algorithm. These algorithms interact with \(A\) only through queries of the form \(v \mapsto Av\) (and possibly \(A^*v\)). The matvec complexity of a problem is the number of such queries needed to solve it; in the worst-case, each query may cost \(O(n^2)\), and so we obtain a substantially faster algorithm than full diagonalization when the query complexity is \(o(n^{\omega-2})\) (or just \(o(n)\) in practice). Important classes of matrix-vector query algorithms are Krylov subspace methods and randomized linear sketching methods, both of which are frequently used to approximate outlying eigenvalues of matrices and solve related spectral problems [11]–[16].
For Hermitian matrices, the well-known stochastic Lanczos quadrature (SLQ) and kernel polynomial method (KPM) both show that \(m\) matrix-vector queries suffice to achieve an earth mover’s approximation error of \(O(1/m)\) to the spectral density of any Hermitian matrix when the matrix dimension \(n\) is sufficiently large [17], [18]. A restricted lower bound has also been shown: if the very small approximation error of \(O(1/n)\) is desired, then \(\Omega(n)\) queries are required [19], [20].
Both SLQ and KPM use the fact that the spectrum of \(A\) lies on the real line, and SLQ additionally uses the property that compressions of Hermitian matrices are themselves Hermitian. Neither of those hold in the general normal setting [21]. Because the spectrum of \(A\) can be distributed throughout a 2D set in \(\mathbb{C}\), one might expect more matrix-vector queries are needed to resolve it. A heuristic based on packing suggests the problem should be quadratically harder: at error \(\varepsilon\), one can fit \(\Theta(1/\varepsilon^2)\) well-separated eigenvalue clusters in a bounded region of the plane, versus \(\Theta(1/\varepsilon)\) on the line. Indeed, such a loss is inherent for related problems like numerical range estimation [22].
Contrary to this intuition, we show that \(m\) matvecs suffice to achieve approximation error \(O(1/m)\) for large normal matrices \(A\). In particular, normal SDE is no harder than Hermitian SDE in the matvec model. We additionally show a substantially broader lower bound as compared to [19], [20]: for any \(\varepsilon>0\), if an algorithm produces an \(\varepsilon\) approximation for large matrices, then we show \(\Omega(1/\varepsilon)\) queries are required.
Before stating our results more formally, we review some standard notation. We use \(\mleft\|\cdot\mright\|\) to denote the \(\ell_2\) (Euclidean) norm of vectors, and the spectral norm of a matrices. \(\mleft\|\cdot\mright\|_p\) denotes the \(\ell_p\) norm of vectors. For a set \(S\), \(\mleft\|f\mright\|_S=\sup_{x\in S}\mleft|f(s)\mright|\) is the sup-norm of a function \(f\) on \(S\). The \(\ell_p\)-Lipschitz constant of \(f:\mathbb{C}^n\to\mathbb{C}\) is \(\sup_{x\neq y}\frac{\mleft|f(x)-f(y)\mright|}{\mleft\|x-y\mright\|_p}\). The Dirac delta distribution \(\delta(z)\) is defined by \(\int f(z)\delta(z-c)\,\textrm dz=f(c)\) for continuous functions \(f\). \(\mathcal{W}_1\) is the earth mover’s distance defined for compactly supported probability measures \(\mu\), \(\nu\) as \(\mathcal{W}_1(\mu,\nu)=\sup_f\mleft(\int f(z)(\text{d}\mu(z)-\,\textrm d\nu(z))\mright)\), where the supremum is over functions \(f\) with \(\ell_2\)-Lipschitz constant bounded by 1. We note that there is an equivalent definition of \(\mathcal{W}_1(\cdot,\cdot)\) in terms of the optimal transport distance between \(\mu\) and \(\nu\). The above characterization, which will be easier to work with, is often known as Kantorovich-Rubinstein duality. We denote the spectral density of a matrix \(A \in \mathbb{C}^{n \times n}\) as \[\begin{align} \label{a0} \mu_A(z) \equiv \dfrac{1}{n} \sum_{j=1}^n \delta(z - \lambda_j), \end{align}\tag{1}\] where \(\lambda_1, \ldots, \lambda_n\) are the eigenvalues of \(A\). Our main algorithmic result is the following.
Theorem 1. Let \(A\in\mathbb{C}^{n\times n}\) be a normal matrix with spectral density \(\mu_A\). Fix \(0 < \delta< 1\). 4 computes an approximating measure \(\widehat \mu_A\) satisfying, with probability \(1-\delta\), \[{\label{a3}\dfrac{1}{{\left\lVert A \right\rVert}}\mathcal{W}_1(\mu_A, \widehat \mu_A) \le O\mleft(\frac{1}{m}+\sqrt{\dfrac{\log(m/\delta) \log(m)}{n}}\mright).}\qquad{(1)}\] The algorithm uses \(O(m)\) matrix-vector products with \(A\) and \(A^*\), \(O(mn)\) storage, and \(O(m^2n)\) additional running time. Moreover, the measure \(\widehat \mu_A\) may be chosen to be supported on just \(O(m^2)\) atoms.
Remark 1. Several variants of 4 are described in 2.5 for various purposes and applications, all still achieving ?? or the appropriate analog.
Remark 2. In addition to applying to a broader class of matrices, 4 improves upon the run-times of the algorithms of [18], [20] for Hermitian SDE. Namely, the \(m\) dependence is improved by avoiding an optimization problem used by those algorithms to enforce strict moment-matching not required in our analysis. Additionally, a variant of 4 described in 2.5 allows us to output \(\widehat\mu_A\) supported on just \(O(m)\) atoms in the Hermitian case, compared to \(O(m^{1.5})\) atoms in the algorithm of [20]. We also note that the arguments of [20] for producing atomic measures do extend to the non-Hermitian / 2D setting, but this only gives a construction for \(O(m^3)\) atoms compared to \(O(m^2)\) of 1.
![Figure 1: Spectral density estimation for two example matrices. The leftmost column is a scatter plot of the eigenvalues of the matrices. The next columns left-to-right are the result of running 4 for m=16, 32, 64, and 128 respectively. The first row is with the 150,000\times150,000 matrix where 100,000 eigenvalues are uniformly distributed in the disk, and 50,000 additional eigenvalues are uniformly distributed in the intersection of the disk with the second and fourth quadrants. The second row is with the 27,078\times27,078 matrix which has has its eigenvalues trace out a version of the “Anscombosaurus” [23], [24]. All plots depict the unit square [-1,1]+i[-1,1]\subset\mathbb{C}. The total time for generating the entire figure (including plotting) on a standard laptop with an i7 core took less than 2 seconds.](/article/2605.31430/pyplot3.png "fig:")
An alternative way to state 1 is that to achieve a target error of \(\varepsilon\mleft\|A\mright\|\) in earth mover’s distance for an \(n \times n\) matrix with \(n\) sufficiently large, 4 need only make \(O(1/\varepsilon)\) queries. We complement this result with a matching lower bound.
Theorem 2. Consider any triplet of parameters \(\varepsilon,\delta,n\) satisfying \(\delta\ge\varepsilon\) and \(n\gg\max\mleft(\frac{1}{\varepsilon^4\delta^2},\mleft(\frac{2}{\varepsilon}\mright)^{1/\delta}\mright)\). If a randomized, adaptive, matrix-vector query algorithm, on input any \(n\times n\) real symmetric matrix \(A\), outputs a distribution that is within \(\varepsilon\mleft\|A\mright\|\) in earth mover’s distance of the spectral density \(A\) with probability more than \(2\delta\), then it must make \(\Omega(1/\varepsilon)\) matrix-vector queries.
Remark 3. The lower bound of 2 is against algorithms which need only succeed for real symmetric inputs. Thinking of \(\varepsilon\) and \(\delta\) as fixed constants, this matches our upper bound in the general normal case, establishing that there is no gap between the complexities of the normal and real symmetric spectral density estimation problems.
Remark 4. A weaker lower bound changes the order of the quantifiers on the parameters: suppose one wants an algorithm which for fixed \(n\) can produce an \(\varepsilon\) approximation for all \(\varepsilon>0\), including \(\varepsilon\) going to 0. In this case, there are several existing lower bounds, e.g. [20], [25], [26] all prove the following: for each \(n\), there exists an \(\varepsilon\) such that approximating the spectral density to \(\varepsilon\) earth mover’s error requires \(O(1/\varepsilon)\) queries. In particular, their constructions are for \(\varepsilon=\Theta(1/n)\). Stated differently, they show that to achieve \(O(1/n)\) error, one must learn a constant fraction of the matrix. However, this is not typically the setting in which SDE is interesting. The goal is to obtain fast algorithms which require very few queries relative to the matrix size and output a relatively coarse approximation.
In the practical regime where \(\varepsilon\) is a small, fixed constant (or tends to 0 slowly), and \(n\) is large, the previous lower bounds do not apply. 2 addresses this gap: it applies whenever \(\varepsilon=\Omega(1/n^{1/4})\). This resolves the complexity of spectral density estimation outside the window \(1/n\ll \varepsilon\ll1/n^{1/4}\).
2 follows from a more general statement, 6, which isolates the information-theoretic bottleneck at the core of spectral density estimation using the ideas of [27]. We concretely find that any algorithm making only \(m\) matrix vector queries to real symmetric \(A\) does not learn much unitarily-invariant information about \(A\) beyond the values of \(\mathop{\mathrm{tr}}(A^j)\) for \(j=1,\ldots,4m\), which are just the first \(4m\) moments of the spectral density of \(A\). These moments are both necessary and sufficient to approximate the spectral density of \(A\).
2 Algorithm and its analysis↩︎
In this section, we derive 4 for approximating the spectral density of a normal matrix and discuss its theoretical properties, most notably by proving 1. We begin by introducing the algorithm itself: given matrix-vector access to \(A\) and \(A^*\), we explain how to compute mixed Chebyshev moments of a random surrogate spectral measure and then use these moments to form an approximation to \(\mu_A\). After presenting the algorithm, we turn to its analysis via a multivariate version of the kernel polynomial method (KPM). Although the KPM is well known in the Hermitian setting [5], rigorous error bounds for approximating probability distributions are more recent [18], and the multivariate version we require here has received comparatively little attention, excepting the appendix of [20]. For this reason, we include a discussion of Jackson’s theorems, their multivariate extension, and how they are deployed our setting. We conclude the section with some algorithmic variants and numerical examples, the results of which are depicted in 1.
Weighted spectral density.
A key object appearing in our analysis is the weighted spectral density of a matrix \(A\). If we have the eigendecomposition \(A=\sum_{j=1}^n\lambda_jv_jv_j^*\), then the weighted spectral density with respect to a unit vector \(b\in\mathbb{C}^n\) is the measure \[\label{a9} \mu_{A,b}(z)\equiv\sum_{j=1}^n\delta(z-\lambda_j)\mleft|b^*{v_j}\mright|^2.\tag{2}\] When \(b=\frac{1}{\sqrt n}\mleft(v_1+\cdots+v_n\mright)\), this coincides with the unweighted spectral measure \(\mu_{A,b}=\mu_A\) defined in 1 . The analysis of our algorithm consists of two steps: (1) show that we can compute an approximation to \(\mu_{A,b}\) for any \(b\) of our choice (without ever explicitly forming \(\mu_{A,b}\)), and (2) show that \(\mu_{A,b}\) is an approximation to \(\mu_A\) with high probability over a Haar random selection of \(b\).
Chebyshev polynomials and moments.
Playing a role both steps mentioned above are the Chebyshev polynomials and Chebyshev moments. We let \(T_j(x)\) denote the \(j^{\textrm{th}}\) Chebyshev polynomial of the first kind, defined by the three-term recurrence, \[\label{a10}\begin{align} T_0(x) &= 1, \quad T_1(x) = x, \\ T_{n+1}(x) &= 2x\,T_n(x) - T_{n-1}(x), \quad n \ge 1.\end{align}\tag{3}\] The Chebyshev polynomials are orthogonal on the interval \([-1, 1]\) with respect to the weight function \(w(x) = (1-x^2)^{-1/2}\): \[\label{a11} \int_{-1}^1 \dfrac{T_j(x)T_k(x)}{\sqrt{1-x^2}}\,\textrm dx = \begin{cases} 0 & j \neq k \\ \pi & j = k = 0 \\ \pi/2 & j = k \geq 1 \end{cases}.\tag{4}\] We write \(\widetilde{T}_0(x) = 1/\sqrt{\pi}\) and \(\widetilde{T}_j(x) = \sqrt{2/\pi} \cdot T_j(x)\) for the normalized Chebyshev polynomials, which are an orthonormal family. For a measure \(\mu\) supported on the square \([-1,1]^2\), define the \((j,k)\) mixed Chebyshev moment as \[\Gamma_{jk} = \int_{-1}^1 \int_{-1}^1 \widetilde{T}_j(x)\widetilde{T}_k(y) \,\textrm d\mu(x, y).\] That is, \(\Gamma_{jk}\) is obtained by integrating the bivariate polynomial \(\widetilde{T}_j(x)\widetilde{T}_k(y)\) over the unit square with respect to \(\mu\). By identifying \(z=x+iy\in\mathbb{C}\), this definition applies to measures \(\mu\) on the unit square \([-1,1]+i[-1,1]\) in \(\mathbb{C}\) as well.
In our application to spectral density estimation, we assume \({\left\lVert A \right\rVert}_2 \leq 1\) so that the real and imaginary parts of the eigenvalues of \(A\) each lie in the interval \([-1,1]\). This is without loss of generality, as one may apply the algorithm to a shifted and rescaled matrix \(\widetilde{A} = c_1(A-c_2)\) and then shift and rescale the output \(\widehat\mu\) to \(\widehat \mu(z/c_1+c_2)\).1
2.1 Chebyshev moments of weighted spectral measure↩︎
Our first goal is to compute the mixed Chebyshev moments \(\Gamma_{jk}\) for the weighted spectral measure \(\mu_{A,b}\). The key observation is that that the Toeplitz decomposition \(A = H + iK\), i.e. \(H = (A+A^*)/2\) and \(K = (A-A^*)/2i\) of a normal matrix allows us to compute these moments without actually forming \(\mu_{A,b}\).
Lemma 1. Let \(A = H+iK\) be the Toeplitz decomposition of the normal matrix \(A\) so that \(H = (A+A^*)/2\) and \(K = (A-A^*)/2i\). Then for any fixed polynomials \(p, q\) it holds that \[\int_{-1}^1\int_{-1}^1\overline{p}(x)q(y)\,\textrm d\mu_{A,b}(x,y) = b^* p(H)^*q(K)b\] where \(\bar p\) denotes the polynomial with coefficients that are the conjugates of the coefficients of \(p\).
Proof. Since \(A\) is normal, by the spectral theorem we may write its eigendecomposition as \(A=U\Lambda U^*\) where \(U\) is unitary. From the definitions \(H=(A+A^*)/2\) and \(K=(A-A^*)/2i\), we have that \(H\) and \(K\) have the corresponding eigendecompositions \[H = U \mathrm{Re}(\Lambda)U^*, \quad K = U \mathrm{Im}(\Lambda)U^*.\] Then we observe that the quadratic form simplifies considerably: \[\begin{align} b^* \overline{p}(H)q(K)b & = b^*U\overline{p}(\operatorname{Re}(\Lambda))U^*Uq(\operatorname{Im}(\Lambda))U^*b \\ & = (p(\operatorname{Re}(\Lambda))U^*b)^*(q(\operatorname{Im}(\Lambda))U^*b) \\ & = \sum_{j=1}^n\overline{p}(\operatorname{Re}(\lambda_j))q(\operatorname{Im}(\lambda_j)) |(U^*b)_j|^2 \\ & = \int_{-1}^1 \int_{-1}^1 \overline{p}(x)q(y) \,\textrm d\mu_{A,b}(x, y). \end{align}\] The first equality follows from the definition of matrix functions that if \(B = VDV^{-1}\) then \(f(B) = Vf(D)V^{-1}\), the second and third equalities rewrite the result as an inner product, and the final equality rewrites the discrete sum as an integral over the weighted spectral measure \(\mu_{A,b}\) as defined in 2 . ◻
If we collect the mixed Chebyshev moments of \(\mu_{A,b}\) up to degree \(m\) into an \((m+1)\)-by-\((m+1)\) matrix \(\Gamma\), then 1 provides the factorization \[\label{a13} \Gamma = X^*Y\tag{5}\] where \[X = \begin{bmatrix}\widetilde{T}_0(H)b& \widetilde{T}_1(H)b&\cdots& \widetilde{T}_m(H)b\end{bmatrix}\] and \[Y= \begin{bmatrix}\widetilde{T}_0(K)b&\widetilde{T}_1(K)b&\cdots&\widetilde{T}_m(K)b\end{bmatrix}.\] The matrices \(X, Y\) can be computed using the three-term recurrence of the Chebyshev polynomials, costing only \(m\) matrix-vector products with each of \(H\) and \(K\). If \(A\) and \(A^*\) are only available implicitly, then one can simulate the application of \(H\) via \(Hx = (Ax+A^*x)/2\), and similarly for \(K\). The resulting procedure is summarized in 2, which as a result requires \(O(m)\) matrix-vector products with \(H\) and \(K\) or \(O(m)\) matrix-vector products with \(A\) and \(A^*\). It also requires \(O(mn)\) storage for \(X\) and \(Y\) and \(O(m^2n)\) additional operations for the formation of \(\Gamma = X^*Y\).
2.2 Approximating the weighted spectral measure.↩︎
The previous subsection gave an algorithm for computing the mixed Chebyshev moments \(\Gamma_{jk}\) of \(\mu_{A,b}\). We now show how to use a multivariate version of the kernel polynomial method (KPM) [5] to convert these moments to an approximation to \(\mu_{A,b}\).2 Our approximation will be \[\label{a15} \,\textrm d\mu_{A,b}(x,y) \approx \,\textrm d\mu_{\textrm{KPM}}(x,y) \equiv \dfrac{p_m(x,y)}{\sqrt{1-x^2}\sqrt{1-y^2}}\, \,\textrm dx\, \,\textrm dy,\tag{6}\] where \[p_m(x,y) = \sum_{j=0}^m\sum_{k=0}^m \Gamma_{jk} \rho_j\rho_k \widetilde{T}_j(x)\widetilde{T}_k(y).\] The \(\widetilde{T}_j\) are the normalized Chebyshev polynomials and the \(\rho_j\) are fixed constants, sometimes called damping coefficients, which we specify momentarily. We summarize the resulting procedure in 3, which uses 2 as a subroutine for computing the Chebyshev moments. The central theoretical result of this subsection is the following theorem, whose proof is found at the end of this subsection:
Theorem 3. Let \(\mu\) be a probability distribution on \([-1, 1]^2\) and let \(\mu_{\textrm{KPM}}\) be defined as in 6 . Then \[\mathcal{W}_1(\mu, \mu_{\textrm{KPM}}) \le \frac{12}{m}.\]
To understand \(\mu_{\textrm{KPM}}\), we first present the usual single-variable KPM formulation and discuss its generalization to multiple variables. The central object is the polynomial kernel \(K_m\), \[\label{a18} K_m(x,r) \equiv \sum_{j=0}^m\rho_j \widetilde{T}_j(x) \widetilde{T}_j(r).\tag{7}\] One associates with \(K_m\) a linear operator (which we will also denote \(K_m\)) mapping bounded functions \(f\) supported on \([-1, 1]\) to polynomials of degree at most \(m\). Application of this operator is performed via the formula \[\label{a19} (K_mf)(x) = \int_{-1}^1 \dfrac{K_m(x,r)f(r)}{\sqrt{1-r^2}}\,\textrm dr.\tag{8}\] If \(g(x)\) is a probability density function on the interval \([-1, 1]\) then the KPM approximation to \(g(x)\) is the measure with density \[{\label{a20}\widetilde{g}(x)= \frac{1}{\sqrt{1-x^2}}\int_{-1}^1K_m(x,r)g(r)\,\textrm dr}\tag{9}\] which can be interpreted as the composition of three operations on \(g\): pointwise multiplication by \(\sqrt{1-x^2}\), application of the kernel \(K_m\), and division by \(\sqrt{1-x^2}\). Importantly, this multiplication by the weight function \(\sqrt{1-x^2}\) exactly cancels the weight function from the definition of the kernel application 8 . By expanding the kernel definition in the resulting approximation, one sees that \(\widetilde{g}\) can be computed from just the knowledge of the Chebyshev moments of \(g\), \(\int_{-1}^1T_j(r)g(r)\,\textrm dr.\)
In order to ensure that \(\widetilde{g}\) approximates \(g\) and that \(\widetilde{g}\) remains a probability density, the damping coefficients \(\rho_j\) must be selected carefully. One choice is to pick \(\rho_j\) to be the optimum of the following linear program: minimize a quantity called the squared kernel resolution \[Q(K_m) \equiv \int_{-1}^1 \int_{-1}^1 \dfrac{K_m(x,r)(x-r)^2}{\sqrt{1-x^2}\sqrt{1-r^2}}\,\textrm dx \,\textrm dr,\] which rewards \(K_m\) for being close to the identity operation, subject to \[\label{a21} \int_{-1}^1\dfrac{K_m(x,r)}{\sqrt{1-r^2}} \,\textrm dr = 1, \quad K_m(x,r) \geq 0 \textrm{ for } x,r \in [-1, 1],\tag{10}\] which ensure \(\widetilde{g}\) is nonnegative and integrates to 1, i.e. is a probability density function. Those optimal values of \(\rho_j\) are obtained in, e.g., [5], [28] and are given by \[\label{a22} \rho_j = \frac{(m+2-j) \cos\mleft(\dfrac{j\pi}{m+2}\mright) + \sin\mleft(\dfrac{j \pi}{m+2}\mright)\cot \mleft(\dfrac{\pi}{m+2}\mright)}{m+2}, \quad j = 0, 1, \cdots, m.\tag{11}\] A different derivation, which emphasizes trigonometric polynomials rather than algebraic polynomials, for these same coefficients was also obtained by Rivlin [29].
KPM posits that \(\widetilde{g}\) is close to \(g\) in the earth mover’s distance. To understand why, we summarize several results from [29], which control how \(K_m\) acts on Lipschitz test functions. When reading that reference, one should note that the coefficients in [29] are not presented in the form of 11 , but they are equivalent.3
Lemma 2. Suppose \(f\) is an \(\ell\)-Lipschitz function on \([-1, 1]\), and let the damping coefficients \(\rho_0\), \(\rho_1\), \(\ldots\), \(\rho_m\) be selected such that \(K_m\) is positive and normalized. Then \[{\left\lVert f -K_mf \right\rVert}_{[-1, 1]} \leq \dfrac{\ell}{m}\left( 1+\dfrac{m \pi}{\sqrt{2}}\sqrt{1-\rho_1}\right).\] In particular, with the choice in 11 , one obtains \[{\label{a24} {\left\lVert f - K_m f \right\rVert}_{[-1, 1]} \leq \dfrac{6\ell}{m}.}\qquad{(2)}\]
Remark 5. The constant \(6\) in ?? may be improved, especially for large \(m\); however, a lower bound on the constant of \(\pi/2\) is known [30].
Remark 6. Results similar to 2 of this kind are often called “Jackson-type theorems,” after Dunham Jackson, who showed that algebraic or trigonometric polynomials can be used to approximate Lipschitz functions with a sup-norm error of \(O(1/m)\). Although the kernel \(K_m\) that we use here arising from the choice 11 is sometimes called “Jackson’s kernel” in the literature, we note that Jackson’s original construction as given in [30] is not exactly equivalent. Still, Jackson proves that the kernel he introduced in [30] obtains an analogous result. The kernel of [30] is also analyzed in [18]. For practical purposes, differences between the kernels seem to be negligible.
2 directly allows one to control the earth mover’s distance between the measures with densities \(g\) and \(\widetilde{g}\): fix a 1-Lipschitz test function \(f\) and observe \[{\label{a25}\int_{-1}^1f(x)\widetilde{g}(x)\,\textrm dx = \int_{-1}^1\int_{-1}^1\frac{f(x)K_m(x,r)g(r)}{\sqrt{1-x^2}}\,\textrm dr\,\textrm dx =\int_{-1}^1 (K_mf)(r)g(r)\,\textrm dr.}\tag{12}\] Since 2 implies \(f\) and \(K_mf\) are close, this means that the integrations of \(f\) against \(g\) and \(\widetilde{g}\) are close, which is precisely what is required to show that \(g\) and \(\widetilde{g}\) are close in the earth mover’s distance. From 12 we can directly bound: \[\begin{align} \mathcal{W}_1(g,\widetilde{g}) = \sup_{{\left\lVert f \right\rVert}_{\text{Lip}} \le 1} \int_{-1}^1 f(z)(g(z) - \widetilde{g}(z)) d x = \sup_{{\left\lVert f \right\rVert}_{\text{Lip}} \le 1} \int_{-1}^1 (f(z) - K_m f(z))g(z) d x \le \frac{6}{m}, \end{align}\] where the last inequality follows from ?? .
One final observation of single-variable KPM before generalizing to multiple variables is that \(K_m\) is a contraction in the following sense.
Lemma 3. Let \(K_m\) be the polynomial kernel defined by 7 , 8 , and 11 . \(K_m\) has the following contraction property: \[{\left\lVert K_m f \right\rVert}_{[-1, 1]} \leq {\left\lVert f \right\rVert}_{[-1,1]}.\]
Proof. Simply expand the definition to see that \[\begin{align} {\left\lVert K_m f \right\rVert}_{[-1,1]} & =\max_{x \in [-1, 1]}\int \dfrac{|K_m(x,r)f(r)|}{\sqrt{1-r^2}} \,\textrm dr \\ & \leq \max_{r \in[-1,1]} |f(y)| \max_{x \in[-1, 1]}\int_{-1}^1 \dfrac{K_m(x,r)}{\sqrt{1-r^2}} \,\textrm dr \\ & = {\left\lVert f \right\rVert}_{[-1,1]}. \end{align}\] ◻
The preceding discussion for the univariate KPM generalizes in a straightforward way to multivariate problems in \(d \geq 1\) variables on the hypercube \([-1,1]^d\). In the following we state multivariate analogues of 2 and 12 , and we consequently bound the KPM approximation error in earth mover’s distance in 3.4 Define for a multivariate function \(f(x_1, x_2, \ldots, x_d)\) on \([-1,1]^d\) the component Lipchitz constants \[\mathrm{Lip}_j(f) \equiv \sup_{\substack{x_k\in[-1,1]\\ k\neq j}} \;\sup_{\substack{u,v\in[-1,1]\\ u\neq v}} \dfrac{\bigl|f(x_1,\dots,x_{j-1},u,x_{j+1}\dots,x_d)-f(x_1,\dots,x_{j-1},v,x_{j+1},\dots,x_d)\bigr|}{|u-v|}.\] \(\mathrm{Lip}_j(f)\) is the maximal Lipschitz constant of \(f\) when viewed as a function of \(x_j\) with all other variables fixed.
The one-dimensional kernel \(K_m\) defined by 8 extends to a \(d\) dimensional analog via the tensor product, i.e., let \(\otimes^d K_m = K_m \otimes K_m \otimes \cdots \otimes K_m\) be the integral operator defined by applying the one-dimensional Jackson operator separately in each coordinate: \[{\label{a28}(\otimes^d K_m f)(x_1, \cdots, x_d) = \int_{[-1,1]^d} f(y_1,\cdots, y_d)\prod_{j=1}^d\left( \dfrac{K_m(x_j, y_j)}{\sqrt{1-y_j^2}}\right) \,\textrm dy_1 \cdots \,\textrm dy_d.}\tag{13}\] Then we have the following multivariate extension of 2.
Theorem 4 (Multivariate Jackson). Let \(f(x_1, x_2, \cdots, x_d)\) be a continuous function on \([-1, 1]^d\) and suppose there exist positive constants \(\ell_1, \cdots, \ell_d\) such that \(\mathrm{Lip}_{j}(f) \leq \ell_j\) for \(j = 1, 2, \cdots, d\). Then \[{\left\lVert f - (\otimes^d K_m)f \right\rVert}_{[-1,1]^d} \leq \frac{6}{m}\sum_{j=1}^d \ell_j,\]
Proof. We prove first \(d = 2\) for simplicity, and then describe the general proof afterward. We use \((I \otimes K_m)\) to indicate the linear integral operator that applies Jackson’s one-dimensional kernel on the second variable: \[((I \otimes K_m)f)(x,y) = \int_{-1}^1f(x, s) \dfrac{K_m(y,s)}{\sqrt{1-s^2}}\,\textrm ds .\] The operator \((K_m \otimes I)\) is defined similarly as operating on the first variable. Under this notation, we may factor the iterated integrals via \((K_m \otimes K_m) = (I \otimes K_m)(K_m \otimes I) = (K_m \otimes I)(I \otimes K_m)\). Then observe \[\begin{align} {\left\lVert f - (K_m \otimes K_m)f \right\rVert}_{[-1,1]^2} & = {\left\lVert f - (I \otimes K_m)f + (I \otimes K_m)(f - (K_m \otimes I)f) \right\rVert}_{[-1,1]^2} \\ & \leq {\left\lVert f - (I \otimes K_m)f \right\rVert}_{[-1,1]^2} + {\left\lVert (I \otimes K_m)(f- (K_m \otimes I)f) \right\rVert}_{[-1,1]^2} \\ & \leq {\left\lVert f - (I \otimes K_m)f \right\rVert}_{[-1,1]^2} + {\left\lVert f- (K_m \otimes I)f \right\rVert}_{[-1,1]^2} \\ & \leq \dfrac{6 \ell_2}{m} + \dfrac{6 \ell_1}{m}. \end{align}\] The first equality comes from adding and subtracting \((I \otimes K_m)f\). The inequality on the second line is an application of the triangle inequality. For the third line, let \(x\) be fixed and define the one-dimensional slice \(y \mapsto h_x(y) = (f - (K_m \otimes I)f)(x, y)\), so that we can apply 3 by observing \[{\left\lVert (I \otimes K_m)(f- (K_m \otimes I)f) \right\rVert}_{[-1,1]^2} = \sup_{x \in [-1, 1]} {\left\lVert K_mh_x \right\rVert} \leq \sup_{x \in [-1,1]} {\left\lVert h_x \right\rVert}_{[-1,1]} = {\left\lVert f - (K_m \otimes I)f \right\rVert}_{[-1,1]^{2}}.\] Finally, the fourth line is a consequence of applying the one-dimensional 2 to the one-dimensional slices analogous to the argument applied to \(h_x\) above. More precisely, for fixed \(x\), observe that \(y \mapsto f_x(y) \equiv f(x,y)\) is an \(\ell_2\)-Lipschitz function of \(y\), so that \[{\left\lVert f - (I \otimes K_m)f \right\rVert}_{[-1,1]^2} = \sup_{x \in [-1,1]} {\left\lVert f_x-K_mf_x \right\rVert}_{[-1,1]} \leq \dfrac{6\ell_2}{m}.\] The same argument may be applied to the analogously defined \(f_y\) to control the second term in the third line.
The \(d > 2\) case follows from a similar argument. If we let \(P_j = (\otimes^{j}K_m)\otimes(\otimes^{d-j}I)\) be the operator which applies \(K_m\) in the first \(j\) coordinates, we apply the triangle inequality to \[f-(\otimes^{d}K_m)f = \sum_{j=1}^{d} P_{j-1}f - P_{j}f,\] whose \(d\) terms lead to the factor of \(\sum_{j=1}^d \ell_j\) in the result. ◻
Remark 7. If \(f\) is \(\ell\)-Lipschitz with respect to any \(p\)-norm then \(\sum_j \mathrm{Lip}_j(f) \leq d \ell\).
Remark 8. Analogous results can be obtained by using a different kernel \(K_m\) corresponding to a different choice of the \(\rho_j\), such as Jackson’s original construction [30].
The calculation in 12 immediately generalizes to the multivariate setting. For simplicity, we only state it for \(d=2\) variables, which is all we need for complex measures.
Lemma 4. Let \(\mu\) be a probability distribution on \([-1,1]^2\). Define further the approximation \[\,\textrm d\mu_\textrm{KPM}(x, y) \equiv \dfrac{p_m(x,y)}{\sqrt{1-x^2} \sqrt{1-y^2}} \,\textrm dx \,\textrm dy\] where \[p_m(x, y) = \sum_{i=0}^m \sum_{j=0}^m \rho_i \rho_j\Gamma_{ij} \widetilde{T}_i(x) \widetilde{T}_j(y), \quad \Gamma_{ij} = \int_{-1}^1 \int_{-1}^1 \widetilde{T}_i(x) \widetilde{T}_j(y) \,\textrm d\mu(x,y).\] Then \(\mu_{\textrm{KPM}}\) is a probability distribution and \[\int_{[-1,1]^2} f(x, y) \,\textrm d\mu_{\textrm{KPM}}(x, y)= \int_{[-1,1]^2}((K_m\otimes K_m)f)(x,y))\,\textrm d\mu(x, y).\]
Remark 9. We observe that 4 shows that the ansatz 6 constructs a measure \(\mu_{\textrm{KPM}}\) whose action on a test function \(f\) reproduces the action of the true measure \(\mu\) on the test function after polynomial smoothing, \(K_mf\). This is independent of the choice of the \(\rho_j\).
We may now prove 3.
Proof of 3. Let \(\mathcal{F}\) denote the class of 1-Lipschitz functions with respect to the Euclidean norm on \([-1, 1]^2\). Then applying 4 and 4 we obtain \[\begin{align} \mathcal{W}_1(\mu, \mu_{\textrm{KPM}}) & = \sup_{f \in \mathcal{F}} \left( \int_{-1}^1\int_{-1}^1f(x,y)\,\textrm d\mu(x,y) - \int_{-1}^1\int_{-1}^1f(x,y)\,\textrm d\mu_{\textrm{KPM}}(x,y)\right) \\ & = \sup_{f \in \mathcal{F}} \left( \int_{-1}^1\int_{-1}^1\left(f - (K_m \otimes K_m)f\right)(x,y)\,\textrm d\mu(x,y)\right) \\ & \leq \dfrac{12}{m}. \end{align}\] ◻
To be explicit about the algorithmic importance of 3, it shows that one may compute an approximation to a measure with just knowledge of its Chebyshev moments. The previous section showed how to compute those moments of \(\mu_{A,b}\), and so we have an algorithm, 3, for computing an approximation \(\widehat\mu_{A,b}\) of it.
2.3 Selection of weight vector.↩︎
For any fixed vector \(b\), the previous two subsections show how to compute an approximation of the weighted spectral density \(\mu_{A,b}\). Our goal, however, is to approximate the unweighted spectral density \(\mu_A\). Observe that if \(b\) is a random vector sampled from an isotropic distribution, then \[\mathbb{E}\mu_{A,b}=\mu_A.\] This section argues that \(\mu_{A,b}\) exhibits sufficient concentration when \(b\) is sampled uniformly at random from the unit sphere. The final algorithm will then be 4.
Our analysis is in two parts, both appealing to known results. Claim I is that the Chebyshev moments of \(\mu_{A,b}\) are highly concentrated around the corresponding moments of \(\mu_A\). Claim II is that this suffices to imply \(\mu_{A,b}\) is close to \(\mu_A\) in earth mover’s distance with high probability.
Claim I is a consequence of the following concentration result for Lipschitz functions on the sphere. We identify the unit sphere living in \(\mathbb{C}^n\) with the unit sphere \(\mathcal{S}^{2n-1}\) living in \(\mathbb{R}^{2n}\).
Lemma 5 ([33]). Suppose \(g:\mathcal{S}^{2n-1} \to \mathbb{R}\) has \(\ell_2\) Lipschitz constant bounded by \(L\). Let \(b \sim \mathbf{Unif}(\mathcal{S}^{2n-1})\). Then \[\Pr\mleft(\mleft|g(b) - \mathbb{E}[g(b)]\mright| \geq r\mright) \le 2e^{-(2n-1)r^2/(2L^2)}.\]
We can use this result to bound the fluctuations in the Chebyshev moments of \(\mu_{A,b}\) around its expectation \(\mu_A\).
Proposition 5. Let \(\Gamma_{jk}\) and \(\Gamma_{jk}'\) be the mixed Chebyshev moments of \(\mu_A\) and \(\mu_{A,b}\), respectively, for \(0 \leq j,k \leq m\). Let \(b \sim \mathbf{Unif}(\mathcal{S}^{2n-1})\). Then for \(0 < \delta < 1\), \[\mathbb{P}\left(\max_{(j,k) \neq 0}|\Gamma_{jk}-\Gamma_{jk}'| \geq \sqrt{\dfrac{8}{2n-1} \log\left(\dfrac{2(m+1)^2-2}{\delta}\right)} \right) \leq \delta.\]
Proof. For any \(f : \mathbb{C}\to \mathbb{R}\), and \(b \sim \mathbf{Unif}(\mathcal{S}^{2n-1})\), we have \[\int f(z) \,\textrm d\mu_{A}(z) = \dfrac{1}{n}\mathop{\mathrm{tr}}(f(A)), \quad \int f(z) \,\textrm d\mu_{A,b}(z) = b^*f(A)b.\] Moreover, \(\mathbb{E}bb^*=\frac{1}{n}I\) so \[\mathbb{E}[b^*f(A)b] = \mathbb{E}[\mathop{\mathrm{tr}}(f(A)bb^*)] = \dfrac{1}{n}\mathop{\mathrm{tr}}(f(A)).\] Observe also that \(g(b)\) defined by \(b \mapsto g(b)= b^*f(A)b\in\mathbb{R}\) is a Lipschitz function with Lipschitz constant bounded by \(2 {\left\lVert f(A) \right\rVert}_2\) since \[|b^*f(A)b-c^*f(A)c| = |b^*f(A)(b-c)-(c-b)^*f(A)c| \leq 2 {\left\lVert f(A) \right\rVert}{\left\lVert b-c \right\rVert}.\] Thus, according to 5, \[\label{a32} \mathbb{P}\left(\left|b^*f(A)b - \dfrac{1}{n}\mathop{\mathrm{tr}}(f(A))\right| \geq t\right) \leq 2 e^{-(2n-1)t^2/(8 {\left\lVert f(A) \right\rVert}^2)}.\tag{14}\] Let us now specialize to the case \(f_{ij}(z) = \widetilde{T}_{i}(x)\widetilde{T}_{j}(y)\) so that \[\Gamma_{ij}-\Gamma_{ij}' = \dfrac{1}{n}\mathop{\mathrm{tr}}(f_{ij}(A))-b^*f_{ij}(A)b,\] to which we may apply 14 with a union bound over the \((m+1)^2-1\) polynomials. Note also that \({\left\lVert f_{ij}(A) \right\rVert} \leq 1\), so for each \(0 < \delta < 1\) we may set \(t=\sqrt{\dfrac{8}{2n-1} \log\left(\dfrac{2(m+1)^2-2}{\delta}\right)}\) for the final result. ◻
Claim II is proven in [20]. They show that if two measures have close Chebyshev moments, they must be close in earth mover’s distance. Although the result is stated for \(d\)-dimensions in [20], we state it only for 2-dimensions here.
Lemma 6 ([20]). Suppose \(\mu, \nu\) are probability distributions on \([-1, 1]^2\) with normalized Chebyshev moments \(\Gamma_{kj}\) and \(\Gamma_{kj}'\), respectively. If the distributions’ normalized Chebyshev moments satisfy the bound \[\label{a34} \sum_{\substack{0 \leq k,j \leq m \\ (k,j) \neq (0, 0)}}\dfrac{1}{k^2+j^2} (\Gamma_{kj}-\Gamma_{kj}')^2 \leq K^2,\qquad{(3)}\] for some constant \(K > 0\), then one may bound the earth mover’s distance \[\mathcal{W}_1(\mu, \nu) \leq \dfrac{4C}{m} + \pi K,\] where \(\pi/2 \leq C \leq 3\).5
In order to combine [a31,a33] to imply a bound on \(\mathcal{W}_1(\mu_A,\mu_{A,b})\), we must bound the sum over the terms of the form \((k^2+j^2)^{-1}\). For this, we introduce the following lemma.
Lemma 7. One may bound \[\sum_{\substack{0 \leq k,j \leq m \\ (k,j) \neq (0, 0)}}\dfrac{1}{k^2+j^2} \leq 5 + \frac{\pi}{2} \log m.\]
Proof. We have \[\sum_{\substack{0 \leq k,j \leq m \\ (k,j) \neq (0, 0)}}\dfrac{1}{k^2+j^2} =2 \sum_{k=1}^m\dfrac{1}{k^2} + \sum_{i=1}^m \sum_{j=1}^m \dfrac{1}{i^2+j^2}.\] The first term is bounded by \(2\zeta(2)=\pi^2/3\) and the second term is bounded by \[\begin{align} \sum_{k=1}^m \sum_{j=1}^m \dfrac{1}{k^2+j^2} \le \sum_{k=1}^m \int_{0}^\infty \frac{1}{k^2+x^2} d x = \sum_{k=1}^m \frac{\pi}{2k} \le \frac{\pi}{2}(1+\log m). \end{align}\] Lastly, \(\pi^2/3+\pi/2<5\). ◻
2.4 Proof of 1↩︎
We now collect the results of the previous two subsections into a proof of our main result 1, which bounds the error in 4.
Proof of 1. Denote the output of 4 by \(\widehat \mu_A\), which is the KPM approximation to the random measure \(\mu_{A,b}\). Then observe that with probability \(1-\delta\) we have that \[\mathcal{W}_1(\mu_A, \widehat\mu_{A}) \leq \mathcal{W}_1(\mu_{A}, \mu_{A,b}) + \mathcal{W}_1(\mu_{A,b}, \widehat\mu_A) \leq O\left(\dfrac{1}{m}\right)+O\left(\sqrt{\dfrac{\log(m) \log(m/\delta)}{n}}\right).\] The first term \(\mathcal{W}_1(\mu_A, \mu_{A,b})\) is bounded via 6, within which the left-hand side of ?? is bounded with high probability using Hölder’s inequality together with 7 and 5. The second term \(\mathcal{W}_1(\mu_{A,b}, \widehat \mu_A)\) is bounded by 3.
Moreover, the measure \(\widehat\mu_A\) may be (optionally) discretized using Gauss-Chebyshev quadrature on \((m+1)^2\) points. That is, let \(x_1, \cdots, x_{m+1}\) and \(w_1, \cdots, w_{m+1}\) be the nodes and weights of the Gauss quadrature associated with the Chebyshev measure. Recall that Gauss quadrature on \(k\) nodes is exact for polynomials of degree \(2k-1\). Then consider the discrete measure \[\widehat \mu_{A, \textrm{discrete}} \mathrel{\vcenter{:}}= \sum_{i=1}^{m+1}\sum_{j=1}^{m+1} w_i w_j p_m(x_i, x_j) \delta_{(x_i, x_j)},\] where \(p_m\) is as in 6 . The exactness of Gauss quadrature ensures \[\int_{[-1,1]^2} T_i(x)T_j(y) (\,\textrm d\widehat \mu_A - \,\textrm d\widehat \mu_{A, \textrm{discrete}}) = 0\] for all \(i, j \leq m\). When combined with 6, this implies that \(\mathcal{W}_1(\widehat \mu_A, \widehat \mu_{A, \textrm{discrete}}) \leq 4C/m\), which merely alters the constant in the bound. ◻
2.5 Algorithmic variants↩︎
We include here some comments about algorithmic variants that may be useful.
It is possible to consider an adjoint-free version of 4 that avoids the use of the adjoint \(A^*\), which is an important feature in some applications, particularly PDEs [34]. One way to accomplish this is to use \(m-1\) matrix-vector products to form the explicit Krylov matrix \(\mathcal{K}_m = [b, Ab, \cdots, A^{m-1}b].\) In analogy to 5 , the Gram matrix \(\mathcal{K}_m^*\mathcal{K}_m\) contains the “polyanalytic” mixed moments corresponding to polynomials in \(z\) and \(\bar z\). Since there is an invertible linear change of coordinates from \(z=x+iy, \bar{z}=x-iy\) into \(x, y\), one could in principle pre-compute the change of basis coefficients from \(z^k\bar z^l\) into the Chebyshev polynomials \(\widetilde{T}_i(x) \widetilde{T}_j(y)\), though this may require high precision. One would likely use the Arnoldi procedure in practice rather than forming the ill-conditioned explicit Krylov matrix \(\mathcal{K}_m\).
It is possible to consider a low-memory variant of 4. We could compute \(X^*Y\) in 2 one row at a time, which requires \(O(m^2)\) matvecs but only \(O(n)\) storage.
Say the spectrum of \(A\) is somewhat oblong, e.g., the range of real parts is decently larger than the range of imaginary parts. Then rather than computing \(m\) matvecs in both the \(x\) and \(y\) coordinates, we could decouple the \(x\) and \(y\) directions into \(m_1\) and \(m_2\) matvecs, so that \(X^*Y \in \mathbb{R}^{m_1 \times m_2}\).
If the \(\widetilde{O}(n^{-1/2})\) error from 1 is too high, one could apply 4 to the block matrix \(A \otimes I_s \in \mathbb{C}^{ns \times ns}\), where \(I_s\) is the \(s\times s\) identity matrix and \(\otimes\) is the Kronecker product. This would decrease the error due to randomness to \(O((ns)^{-1/2})\). Note that this is equivalent to replacing \(\mu_{A,b}\) with a weighted average of \(s\) copies of \(\mu_{A,b}\) with independent \(b\). By picking \(s=\lceil m^{2/3}/n^{1/3}\rceil\), this improves the overall error to \(\widetilde{O}(1/(nm)^{1/3})\) when \(m=\Omega(\sqrt n)\).
Because of the decoupling of the \(x\)- and \(y\)-coordinate directions, the resulting measure seems to exhibit axis-aligned artifacts. One can remove these artifacts by “rotating” the input, i.e. applying the algorithm to \(e^{i\theta }A\) for a few angles \(\theta\), rotating back the returned measures, and averaging. That is, use the estimate \(\frac{1}{\mleft|\Theta\mright|}\sum_{\theta\in\Theta}\widehat\mu_{e^{i\theta}A}\circ(z\mapsto e^{i\theta}z)\) for some set of angles \(\Theta\). It may also be possible to work with other polynomial bases, such as disk polynomials [35].
3 Lower bounds for spectral density estimation↩︎
In this section, we prove a lower bound for all matrix-vector query algorithms. That is, we consider any algorithm interacts with the input matrix \(A\) only by constructing query vectors \(x\) and receiving \(Ax\) and \(A^*x\). We show that for large matrices, 4 makes the asymptotically lowest possible number of queries to achieve its error guarantee stated in 1. In fact, we show the lower bound holds even for algorithms only required to succeed when \(A\) is real symmetric.
The main result for this purpose, 6, is a characterization of all adaptive matrix-vector query algorithms which make a small number of queries relative to the matrix size. In particular, closely following ideas of Chewi et al [27], we show that any adaptive algorithm which makes at most \(m\) matvec queries to large \(A\) learns only a small amount of information about \(A\) beyond possibly the values of \(\mathop{\mathrm{tr}}(A^j)\) for \(j=1,\ldots,4m\). We believe this result may be independently useful for producing related lower bounds.
The second step of our lower bound, 8, is to show the values of \(\mathop{\mathrm{tr}}(A^j)\) for \(j=1,\ldots,4m\) is not enough information about \(A\) to reconstruct its spectral density to \(o(1/m)\) error. In particular, we produce a family of matrices with far spectral densities, but matching values of \(\mathop{\mathrm{tr}}(A^j)\).
In order to argue that matvec algorithms do not learn much information beyond \(\mathop{\mathrm{tr}}(A^j)\) for \(j=1,\ldots,4m\), we show that the output of any matvec algorithm can be simulated with high probability by an algorithm which takes as input just the values of \(\mathop{\mathrm{tr}}(A^j)\) for \(j\le 4m\). To make this formal, we need the notion of total variation distance, \[{\operatorname{dist}}_{\textrm{TV}}(X,Y)=\inf_{(X,Y)\sim\gamma}\Pr\mleft(X\neq Y\mright),\]where the infimum is over all joint distributions \(\gamma\) whose marginals are the laws of \(X\) and \(Y\). From this formulation of total variation distance, it is clear for every \(f\) that\[{\label{a38}{\operatorname{dist}}_{\textrm{TV}}(f(X),f(Y))\le{\operatorname{dist}}_{\textrm{TV}}(X,Y).}\tag{15}\] The notation \(X\overset d=Y\) means \(X\) and \(Y\) are equal in distribution, or equivalently \({\operatorname{dist}}_{\textrm{TV}}(X,Y)=0\). In this section, \(\otimes\) denotes the Kronecker product, \(\mathbf{1}_p\in\mathbb{R}^p\) the column vector of 1s, \(I_m\) the \(m\times m\) identity, and \(e_j\) the \(j\)th elementary basis column vector of the appropriate dimension.
Proposition 6. Let \(v^{\textrm{alg}}_0\in\mathbb{R}^n\) be any fixed vector and let \(v^{\textrm{alg}}_j:(\mathbb{R}^n)^j\to\mathbb{R}^n\) be measurable functions. Fix any real symmetric matrix \(B\in\mathbb{R}^{p\times p}\) and let \(A=U(I_{n/p}\otimes B)U^\top\) be the random matrix where \(U\) is a Haar-distributed random orthogonal matrix. If \[v_0=v^{\textrm{alg}}_0,\quad v_j=v^{\textrm{alg}}_j\mleft(Av_0,Av_1,\ldots,Av_{j-1}\mright).\] then for every \(f_{\textrm{alg}}:(\mathbb{R}^n)^m\to\mathcal{X}\) (for any output space \(\mathcal{X}\)), there is a measurable function \(f_{\textrm{sim}}:\mathbb{R}^{4m}\to\mathcal{X}\) such that \[{\operatorname{dist}}_{\textrm{TV}}\mleft( f_{\textrm{alg}}\mleft( Av_0,\ldots,Av_{m-1} \mright), f_{\textrm{sim}}(\mathop{\mathrm{tr}}(A),\mathop{\mathrm{tr}}(A^2),\cdots,\mathop{\mathrm{tr}}(A^{4m}))\mright)=O\mleft(\frac{m^{3/2}p^{1/2}}{n^{1/2}}\mright).\]
3.1 Proof of 2 assuming 6↩︎
In the following lemma we give a family of matrices with matching traces of powers but separated spectral densities.
Lemma 8. Let \(B_m\in\mathbb{R}^{m\times m}\) be the diagonal matrix with \((B_m)_{jj}=\cos\mleft(\frac{2j-1}{2m}\pi\mright)\). Let \(k\), \(k'\), \(m\), \(m'\), \(n\) be such that \(km=k'm' = n\). Then \(\mathop{\mathrm{tr}}((B_{m}\otimes I_{k})^i)=\mathop{\mathrm{tr}}((B_{m'}\otimes I_{k'})^i)\) for \(i\le\min(m,m')\) but if \(m\neq m'\) then \[\mathcal{W}_1(\mu_{B_{m}\otimes I_{k}},\mu_{B_{m'}\otimes I_{k'}})=\Omega\mleft(\frac{1}{\min(m,m')}\mright).\]
Proof. By Chebyshev quadrature, one has \(\frac{1}{m}\mathop{\mathrm{tr}}\mleft(B_m^i\mright)=\frac{1}{\pi}\int_{-1}^1\frac{x^i}{\sqrt{1-x^2}} \,\textrm dx\) for \(i\le m\). In particular, \[\mathop{\mathrm{tr}}\mleft(B_{km}^i\mright)=\frac{km}{\pi}\int_{-1}^1\frac{x^i}{\sqrt{1-x^2}}\,\textrm dx =k\mathop{\mathrm{tr}}\mleft(B_m^i\mright)=\mathop{\mathrm{tr}}\mleft((B_m\otimes I_k)^i\mright).\] Since \(km=k'm'\), we thus have \(\mathop{\mathrm{tr}}\mleft((B_m\otimes I_k)^i\mright)=\mathop{\mathrm{tr}}\mleft((B_{m'}\otimes I_{k'})^i\mright)\) for \(i\le\min(m,m')\). On the other hand, the earth mover’s distance between the spectral measures of \(B_{m}\otimes I_{k}\) and \(B_{m'}\otimes I_{k'}\) is the minimum \(\ell_1\) cost of a coupling between the diagonal entries of those two matrices. Specifically, denoting \(q_r(x)=\frac{\mleft\lfloor rx\mright\rfloor+\frac{1}{2}}{r}\) and replacing \(x=j/n\), we have \[\label{a39}\begin{align} W\equiv\mathcal{W}_1(\mu_{B_{m}\otimes I_{k}},\mu_{B_{m'}\otimes I_{k'}}) &= \frac{1}{n}\sum_{j=1}^{n}\mleft|(B_m\otimes I_k)_{jj}-(B_{m'}\otimes I_{k'})_{jj}\mright| \\&=\frac{1}{n}\sum_{j=0}^{n-1}\mleft|\cos\mleft(\frac{2\mleft\lfloor j/k\mright\rfloor+1}{2m}\pi\mright)-\cos\mleft(\frac{2\mleft\lfloor j/k'\mright\rfloor+1}{2m'}\pi\mright)\mright| \\&=\int_0^1\mleft|\cos(\pi q_m(x))-\cos(\pi q_{m'}(x))\mright|\,\textrm dx \\&=\int_0^1\mleft|\int_{q_{m'}(x)}^{q_m(x)}\sin(\pi t)\,\textrm dt\mright|\,\textrm dx \end{align}\tag{16}\] Let \(S\) be the set \[S={\left\{{(x,t)\in[0,1]^2: q_m(x)\le t\le q_{m'}(x)\text{ or } q_{m'}(x)\le t\le q_{m}(x)}\right\}}\] so that \(W\) is simply the integral of \(\sin(\pi t)\) over \(S\). Let \[w(t)=\max{\left\{{x:(x,t)\in S}\right\}}-\min{\left\{{x:(x,t)\in S}\right\}}.\] so that 16 reduces to \[{\label{a40}W=\int_0^1w(t)\sin(\pi t)\,\textrm dt.}\tag{17}\] Since \(\mleft|q_r(x)-x\mright|\le\frac{1}{2r}\), we have \(w(t)\le\frac{1}{m}+\frac{1}{m'}\). Suppose \(\int_0^1w(t)\ge c\mleft(\frac{1}{m}+\frac{1}{m'}\mright)\) for some positive constant \(c\). Then 17 would be minimized for \(w(t)=\frac{1}{m}+\frac{1}{m'}\) for \(t\le c\) and \(w(t)=0\) for \(t>c\). In particular, \[W \ge\mleft(\frac{1}{m}+\frac{1}{m'}\mright)\int_0^c\sin(\pi t)\,\textrm dt \ge\frac{1}{\pi}\mleft(\frac{1}{m}+\frac{1}{m'}\mright)\mleft(1-\cos(\pi c)\mright).\] It therefore suffices to prove this supposition. To this end, note \[\begin{align} \int_0^1w(t)\,\textrm dt &=\int_0^1\mleft|q_m(x)-q_{m'}(x)\mright|\,\textrm dx \\&=\frac{1}{n}\int_0^1\mleft|k\mleft\lfloor mx\mright\rfloor-k'\mleft\lfloor m'x\mright\rfloor+\frac{k-k'}{2}\mright|\,\textrm dx \\&=\frac{1}{n^2}\sum_{j=0}^{n-1}\mleft|k\mleft\lfloor j/k\mright\rfloor-k'\mleft\lfloor j/k'\mright\rfloor+\frac{k-k'}{2}\mright|\,\textrm dx \\&=\frac{1}{n^2}\sum_{j=0}^{n-1}\mleft|(j\bmod k)-(j\bmod k')+\frac{k-k'}{2}\mright|\,\textrm dx \end{align}\] Let \(d=\gcd(k,k')\). Then we group the terms in the sum based on the residue of \(j=id+r\) modulo \(d\). \[\begin{align} \int_0^1w(t)\,\textrm dt &=\frac{1}{n^2}\sum_{r=0}^{d-1}\sum_{i=0}^{n/d-1} \mleft|(j\bmod k)-(j\bmod k')+\frac{k-k'}{2}\mright| \\&=\frac{d^2}{n^2}\sum_{i=0}^{n/d-1} \mleft|(i\bmod k/d)-(i\bmod k'/d)+\frac{k-k'}{2d}\mright|.\end{align}\] Notice that the map \(i\mapsto(i\bmod k/d,i\bmod k'/d)\) is a \(n/\mathop{\mathrm{lcm}}(k,k')\)-to-1 map from \({\left\{{0,\ldots,n/d-1}\right\}}\) to \({\left\{{0,\ldots,k/d-1}\right\}}\times{\left\{{0,\ldots,k'/d-1}\right\}}\). That is, \[\int_0^1w(t)\,\textrm dt =\frac{d^2}{n^2}\frac{n}{\mathop{\mathrm{lcm}}(k,k')}\sum_{x=0}^{k/d-1}\sum_{y=0}^{k'/d-1}\mleft|x-y+\frac{k-k'}{2d}\mright| =\Theta\mleft(\frac{d^2}{n^2}\cdot\frac{n}{\mathop{\mathrm{lcm}}(k,k')}\cdot\frac{kk'}{d^2}\frac{k+k'}{d}\mright).\] Which, since \(kk'/d=\mathop{\mathrm{lcm}}(k,k')\), simplifies to \(\Theta\mleft( \frac{1}{m}+\frac{1}{m'} \mright)\) as desired. ◻
Proof of 2. Let \(B_{(\cdot)}\) be as in 8. Let \(m_1<\cdots<m_s\) be some divisors of \(n\). Set \[A_j = B_{m_j}\otimes I_{n/m_j}\] so that \(A_j\) is an \(n\times n\) matrix with a spectral measure matching that of \(B_{m_j}\). By 8, \[\mathop{\mathrm{tr}}(A_j^i)= \mathop{\mathrm{tr}}(A_{j'}^i)\] for all \(i\le m_1\) and \[\mathcal{W}_1(\mu_{A_j},\mu_{A_{j'}})>\frac{1}{Cm_{s-1}}.\] Let \(X_j\) be the result of some algorithm making at most \(m_1/4\) matrix vector queries on input \(UA_jU^\top\). By 6, each \(X_j\) are all within \(O\mleft(\frac{m_1^{3/2}m_s^{1/2}}{n^{1/2}}\mright)\) in total variation distance to the same random variable (produced by \(f_{\textrm{sim}}\)). Thus if the target accuracy \(\varepsilon\) is less than \(\frac{1}{2Cm_{s-1}}\), then with probability \(1-O\mleft(\frac{m_1^{3/2}m_s^{1/2}}{n^{1/2}}\mright)\) over the randomness of \(U\), the algorithm must fail for all but one of the \(s\) inputs. The algorithm therefore cannot succeed with overall probability more than \[\frac{1}{s}+O\mleft(\frac{m_1^{3/2}m_s^{1/2}}{n^{1/2}}\mright).\] Now, simply set \(m_j=\frac{1}{\varepsilon}+j-1\), \(s=1/\delta\), and \(n = \prod_{j=1}^s m_j\). ◻
3.2 Proof of 6↩︎
We now show how to simulate any matvec algorithm with an algorithm taking as input the traces of powers of \(A\).
Lemma 9 ([36]). Let \(v^{\textrm{alg}}_j\), \(v_j\), and \(A\) be as in 6. Let \(z_j\sim{\mathcal{N}}(0,I)\). Then there is a function \(g\) depending only on \(v^{\textrm{alg}}\) such that \[g({\left\{{A^iz_j}\right\}}_{i+j\le m})\overset d={\left\{{Av_j}\right\}}_{j<m}.\]
Lemma 10. Let \(A\) be as in 6 and \(z_j\sim{\mathcal{N}}(0,I)\). Then there is a function \(h\) such that \[h({\left\{{z_{j'}^TA^{i}z_j}\right\}}_{(i,j,j')\in H})\overset d={\left\{{A^iz_j}\right\}}_{i+j\le m}\] where \[H={\left\{{(i,j,j'):i+j+j'\le 2m,j,j'\le m}\right\}}.\]
Proof. Since \(A\) is real symmetric, \({\left\{{z_j^TA^iz_{j'}}\right\}}_{(i,j,j')\in H}\) is equivalent (up to some duplicated entries) to the set of inner products \(\label{a43}{\left\{{(A^{i'}z_{j'})^\top(A^{i}z_j)}\right\}}_{i+j\le m,i'+j'\le m}\) which \(h\) may compute a factorization of. For two collections of vectors arranged into the columns of matrices \(X\) and \(Y\), if all the inner products agree, i.e. \(X^*X=Y^*Y\), then there exists a orthogonal \(U\) such that \(X=UY\). Finally note by construction of \(A\) that\({\left\{{UA^iz_j}\right\}}_{i+j\le m}\overset{d}{=}{\left\{{A^iz_j}\right\}}_{i+j\le m}\), so applying a random orthogonal matrix to the factorization found by \(h\) achieves the desired distribution. ◻
The following 11 is a non-asymptotic version of a lemma appearing several places, such as [37], [38], and [39]. This non-asymptotic version was first stated without a proof in [27], and then again in a concurrent preprint [40]. A non-trivial refinement of the existing asymptotic proofs verifies the non-asymptotic version.6
11 will be used to argue that the bilinear forms \(z_{j'}^\top A^iz_j\) appearing in 10 look like Gaussian fluctuations around their means.
Lemma 11. Let \(Z\in\mathbb{R}^{t\times m}\) and \(X\in\mathbb{R}^{m\times m}\) have i.i.d. \({\mathcal{N}}(0,1)\) entries. Then \[{\operatorname{dist}}_{\textrm{TV}}\mleft(Z^\top Z, tI+\sqrt{\frac{t}{2}}\mleft(X+X^\top\mright)\mright)\le O\mleft(\frac{m^{3/2}}{t^{1/2}}\mright).\] The matrix \(Z^\top Z\) is called a Wishart matrix.
Lemma 12. Let \(Z\in\mathbb{R}^{tp\times m}\) have i.i.d. \({\mathcal{N}}(0,1)\) entries. Fix any symmetric \(B\in\mathbb{R}^{p\times p}\) and let \(A=B\otimes I_t\). Then there is a Gaussian matrix \(G\) with correlated entries entirely characterized by \({\left\{{\mathop{\mathrm{tr}}(A^i)}\right\}}_{i\le\ell}\) and \({\left\{{\mathop{\mathrm{tr}}(A^{2i})}\right\}}_{i\le\ell}\) such that \[{\operatorname{dist}}_{\textrm{TV}}\mleft(\begin{pmatrix}Z^\top A Z&\cdots&Z^\top A^\ell Z\end{pmatrix},G\mright)\le O\mleft(\frac{m^{3/2}}{t^{1/2}}\mright).\]
Proof. Let \(\lambda_1,\ldots,\lambda_p\) be the eigenvalues of \(B\) and denote the Vandermonde matrix \[V=\begin{bmatrix} \lambda_1&\cdots&\lambda_1^\ell\\\vdots&&\vdots\\ \lambda_p&\cdots&\lambda_p^\ell \end{bmatrix}.\] Let \(B=UDU^\top\) be the diagonalization of \(B\). Set \(\widetilde{Z}=(U^\top\otimes I_t)Z\). \[Z^\top A^iZ =Z^\top(B^i\otimes I_t)Z =\widetilde{Z}^\top(U^\top\otimes I_t)(B^i\otimes I_t)(U\otimes I_t)\widetilde{Z} =\widetilde{Z}^\top(D^i\otimes I_t)\widetilde{Z}\] If one breaks \(\widetilde{Z}\) into a column of \(p\) many \(t\times m\) blocks, then this becomes simply \[Z^\top A^iZ =\sum_{j=1}^p\lambda_j^i\widetilde{Z}_j^\top\widetilde{Z}_j.\] Rearranging this equality gives \[\begin{pmatrix}Z^\top AZ&\cdots&Z^\top A^\ell Z\end{pmatrix} \overset d= \begin{pmatrix}\widetilde{Z}_1^\top\widetilde{Z}_1&\cdots&\widetilde{Z}_p^\top \widetilde{Z}_p\end{pmatrix}\mleft(V\otimes I_m\mright).\] Notably, by rotational invariance of \(Z\), the matrices \(\widetilde{Z}_j^\top\widetilde{Z}_j\) on the right are all independent Wishart matrices. Set \[G={\left[{\mleft(t\mathbf{1}_p^\top\otimes I_m\mright)+\sqrt{t/2}\begin{pmatrix}X_1+X_1^\top&\cdots&X_p+X_p^\top\end{pmatrix}}\right]}\mleft(V\otimes I_m\mright)\] where \(X_j\) has i.i.d. \({\mathcal{N}}(0,1)\) entries and note that the entries of \(G\) are indeed Gaussian random variables as they are affine functions of the entries of the \(X_j\). By 15 and 11, we have \[\begin{align} &{\operatorname{dist}}_{\textrm{TV}}\mleft( \begin{pmatrix}Z^\top AZ&\cdots&Z^\top A^\ell Z\end{pmatrix},G\mright) \\&\le {\operatorname{dist}}_{\textrm{TV}}\mleft( \begin{pmatrix}\widetilde{Z}_1^\top\widetilde{Z}_1&\cdots&\widetilde{Z}_p^\top\widetilde{Z}_p\end{pmatrix},\mleft(t\mathbf{1}^\top\otimes I\mright)+\sqrt{t/2}\begin{pmatrix}X_1+X_1^\top&\cdots&X_p+X_p^\top\end{pmatrix} \mright) \\&\le O\mleft(\frac{m^{3/2}}{t^{1/2}}\mright). \end{align}\] Let’s now compute the first two moments of \(G\). Since \(\mathbb{E}(X)=0\), the first moment is simply \[\mathbb{E}(G) = (t\mathbf{1}_p^\top\otimes I_m)(V\otimes I_m) = (t\mathbf{1}_p^\top V\otimes I_m) = \mleft(\begin{pmatrix}\mathop{\mathrm{tr}}(A)&\cdots&\mathop{\mathrm{tr}}(A^\ell)\end{pmatrix}\otimes I_m\mright)\] so is determined by \({\left\{{\mathop{\mathrm{tr}}(A^i)}\right\}}_{i\le\ell}\). The covariance of the entries is given by \[\begin{align} \mathop{\mathrm{Cov}}\mleft( (1\otimes e_b^\top)G(e_i\otimes e_a),\, (1\otimes e_{b'}^\top)G(e_{i'}\otimes e_{a'}) \mright) &=\frac{t}{2} \sum_{j=1}^p \sum_{j'=1}^p \lambda_j^i \lambda_{j'}^i\mathbb{E}\mleft(e_b^\top(X_j+X_j^\top)e_ae_{b'}^\top(Y_{j'}+Y_{j'}^\top)e_{a'}\mright) \\&=\frac{t}{2} \sum_{j=1}^p \lambda_j^{2i} \mathbb{E}\mleft(e_b^\top(X_j+X_j^\top)e_ae_{b'}^\top(Y_{j}+Y_{j}^\top)e_{a'}\mright) \\&=\frac{1}{2}\mathop{\mathrm{tr}}(A^{2i})\mathbb{E}\mleft(e_b^\top(X_1+X_1^\top)e_ae_{b'}^\top(X_1+X_1^\top)e_{a'}\mright)\end{align}\] so is determined by \({\left\{{\mathop{\mathrm{tr}}(A^{2i})}\right\}}_{i\le\ell}\). ◻
Proof of 6. We describe the operations performed by \(f_{\textrm{sim}}\). First, because it is given \(\mathop{\mathrm{tr}}(A^i)\) for \(i\le 4m\), it is able to sample a Gaussian matrix \(G\) which, by 12 for \(\ell=2m\) and \(t=n/p\) and the definition of total variation distance, can be written as the mixture \[G=\begin{cases} \begin{pmatrix}Z^\top AZ&\cdots&Z^\top A^{2m} Z\end{pmatrix} & \text{w.p. } 1-O\mleft(\frac{m^{3/2}}{(n/p)^{1/2}}\mright)\\ \star & \text{o.w.} \end{cases}\] where \(Z\in\mathbb{R}^{n\times m}\) has i.i.d. \({\mathcal{N}}(0,1)\) entries and \(\star\) is some unspecified random variable. The blocks \(Z^\top A^iZ\) contain the inner products \({\left\{{z_j^\top A^iz_{j'}}\right\}}_{(i,j,j')\in H}\) where \(H\) is the index set in 10. By 9 and 10 respectively, there are \(g\) and \(h\) such that \[g\mleft(h\mleft({\left\{{z_j^\top A^iz_{j'}}\right\}}_{(i,j,j')\in H}\mright)\mright) \overset d=g\mleft({\left\{{A^iz_j}\right\}}_{i+j\le m}\mright) \overset d={\left\{{Av_j}\right\}}_{j<m}.\] From there, we may simply apply \(f_{\textrm{alg}}\). The result follows by 15 . ◻
4 Conclusion↩︎
In this work, we studied the problem of estimating the spectral density of a large normal matrix using only matrix-vector products with \(A\) and \(A^*\). Our main conclusion is that the spectral density of a normal matrix can be estimated with roughly the same complexity as a Hermitian matrix More specifically, we showed that when \(n \gg 1/\varepsilon^2\), for any normal matrix \(A\) with \({\left\lVert A \right\rVert}_2\leq 1\), it is sufficient to use \(O(1/\varepsilon)\) matrix-vector products in order to return with high probability an \(O(\varepsilon)\)-error spectral density estimation in earth mover’s distance. Moreover, we proved that this rate is asymptotically optimal, even for Hermitian input matrices, in that any algorithm with such a guarantee must in general make \(\Omega(1/\varepsilon)\) matrix-vector queries.
An important direction for future research is to determine whether our arguments extend beyond the normal matrix setting. Our algorithm and analysis relies crucially on the spectral theorem and on the commutativity of \(A\) and \(A^*\) or the Hermitian and skew-Hermitian parts of \(A\). For non-normal matrices, the matrix will not be unitarily diagonalizable and these commutativity relations break down. Moreover, the eigenvalues of \(A\) may not be well-conditioned. It would therefore be interesting either to identify useful classes of nearly normal matrices for which analogues of our results still hold or to prove hardness results showing that the normal case is a genuine boundary of tractability. In our view, understanding which parts of non-Hermitian spectral computation are rendered intractable by non-normality remains one of the central problems of numerical linear algebra.
Acknowledgments↩︎
This material is based upon work supported by the National Science Foundation under grant nos. DMS-2513687, AF-2427362, AF-2427363, AF-2046235, and AF-2045590. Nicholas thanks MathWorks for its support through the MathWorks Fellowship. We would like to thank Rajarshi Bhattacharjee for helpful conversations.
References↩︎
The desired values of \(c_1\) and \(c_2\) can be computed efficiently. For example, one may take \(c_2 = (x^*Ax)/(x^*x)\) with \(x\) sampled randomly; in fact, any element of the numerical range of \(A\) will suffice. Then since a constant-factor estimate of \({\left\lVert A-c_2I \right\rVert}_2\) suffices one may use a few steps of power iteration in order to estimate \(c_1 \leq 1/{\left\lVert A-c_2 I \right\rVert}_2\).↩︎
This multivariate generalization of the KPM was described briefly but not analyzed in [5].↩︎
Make the substitutions \(x = \cos \theta\) and \(r = \cos \phi\) into our definition of \(K_m\) and use the property of Chebyshev polynomials that \(T_j(\cos \theta) = \cos (j \theta)\).↩︎
While multivariate extensions of Jackson’s theorem have been known for a long time (c.f. [31], [32]), they do not seem to be well-known and we provide a simple proof of 4 using known properties of the one-dimensional \(K_m\).↩︎
We have modified the result of [20] by using \(\pi/2 \leq C \leq 3\), as derived by Jackson in [30], rather than the number 18 cited in Fact 12 of [20]. This improves the constant in [20] in the univariate case from \(c \leq36\) to \(c \leq 6\). A similar gain may be obtained by replacing the 18 in Fact 3.2 of [18], since both of these works rely on Jackson’s construction in [30].↩︎
A full proof was provided in private correspondence with the authors of [40].↩︎