# Efficiently populate a Sparsearray for a set of rules for a constrained basis

I want to populate quite a large SparseArray(10^6 x10^6) efficiently. It is basically a spin system Hamiltonian with a constrained Hilbert space. Unlike the examples I have looked at in this forum this introduces a problem. On use of BitXor to emulate terms like $$\sigma_x$$ , in many cases it would yield a value which is not in the Hilbert space, which means I have to do a search in the basis to find the position it corresponds to before storing it in the sparsearray. Let me first show the code which I have written.

First to generate the constrained Hilbert space. I am not going to explain the constraint in detail because that may be confusing, but the following code generates it efficiently enough for me so one can assume this is the Hilbert space.

N1 =24;
X1 = Select[Range[0, 2^N1 - 1], BitAnd[#1, BitShiftLeft[#1, 1]] == 0 & ];
X1 = Select[X1, BitGet[#1, 0]*BitGet[#1, N1 - 1] == 0 & ];
len=Length[X1];


X1 stores the allowed basis vectors and len stores the length of X1. Now the Hamiltonian term which I want to calculate is $$\sigma_x$$ or in terms of creation and annihilation operators is $$\tau \sum_l^{N_1} (d_l^{\dagger}+d_l)$$. To generate this I wrote the following code,

SetAttributes[f, Listable];
Needs["Combinatorica"];
g[j_,k_] := BinarySearch[X1, BitXor[X1[[j]], 2^k]];
f[i_] := (s1 = Mod[i, len] + 1; s2 = IntegerPart[i/len]; k1 = g[s1, s2];
If[IntegerQ[k1], {s1, k1} -> \[Tau]]);
Print[AbsoluteTiming[Q1 = SparseArray[DeleteCases[f[Range[0, len*N1 - 1]],
Null]]; ]];


This takes approximately 217 seconds in my i7 machine. However a similar code written in normal loop constructs in FORTRAN takes a second to evaluate. Now my experience tells me while Mathematica will take more than a second of course because it is not a low level programming language, but the time taken being almost two orders of magnitude higher indicates there is a faster way of doing this. Can anyone help me find it?

I actually need to go to sizes an order of magnitude higher than in the example, hence I need a speedup, plus while I can generate this data via FORTRAN and import to MATHEMATICA from a file, I want to do that as a last resort after exhausting all possibilities of speedup here.

Edit:- In the slow code snippet I want to do the following, I take a basis vector, then operate a BitXor gate one by one on all the bit positions, i.e. from 0 to N1-1. Each would give me a new integer, which may or may not be present in the Hilbert space. I check whether it is present and also its position in the basis vector array by using BinarySearch, once I do that I create a list like {x1,x2}->\tau where x1 is position the basis vector which I started with and x2 is position of the basis vector I get after the operation. If it is not present a "Null" is returned which I delete from the SparseArray. This operation is tried to be made listable using the function "f". Where using "Range" I try to generate both the basis vector positions as well as the bit position using one value to vectorized the code. The "Mod" and "IntegerPart" operation then splits it up to basis vector position and bit position inside the function f.

• Perhaps you could describe what operations the slow code is trying to carry out, without any reference to the underlying physical problem? I'm not sure that I follow what you are doing so it may be difficult to help. – MarcoB Mar 2 at 14:19
• Sure, I shall add it in the main question. – Roopayan Ghosh Mar 2 at 14:22
• @MarcoB I have tried to explain it, is it still confusing? – Roopayan Ghosh Mar 2 at 14:29

There are a couple of things that can be improved:

Binary search can simply be replaced by a lookup table that is very cheap to build.

Moreover, IntegerPart[i/len] leads to rational numbers (so no machine numbers), while Quotient can to the thing in vectorized way.

Avoiding the use of a list of rules allows one to work with packed arrays.

Using Pick instead of Select.

Exploiting that bitwise operations are vectorized.

Instead of using DeleteCases (which involves slow pattern mathing), we rephrase thr problem so that the undocumented function RandomPrivatePositionsOf can be used.

This is my version of the code:

N1 = 24;
X1 = Range[0, 2^N1 - 1];
X1 = Pick[X1, BitAnd[X1, BitShiftLeft[X1, 1]], 0];
X1 = Pick[X1, BitGet[X1, 0] BitGet[X1, N1 - 1], 0];
len = Length[X1];

i = Range[0, len N1 - 1];
s2 = Quotient[i, len];
s1 = Subtract[i, len s2] + 1;
labels = X1 + 1;
keys = BitXor[X1[[s1]], Power[2, s2]] + 1;
lookuptable = Normal[SparseArray[ Partition[labels, 1] -> Range[len], {2^N1 - 1}]];
k1 = lookuptable[[keys]];
pos = RandomPrivatePositionsOf[Unitize[k1], 1];

Q1 = SparseArray[Transpose[{s1[[pos]], k1[[pos]]}] -> \[Tau]];


It takes about 0.55 second on my machine while the original code took about 200 seconds.

This could be faster if QuotientRemainder were properly vectorized.

Not applying Normal and using

 lookuptable = SparseArray[Partition[labels, 1] -> Range[len], {2^N1 - 1}];


instead saves a lot of memory, but it takes also ten times longer, because SparseArray has to look up keys in a sorted list of labels (namely in the list lookuptable["ColumnIndices"]). And even though it can be done by a binary search, it is still substantially faster than indexing into a dense array.

Edit:

Because labels is sorted, you can also use binary search:

k1 = cBinaryLookUp[labels, keys];


with the compiled function cBinaryLookUp  below. This is even slower than indexing into the SparseArray, but it should use less memory.

cBinaryLookUp = Compile[{{labels, _Integer, 1}, {key, _Integer}},
Block[{i, j, k, a, b, c},
j = Length[labels];
b = CompileGetElement[labels, j];
If[key > b,
0,
(
i = 1;
a = CompileGetElement[labels, i];

k = i + Quotient[j - i, 2];
c = CompileGetElement[labels, k];
While[j - i > 1,
If[key <= c,
j = k;
b = c;
,
i = k;
a = c;
];
k = i + Quotient[j - i, 2];
c = CompileGetElement[labels, k];
];
If[c < key,
If[b == key, j, 0],
If[c == key, k, 0]
]
)
]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];


Please check this function for correctness. I can give no guarantees for it.

Edit 2:

And this is a compiled function that should generate X1, not using more memory that the final X1 requires. It is a bit slow, though. But it can certainly be parallelized.

getX1 = Compile[{{N1, _Integer}},
Block[{bag},
bag = InternalBag[Most[{0}]];
Do[
If[BitAnd[i, 2 i] == 0,
InternalStuffBag[bag, i];
]
, {i, 0, 2^(N1 - 1)}];
Do[
If[BitAnd[i, 2 i] == 0,
InternalStuffBag[bag, i];
]
, {i, 2^(N1 - 1) + 2, 2^N1 - 1, 2}];
InternalBagPart[bag, All]
],
CompilationTarget -> "C"
];


Have fun finding out how it works! =)

It still requires about 2^N1` work, so it is still far from optimal.

• This is amazing. Thank you so much. Just wanted to add I have been following your answers on linear algebra for a long time and it has helped me a lot. There is so much to learn in Matematica. I have accepted this answer as the correct one. – Roopayan Ghosh Mar 3 at 5:18
• Thank you. You're welcome. – Henrik Schumacher Mar 3 at 7:38
• Sorry its been 2 months since this reply but I had a query, why did you use Normal[SparseArray[]] in lookuptable, I mean if I am converting it to normal form, whats the use of generating a SparseArray? – Roopayan Ghosh May 23 at 8:16
• Okay thanks, I will implement that for memory intensive problems then. – Roopayan Ghosh May 25 at 8:32
• I implemented it and the speeds are not really that far off, maybe just 2x slower. My current code has an even more constrained Hilbert space, so the previous method had too much of a bottleneck. The memories now are 20x less almost. Amazing. Also yes that is what is needed in X1. – Roopayan Ghosh May 25 at 12:59