What is the difference between a covariance matrix created by an RBF kernel and a covariance matrix created by I can't explain something simple to myself and it is probably a matter of vocabulary, I am not sure...
If I create and random normal $Z \in \mathbb{R}^{3\times5}$, each row and column has a mean of 0. If I then wanted to created a covariance matrix of the first dimension, I could either calculate
$$
\frac{XX^\top}{n-1}
$$
which would give the average squared distance from the mean of 0 between each of the rows in $X$, which is the covariance. Alternatively, I could use the RBF kernel and calculate $RBF(x)$ which would also be called the covariance matrix between the rows of $X$. Unless I have gotten the approach wrong, these two things are not equal as the code demonstrates below. I would like to gain a deeper understanding of the difference between these two things and why they are both called covariance matrices even though they are different.
For example, could I use the $XX^\top$ covariance matrix for Gaussian process regression if it is a valid covariance matrix, or is there something that prevents it from working?
import gpytorch  # type: ignore
import torch  # type: ignore

RBF_KERN_FUNC = gpytorch.kernels.RBFKernel()


def run():
    x_tr = torch.randn(3, 5)
    cov_rbf = RBF_KERN_FUNC(x_tr).evaluate()
    print(f"rbf size: {cov_rbf.size()}")

    cov = x_tr @ x_tr.t()
    print(f"cov size: {cov.size()}")

    print(f"rbf: {cov_rbf}")
    print(f"cov: {torch.exp(-cov / 4)}")

if __name__ == '__main__':
    run()

Output:
rbf size: torch.Size([3, 3])
cov size: torch.Size([3, 3])
rbf: tensor([[1.0000, 0.0970, 0.0107],
        [0.0970, 1.0000, 0.0320],
        [0.0107, 0.0320, 1.0000]], grad_fn=<RBFCovarianceBackward>)
cov: tensor([[0.1803, 0.4715, 0.8288],
        [0.4715, 0.3840, 0.9206],
        [0.8288, 0.9206, 0.3947]])

 A: It appears that the matrix from rbf is a correlation matrix, with ones on the diagonal.  It's also a real-square symmetric matrix, since the upper triangular (off-diagonal elements) are a mirror image of the lower triangular.
The cov matrix looks like a covariance matrix, with diagonal elements that are the variance of each column, and off-diagonals are covariances.  You can check if the following are true:
$\sigma_{1}^2 = \frac{\sum_{i=1}^3 (x_{i1} - \bar{x}_1 )^2}{n-1} = 0.1803    \rightarrow\textrm{variance of feature 1}$
$\sigma_{12} = \frac{\sum_{i=1}^3 (x_{i1} - \bar{x}_1 )(x_{i2} - \bar{x}_2 )}{n-1} = 0.4715   \rightarrow\textrm{covariance between feature 1 and 2}$
If so, the cov is the covariance matrix.
I would not commingle kernel methods with the covariance matrix, mostly because kernel tricks could throw you off.  If you stick with statistical notation and calculations for obtaining the covariance matrix:
\begin{equation}
\boldsymbol{\sigma} = 
\begin{bmatrix}
\sigma_1^2       & \sigma_{12}    & \cdots  & \sigma_{14} \\
\sigma_{21}   &  \sigma_2^2    & \cdots  &  \sigma_{24} \\
\vdots  & \vdots  & \ddots  & \vdots  \\
\sigma_{41} & \sigma_{42}  & \cdots  &  \sigma_4^2  \\
\end{bmatrix},
\end{equation}
then you can't go wrong.   The correlation matrix is simply:
\begin{equation}
\boldsymbol{\rho}  = 
\begin{bmatrix}
1       & \rho_{12}    & \cdots  & \rho_{14} \\
\rho_{21}   &  1    & \cdots  &  \rho_{24} \\
\vdots  & \vdots  & \ddots  & \vdots  \\
\rho_{41} & \rho_{42}  & \cdots  &  1 \\
\end{bmatrix},
\end{equation}
where all diagonal elements are ones, and elements in the off-diagonal are calculated as
$\rho_{jk} = \frac{\sigma_{jk}}{\sigma_j \sigma_j}$,
where $\sigma_j $ = $\sqrt{\sigma_j^2}$, which are on the diagonal of the covariance matrix.
Now, more to the answer for what you did.  You simulated standard normal variates, $z_{ij}$, which have mean zero and variance unity, i.e. $N(0,1)$.  The correlation between feature $j$ and feature $k$ (in columns) is defined as
$\rho_{jk} = \frac{\sum_{i=1}^n \left( \frac{x_{ij} - \bar{x}_j}{\sigma_j} \right)  \left( \frac{x_{ik} - \bar{x}_k}{\sigma_k} \right)   }{n-1} $,
which is simply equal to the product of two z-variates, divided by $n-1$:
$\rho_{jk} = \frac{\sum_i^n z_{ij} z_{ik}}{n-1} $,
since
$z_{ij} =  \frac{x_{ij} - \bar{x}_j}{\sigma_j}$.
So your notation should be:
$\boldsymbol{\rho} = \frac{\mathbf{Z}^\top \mathbf{Z}}{n-1}$.
With $n$ rows and $p$ columns in a $\mathbf{Z}$ matrix, plugging in the dimensions, gives
$\underset{p \times p}{\boldsymbol{\rho}} = \frac{\underset{p \times n}{\mathbf{Z}^\top} \underset{n \times p}{\mathbf{Z}}}{n-1}$.
So for column-based correlation matrices, the tranposed matrix, $\mathbf{Z}^\top$, is on the left, and for column-based covariance matrices, $\mathbf{X}^\top$, is also on the left.
It does look like you simulated a $(3 \times 5)$ matrix, and since your matrices are ($3 \times 3)$, you ran row-based correlation and covariance -- which yields ($n \times n)$ matrices.  In statistics, it's more common to run correlation and covariance on columns to yield ($p \times p)$ matrices, where $j=1,2,\ldots,p$ represents the features (in columns).
