Specify positions in a SparseArray using patterns with functions

Question

I spent some time trying to figure out what was what was making so slow the calculation of the following array(actually this is only a SUPERsimple toy model):

Input to the array:

Nx=2;
Ny=1;
nn =Nx*Ny; 
Nband = 2;
Nstates = 2*nn*Nband;

eigenvectors = Table[Range[Nstates] + i, {i, 1, Nstates}];

InverseFlatten[l_, dimensions_] := Fold[Partition, Flatten@l,Most[Reverse[dimensions]]];
(*************)
N1[i_] := Module[{ix, iy, n1x},
ix = Mod[i, Nx, 1];
iy = Quotient[i, Nx, 1] + 1;
If[ix + 1 > Nx, n1x = 1, n1x = ix + 1];
{n1x + (iy - 1)*Nx}];
(*************)
uv = InverseFlatten[eigenvectors, {Nstates, Nband, 2, nn}];
u = uv[[1 ;; Nstates, 1 ;; Nband, 1]];
v = uv[[1 ;; Nstates, 1 ;; Nband, 2]];
f = Range[Nstates];
V = Table[Which[MemberQ[N1[i], j] == True, 2], 
{l, 1, Nband}, {m, 1, Nband}, {s, 1, Nband},{q, 1, Nband}, {i, 1, nn}, {j, 1, nn}];

Array:

delta = 
Table[
Total[Flatten[V[[μ, ν, ;; , ;; , ;; , ;;]]*
  Table[
     Table[
        Which[
            MemberQ[N1[i], j] == True,
            f.(u[[;; , q, i]]*Conjugate[v[[;; , s, j]]])    ], 
          {i, 1, nn}, {j, 1, nn}],
      {q, 1, Nband}, {s, 1, Nband}], 1]], 
 {μ, 1, Nband}, {ν, 1, Nband}];

After having tried several things like using Packed arrays, using Compile...and seen that nothing really helped, I thing that what it actually slows down the thing is Which. I think, that the process should be a lot faster if instead of using Which, which gives a condition for which matrix elements to calculate, I directly specify the positions for which I would like to calculate the matrix element; I think so because the number of elements to be calculated is really small compared to the number of elements(for big nn).

Then looking at the SparseArray list manipulating tutorial I saw that one can do this using patterns, for example:

Normal[SparseArray[{{i_, i_} :> 1}, {nn, nn}]]

Now, what specifies the second index of my matrix is the function N1 defined above, therefore, I would like to use something like:

 A=Normal[SparseArray[{{i_, N1[i_][[1]]} :> 1}, {nn, nn}]]  
 SparseArray::posd: The left-hand side of {i_,n1x$671+2 Quotient[i_,2,1]}:>1  in {{i_,n1x$671+2 Quotient[i_,2,1]}:>1} is not a position or a pattern that will match the position of an element in an array with dimensions {2,2}. >>

but it complaints. I expect to get:

 A={{0,1},{1,0}}

since N1[1][[1]]=2 and N1[2][[1]]=1.

Does anybody now how could I make a "pattern SparseArray" having this function N1 in the pattern?

Thanks

Try Condition: A = Normal[SparseArray[{{i_, j_} /; j == N1[i][[1]]} :> 1, {nn, nn}]] — Michael E2
– Michael E2, Commented Dec 16, 2013 at 15:17

Michael E2 · Accepted Answer · 2013-12-16 15:50:14Z

If you want a fast way, then you don't really want to pursue SparseArrays with patterns and condition. Built-in patterns like Band are optimized and work relatively fast.

Slow ways

With the parameters,

Nx = 10;
Ny = 20;
nn = Nx*Ny;
Nband = 5;
(* etc. *)

the following took nearly 90 sec., while I went to get some coffee.

SparseArray[{l_, m_, s_, q_, i_, j_} /; 
      Mod[i, Nx] + Quotient[i, Nx, 1]*Nx + 1 == j -> 2,
   {Nband, Nband, Nband, Nband, nn, nn}]; // AbsoluteTiming

(* {87.753965, Null} *)

The OP's Table method took too long for me to wait. A compiled version took almost 6 seconds:

Vc = With[{Nband = Nband, nn = nn, Nx = Nx, Ny = Ny},
     Compile[{},
      Table[
       If[Mod[i, Nx] + Quotient[i, Nx, 1]*Nx + 1 == j, 2, 0],
       {l, 1, Nband}, {m, 1, Nband}, {s, 1, Nband}, {q, 1, Nband}, {i, 1, nn}, {j, 1, nn}]]
      ][]; // AbsoluteTiming

(* {5.772137, Null} *)

Much faster

This SparseArray approach was several hundred times faster:

spV = SparseArray@ConstantArray[
     SparseArray[
      Table[Band[{i, Mod[i + 1, Nx, 1]}, Automatic, {Nx, Nx}] -> 2, {i, Nx}],
      {nn, nn}],
     {Nband, Nband, Nband, Nband}]; // AbsoluteTiming

(* {0.006551, Null} *)

Check:

 Vc == spV
 (* True *)

P.S. I replaced Null with 0. That seemed to be acceptable in one of the OP's earlier posts.

thanks so much for answering(always there lately :)), I will see if I can understand what you did and apply it to calculate delta in a fast way, which is my problem because, although the computation of V is also slow, V is not in an iterating convergence process, whereas delta is, so I really need delta to be fast!!! — Mencia
– Mencia, Commented Dec 16, 2013 at 15:39

Stack Exchange Network

Specify positions in a SparseArray using patterns with functions

1 Answer 1

Hot Network Questions

Specify positions in a SparseArray using patterns with functions

1 Answer 1

Related

Hot Network Questions