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.
CholeskyDecomposition
? Use something likeCholeskyDecomposition[(m + Transpose[m])/2]
instead ofCholeskyDecomposition[m]
. $\endgroup$CholeskyDecomposition[]
that only needs the upper triangle of the symmetric positive definite matrix?" Because, hell, as already said, Cholesky's thing is really only intended for Hermitians... $\endgroup$