23
$\begingroup$

I'm struggling with the following problem. I have $48$ square matrices (full, filled with real machine precision numbers, thus are packed, all different) of size $128$. I wolud like to place them on a diagonal of sparse array of dimension $48\times 128=6144$. The method(1)

SparseArray @ ArrayFlatten @ ReleaseHold @ DiagonalMatrix[Hold /@ matrices]
(* matrices is a list of 48 matices 128 x 128,
   e.g. matrices = RandomReal[{}, {128, 128}] & /@ Range[48] *)

is too slow (it takes ~6s on my laptop). I'm suspecting that the problem is with the ArrayFlatten function, since this produces huge matrix $6144\times 6144$ filled moslty with zeros (in some sense it unpacks sparse array). Is there any way to do the same but much faster (more efficient)? In a fraction of a second (I'm optimistic)? I've looked at "SparseArray`" context but without much success (SparseArray`VectorToDiagonalSparseArray seems to be equivalent to DiagonalMatrix so accepts only vectiors/lists). (Specific numbers given here are just for tests, in the end I would like to increase size of my problem, but then it of course gets even worse.)

After posting this question I've found the code on MathWorld which gives me the result in ~3.63s. Code by ybeltukov SparseArray[Band@{1, 1} -> matrices] is even faster (~2.48s) but still far from being ideal.

Update: I've checked that asymptotically execution time scales as (based on AbsoluteTiming):

  • $m^{2}n^{2}$ for BlockDiagonalMatrix
  • $m^{2}n^{1}$ for recent version of blockArray by ybeltukov

where: $n$ is a number of matrices/blocks and $m$ is a size of a single matix/block.

$\endgroup$
26
$\begingroup$

You are right, it can be done in a fraction of second. One can explicitly construct an array of indexes

blockArray[mat_] := SparseArray[
  Tuples[Range@# - {1, 0, 0}].{Rest@#, {1, 0}, {0, 1}} &@Dimensions@mat -> 
   Flatten@mat]

Timings:

matrices = RandomReal[1, {48, 128, 128}];

s1 = 
   SparseArray@
    ArrayFlatten@ReleaseHold@DiagonalMatrix[Hold /@ matrices]; // RepeatedTiming
(* {7.56, Null} *)

s2 = SparseArray[Band@{1, 1} -> matrices]; // RepeatedTiming
(* {4.03, Null} *)

s3 = blockArray[matrices]; // RepeatedTiming
(* {0.097, Null} *)

TrueQ[s1 == s2 == s3]
(* True *)

For further acceleration you can create the internal structure of the SparseArray directly

c = Compile[{{b, _Integer}, {h, _Integer}, {w, _Integer}},
   Partition[Flatten@Table[k + i w, {i, 0, b - 1}, {j, h}, {k, w}], 1], 
   CompilationTarget -> "C", RuntimeOptions -> "Speed"];

blockArray2[mat_] := 
 SparseArray @@ {Automatic, # {##2}, 
     0, {1, {Range[0, 1 ##, #3], c@##}, Flatten@mat}} & @@ Dimensions@mat

s4 = blockArray2[matrices]; // RepeatedTiming
(* {0.031, Null} *)

s3 == s4
(* True *)
$\endgroup$
  • $\begingroup$ Many thanks! This is great. I was going to ask whetere it can be further improved, but since you are still improving timings... $\endgroup$ – mmal Nov 14 '15 at 19:11
  • $\begingroup$ @mmal Yes, I know one more hack :) $\endgroup$ – ybeltukov Nov 14 '15 at 19:15
  • $\begingroup$ I would be grateful. $\endgroup$ – mmal Nov 14 '15 at 19:25
  • $\begingroup$ @mmal Please check the update. Note, that you can further accelerate SparseArray construction with LibraryLink and high-efficient C code. $\endgroup$ – ybeltukov Nov 14 '15 at 20:50
  • $\begingroup$ Compiled generation of sparse array structure is very promising. Unfortunately I'm not aware of LibraryLink technology. I'm still wondering whether any undocumented function from SparseArray could do the trick. $\endgroup$ – mmal Nov 15 '15 at 10:32
10
$\begingroup$

There is an undocumented built-in solution:

rules = {#, #} -> RandomReal[{}, {128, 128}] & /@ Range[48];
SparseArray`SparseBlockMatrix[rules]; // RepeatedTiming

(* {0.042, Null} *)
$\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.