8
$\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$
0

3 Answers 3

10
$\begingroup$

Maybe this ?

SparseArray[{Band[{1, 1}, Dimensions[tmp] {m, m}] -> {tmp}}, Dimensions[tmp] {m, m}]
$\endgroup$
1
  • $\begingroup$ Thanks a lot, this seems to solve the issue! $\endgroup$
    – Jesper
    Jan 22, 2013 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$
4
  • $\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, 2017 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, 2017 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, 2017 at 18:41
  • $\begingroup$ Agreed, that is nicer! $\endgroup$
    – Carl Woll
    Jun 11, 2017 at 18:43
2
$\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 and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.