7
$\begingroup$

I need to generate a very large sparse block matrix, with blocks consisting only of ones along the diagonal. I have tried several ways of doing this, but I seem to always run out of memory.

The fastest way of doing this that I've come up with so far is as follows:

(typically, I will need n to be at least 2500 and m of the order 50).

tmp= SparseArray[{}, {n,n}, 1];

SparseArray@
 ArrayFlatten@
  Table[If[i == j, tmp, 0], {i, m}, {j, m}]

Example when n=2, m=4:

Matrix

The problem with this construction is that ArrayFlatten for some reason converts the result to a normal matrix, and I run out of memory. That is, when it works, this code computes the end result very quickly and the result does not take up much memory. At some specific number however, it suddenly crashes as the intermediate ArrayFlatten step clogs up the memory.

Any help will be greatly appreciated!

$\endgroup$
9
$\begingroup$

Maybe this ?

SparseArray[{Band[{1, 1}, Dimensions[tmp] {m, m}] -> {tmp}}, Dimensions[tmp] {m, m}]
$\endgroup$
  • $\begingroup$ Thanks a lot, this seems to solve the issue! $\endgroup$ – Jesper Jan 22 '13 at 11:13
3
$\begingroup$

I think KroneckerProduct is the right tool for this question:

kp[n_,m_] := KroneckerProduct[
    SparseArray[IdentityMatrix[m]],
    SparseArray[SparseArray[{}, {n,n}, 1], Automatic, 0]
]

Note that I convert the SparseArray in the second argument into a SparseArray that has 0 as the default element. Otherwise, the output of KroneckerProduct is no longer a SparseArray object. Here is Mr.Wizard's last suggestion:

fn[n_,m_] := SparseArray[
    Tuples[Range@# - {1,0,0}] . {Rest@#, {1,0},{0,1}}&@{m,n,n} -> 1
]

And a comparison:

r1 = fn[300, 20]; //MaxMemoryUsed //AbsoluteTiming
r2 = kp[300, 20]; //MaxMemoryUsed //AbsoluteTiming
r1 === r2

{0.959416, 379069184}

{0.023961, 59095360}

True

$\endgroup$
  • $\begingroup$ I bow to the the master. :-) I didn't even realize KroneckerProduct could be used this way, much less how to optimize it as you did. Thank you, again, for lending your expertise to this community. Are you able to provide any additional insight into why this outperforms SparseArray[Band[{1, 1}] -> SparseArray[{}, {m, n, n}, 1]];] to such a great degree? $\endgroup$ – Mr.Wizard Jun 11 '17 at 18:08
  • $\begingroup$ @Mr.Wizard My guess is that "unsparsing" is occurring, i.e., the SparseArray in the Band is getting converted to a normal matrix under the hood. $\endgroup$ – Carl Woll Jun 11 '17 at 18:24
  • 3
    $\begingroup$ For what it's worth I think IdentityMatrix[m, SparseArray] is superior to SparseArray[IdentityMatrix[m]] $\endgroup$ – Mr.Wizard Jun 11 '17 at 18:41
  • $\begingroup$ Agreed, that is nicer! $\endgroup$ – Carl Woll Jun 11 '17 at 18:43
1
$\begingroup$

I find that this performs equivalently to b.gatessucks' code while being significantly cleaner.

n = 3;
m = 2;

SparseArray[Band[{1, 1}] -> SparseArray[{}, {m, n, n}, 1]]

$\left( \begin{array}{cccccc} 1 & 1 & 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 1 \\ 0 & 0 & 0 & 1 & 1 & 1 \\ 0 & 0 & 0 & 1 & 1 & 1 \\ \end{array} \right)$

However ArrayFlatten is still worth your consideration. Using a modified form of a method I posted for this slightly more general Q&A How to form a block-diagonal Matrix from a list of matrices?:

IdentityMatrix[m] /. {1 -> SparseArray[{}, {n, n}, 1]} // ArrayFlatten

(* same output as above, but not a SparseArray *)

In cases were n >> m (as described in the question) the peak memory consumption is actually not much more, it is easily an order of magnitude faster, and conversion into a sparse array is very fast:

{n, m} = {300, 20};

MaxMemoryUsed[
  r1 = SparseArray[Band[{1, 1}] -> SparseArray[{}, {m, n, n}, 1]];] // AbsoluteTiming

MaxMemoryUsed[
  r2 = IdentityMatrix[m] /. {1 -> SparseArray[{}, {n, n}, 1]} // 
     ArrayFlatten;] // AbsoluteTiming

r1 == r2

SparseArray[r2]; // RepeatedTiming
{3.37077, 459439920}

{0.13651, 576105096}

True

{0.091, Null}

For optimal memory performance a variation of ybeltukov's blockArray from Speeding up generation of block diagonal matrix surpasses both methods:

fn[n_, m_] := 
  SparseArray[Tuples[Range@# - {1, 0, 0}].{Rest@#, {1, 0}, {0, 1}} &@{m, n, n} -> 1]

MaxMemoryUsed[r3 = fn[300, 20];] // RepeatedTiming

r1 == r3
{0.769, 379069184}

True
$\endgroup$

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.