Approximate the singular values of a certain random dot-product kernel matrix (in the sense of El Karoui, Cheng-Singer, etc.) 
*

*Let $g:\mathbb R \to \mathbb R $ be a continuous function which is

*

*"sufficiently smooth" (e.g $\mathcal C^3$) around $0$, and

*"sufficiently integrable" (e.g integrable w.r.t $N(0,1)$).



*Let $d'$ and $d$ be positive integers tending to infinity such that $d'/d \to \rho \in (0,\infty)$.


*Let $X$ be an $d' \times d$ random matrix with iid rows $x_1,\ldots,x_D$ from $N(0,(1/d)I_d)$, let $\widetilde{X}$ be the $d'\times d$ random matrix with rows $\widetilde{x}_i := x_i/\|x_i\|$. Thus each $\widetilde{x}_i$ is uniformly distributed on the unit-sphere in $\mathbb R^{d}$.


*For $x \sim N(0,(1/d)I_d)$ independent of $X$, define an $d' \times d'$ random psd matrices $B$ and $\widetilde{B}$ by
$$
\begin{split}
b_{ij} &:= (x_i^\top x_j)\mathbb E_x[g(x^\top x_i)g(x^\top x_j)],\\
\widetilde{b}_{ij} &:= (\widetilde{x}_i^\top \widetilde{x}_j)\mathbb E_x[g(\widetilde{x}^\top \widetilde{x}_i)g(\widetilde{x}^\top \widetilde{x}_j)].
\end{split}
$$

Question. Are the following true ?
$$
\begin{split}
\|B-g(0)^2 XX^\top\|_{op} &\overset{p}{\to} 0,\\
\|\widetilde{B}-g(0)^2 \widetilde{X}\widetilde{X}^\top\|_{op} &\overset{p}{\to} 0.
\end{split}
$$

Empirical evidence
Here are results for some experiments with $d'=200$ and $d=300$, and different choices for $g$, both smooth and rough. Judgin from these observations, it seems the above question is affirmative!


An sloppy attempt

Let us give a heuristic (and incorrect!) argument why one would expect
$$
\|\widetilde{B}-g(0)^2\widetilde{X}\widetilde{X}^\top\|_{op} \overset{p}{\to} 0
$$
to hold.

For simplicity of notation, let $w_i:=\widetilde{x}_i$ henceforth. Let the matrix $\widetilde{U}$ be the $n \times n$ random psd matrix be defined by $\widetilde{u}_{ij} := \mathbb E_z[g(z^\top w_i)g(z^\top w_j)]$, where $z=(z_1,\ldots,z_d)$ is uniform on the unit-sphere on $\mathbb R^d$, and independent of the $w_i$'s.

Note that $\widetilde{B}$ the Hadamard product of $\widetilde{U}$ with $\widetilde{X}\widetilde{X}^\top$.

Now, due to rotational-invariance, we may write $\widetilde{u}_{ij} = u_d(z_i^\top z_j)$, where $u_d:[-1,1] \to \mathbb R$ is defined by
$$
u_d(t) := \mathbb E_z[g(z_1)g(tz_1+(1-t^2)^{1/2} z_2)].
$$
Thus, $\widetilde{B}$ is a dot-product kernel matrix, via an envelop function $u_d$ which varies with the dimension $d$. This dependence on $d$ is bad news for us.

Let us ignore the dependence of $u_d$ on $d$. Of course, this is incorrect (and will be the only sloppy part of our arguments), but let's do it anyways.

We can then invoke Theorem 2.3 of El Karoui '10 to get that $\|\widetilde{U}-A\|_{op} \overset{p}{\to} 0$, where
$$
A := u_d(0) 1_{d'}1_{d'}^\top + u_d'(0)\widetilde{X}\widetilde{X}^\top + \gamma I_{d'},
$$
with $\gamma := u_d(1)-u_d(0)-u_d'(0)$.
Finally, using exponential-in-$d$ concentration of $z_1^2 + z_2^2$, it is to see that

Fact. If $G:\mathbb R^2 \to \mathbb R$ is  a continuous function, then $\mathbb E_z[G(z_1,z_2)] \to G(0,0)$.

We deduce that if our $g$ is $\mathcal C^1$ at $0$, then

$$
\begin{split}
u_d(1) &= \mathbb E_z[g(z_1)^2] \to g(0)^2,\\
u_d(0) &=\mathbb E_z[g(z_1)g(z_2)] \to g(0)^2,\\
u_d'(0) &= \mathbb E_z[z_1 g(z_1)g'(z_2)] \to 0,\\
\gamma &= u_d(1)-u_d(0)-u_d'(0) \to 0,
\end{split}
$$

from which it would follow that $\|\widetilde{U}-g(0)^2 1_{d'}1_{d'}\|_{op} \overset{p}{\to} 0$, and so $\|\widetilde{B}-g(0)^2 \widetilde{X}\widetilde{X}^\top\|_{op} \overset{p}{\to} 0$, as claimed.
 A: 
Claim (Nonasymptotic result under smoothenss condition). Suppose $g$ is $\mathcal C^5$ at $0$ and that $d'$ and $d$ are sufficiently large with $c_1 \le n'/d \le c_2$ for some absolute constants $0 < c_2 \le c_2 < \infty$. Then w.p $1-e^{-Cd}$, it holds that
$$
\|\widetilde{B}-g(0)^2\widetilde{X}\widetilde{X}^\top\|_{op} = o_d(1).
$$

Of course, ideally, one would want to work with much less regularity (e.g think of $g(z) := 1 / (|z| + 1)$ in the experiments above).

Proof. Indeed, let $g(x) = a_0 + a_1 x + a_2x^2 + a_3 x^3 + a_4x^4 + a_5x^5 +  \mathcal O(x^6)$ be its Maclaurin approximation of order $5$. Fix $i,j \in [n]$, and let $a := \widetilde{x}_i = x_i/\|x\|$, $b := \widetilde{x}_j = x_j/\|x\|$, and $t := a^\top b$. Note that $|t| \le 1$ unconditionally.
By linearity of expectation, we have $\widetilde{u}_{ij} := \mathbb E_w[g(w^\top a)g(w^\top b)] = c_{ij} + \mathcal O(1/d^3)$, where
$$
c_{ij} := \sum_{n=0}^5\sum_{m=0}^5 a_0 a_m \mathbb E[(w^\top a)^n(w^\top b)^m].
$$
On the other hand, for any $n,m \ge 0$, we have

Fact. $\mathbb E_w[(w^\top a)^n(w^\top b)^m]=0$ if $n$ and $m$ have different parities, and otherwise
$$
\mathbb E_w[(w^\top a)^n(w^\top b)^m] = \frac{n!m!\Gamma(d/2)}{2^{n/2+m/2}\Gamma(d/2+n/2+m/2)}\sum_k \frac{t^k}{(n/2-k/2)!(m/2-k/2)!},
$$
where the sum is over all $k \in \{0,\ldots,\min(n,m)\}$ which have the same parity as $n$ and $m$.

Thus, $c_{ij}$ is a polynomial in $t$ of degree at most 5. The constant term (i.e the coefficient of $t^0$) in this polynomial is
$$
c_{ij}[t^0] = a_0^2 + \sum_{n,m}a_na_m\frac{n!m!}{2^{n/2+m/2}(n/2)!(m/2)!} \cdot \frac{\Gamma(d/2)}{\Gamma(d/2+n/2+m/2)}
$$
where the sum is over even $n,m \in \{0,\ldots,5\}$ not both equal to $0$.
For any such $n$ and $m$, we have $n/2 + m/2 \ge 1$ and so
$$
\dfrac{\Gamma(d/2)}{\Gamma(d/2+n/2+m/2)} \le \dfrac{\Gamma(d/2)}{\Gamma(d/2+1)} = \frac{2}{d}.
$$
We conclude that $c_{ij}[t^0] = a_0^2 + \mathcal O(1/d)$ for large $d$.
Similarly, for any $k \in \{1,\ldots,5\}$, we have that the coefficient of $t^k$  in $c_{ij}$ is
$$
c_{ij}[t^k] = \sum_{n,m}a_na_m\frac{n!m!}{2^{n/2+m/2}(n/2-k/2)!(m/2-k/2)!} \cdot \frac{\Gamma(d/2)}{\Gamma(d/2+n/2+m/2)}
$$
where the sum is over $n,m \in \{k,\ldots,5\}$ which have same parity as $k$. Note that for any such $n,m$ we have $n/2+m/2 \ge k$, and so
$$
\dfrac{\Gamma(d/2)}{\Gamma(d/2+n/2+m/2)} \le \dfrac{\Gamma(d/2)}{\Gamma(d/2+k)}=\dfrac{1}{(d+2(k-1))(d+2(k-1))\ldots d} \le \frac{1}{d^k},
$$
We deduce that $u_{ij}[t^k] = \mathcal O(1/d^k)$ for large $d$, and so because $\|\widetilde{X}\widetilde{X}^\top \|_{op} = \mathcal O(1)$ w.p $1-e^{-Cd}$ (by way of classical nonasymptotic RMT), we deduce that
$$
\|\widetilde{U}-(a_0^2+\mathcal O_d(1/d)) 1_{d'}1_{d'}^\top\|_{op}=o_{d}(1).
$$
Finally, since $(1_{d'}1_{d'}^\top) \circ (\widetilde{X}\widetilde{X}^\top) = \widetilde{X}\widetilde{X}^\top$, the claim follows.  $\quad\quad\quad\quad\Box$
