Solving the matrix with assumption?

Question

I want to solve equations analytically and found my equation can be composed as Hermitian matrix. So I found on MMA that there is a Cholesky method that I can use to solve this matrix. However, since it is not the numerical solution, I am having little trouble solving my equations. Here is the simple example what I am struggling with

LinearSolve[{{a, b}, {Conjugate[b], c}}, {d, e}, Method -> "Cholesky"]

LinearSolve::herm: The matrix {{a,b},{Conjugate[b],c}} is not 
Hermitian or real and symmetric.

If I do not use the methods, then I will get the answer.

LinearSolve[{{a, b}, {Conjugate[b], c}}, {d, e}]

$\left\{\frac{c d-b e}{a c-b b^*},\frac{a e-d b^*}{a c-b b^*}\right\}$

However, I want to use the method to boost speed of my calculation. Do you have any suggestion?

Answer 1

The fact that Mathematica does not provide a CholeskyDecomposition for symbolic matrices suggests that it may not be very useful. However, it is possible to find a decomposition manually. All the assumptions made in the following are consequences of the matrix M being positive-definite Hermitian.

Given a matrix

M = {{a, b}, {Conjugate[b], c}};

Find a matrix

L = {{u, 0}, {v, w}};

satisfying

Assuming[{u, w} ∈ Reals, Refine[ConjugateTranspose[L].L]] == M // LogicalExpand
(* w^2 == c && Conjugate[b] == v w && w Conjugate[v] == b && 
 u^2 + v Conjugate[v] == a *)

Solve this one component at a time, taking the most convenient solution when more than one is offered

Last[Solve[First[%], w]]
(* {w -> Sqrt[c]} *)

%% /. %
(* Conjugate[b] == Sqrt[c] v && Sqrt[c] Conjugate[v] == b && 
 u^2 + v Conjugate[v] == a *)

First[Solve[First[%], v]]
(* {v -> Conjugate[b]/Sqrt[c]} *)

Assuming[c > 0, FullSimplify[%% /. %]]
(* a == u^2 + (b Conjugate[b])/c *)

Last[Solve[%, u]]
(* {u -> Sqrt[a c - b Conjugate[b]]/Sqrt[c]} *)

This gives the decomposition as

L /. % /. %%% /. %%%%%
(* {{Sqrt[a c - b Conjugate[b]]/Sqrt[c], 0}, {Conjugate[b]/Sqrt[c], Sqrt[c]}} *)

This satisfies the required equality

Assuming[c > 0 && a > 0 && a c > b Conjugate[b] && a c > Abs[b]^2, 
 M == ConjugateTranspose[%].% // FullSimplify]
(* True *)

Extending this to higher dimensions is left as an exercise to the determined reader.

Answer 2

Here is a method that is independent of the dimensions of the matrix. This method is just a way of telling MMA that our matrix is indeed hermitian. First, define only the upper triangle of the matrix as

mat = {{a, b}, {0, c}};
v = {d, e};

Explicitly take the real part of the diagonal elements and the hermitian transpose of the upper triangle. Put these pieces together to obtain a matrix that is explicitly hermitian.

diag = DiagonalMatrix[Re[Diagonal[mat]]];
upper = UpperTriangularize[mat, 1];
lower = Transpose[Conjugate /@ upper];
m = diag + lower + upper;
HermitianMatrixQ[m]

(* True *)

Apply the Cholesky method, but the result is rather ugly, so we don't show it.

s = LinearSolve[m, v, Method -> "Cholesky"];

We do not really want all of the expressions like Re[a] that we introduced, so we create a list of rules and apply that list to the solution s to obtain a cleaner solution t, which is still too ugly to show. Here again we are using the fact that our diagonal elements are real numbers.

r = Thread[Diagonal[diag] -> Diagonal[mat]];
t = s /. r;

We would like to understand how that expression could be the same as the much simpler expression we get when we use the default method instead of Cholesky. If we look at the special case in which a>0, the result is still not a bit more complicated than the default solution.

Assuming[a > 0,
  Refine[t]
  ] // FullSimplify

We end up an expression that has

a Sqrt[c - (b Conjugate[b])/a] Conjugate[Sqrt[c - (b Conjugate[b])/a]]

in the denominator. After studying it, we can see that it is a c - b Conjugate[b], which is what we get using the default method.

Edit:

More succinctly, if we know all of our symbols are real numbers, we can wrap them with Re before sending them to Cholesky. Then, when get the results, we remove the wrappers. Still, the result too ugly to show. These commands do the trick:

m = {{a, bx - I by}, {bx + I by, c}};
symb = Cases[Flatten[Apply[List, m, Infinity]], _Symbol] // Union;
herm = m /. Thread[symb -> Re[symb]]
LinearSolve[herm, {d, e}, Method -> "Cholesky"] /. 
 Thread[Re[symb] -> symb]