on eigenvalue estimate of kernel

🏷️Introduction

Given a positive definite kernel function $\mathcal{G}(x,y)$, it is natural to sample some points $(x_i,x_j)$ (so it creates a matrix $G$) as an approximation to the kernel function. It is a classical topic to estimate the eigenvalues of $G$ and compare with $\mathcal{G}$‘s eigenvalues.

🌊Weyl’s inequality

The simplest idea to quantify the difference is probably Weyl’s inequality for self-adjoint compact operators.

Estimate by Weyl's inequality

Let $\Omega =[0,1{]}^n$ be a unit cube and $\mathcal{G}(x,y)\in C^1(\Omega \times \Omega )$ is a positive definite kernel. The matrix $G\in {\mathbb{R}}^{m\times m}$ has entries as $\mathcal{G}(x_i,x_j)$ for $\{x_i{\}}_{i=1}^m$ be the regular “lattice points” with $\sqrt[n]{m}$ points along each axis. Then

$$|{\mu }_i(\mathcal{G})-\frac{1}{m}{\mu }_i(G)|=O\left(\frac{1}{\sqrt[n]{m}}\right)$$

where ${\mu }_i$ denotes the $i$th eigenvalue in descending order.

Proof

The lattice points automatically generate a partition of $\Omega$, we denote $C(x_i)$ as the small cube centered at $x_i$ which is disjoint from other small cubes.

The matrix $G$ constructs an approximated kernel function on each small cube around the lattice points. We define ${\mathcal{G}}^{\ast }(x,y)$ as

$${\mathcal{G}}^{\ast }(x,y)=\mathcal{G}(x_i,y_j)(x,y)\in C(x_i)\times C(y_j).$$

It is straightforward to see ${\mu }_i({\mathcal{G}}^{\ast })=\frac{1}{m}{\mu }_i(G)$. Then by Weyl’s inequality,

$$\begin{array}{rl}|{\mu }_i(\mathcal{G})-{\mu }_i({\mathcal{G}}^{\ast })|& \le \Vert \mathcal{G}-{\mathcal{G}}^{\ast }{\Vert }_{op}\le \sqrt{{\int }_{\Omega \times \Omega }|\mathcal{G}(x,y)-{\mathcal{G}}^{\ast }(x,y){|}^{2}dxdy}\\ & =\sqrt{\sum _{i=1}^m\sum _{j=1}^m{\int }_{C(x_i)\times C(x_j)}|\mathcal{G}(x,y)-{\mathcal{G}}^{\ast }(x,y){|}^{2}dxdy}\\ & =O\left(\frac{\Vert \nabla \mathcal{G}{\Vert }_{\infty }}{\sqrt[n]{m}}\right).\end{array}$$

🌵Ostrowski Theorem

Intuitively, the insufficiency of samples will only affect smaller eigenvalues (more oscillatory eigenfunctions), thus isolating the “safer” eigenfunctions is natural.

Let the kernel $\mathcal{G}$ be expanded into its eigenfunctions:

$$\mathcal{G}(x,y)=\sum _{k=1}^{\infty }{\mu }_k(\mathcal{G}){\psi }_k(x){\psi }_k(y).$$

And we pick a truncation ${\mathcal{G}}_l:={\sum }_{k=1}^l{\lambda }_k{\psi }_k(x){\psi }_k(y)$. Then the matrix $G$ can be split into two parts:

$$G_{i,j}={\mathcal{G}}_l(x_i,x_j)+(\mathcal{G}(x_i,x_j)-{\mathcal{G}}_l(x_i,x_j)).$$

The first part should be less affected by the sampling, thus we should expect the correlation

$$\frac{1}{m}\sum _{i=1}^m{\psi }_k(x_i){\varphi }_r(x_i)\approx {\delta }_{k,r}$$

The common tool we use is Ostrowski’s theorem(Bhatia, 2013).

Ostrowski theorem

If $A\in {\mathbb{C}}^{n\times n}$ is Hermitian and $X\in {\mathbb{C}}^{n\times n}$, then

$${\mu }_i(X^{\ast }AX)={\theta }_i{\mu }_i(A),i\in [n].$$

where ${\mu }_n(X^{\ast }X)\le {\theta }_i\le {\mu }_1(X^{\ast }X)$.

Let $({\Phi }_l{)}_{i,j}={\psi }_j(x_i)$ and ${\Lambda }_l=\mathrm{diag}({\mu }_i(\mathcal{G}){)}_{i=1}^l$, then the theorem implies

$$\left|\frac{1}{m}{\mu }_i({\Phi }_l{\Lambda }_l{\Phi }_l^{\ast })-{\mu }_i(\mathcal{G})\right|\le {\mu }_i(\mathcal{G}){\Vert \frac{1}{m}{\Phi }_l^T{\Phi }_l-I{d}_l\Vert }_{op}.$$

Then, we obtain the bound (Ostrowski theorem and Weyl’s theorem, respectively)

$$\begin{array}{rl}|{\mu }_i(\mathcal{G})-\frac{1}{m}{\mu }_i(G)|& \le |{\mu }_i(\mathcal{G})-\frac{1}{m}{\mu }_i({\Phi }_l{\Lambda }_l{\Phi }_l^T)|+|\frac{1}{m}{\mu }_i({\Phi }_l{\Lambda }_l{\Phi }_l^T)-\frac{1}{m}{\mu }_i(G)|\\ & \le {\mu }_i(\mathcal{G}){\Vert \frac{1}{m}{\Phi }_l^T{\Phi }_l-I{d}_l\Vert }_{op}+\underset{\text{by Weyl’s inequality}}{\underbrace{\frac{1}{m}\Vert G-{\Phi }_l{\Lambda }_l{\Phi }_l^{\ast }{\Vert }_{op}}}.\end{array}$$

🌀Conservative estimates

The remainder term $G-{\Phi }_l{\Lambda }_l{\Phi }_l^{\ast }={\sum }_{k>l}{\mu }_k{\psi }_k(x_i){\psi }_k(x_j)$, a naive bound of the 2nd term will be $C{\sum }_{k>l}{\mu }_k$ if the eigenfunctions are uniformly bounded, otherwise the growth needs to be considered. Of course, this bound is quite conservative. Few possible directions are viable.

• The $\{{\psi }_k(x_i){\}}_{i=1}^m$ can be viewed as some random vector for large $k$, then the central limit theorem can create a more aggressive estimate, see related work(Chun, Chung & Lee, 2025; Kong & Valiant, 2017).
• If ${\psi }_k$ is close to a standard Fourier mode (up to a fixed phase shift), then on an equispaced lattice of $\{x_i{\}}_{i=1}^m$, it is still expecting the (somewhat) orthogonality between the discrete vectors ${\Phi }_l$.

Random vector perspective

Consider $\{{\psi }_k(x_i){\}}_{i=1}^n$ as a random vector with respect to i.i.d random variables $x_i$. We set the truncated kernel (still positive definite) $g(x,y)=\mathcal{G}(x,y)-{\mathcal{G}}_l(x,y)$, then we can use the spectrum estimator idea, take $\hat{m}(n)$ as a moment estimator (biased), then

$$\begin{array}{r}{\mathbb{E}}_x(\hat{m}(n))=\mathbb{E}\left[\mathrm{tr}\left((\frac{1}{m}G-\frac{1}{m}{\Phi }_l{\Lambda }_l{\Phi }_l^{\ast }{)}^n\right)\right]\end{array}$$

For instance, $n=2$,

$$\begin{array}{rl}\mathbb{E}\left[\mathrm{tr}\left((\frac{1}{m}G-\frac{1}{m}{\Phi }_l{\Lambda }_l{\Phi }_l^{\ast }{)}^2\right)\right]& =\frac{1}{m^2}\sum _{i,j=1}^m\mathbb{E}\left[g(x_i,x_j)g(x_j,x_i)\right]\\ & =\frac{1}{m^2}\sum _{k,s>l}\sum _{i,j\ne i}^m{\mu }_k{\mu }_s\mathbb{E}\left[{\psi }_k(x_i){\psi }_s(x_i)\right]\mathbb{E}[{\psi }_k(x_j){\psi }_s(x_j)]\\ & +\frac{1}{m^2}\sum _{k,s>l}\sum _{i=1}^m{\mu }_k{\mu }_s\mathbb{E}\left[{\psi }_k(x_i){\psi }_s(x_i){\psi }_k(x_i){\psi }_s(x_i)\right]\\ & =\frac{m-1}{m}\sum _{k>l}{\mu }_k^2+\frac{1}{m}{\int }_{D}g(x,x{)}^{2}dx.\end{array}$$

The last term is bounded by $m^{-1}O({\sum }_{k>l}{\mu }_k{)}^2$. Once the decay rate of ${\mu }_k=O(k^{-\beta })$, we can expect a bound of $O(l^{1-2\beta }+m^{-1}{l}^{2-2\beta })$.

The variance can be estimated similarly. We only need

$$\begin{array}{rl}\mathbb{E}\left[\mathrm{tr}{\left((\frac{1}{m}G-\frac{1}{m}{\Phi }_l{\Lambda }_l{\Phi }_l^{\ast }{)}^2\right)}^2\right]& =\mathbb{E}\left[\frac{1}{m^4}\sum _{i,j,s,t=1}^{m}g(x_i,x_j{)}^{2}g(x_s,x_t{)}^2\right]\\ & =\mathbb{E}\left[\frac{1}{m^4}\sum _{\text{distinct},i,j,s,t=1}^{m}g(x_i,x_j{)}^{2}g(x_s,x_t{)}^2\right]+O\left(\frac{1}{m}\right)(\sum _{k>l}{\mu }_k{)}^4\\ & =(1-O\left(m^{-1}\right)){\left[{\int }_{D\times D}|g(x,y){|}^{2}dxdy\right]}^2+O\left(\frac{1}{m}\right)(\sum _{k>l}{\mu }_k{)}^4\\ & =(1-O\left(m^{-1}\right))(\sum _{k>l}{\mu }_k^2{)}^2+O\left(\frac{1}{m}\right)(\sum _{k>l}{\mu }_k{)}^4\end{array}$$

Therefore, the variance is bounded by $O(m^{-1})(({\sum }_{k>l}{\mu }_k^2{)}^2+({\sum }_{k>l}{\mu }_k{)}^4)=O(m^{-1})l^{4-4\beta }$ and

$$\Vert \frac{1}{m}G-\frac{1}{m}{\Phi }_l{\Lambda }_l{\Phi }_l^{\ast }{\Vert }_{op}^2\le \hat{m}(2)\le O_p(l^{1-2\beta }+m^{-1/2}{l}^{2-2\beta })).$$

Extension to higher moments

When a higher moment $m(n)$ is involved, we should still expect something like

$$O(1)\sum _{k>l}{\mu }_k^n+O_p(m^{-1/2})(\sum _{k>l}{\mu }_k{)}^n=O_p(l^{1-n\beta }+m^{-1/2}{l}^{n-n\beta })$$

as the upper bound. As we can see later, $n=2$ probably is the best we can do for the condition number estimate (since the first term dominates for $n=1$ and $l=m$).

〰️Concentration argument

The first term can be quantified by viewing $x_i$ as a certain quadrature rule, or ready for Hoeffding’s inequality. We should expect the stochastic bound $O(m^{-1/2}\sqrt{\log l^2/p})$ for each entry with probability $1-\frac{p}{l^2}$, therefore, with probability $1-O(p)$,

$${\Vert \frac{1}{m}{\Phi }_l^T{\Phi }_l-I{d}_l\Vert }_{op}=O\left(m^{-1/2}l\sqrt{\log l^2/p}\right).$$

The final bound will be a compromise between these two terms (in high probability). Under the assumption that ${\mu }_k(\mathcal{G})=O(k^{-\beta })$ that $\beta >1$,

$$|{\mu }_i(\mathcal{G})-\frac{1}{m}{\mu }_i(G)|\le {\mu }_i(\mathcal{G})O_p(m^{-1/2}l\sqrt{\log l^2/p})+O_p\left(l^{1/2-\beta }+m^{-1/4}{l}^{1-\beta }\right).$$

Then, we should truncate at $i^{-\beta }{m}^{-1/2}l\approx l^{-\beta +1/2}$ if $l\ge \sqrt{m}$, that is, $l\approx m^{1/(2\beta +1)}{i}^{\beta /(\beta +1/2)}$ if it is not exceeding $\sqrt{m}$. Otherwise, $i^{-\beta }{m}^{-1/2}l\approx m^{-1/4}{l}^{1-\beta }$, that is, $l\approx i\cdot m^{1/2\beta }$ if not exceeding $m$. The condition number of the sample matrix can be estimated from below.

$$\kappa =\frac{{\mu }_1(G)}{{\mu }_m(G)}={\Omega }_p(m^{\beta -3/4}).$$

Since ${\mu }_1(G)=\Theta (m)$, and ${\mu }_m(G)=O(m^{-\beta +3/4})$ in high probability sense.

🔦 Examples

Two-layer ReLU neural network with random sampling

Consider the two-layer ReLU network’s NTK regime (considered in on two-layer ReLU networks), the kernel permits a decay rate ${\mu }_k\sim O(k^{-(n+3)/n})$, then we can expect the condition number is bounded below by $O(m^{1/4+3/n})$ for a total number of $m$ samples.

Two-layer ReLU neural network with equispaced sampling

In this case, the estimate can be more aggressive since the eigenfunctions are almost “Fourier” modes. Because the equispaced lattice points $\{x_i{\}}_{i=1}^m$ still sustain the orthogonality relation for the sampled eigenfunctions. In this scenario, we can expect ${\mu }_m(G)=O(m^{-\beta })$.

References

🐻  Bhatia, R. 2013. Matrix Analysis, Springer Science & Business Media,p.

🐻  Chun, C., Chung, S. & Lee, D.D. 2025. Estimating the Spectral Moments of the Kernel Integral Operator from Finite Sample Matrices.

🐻  Kong, W. & Valiant, G. 2017. Spectrum Estimation from Samples.