7
$\begingroup$

While doing a computation, I needed to take numerical derivatives of eigenvectors of a 4x4 hermitian matrix with respect to a parameter. I ran into the issue of phase jumps -- a random phase unpredictability added to numerically solved for eigenvectors -- and played with solutions like forcing the phase of the first element to be zero, setting the weighted average of phases zero, etc.

This is a 4x4 matrix though, and I could avoid jumps by solving for an analytic solution of the full eigensystem for an arbitrary hermitian matrix. I could then simply plug in my matrix values to get eigenvectors. I can do this twice and use finite difference to numerically compute the derivative. Here's what I did:

 x = {{aa/2, ab, ac, ad}, {0, bb/2, bc, bd}, {0, 0, cc/2, cd}, {0, 0,    0, dd/2}}; m = (x + ConjugateTranspose[x]); mEigensystem = Eigensystem[m, Cubics -> True, Quartics -> True]; give4x4Eigensystem[matrix_] :=   mEigensystem /. {aa -> matrix[[1, 1]], ab -> matrix[[1, 2]],     ac -> matrix[[1, 3]], ad -> matrix[[1, 4]], bb -> matrix[[2, 2]],     bc -> matrix[[2, 3]], bd -> matrix[[2, 4]], cc -> matrix[[3, 3]],     cd -> matrix[[3, 4]], dd -> matrix[[4, 4]]} 

Not the most elegant solution, but it gave me the right answers. Unfortunately, the calculation took surprisingly long.

I understand that the formulas I'm plugging elements into are very large, but they're stored in memory, and all mathematica has to do is plug numbers into these formulas and evaluate.

I'm not great at Mathematica so I'm wondering if there something obvious I'm missing that slows down this approach? Am I implicitly asking it to symbolically simplify the huge expressions involved? How could I speed up this approach?

Thanks!

Improve this question
$\endgroup$
6
  • $\begingroup$ You can gain a factor of 25 in speed by setting Cubics -> False, Quartics -> False. $\endgroup$ Commented Mar 14, 2020 at 12:12
  • $\begingroup$ Finite differences are usually a bad idea. Instead, matrix perturbation theory might be precisely what you need. $\endgroup$ Commented Mar 14, 2020 at 12:18
  • 1
    $\begingroup$ related: mathematica.stackexchange.com/q/100647/34893 $\endgroup$ Commented Mar 14, 2020 at 17:29
  • 1
    $\begingroup$ Thanks for the response -- why are finite differences a bad idea? $\endgroup$ Commented Mar 14, 2020 at 18:43
  • $\begingroup$ Eigenvectors are only defined up to scalar multiple. What exactly is it that you wish to do here? A full example would be helpful. $\endgroup$ Commented Mar 14, 2020 at 21:10

1 Answer 1

Reset to default
9
$\begingroup$

You can try to autocompile this:

x = {{aa/2, ab, ac, ad}, {0, bb/2, bc, bd}, {0, 0, cc/2, cd}, {0, 0, 0, dd/2}}; m = (x + ConjugateTranspose[x]); mEigensystem = Eigensystem[m, Cubics -> True, Quartics -> True]; evalMX = With[{mm = mEigensystem[[2]] // N}, Hold[compile[ {{mx, _Complex, 2}}, block[ { aa = mx[[1, 1]], ab = mx[[1, 2]], ac = mx[[1, 3]], ad = mx[[1, 4]], bb = mx[[2, 2]], bc = mx[[2, 3]], bd = mx[[2, 4]], cc = mx[[3, 3]], cd = mx[[3, 4]], dd = mx[[4, 4]] }, mm ], CompilationTarget -> "C" ] ] /. {compile -> Compile, block -> Block} // ReleaseHold ]; evalMX[{{1, 2, 3, 4}, {2, 3, 4, 1}, {3, 4, 1, 2}, {4, 1, 2, 3}}] // Re // RepeatedTiming {7.0*10^-6, {{-2.41421, -1., 2.41421, 1.}, {-1., 1., -1., 1.}, {0.414214, -1., -0.414214, 1.}, {1., 1., 1., 1.}}}
Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ This works incredibly well. I didn't even know this was possible in Mathematica. Thanks so much! $\endgroup$ Commented Mar 14, 2020 at 18:44
  • 2
    $\begingroup$ One of Mathematica's great uses is for automated code generation. Given how simple this kind of stuff is, you could imagine also marshalling down the expressions to something that could be used directly in another language if you primarily work in a different one. $\endgroup$ Commented Mar 14, 2020 at 20:46

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.