Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a module that I need to call 1-10 million times in my program. Currently, it is taking several hours to run so I am hoping that I can cut down some runtime with your help.

r = RandomReal[NormalDistribution[0., 1./2.], 6];

es = Eigensystem[H0[\[Omega]0, r[[1]], r[[2]], r[[3]], r[[4]], r[[5]], r[[6]] ];

\[Epsilon] = es[[1]];
v = es[[2]];
vS = Conjugate[v];

(*elements of v and vS are called later; v[[1]], vS[[1]] etc...*)

H0 is a compiled function which sped things up a little. It looks like this:

enter image description here.

In copy-paste form, H0[] is

 {{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)}}

Is there anything else that can be optimized here?

share|improve this question
2  
Without you saying more on H0[], there's really nothing else to say, except that the computation of an eigensystem can be simplified a bit, if the matrix given to it has some structure (e.g. symmetry, sparseness). –  J. M. Oct 22 '12 at 2:24
    
Just to check - is whatever loop you have running 10 mil times parallellized? –  VF1 Oct 22 '12 at 2:43
    
@J. M. H0[] is 4x4 Hermitian. Four elements of the matrix are 0. –  BeauGeste Oct 22 '12 at 3:06
3  
"Four elements of the matrix are 0." - which ones? You'll have to be more explicit than that for answerers to get anywhere... –  J. M. Oct 22 '12 at 3:53
2  
@BeauGeste Are you calculating the eigensystem for each choice of random parameters (rather than doing it once symbolically and then substitute specific values) ? –  b.gatessucks Oct 22 '12 at 7:43
show 6 more comments

2 Answers

up vote 8 down vote accepted

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)}}];
share|improve this answer
    
thanks for you analysis. The times look very promising. I need the eigenvectors to be normalized. Is it possible to get normalized vectors out of Eigensystem? –  BeauGeste Oct 23 '12 at 3:59
1  
@Beau, look up Normalize[] and Orthogonalize[]. –  J. M. Oct 23 '12 at 4:05
add comment

If Length[r]==6 the following will be faster.

es = Apply[Eigensystem[H0[\[Omega]0,##]]&,r];

Well I may not have the syntax right because I am doing this with my iPhone. My point is use one call to Apply instead of six calls to Part.

If Length[r]>6, then change ## to (#1,#2,#3,#4,#5,#6).

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.