CholeskyDecomposition and Simplify

Question

I have a problem when trying to simplify a $6\times 6$ matrix after Cholesky decomposition. I tried all the regular operations such as

FullSimplify[Hnew[z], Element[z, Reals]]

or

$Assumptions = z ∈ Reals 

where my matrix is

Hnew[z_] := 
  {{1, 1 + z, 1 + z, 0, -3 z, -3 z}, {1 + z, 1, 1 + z, -3 z, 0, -3 z}, 
   {1 + z, 1 + z, 1, -3 z, -3 z, 0}, {0, -3 z, -3 z, 1, 1 + z, 1 + z}, 
   {-3 z, 0, -3 z, 1 + z, 1, 1 + z}, {-3 z, -3 z, 0, 1 + z, 1 + z, 1}}

but Mathematica evaluates indefinitely when I gave it

FullSimplify[CholeskyDecomposition[Hnew[z]], z > 0]

and it ignores assumptions. Also I tried Refine, Simplify and Assuming, but nothing makes Mathematica delete all the conjugates are reals. It just calculates so long, that I need to abort the calculation. Does anybody has experience with CholeskyDecomposition who is willing to help me out?

P.S. I'm new here.

Answer 1

The documentation for CholeskyDecomposition tells us the function argument must be a positive definite matrix. We can prove, however, that your matrix is not positive definite. Here's how:

For a matrix to be positive definite, all of its eigenvalues must be positive real numbers. So, we look at its eigenvalues, like this:

Eigenvalues[ { {1, 1 + z, 1 + z, 0, -3 z, -3 z}, 
    {1 + z, 1, 1 + z, -3 z,0, -3 z}, 
    {1 + z, 1 + z, 1, -3 z, -3 z, 0}, 
    {0, -3 z, -3 z, 1, 1 + z, 1 + z}, 
    {-3 z, 0, -3 z, 1 + z, 1, 1 + z}, 
    {-3 z, -3 z, 0, 1 + z, 1 + z, 1} } ]

(*  {3 - 4 z, -4 z, -4 z, 2 z, 2 z, 3 + 8 z}  *)

We quickly see there is no real value of $z$ that gives all positive eigenvalues, so CholeskyDecomposition should not be used.

Answer 2

Alternatively, one can use the $\mathbf L\mathbf D\mathbf L^\top$ decomposition to avoid the square roots needed by Cholesky. Using the routine in this answer, we get the diagonal factor $\mathbf D$ and check for conditions such that all of them are positive:

LDLT[mat_?SymmetricMatrixQ] := 
     Module[{n = Length[mat], mt = mat, v, w},
            Do[
               If[j > 1,
                  w = mt[[j, ;; j - 1]]; v = w Take[Diagonal[mt], j - 1];
                  mt[[j, j]] -= w.v;
                  If[j < n,
                     mt[[j + 1 ;;, j]] -= mt[[j + 1 ;;, ;; j - 1]].v]];
               mt[[j + 1 ;;, j]] /= mt[[j, j]],
               {j, n}];
            {LowerTriangularize[mt, -1] + IdentityMatrix[n], Diagonal[mt]}]

LDLT[{{1, 1 + z, 1 + z, 0, -3 z, -3 z},
      {1 + z, 1, 1 + z, -3 z, 0, -3 z},
      {1 + z, 1 + z, 1, -3 z, -3 z, 0},
      {0, -3 z, -3 z, 1, 1 + z, 1 + z},
      {-3 z, 0, -3 z, 1 + z, 1, 1 + z},
      {-3 z, -3 z, 0, 1 + z, 1 + z, 1}}] // Last // Simplify
   {1, -z (2 + z), -((z (3 + 2 z))/(2 + z)), (3 + 20 z)/(3 + 2 z),
    -((16 z (-3 - 8 z + 8 z^2))/(3 + 20 z)), (4 z (-9 - 12 z + 32 z^2))/(-3 - 8 z + 8 z^2)}

Reduce[And @@ Thread[% > 0], z]
   False

and thus, we come to the same conclusion as in Louis's answer.