1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Remove the Cholesky method and it works fine. What is the problem? $\endgroup$ – bill s Oct 21 '16 at 13:56
  • $\begingroup$ I want to use Cholesky method because it will boost my calculation time. I edited little bit to make the question clear. Thank you $\endgroup$ – Saesun Kim Oct 21 '16 at 14:10
  • $\begingroup$ Cholesky method is for Hermitian matrices (or symmetric in the case of real matrices). You matrix is not, why would you expect an improvement in time? The algorithm simply does not apply to your matrix. $\endgroup$ – anderstood Oct 21 '16 at 19:03
1
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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]
$\endgroup$
  • $\begingroup$ Thank you for your work. I am sorry that I did not accept you, but your answer is also good as accepted answer. I had hard time to choose one. Thank you! $\endgroup$ – Saesun Kim Oct 24 '16 at 13:36

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.