What is the difference between a covariance matrix created by an RBF kernel and a covariance matrix created by

Question

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]])

Answer 1

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).