How to solve an eigensystem faster?
This is too long for a comment and honestly, to give a real answer, there is more information required in your question. Isn't it possible, that you give a working example, so that we see what takes long and how you implemented it?
If you are calling Eigensystem
for many different input values which are know, there is still some place for speed-up. Since your expressions are very lengthy, please find the initialization in an extra section.
First we measure how long it takes to calculate the Eigenvectors of H0
for 1 Million random values
data = RandomReal[{-1, 1}, {1000000, 7}];
First@AbsoluteTiming[Eigenvectors[H0[#]] & /@ data]
This took 44.4 sec here.
The next thing you can try is to distribute H0
over parallel kernels and use ParallelMap
DistributeDefinitions[H0];
First@AbsoluteTiming[ParallelMap[Eigenvectors[H0[#]] &, data]]
This took 25.3 sec with 4 subkernels. Let's test the compiled code. First when we apply it non-parallel
First@AbsoluteTiming[evectors @@@ data]
This took 2.4 sec which is almost 20 times faster then the initial version. Let's see what we can get if we call it parallel
First@AbsoluteTiming[evectors @@ Transpose[data]]
This took only 0.24 sec. If this scales well, that it means I can run 10 million samples in about 2.5 seconds. An indeed, a test with $10^7$ runs required 2.75 sec.
Now you might ask, whaat??, why is evectors @@@ data
a serial call while evectors @@ Transpose[data]
is parallel? It's because of the Listable
attribute in the Compile
-call and since we turn on Parallelization
. Sure is
evectors @@@ data == evectors @@ Transpose[data]
(* Out[21]= True *)
Initialization
Compiled parallel "C" versions of Eigenvalues and Eigenvectors
{evalues, evectors} =
Compile[{{ω0, _Complex, 0}, {ωx1, _Complex,
0}, {ωx2, _Complex, 0}, {ωy1, _Complex,
0}, {ωy2, _Complex, 0}, {ωz1, _Complex,
0}, {ωz2, _Complex, 0}}, #, Parallelization -> True,
CompilationTarget -> "C", RuntimeAttributes -> {Listable}] & /@
Eigensystem[{{0, (ωz1 - ωz2)/
2, (-ωx1 + ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2]), (ωx1 - ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2])}, {(ωz1 - ωz2)/2,
0, (ωx1 + ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2]), (ωx1 + ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2])}, {(-ωx1 + ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2]), (ωx1 + ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2]), ω0 + (ωz1 + ωz2)/2,
0}, {(ωx1 - ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2]), (ωx1 + ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2]), 0, -ω0 + 1/2 (-ωz1 - ωz2)}}];
Furthermore, I try to copy your approach by defining H0
H0[{ω0_, ωx1_, ωx2_, ωy1_, ωy2_,
ωz1_, ωz2_}] =
N[{{0, (ωz1 - ωz2)/
2, (-ωx1 + ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2]), (ωx1 - ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2])}, {(ωz1 - ωz2)/2,
0, (ωx1 + ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2]), (ωx1 + ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2])}, {(-ωx1 + ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2]), (ωx1 + ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2]), ω0 + (ωz1 + ωz2)/2,
0}, {(ωx1 - ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2]), (ωx1 + ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2]), 0, -ω0 + 1/2 (-ωz1 - ωz2)}}];
You may want to consider simplifying your matrix by substituting the [Omega] pairs that occur often with a dummy. This will simplify your matrix as there are only very few variables left. All you need to do is design a pattern to replace the most common [Omega] pairs (in that order). After the Eigensystem calculation, you can substitute them back.
I did something similar where several component were frequently occuring however were nested in divisions on several levels. Standard Simplify did not work, so I wrote my own which substantially increased performance of the simplified matrix.
The pattern I used was the simple Power[x_,y__] How you identify what combination to substitute depends on what's underlying your system, I am not really up to speed on Fourier transformations but this seems to be underlying what you are analysing here.
Below example gives your a transformastionRules fir the matrix and a reverseTransformation after the Eigensystem is calculated.
dummycount=1;
Position[HO,(*pattern*)];
H0[[#/.List->Sequence]]&/@Position[H0,(*pattern*)];
Sort[Cases[Tally[%],Except[{_,1}]],#1[[2]]>#2[[2]]&];
transformationRules=Rule[#[[1]],dummy[dummycount++]]&/@%;
reverseTransformation=transformationRules/.Rule[x_,y_]->Rule[y,x];
@halirutan : Thanks for the excellent example, I achieved several orders of magnitude speed-up for a program I was working on; however when I tested the execution of the compiled code without CompilationTarget->"C" something changes:
evectors@@@data==Eigenvectors[H0[#]]&/@data
(* False *)
Eigenvectors[H0[#]]&/@data==ParallelMap[Eigenvectors[H0[#]]&,data]
(* True *)
ParallelMap[Eigenvectors[H0[#]]&,data]==evectors@@Transpose[data]
(* False *)
evectors@@@data==evectors@@Transpose[data]
(* True *)
It seems to me the compiled functions behave differently from the non-compiled but I cannot discover where the difference is. I did create a trivial (2 dimension) example which simply swapped output and hence was easily corrected by
Reverse/@(evectors@@Transpose[data])
but for higher dimensions it is unclear what causes it to change its output.
Also, removing the RuntimeAttributes -> {Listable} attribute did not change anything. The problem seems to be that parameters are passed on to EigenSystem in listable manner and Eigensystem seems not to work in Listable mode; Eigensystem::matsq : "Argument {…} is not a non-empty square matrix.
However rewriting the code in what I believe to be functional identical by:
cpfnc = Compile[{{ω0,ωx1,ωx2,ωy1,ωy2,ωz1,ωz2},
Module[{eval,evec},{eval,evec} = Eigensystem[{{0, (ωz1 - ωz2)/
2, (-ωx1 + ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2]), (ωx1 - ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2])}, {(ωz1 - ωz2)/2,
0, (ωx1 + ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2]), (ωx1 + ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2])}, {(-ωx1 + ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2]), (ωx1 + ωx2)/(2 Sqrt[
2]) - (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2]), ω0 + (ωz1 + ωz2)/2,
0}, {(ωx1 - ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) - (I ωy2)/(2 Sqrt[
2]), (ωx1 + ωx2)/(2 Sqrt[
2]) + (I ωy1)/(2 Sqrt[2]) + (I ωy2)/(2 Sqrt[
2]), 0, -ω0 + 1/2 (-ωz1 - ωz2)}}];
{eval,evec}],
RuntimeAttributes -> {Listable}],Parallelization -> True];
does give the correct output, however again Listable provides problems for EigenSystem[] once parallel input is provided.
Love to have your views!