Numerical stability regarding posterior inference

[Reimilia]

Problem Description

In phi_opt function from sparse_gp.py (matrix inversion)[https://github.com/mpnguyen2/motion_code/blob/26b7d279d64b714dd1a183dac33a87143ce82a87/sparse_gp.py#L201] uses np.linalg.inv, which is numerically stable and cause the symmetric matrix no longer symmetric. The inference covariance thus is not symmetric.

Proposed solution

@jit
def safe_symmetric_matrix_inv(A):
    """
        Use cholesky factorization to safely inverse a symmetric matrix
    """ 
    L = jsp.linalg.cho_factor(A)
    Ainv = jsp.linalg.cho_solve(L, jnp.eye(A.shape[0]))
    return Ainv

Note that gradient propagation in cho_solve function has some potential problems, but for inference it is better compared with linalg.inv. There are some marginal improvement on several datasets after replacing the matrix inversion function in phi_opt.