Can the CholeskyDecomposition function in Mathematica be made to work on non-symmetric matrices? - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-11-12T20:51:03Z https://mathematica.stackexchange.com/feeds/question/13808 https://creativecommons.org/licenses/by-sa/4.0/rdf https://mathematica.stackexchange.com/q/13808 6 Can the CholeskyDecomposition function in Mathematica be made to work on non-symmetric matrices? asim https://mathematica.stackexchange.com/users/210 2012-10-28T18:12:39Z 2012-10-29T00:03:52Z <p>The <code>CholeskyDecomposition[m]</code> function in <em>Mathematica</em> requires a symmetric and positive definite matrix <code>m</code>.</p> <p>For instance, the following works:</p> <pre><code>In:= m = {{1, 0.5, 0.5}, {0.5, 1, 0.5}, {0.5, 0.5, 1}}; In:= CholeskyDecomposition[m] Out= {{1., 0.5, 0.5}, {0., 0.866025, 0.288675}, {0., 0., 0.816497}} </code></pre> <p>However, consider another matrix <code>m1</code></p> <pre><code>In:= m1 = m; m1[[3, 2]] = 0.5000000000001 Out= 0.5 In:= m1 Out= {{1, 0.5, 0.5}, {0.5, 1, 0.5}, {0.5, 0.5, 1}} </code></pre> <p>The Cholesky decomposition does not work here:</p> <pre><code>In:= CholeskyDecomposition[m1] During evaluation of In:= 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. &gt;&gt; Out= CholeskyDecomposition[{{1, 0.5, 0.5}, {0.5, 1, 0.5}, {0.5, 0.5, 1}}] </code></pre> <p>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 <code>chol</code> function and it is possible to use either the upper or lower triangular portion. <em>Mathematica</em> does not offer this functionality. </p> <p>For example, in R we can have</p> <pre><code>&gt; m3&lt;-matrix(c(1.0, 0.5, 0.5, 0.5, 1.0, 0.5000000000001, 0.5, 0.5, 1.0), 3,3) &gt; chol(m3) [,1] [,2] [,3] [1,] 1 0.5000000 0.5000000 [2,] 0 0.8660254 0.2886751 [3,] 0 0.0000000 0.8164966 &gt; </code></pre> <p>Why is this important? </p> <p>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. </p> <p>Is there any way in which one can modify the built in <code>CholeskyDecomposition</code> 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 <em>Mathematica</em> could be using but will not enforce symmetry.</p> https://mathematica.stackexchange.com/questions/13808/-/13809#13809 4 Answer by bill s for Can the CholeskyDecomposition function in Mathematica be made to work on non-symmetric matrices? bill s https://mathematica.stackexchange.com/users/1783 2012-10-28T18:32:47Z 2012-10-28T23:21:05Z <p>The <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition" rel="nofollow">Cholesky decomposition</a> 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</p> <pre><code>a+Transpose[a] </code></pre> <p>in <em>Mathematica</em>. Indeed, this is what MATLAB (and most likely the others) are doing, as can be seen in the <a href="http://www.mathworks.com/help/matlab/ref/chol.html" rel="nofollow"><code>chol</code> help file</a>.</p> https://mathematica.stackexchange.com/questions/13808/-/13823#13823 3 Answer by J. M. will be back soon for Can the CholeskyDecomposition function in Mathematica be made to work on non-symmetric matrices? J. M. will be back soon https://mathematica.stackexchange.com/users/50 2012-10-29T00:03:52Z 2012-10-29T00:03:52Z <p>If all you really want is to construct a Hermitian positive definite matrix that is suitable for <code>CholeskyDecomposition[]</code>'s use from the upper triangle of your input matrix, here are a few preprocessing routines you can use:</p> <pre><code>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_} -&gt; _) /; i &lt;= j]; SparseArray[Join[ru, Cases[ru, (({i_, j_} -&gt; r_) /; i &lt; j) :&gt; ({j, i} -&gt; Conjugate[r])]], Dimensions[mat]]] </code></pre> <p>For instance:</p> <pre><code>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}} </code></pre>