Why is the posterior of a neural network gaussian process equal to the posterior of a neural network in the limit of infinite width layers? Its noted in this paper, Deep Neural Networks as Gaussian Processes, that "a single-layer fully-connected neural network with an i.i.d. prior over its parameters is equivalent to a Gaussian process (GP), in the limit of infinite network width".
I'm wondering, if that's the case, why so? And, how to prove it? And, under which conditions is this true for a particular model $f$?
My thoughts so far have been, if I get the kernel of a neural network $f_\theta$, I have
$$k(x,y) = cov(f_\theta(x), f_\theta(y)) = E_\theta[f_\theta(x)f_\theta(y)^T] = \Phi(x)\Phi(y)^T$$
for some $\Phi$.
If I then did some learning process over the linear learning problem $w\Phi(x)^T + \lambda ||w||_2$, and the equivalence claim above is true, then the posterior over $w$ in the linear problem would lead to the same model as the posterior over $\theta$ in $f_\theta$. This seems somewhat surprising to me? My intuition says that for some $f_\theta$s, some information in the original posterior would be lost when you take the covariance, so why does it hold for the covariance function of a neural network?
Thanks in advance. If there are other resources I can read up on this topic to understand it better, would be curious to learn more too!

See section 3 of the paper Emulating computer models with step-discontinuous outputs using Gaussian processes or section 15.4.5 of "Machine Learning: a probabilistic perspective" by Kevin Murphy.
Write your one-layer neural network of length $L$ and activation function $\sigma$ as
$$f(x) = b + \sum_{j=1}^L v_j \sigma(w_j x + b_j)$$
where the weights $w_j$ and biases $b_j$ are i.i.d. Gaussians with appropriate variances; similarly $v_j$ are i.i.d. Gaussian and $b$ is a mean-zero Gaussian. Denote all the model parameters as $\theta$, then
$$\mathbb{E}_\theta[f(x)] = 0 \quad \text{ and} \quad \text{Cov}(f(x), f(x^\prime)) = \sigma_b^2 + L \sigma_v^2 \mathbb{E}_{w, b}[\sigma(w x + b) \sigma(w x^\prime + b)].$$
Take $\sigma_v^2 = \omega^2/L$ for a constant $\omega^2$ , then by the central limit theorem, since $f(x)$ is the sum of appropriately scaled i.i.d. r.v.s we have that $f(x_1), \ldots, f(x_n)$ is asymptotically (in L) jointly Gaussian with mean zero and covariance function
$$k(x, x^\prime) = \omega^2 \mathbb{E}_{w, b}[\sigma(w x + b) \sigma(w x^\prime + b)]$$
for any finite collection of points $x_1, \ldots, x_n$.
Now intuitively, this means that $f(\cdot)$ is a Gaussian process, but you need some limiting argument to be rigorous, since you're considering convergence as $L \rightarrow \infty$ for all finite-dimensional marginals.
