6
$\begingroup$

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.

|improve this question
$\endgroup$
  • 3
    $\begingroup$ The Cholesky decomposition is defined only for Hermitian PD matrices. Are you sure your R function is not simply computing the SVD and taking the square root of the singular values before piecing it together and calling that as a "Cholesky decomposition"? I know that MATLAB does something like this in certain functions when the matrix is PSD and is a common way to do it. As to fixing roundoff errors in computed matrices, have you seen this question? $\endgroup$ – rm -rf♦ Oct 28 '12 at 18:31
  • $\begingroup$ @rm -rf R is doing cholesky decomposition, but it uses only the upper triangular portion of the original matrix. So do most other numerical systems. $\endgroup$ – asim Oct 28 '12 at 20:00
  • 2
    $\begingroup$ Why not automatically symmetrize all arguments to CholeskyDecomposition? Use something like CholeskyDecomposition[(m + Transpose[m])/2] instead of CholeskyDecomposition[m]. $\endgroup$ – whuber Oct 28 '12 at 20:11
  • 2
    $\begingroup$ I'm not too fond of how the question is phrased. It sounds like being forced to build a skyscraper out of marshmallows... maybe, "is it possible to have a version of 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$ – J. M. will be back soon♦ Oct 28 '12 at 23:30
4
$\begingroup$

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.

|improve this answer
$\endgroup$
  • 2
    $\begingroup$ The chol documentation you linked to states that "If [the matrix] is not [(complex Hermitian) symmetric], chol uses the (complex conjugate) transpose of the upper triangle as the lower triangle." This agrees with @asim's question, and contradicts your assertion that it uses the symmetric part of the matrix. $\endgroup$ – Rahul Oct 28 '12 at 20:46
3
$\begingroup$

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}}
|improve this answer
$\endgroup$
  • $\begingroup$ Thanks. This works, but will probably slow down things, when used on millions of iterations. $\endgroup$ – asim Oct 29 '12 at 19:50
  • $\begingroup$ You think? How big are your matrices? The effort taken for symmetrization is certainly not as much as the effort taken for the actual decomposition... $\endgroup$ – J. M. will be back soon♦ Oct 29 '12 at 21:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.