# Efficiently populate sparse matrix with band structure

I'm trying to efficiently populate elements of a very large (2^20 x 2^20) symmetric matrix with 1s - luckily the matrix is very sparse, <0.1% filling. Further, the matrix has a very well defined periodic banded structure, as shown here:

In reality, this matrix is the result of a series of KroneckerProducts of 2x2 matrices (or in other words, a Kronecker Sum of Pauli X matrices), which is what gives it that characteristic banded structure - I'm hoping to find a way to speed up the generation without using kronecker products, because even with sparse matrices, the computation can take several seconds or minutes depending on the final dimensionality.

My first question relates to creating this sparse matrix efficiently. I've toyed with lots of different ways of generating even simple bands for the sparse array. For simply populating on the diagonal, the quickest method clearly seems to be to use the {i_,i_} notation, as shown here:

dim = 15;

aa = SparseArray[{i_, i_} -> 1, {2^dim, 2^dim}] // RepeatedTiming;
bb = SparseArray[Band[{1, 1}] -> 1, {2^dim, 2^dim}] // RepeatedTiming;
cc = SparseArray[Table[{ii, ii} -> 1, {ii, 2^dim}], {2^dim, 2^dim}] //RepeatedTiming;
dd = SparseArray[Normal[AssociationThread[Table[{ii, ii}, {ii, 2^dim}] -> Table[1, {ii, 2^dim}]]], {2^dim,2^dim}] // RepeatedTiming;

Column[{aa[], bb[], cc[], dd[]}]

aa[] == bb[] == cc[] == dd[]

0.000309
0.031
0.081
0.054

True


However, when we try to do off-diagonal entries, this gets much worse, presumably because the condition has to be continually checked:

dim = 15;

aa = SparseArray[{i_, j_} /; j - i == 1 -> 1., {2^dim, 2^dim}] // RepeatedTiming;
bb = SparseArray[Band[{1, 2}] -> 1, {2^dim, 2^dim}] // RepeatedTiming;
cc = SparseArray[Table[{ii, ii + 1} -> 1, {ii, 2^dim - 1}], {2^dim, 2^dim}] // RepeatedTiming;
dd = SparseArray[Normal[AssociationThread[Table[{ii, ii + 1}, {ii, 2^dim - 1}] -> Table[1, {ii, 2^dim - 1}]]], {2^dim, 2^dim}] // RepeatedTiming;

Column[{aa[], bb[], cc[], dd[]}]

aa[] == bb[] == cc[] == dd[]

0.185
0.031
0.095
0.052

True


From those two examples then it seems like Band is our best choice, but Band is still painfully slow, especially when compared to the {i_,i_} for the diagonal. Further, this is more frustrating, because in MATLAB the same problem can be accomplished an order of magnitude faster (this took ~1.4 ms): But the fact that the original {i_,i_} case for the diagonal was so fast suggests that there is a more efficient way to do this.

So then my first question is: given all of that, is there a more efficient way to populate the bands of our sparse matrix, so that the speed can at least rival the equivalent in MATLAB?

And my second question, a bit predicated on the first: with whatever method is the most efficient, what is the best way to generate the periodic banding structure present in the final matrix (see above). You can accomplish it with Band by manually inserting spaces with 0s, but doing so can't be the most efficient way.

Finally, because of that period-2 banded structure of the final matrix, where each quadrant is a recursive block of ever smaller diagonal matrices with side length smaller by a factor of 2, maybe you could generate all the smaller diagonal blocks, and then just place them in the final matrix - I'm not sure how this would be accomplished however. Of course, remember that the matrix is symmetric, so I would think that would help with efficient generation because really just one triangle has to be generated and then flipped.

Addendum: MATLAB code for generating the plot, as requested in the comments. This takes on the order of milliseconds for my machine, even with N=15.

N = 4;
a_id    = (0:2^N-1)';

dimH        = length(a_id);
AA          = sparse(dimH, dimH);

for i = 1:N
[~,k1,k2] = intersect(a_id, bitxor(a_id,2^(i-1)));
AA        = AA + sparse(k1,k2,ones(length(k1),1)/2,dimH,dimH);
end


Addendum 2: Henrik's answer is very good, and gives what I'm looking for. Still, it's a bit disappointing that the solution is almost an order of magnitude slower than the equivalent in MATLAB, but I'll take it! As a further question, I took a stab at the method briefly mentioned above of manually placing subarrays within the master array. This takes advantage of the very quick generation of diagonal sparse matrices as I showed above. My current implementation is not very efficient, but I'm wondering if such a method has any possibility for efficiency, and if so how? This is more a curiosity than anything as Henrik's answer is already fast enough for my usecase. For n=14 this takes 3 seconds for me.

func[n_] := Module[{
subarrays =
Table[SparseArray[{i_, i_} -> 0.5, {2^dim, 2^dim}], {dim, 0,
n - 1}],
master = SparseArray[{}, {2^n, 2^n}]},
Do[master[[(jj - 1) 2^(ii + 1) + 1 ;; 2^ii (2 jj - 1),
1 - 2^ii + 2^(1 + ii) jj ;; jj 2^(ii + 1)]] =
subarrays[[ii + 1]], {ii, 0, n - 1}, {jj, 1, 2^(n - 1 - ii)}];
master + Transpose[master]
]


Addendum 3: In response to the comment, this is indeed for the purpose of spins on a lattice, and is simply the Kronecker sum of Pauli X matrices. The equivalent of this generation using KroneckerProduct takes 400ms for N=15 (though it's certainly possible my implementation is not the best).

• A related question. Would you happen to have the MATLAB code for generating the matrix corresponding to your first picture? May 30 '20 at 7:01
• I searched for what felt like hours and never ran across that answer... As for the MATLAB code, I've added it as an addendum. May 30 '20 at 7:40
• Would you please add also a couple of exampes of the resulting matrices for small N for testing reasons? May 30 '20 at 8:14
• Could you please provide an example of the final matrices you want to get (spin operators on a lattice? Fermionic operators?) and of the time you need to generate them via KroneckerProduct[]? May 30 '20 at 12:35

For a start, I think this should provide the correct result. Do you confirm that?

cf = Compile[{{a, _Integer, 1}, {i, _Integer}},
Transpose[{a + 1, 1 + BitXor[a, 2^(i - 1)]}],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
];

ClearAll[A];
A[n_] := With[{a = Range[0, 2^n - 1]},
SparseArray[Flatten[cf[a, Range[1, n]], 1] -> 0.5, {2^n, 2^n}, 0.]
]


It takes 13 ms and 1600 ms on my machine to evaluate A and A, respectively. Is that fast enough? If not then I have also this slightly faster version that happens to produce also CRS-conform matrices (i.e., SparseArraySparseArraySortedQ evaluates to True on the output):

cg = Compile[{{i, _Integer}, {n, _Integer}},
Sort[1 + BitXor[i, 2^Range[0, n - 1]]],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
];
B[n_] := SparseArray @@ {Automatic, {2^n, 2^n}, 0., {1, {
Range[0, n 2 ^n, n],
Partition[Flatten[cg[Range[0, 2^n - 1], n]], 1]
},
ConstantArray[0.5, n 2 ^n]
}
};


Evaluating B and B on my machine takes 8 ms and 620 ms, respectively.

• @HenrikSchumacher, just a curiosity: how much do you gain by compiling? if you do not compile (but you parallelize) how much the performance is worse? May 30 '20 at 14:48
• Well, compilation on its own does not really make a difference, because I call BitXor in a listable way. I really use it for its parallelism: cg[Range[0, 2^n - 1], n] is only about 4 times faster than Table[Sort[1 + BitXor[i, 2^Range[0, n - 1]]], {i, 0, 2^n - 1}] on my Quad Core. But the task is also too small to get any benefit from ParallelTable. Since the body of the CompiledFunction consists merely of calls to other libraries, there is also little difference between CompilationTarget -> "C" and CompilationTarget -> "WVM". May 30 '20 at 14:56
• Thank you! My question was motivated because rarely I can get great benefits from compilation in my codes. And I was wondering whether I am losing some important points. Could you explain me why the parallelism obtained by compilation is better than the ParallelTable? May 30 '20 at 14:59
• That's a good question. I do not why this is, but from my experience I can tell that ParallelTable has a tremendous overhead that pays off only if the tasks in its body are big enough and if only very few data has to be exchanged. It works great for scheduling tasks that run over a couple of minutes. IIRC correctly, Compile` uses a lightweight OpenMP loop to do its thing: It just chops the list of tasks into as many chunks as there are threads available, tells the threads where to read from and where to write to and then lets them loose. Scales pretty well already for smaller problems. May 30 '20 at 15:25
• Nah, it is undocumented. It creates a sparse array directly from the row pointers and column indices and maybe also skips some consistency checks. You can find more details in this post (also follow the links posted by user21 ). Jun 3 '20 at 7:41