Can the CholeskyDecomposition function in Mathematica be made to work on non-symmetric matrices?

Question

The CholeskyDecomposition[m] function in Mathematica requires a symmetric and positive definite matrix m.

For instance, the following works:

In[1]:= m = {{1, 0.5, 0.5}, {0.5, 1, 0.5}, {0.5, 0.5, 1}};

In[2]:= CholeskyDecomposition[m]

Out[2]= {{1., 0.5, 0.5}, {0., 0.866025, 0.288675}, {0., 0., 0.816497}}

However, consider another matrix m1

In[3]:= m1 = m; m1[[3, 2]] = 0.5000000000001

Out[3]= 0.5

In[4]:= m1

Out[4]= {{1, 0.5, 0.5}, {0.5, 1, 0.5}, {0.5, 0.5, 1}}

The Cholesky decomposition does not work here:

In[5]:= CholeskyDecomposition[m1]

During evaluation of In[5]:= CholeskyDecomposition::herm: The matrix {{1.,0.5,0.5},{0.5,1.,0.5},{0.5,0.5,1.}} is not Hermitian or real and symmetric. >>

Out[5]= CholeskyDecomposition[{{1, 0.5, 0.5}, {0.5, 1, 0.5}, {0.5, 
   0.5, 1}}]

Most other matrix based systems use either the lower triangular or upper triangular portion of a matrix when computing the Cholesky decomposition. For example, Eigen, LAPACK and R all do this. MATLAB offers many different versions of its chol function and it is possible to use either the upper or lower triangular portion. Mathematica does not offer this functionality.

For example, in R we can have

> m3<-matrix(c(1.0, 0.5, 0.5, 0.5, 1.0, 0.5000000000001, 0.5, 0.5, 1.0), 3,3)
> chol(m3)
     [,1]      [,2]      [,3]
[1,]    1 0.5000000 0.5000000
[2,]    0 0.8660254 0.2886751
[3,]    0 0.0000000 0.8164966
> 

Why is this important?

I find that checking for symmetry can be problematic in numerical computations where because of roundoff errors etc, slight asymmetries appear in matrices in iterative computations such as MCMC simulations. In such cases, the checking for symmetry results in the computation getting aborted midway.

Is there any way in which one can modify the built in CholeskyDecomposition function so that it does not check for symmetry, but aborts if the matrix is not positive definite? Such a modification would still retail all the speed benefits of the underlying LAPACK or MKL functions that Mathematica could be using but will not enforce symmetry.

Answer 1

The Cholesky decomposition only works for symmetric, positive definite matrices, though it can be generalized to complex-valued matrices as well. If you have a matrix $A$ then you can always make it into a symmetric matrix by taking $(A+A^\top)/2$. So if you have a lower triangular matrix, this operation would fill in the upper portion. If $A$ is full, then it just makes it symmetric (to which you can apply Cholesky). In case it's not obvious, this is just

a+Transpose[a]

in Mathematica. Indeed, this is what MATLAB (and most likely the others) are doing, as can be seen in the chol help file.

Answer 2

If all you really want is to construct a Hermitian positive definite matrix that is suitable for CholeskyDecomposition[]'s use from the upper triangle of your input matrix, here are a few preprocessing routines you can use:

toHermitian[mat_List?MatrixQ] := Module[{ma = mat}, 
          Do[ma[[(k + 1) ;;, k]] = Conjugate[ma[[k, (k + 1) ;;]]], {k, Length[ma] - 1}]; ma]

toHermitian[mat_SparseArray?MatrixQ] := Module[{ar = ArrayRules[mat], ru},
          ru = Cases[ar, ({i_, j_} -> _) /; i <= j];
          SparseArray[Join[ru,
                      Cases[ru, (({i_, j_} -> r_) /; i < j) :> ({j, i} -> Conjugate[r])]],
                      Dimensions[mat]]]

For instance:

CholeskyDecomposition[{{1, 0.5, 0.5}, {0.5, 1, 0.5}, {0.5, 0.5 + Sqrt[$MachineEpsilon], 1}}
                      // toHermitian]
   {{1., 0.5, 0.5}, {0., 0.866025, 0.288675}, {0., 0., 0.816497}}