There is a method I have been studying called Spectral Normalised Neural Gaussian Processes which leaves me with a question I cannot answer. In this method, they utilize Random Fourier Features but instead of constructing the kernel via an inner product ...
$$ k(x_, x_j) = z(x_i)^\top z(x_j) $$
they instead use the random features to construct the Hessian of the last latent layer. (equation 9 in the SNGP paper uses $\Phi$ instead of $\mathbf{Z}$, but I use $\mathbf{Z}$ to highlight the connection to the RFF features.
$$ H = \sum_i(I + \hat{p}(1 - \hat{p}) \mathbf{Z}_i\mathbf{Z}_i^\top) $$
So my question is this? How is it justified to use the RFF features in constructing the Hessian when they were derived (in this case) to approximate the Gaussian kernel through inner products? Since the Hessian is the inverse Kernel in the Laplace approximation, it seems we are using the RFF features to construct $H = K^{-1}$ instead. IF it was this easy to construct a kernel inverse, it would be done all the time so I must be missing something important ...
