5
$\begingroup$

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"}][[1]], {x, 0.1,5., 0.1}];][[1]]
blocks[x_] :=ConnectedComponents@AdjacencyGraph@Unitize@h[x]
eigbloc0[x_, bl_] := -Eigenvalues[-h[x][[#, #]], 1, Method -> {"Arnoldi", "Criteria" -> "RealPart"}][[1]] & /@ {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}];][[1]]

Output:

7.98266

11.7006

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

Improve this question
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    Sep 11, 2021 at 7:59
  • $\begingroup$ @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 ? $\endgroup$
    – 2.1
    Sep 11, 2021 at 8:57
  • 1
    $\begingroup$ 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} ] ]; $\endgroup$
    – Henrik Schumacher
    Sep 11, 2021 at 9:09
  • $\begingroup$ @HenrikSchumacher, You are wonderful, Worked! $\endgroup$
    – 2.1
    Sep 11, 2021 at 9:20

1 Answer 1

Reset to default
6
$\begingroup$

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).)

Improve this answer
$\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.