Something like MINRES or GMRES or CGSTAB. Correct me if I'm wrong, but the covariance matrix in GPR should be diagonally dominant, so the diagonal of the matrix can be used as a somewhat effective preconditioner.

With most kernels the matrix is symmetric and positive-definite so I would probably use a conjugate gradient: https://en.wikipedia.org/wiki/Conjugate_gradient_method

yes