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
```