'Solves Ax = B using Cholesky decomposition,
for a general square non-singular matrix.
For a non-square matrix with rows > cols,
the least squares solution is returned.
See Also
lower_triangular_solve
upper_triangular_solve
gauss_jordan_solve
diagonal_solve
LDLsolve
LUsolve
QRsolve
pinv_solve'

def cholesky_solve(self, rhs):

    if self.is_symmetric():
        L = self._cholesky()
    elif (self.rows >= self.cols):
        L = (self.T * self)._cholesky()
        rhs = (self.T * rhs)
    else:
        raise NotImplementedError('Under-determined    System.    Try    M.gauss_jordan_solve(rhs)')
    Y = L._lower_triangular_solve(rhs)
    return L.T._upper_triangular_solve(Y)
