3
$\begingroup$

I have to create the following sparse arrays, which are reshaped in a two dimensional matrix, so that i can use the build in Dot[] function to do a kind of multiplication.

c[theta_Real, xMax_Integer, yMax_Integer] := 
 ArrayReshape[SparseArray[{
    {x_, y_, s_, x_, y_, s_} -> Cos[theta/2.],
    {x_, y_, 1, x_, y_, 2} -> -Sin[theta/2.],
    {x_, y_, 2, x_, y_, 1} -> Sin[theta/2.]
    }, {2 xMax + 1, 2 yMax + 1, 2, 2 xMax + 1, 2 yMax + 1, 2}, 
   0.], {2 (2 xMax + 1) (2 yMax + 1), 2 (2 xMax + 1) (2 yMax + 1)}]

sX[xMax_Integer, yMax_Integer] := ArrayReshape[SparseArray[{
    {x1_, y_, 1, x2_, y_, 1} /; x1 == x2 + 1 -> 1.,
    {x1_, y_, 2, x2_, y_, 2} /; x1 == x2 - 1 -> 1.
    }, {2 xMax + 1, 2 yMax + 1, 2, 2 xMax + 1, 2 yMax + 1, 2}, 
   0.], {2 (2 xMax + 1) (2 yMax + 1), 2 (2 xMax + 1) (2 yMax + 1)}]

sY[xMax_Integer, yMax_Integer] := ArrayReshape[SparseArray[{
    {x_, y1_, 1, x_, y2_, 1} /; y1 == y2 + 1 -> 1.,
    {x_, y1_, 2, x_, y2_, 2} /; y1 == y2 - 1 -> 1.
    }, {2 xMax + 1, 2 yMax + 1, 2, 2 xMax + 1, 2 yMax + 1, 2}, 
   0.], {2 (2 xMax + 1) (2 yMax + 1), 2 (2 xMax + 1) (2 yMax + 1)}]

I'm trying to avoid direct definition of the 2D sparse matrices, because I have to handle a lot of similiar arrays with given multidimensional representation and I realized that reshaping is very fast, compared to the rest of calculation. Now I have to calculate

rF2D[a_, b_, x0_, y0_, theta1_, theta2_, 
  xMax_, yMax_, steps_] := 
 Block[{r = r0[a, b, x0, y0, xMax, yMax], 
   walk = sY[xMax, yMax].c[theta2, xMax, yMax].sX[xMax, 
      yMax].c[theta1, xMax, yMax]},
  Do[r = walk.r.ConjugateTranspose[walk], {steps}]; r]

for large xMax,yMAx, e.g. rF2D[1.,0.,0,0,Pi/2.,Pi/4.,50,50,50] where r0 is a matrix with complex entries:

r0[a_, b_, x0_, y0_, xMax_, yMax_] := 
 ArrayReshape[SparseArray[{
    {x0 + xMax + 1, y0 + yMax + 1, 1, x0 + xMax + 1, y0 + yMax + 1, 
      1} -> Abs[N[a]]^2,
    {x0 + xMax + 1, y0 + yMax + 1, 1, x0 + xMax + 1, y0 + yMax + 1, 
      2} -> N[a]\[Conjugate] N[b],
    {x0 + xMax + 1, y0 + yMax + 1, 2, x0 + xMax + 1, y0 + yMax + 1, 
      1} -> N[a] N[b]\[Conjugate],
    {x0 + xMax + 1, y0 + yMax + 1, 2, x0 + xMax + 1, y0 + yMax + 1, 
      2} -> Abs[N[b]]^2
    }, {2 xMax + 1, 2 yMax + 1, 2, 2 xMax + 1, 2 yMax + 1, 2}, 
   0.], {2 (2 xMax + 1) (2 yMax + 1), 2 (2 xMax + 1) (2 yMax + 1)}] 

At the moment this is impossible, because calculation needs too much time. For much smaller xMax,yMAx, e.g. rF2D[1.,0.,0,0,Pi/2.,Pi/4.,10,10,10] finishes after 5 seconds on may machine.

Is there any possibility to speed up the creation of the sparse arrays, because this seems to eat most of the time and the multiplication is doing quite fast? Thanks for your help!

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – user9660 Mar 27 '15 at 18:12
4
$\begingroup$

Updated with cleaner and easier to adapt method:

The checking of the patterns gets expensive, more so when conditions are attached. Looking at your generators, it's clear that the fulfillment of the patterns and/or conditions is quite sparse.

Better to generate directly the hits, and create the array from those.

Using your 'c' generator with just the creation part for brevity (that's where the time is getting burned):

c[theta_Real, xMax_Integer, yMax_Integer] := 
  SparseArray[{{x_, y_, s_, x_, y_, s_} -> 
     Cos[theta/2.], {x_, y_, 1, x_, y_, 2} -> -Sin[theta/2.], {x_, y_,
       2, x_, y_, 1} -> Sin[theta/2.]}, {2 xMax + 1, 2 yMax + 1, 2, 
    2 xMax + 1, 2 yMax + 1, 2}, 0.];

newc[theta_Real, xMax_Integer, yMax_Integer] :=
 Module[{k1 = Cos[theta/2.], k2 = -Sin[theta/2.], k3 = Sin[theta/2.],
   d1 = 2 xMax + 1, d2 = 2 yMax + 1, d3 = 2, t1, t2, t3},
  t1 = Thread[Rule[Flatten[Table[{x, y, s, x, y, s}, {x, d1}, {y, d2}, {s, d3}],2], k1]];
  t2 = Thread[Rule[Flatten[Table[{x, y, 1, x, y, 2}, {x, d1}, {y, d2}], 1], k2]];
  t3 = Thread[Rule[Flatten[Table[{x, y, 2, x, y, 1}, {x, d1}, {y, d2}], 1], k3]];
  SparseArray[Join[t2, t3,t1], {d1, d2, d3, d1, d2, d3}, 0.]]

Timings:

tOld = First@Timing[r1 = c[2.5, 10, 5];]
tNew = First@Timing[r2 = newc[2.5, 10, 5];]
r1 == r2
tOld/tNew

(*
1.950012
0.015600
True
125.00
*)

And from n=1 to 20 for c[2.5,n,n] and newc[2.5,n,n]:

enter image description here

About 900X faster by 20x20 case, growing with size...

Simply adapt the rule generation to your other cases. You can hit the edit to see my earlier, much clumsier approach.

$\endgroup$
  • $\begingroup$ Oh, that's quite nice!! Thanks. Seems to be that SparseArray is eating to much time evaluating the given pattern for filling the array. I will check this tomorrow on my main machine. Sorry for answering this way. I've got no stackexchange account, so marking your answer as useful as well as commenting seems to be not possible anymore... $\endgroup$ – user27376 Mar 29 '15 at 12:39
  • $\begingroup$ Please register an account. It only takes a few minutes. I am going to convert this "answer" to a comment now. $\endgroup$ – Mr.Wizard Mar 29 '15 at 12:48

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.