# Individual block diagonalisations take more time than single diagonalization of block matrix

Sample Problem of a large block matrix:

I want to find the lowest eigenvalue of a blockdiagonal matrix h[x] for a range of values of x

• To save the time, I do so by finding the lowest eigenvalue of each block individually and then lowest out of them but

Individual block diagonalizations take more time than single diagonalization of block matrix h[x] for a range of values of x.

b0 = RandomReal[{-1, 1}, {16, 16}]; c0 = RandomReal[{-1, 1}, {64, 64}]; d0 = RandomReal[{-1, 1}, {128, 128}]; e0 = RandomReal[{-1, 1}, {32, 32}];
b = b0 + Transpose[b0];
c = c0 + Transpose[c0];
d = d0 + Transpose[d0];
e = e0 + Transpose[e0];
h[x_] := x*SparseArray[Band[{1, 1}] -> {b, c, d, e}];
AbsoluteTiming[ParallelTable[-Eigenvalues[-h[x], 1, Method -> {"Arnoldi", "Criteria" -> "RealPart"}][], {x, 0.1,5., 0.1}];][]
eigbloc0[x_, bl_] := -Eigenvalues[-h[x][[#, #]], 1, Method -> {"Arnoldi", "Criteria" -> "RealPart"}][] & /@ {blocks[x][[bl]]}
eigbloc1[x_] := Table[eigbloc0[x, bl], {bl, 1, 4, 1}]
AbsoluteTiming[ParallelTable[Min[Flatten[{eigbloc1[x]}]], {x, 0.1, 5., 0.1}];][]


Output:

7.98266

11.7006

I was expecting lower time for individual block diagonalization, Could anyone suggest me the resolution ?

• One issue is that you rebuild the full matrix each time you want to compute only a block, because eigbloc0 calls h. Ideally, you know the blocks already independently of the construction of the matrix -- because you know how to construct the matrix. Sep 11 at 7:59
• @HenrikSchumacher Yes, For a single computation, say h[x=1.0], Individual digonalization procedure is faster. Could you please tell me how to avoid the repetitive computations in above program ? Sep 11 at 8:57
• Try this: ClearAll[blocks]; blocks = Range[#1 + 1, #2] & @@@ Partition[Accumulate[{0, 16, 64, 128, 32}], 2, 1]; eig[x_] := Module[{A}, A = -h[x]; Min@Table[ -Eigenvalues[A[[block, block]], 1, Method -> {"Arnoldi", "Criteria" -> "RealPart"}], {block, blocks} ] ]; Sep 11 at 9:09
• @HenrikSchumacher, You are wonderful, Worked! Sep 11 at 9:20

First issue is that you compute the matrix each time you try to compute within a block. Second issue is that you recompute all the blocks themselves.

So this was just not a fair comparison. Have a look at this example. Even without parallelization, computing the eigenvalues of the dense diagonal blocks is two orders of magnitude faster:

n = 100;
diagonal = # + #\[Transpose] &[RandomReal[{-1, 1}, {#, #}]] & /@
RandomInteger[{120, 240}, {n}];
A0 = SparseArray[Band[{1, 1}] -> diagonal];

val1 = -Max[Eigenvalues[-A0, 1,
Method -> {"Arnoldi", "Criteria" -> "RealPart"}]]; //
AbsoluteTiming // First
val2 = -Max[Eigenvalues[-#, 1,
Method -> {"Arnoldi", "Criteria" -> "RealPart"}] & /@
diagonal]; // AbsoluteTiming // First
val1 == val2


5.35384

0.063925

True

I guess the main speedup here comes from the fact that replace a the sparse matrix vector multiplication kernel by many small, highly efficient dense matrix-vector multiplication kernels. (Sparse matrix multiplication typically has quite a "random" memory access pattern, so the compute kernels do not optimize for contiguous memory access (because they cannot).)