speeding up creation of SparseArray

Question

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!

Answer 1

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]:

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.