7
$\begingroup$

Is there a better method to invert a large block tridiagonal Hermitian block matrix, other than treating it as a ordinary matrix?

For example:

Do[Evaluate@ToExpression["a" <> ToString[i] <> ToString[j]] = 
  ConstantArray[0, {3, 3}], {i, 1, 4}, {j, 1, 4}]

a11 = DiagonalMatrix@RandomInteger[{1, 10}, 3];
a22 = DiagonalMatrix@RandomInteger[{1, 10}, 3];
a33 = DiagonalMatrix@RandomInteger[{1, 10}, 3];
a44 = DiagonalMatrix@RandomInteger[{1, 10}, 3];

a12 = a21 = Partition[RandomInteger[{1, 10}, 9], 3];
a23 = a32 = Partition[RandomInteger[{1, 10}, 9], 3];
a34 = a43 = Partition[RandomInteger[{1, 10}, 9], 3];

s = ArrayFlatten[Evaluate@Table[
      ToExpression["a" <> ToString[i] <> ToString[j]], {i, 1, 4}, {j, 1,4}]]

The elements live in tridiagonal blocks as well as along the diagonal line.

enter image description here

Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ If you are solving a linear equation, then you can easily remove the diagonal terms from both side of the equation. Then the block tridiagonal matrix can be inverse block by block. $\endgroup$ – xslittlegrass Mar 23 '13 at 3:01
  • 3
    $\begingroup$ The algorithm for tridiagonal matrices works with blocks too if you pay attention to the order of multiplication (which is not commutative). You can easily write a procedural function that implements it. @xslittlegrass -- What? How? The best algorithm for tridiagonal systems requires a forward-backward scheme. You can't just remove the diagonal from both sides of the equation (at least for there is NO diagonal on the other side of the equation!). $\endgroup$ – Federico Mar 23 '13 at 3:12
  • $\begingroup$ @Federico Yes you are right, there is no easy way to remove the diagonal terms. When I saw this matrix, I was thinking about Shrodinger equation for a system of several energy levels interaction with a linear polarized classical electric field. The matrix is the total Hamiltonian of the the atomic Hamiltonian plus the dipole interaction Hamiltonian. The atom Hamiltonian, ie. the diagonal term will cancel out, only left the tridiagonal blocks come from the dipole coupling of different states. $\endgroup$ – xslittlegrass Mar 23 '13 at 4:12
  • $\begingroup$ Do you really, truly need the inverse? For instance, if you're solving a linear system whose coefficient matrix just happens to be block tridiagonal, then use a sparse solver instead of computing the inverse, which is more often than not much denser than the original matrix... $\endgroup$ – J. M.'s ennui♦ Mar 23 '13 at 13:52
6
$\begingroup$

Basic implementation

Here is a function BlockTridiagonalSolve that takes three lists of blocks (diag, lower and upper) and a list of vector pieces (vec) and solves the corresponding linear system:

BlockTridiagonalSolve[diag_?(ArrayQ[#, 3] &), lower_?(ArrayQ[#, 3] &), upper_?(ArrayQ[#, 3] &), vec_?MatrixQ] := 
 Module[{a, i, n = Length[diag], d = diag, v = vec},
  For[i = 1, i < n, i++,
   a = lower[[i]].Inverse[d[[i]]];
   d[[i + 1]] -= a.upper[[i]];
   v[[i + 1]] -= a.v[[i]]];
  v[[n]] = Inverse[d[[n]]].v[[n]];
  For[i = n - 1, i > 0, i--,
   v[[i]] = Inverse[d[[i]]].(v[[i]] - upper[[i]].v[[i + 1]])];
  v]

To see it in action you can use the following piece of code:

nBlocks = 50;
blockSize = 10;
diag = DiagonalMatrix /@ RandomReal[{-1, 1}, {nBlocks, blockSize}];
lower = RandomReal[1, {nBlocks - 1, blockSize, blockSize}];
upper = RandomReal[-1, {nBlocks - 1, blockSize, blockSize}];
vec = RandomReal[{-1, 1}, {nBlocks, blockSize}];
trimat = SparseArray[{
    Band[{1, 1}] -> diag,
    Band[{blockSize + 1, 1}] -> lower,
    Band[{1, blockSize + 1}] -> upper}];
bts = BlockTridiagonalSolve[diag, lower, upper, vec];
ls = LinearSolve[trimat, Flatten[vec]];
maxError = Max[Abs[Flatten[bts] - ls]]

Note that BlockTridiagonalSolve is slower than LinearSolve for small inputs, and grows comparable only for large blockSize values.

However, the following fact has to be taken into consideration: if you happen to be able to naturally construct the matrices diag, lower, and upper the way BlockTridiagonalSolve requires them, then BlockTridiagonalSolve may be the way to go, because a huge amount of time is spent (i.e. wasted...) contructing trimat, which is required to call LinearSolve. You can experience this for example with nBlocks = 100 and blockSize = 100, in which case trimat = SparseArray[...]; takes roughly 10 seconds on my Intel-i7 PC.

This is an important general observation: the best algorithm may depend on the type of input you are able to provide, because conversion can be time consuming too.

Considerations

Of course you may try and compile BlockTridiagonalSolve or perform other tricks (such as replacing Inverse with LinearSolve at lines 7 and 9). Last but not least, you should consider how BlockTridiagonalSolve is meant to be used: will you sporadically call it with different matrices, or will you repeatedly call it with the same matrix always? This makes a huge difference, because in the latter case you may want to "precompile" a specific solver as you would do by calling LinearSolve with one argument so that it generates a LinearSolveFunction.

Improve this answer
$\endgroup$
8
  • $\begingroup$ I would recommend using LinearSolve[] instead of Inverse, actually... to use one of the lines of your code as an example, v[[i]] = LinearSolve[d[[i]], v[[i]] - upper[[i]].v[[i + 1]]]. $\endgroup$ – J. M.'s ennui♦ Mar 23 '13 at 16:32
  • $\begingroup$ Actually, I tested both (in the compiled version, to be precise) and the one with LinearSolve does not always perform better, depending on the input size. My conjecture is that LinearSolve spends a little amount of time to examine the matrix and figure out which algorithm to use $\endgroup$ – Federico Mar 23 '13 at 16:50
  • $\begingroup$ @Federico Thanks very much for the answer. Sometimes I have zero blocks in place of the upper or lower band, so a inverse of that block will gives error. Do you know how to deal with that? $\endgroup$ – user0501 Mar 27 '13 at 21:07
  • $\begingroup$ Note that the algorithm requires only the inversion of blocks on the main diagonal. Anyway, if you have zero blocks in place of the upper or lower band, your matrix is actually block-triangular, so that you need to perform only the forward or the backward part of the algorithm (depending whether the matrix is lower or upper triangular). $\endgroup$ – Federico Mar 27 '13 at 23:38
  • $\begingroup$ @Federico Thanks, I see you point. Did you use the formula of the Thomas algorithm in the link in your first comment, and change the "times" to "dot" and "divide" to "dot inverse"? I have a difficult time to associate the variables in your code to the variables in that formula. Could you maybe write down the formula used in your code? Thank you very much. $\endgroup$ – user0501 Mar 28 '13 at 1:26

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.