Implementation of predictive variance in Gaussian process regression of scikit-learn

Question

I'm studying the implementation of Gaussian Process Regression in scikit-learn to get a better understanding of the topic. There I've stumbled upon the following snippet:

if return_cov:
    v = cho_solve((self.L_, True), K_trans.T)  # Line 5
    y_cov = self.kernel_(X) - K_trans.dot(v)  # Line 6
    return y_mean, y_cov

It is from the file _gpr.py, around line 340.

The line numbers 5 and 6 in the comments apparently refer to the corresponding lines in algorithm 2.1 of the book Gaussian Processes for Machine Learning:

\ is the matrix inversion operator like in MATLAB.

My question is Why is it not $\mathbf{v}^T\mathbf{v}$ in sklearn but instead ${\mathbf{k}^\ast}^T\mathbf{v}$? I think it is wrong from dimensional analysis: $L$ has units of a standard deviation and $\mathbf{k}_\ast$ of a variance, so $\mathbf{v}$ has those of a standard deviation too. But then ${\mathbf{k}^\ast}^T\mathbf{v}$ has units of a standard deviation to the power of $1.5$ and not units of a variance.

Answer 1

This has to do with the fact, that cho_solve is different from \, and so sklearn's defintion of $v$ somewhat confusingly ends up not being the same as in the book.

In sklearn v_sk = cho_solve((self.L_, True), K_trans.T) would correspond to solving the system $K v_{sk} = \mathbf{k}_*^\top$ with $LL^\top=K$ being the Cholesky decomposition of training data covariance, self.L_$=L$, K_trans$=\mathbf{k}_*$ and K=$K$.

However line 5 in the book $\mathbf{v} := L \setminus \mathbf{k}_*$ corresponds to solving $ L \mathbf{v} = \mathbf{k}_*$. Substitute the corresponding definitions in the next line and you will end up with the same posterior covariance matrices.

For a better understanding of how cho_solve corresponds to forward/backward substitution compare how line 3 is implemented in _gpr.py :

$\mathbf{\alpha} := L^\top \setminus ( L \setminus \mathbf{y} )$

self.alpha_ = cho_solve((self.L_, True), self.y_train_)