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.