Hello.
I am running some Gaussian processes models. When computing the posterior mean I am facing some numerical errors in a noise free scenario even using float64 precision. I would like to now if pytorch does any approximation when calling the cholesky solver method.
Basically if X are my training locations and Y my observations. The posterior mean at locations X*=X should be equal to Y in a noise free scenario.
This posterior mean is given by: pos_mean = K(x*,x) @ K(x,x)^{-1}Y for a zero mean GP. In pytorch this is coded by:
# X represent training points and X* points where the posterior wants to be computed.
# This post asks for the case in which X*=X but the code presented is general.
# K_xsx_t represents the kernel K(X*,X).
# K_xx represents the kernel K(X,X)
Lxx_t = psd_safe_cholesky(Kxx.transpose(1,2),upper = False)
K_xsx_t = K_xsx.transpose(1,2)
# compute the solver and transpose
lhs = torch.cholesky_solve(Kxsx_t, Lxx_t, upper = False).transpose(1,2)
pos_mean = torch.matmul(lhs,Y)
torch.sqrt((pos_mean-Y)**2).sum() is different from zero or from a small value