How to symbolically do matrix "Block Inversion"?

Question

Consider a block (partitioned) matrix

matrix = ArrayFlatten[{{a, b}, {c, d}}]

where, a, b, c and d are each matrices themselves. Say, for example,

a = {{a11, a12}, {a21, a22}}
b = {{b11, b12}, {b21, b22}}
c = {{0, 0}, {0, 0}}
d = {{d11, d12}, {d21, d22}}

How can you find the block inverse of this matrix? A desired solution is, using the example above

{{Inverse[a] , -Inverse[a].b.Inverse[d]},{0,Inverse[d]}}

which is easily verified using

Simplify[Inverse[ArrayFlatten[{{a, b}, {c, d}}]] == 
ArrayFlatten[{{Inverse[a], -Inverse[a].b.Inverse[d]}, {0, 
 Inverse[d]}}]]

which yields True.

How can you solve the block inverse problem for arbitrary submatrices, and for block matrices of larger sizes (i.e. 3x3, 4x4, etc)?

Answer 1

Mathematica does not support this directly. You can do things of this sort using an external package called NCAlgebra.

http://math.ucsd.edu/~ncalg/

The relevant documentation may be found at

http://math.ucsd.edu/~ncalg/DOWNLOAD2010/DOCUMENTATION/html/NCBIGDOCch4.html#x8-510004.4

In particular have a look at "4.4.8 NCLDUDecomposition[aMatrix, Options]"

Using this package, you would find the block inverse of the example matrix using:

c=0;    
inverse = NCInverse[matrix]  
(* Out[] = {{inv[a], -inv[a] ** b ** inv[d]}, {0, inv[d]}} *)

Here inv[a] represents the general inverse of the a block of the matrix and the ** represents non-commutative (i.e. matrix) multiplication. This approach works for larger (3x3, 4x4, etc) square block matrices as well.

Answer 2

For a general square matrix m and arbitrary partition of it into conformable parts m={{a,b},{c,d}} (i.e., a and d are square matrices, and b and c have appropriate dimensions),

the formula for the inverse (which can be found in, for example, Review of Matrix Algebra) is

 m={{a,b},{c,d}};
 e = d - c.Inverse[a].b; 
 minv={{Inverse[a] + Inverse[a].b.Inverse[e].c.Inverse[a],
          -Inverse[a].b.Inverse[e]}, 
       {-Inverse[e].c.Inverse[a], Inverse[e]}};

This can be used to define a function that takes a matrix and a partition specification and does the needed type and conformability checks before applying the formula. You can use @Mike's formulation (and, off course, heed his warnings) as a first step:

Clear[sInverse2];
sInverse2[{{a_, b_}, {c_, d_}}] := 
With[{e = d - c.Inverse[a].b}, 
     {{Inverse[a] + Inverse[a].b.Inverse[e].c.Inverse[a], 
         -Inverse[a].b.Inverse[e]}, 
      {-Inverse[e].c.Inverse[a], Inverse[e]}}]

Answer 3

In general this is not an easy thing to do, and a package as Daniel Lichtblau suggested may be your best bet.

However in the specialized case of a 2^n x 2^n matrix, the inversion is very well known.

The case mentioned in your post is:

Clear[sInverse]; (* s for symbolic *)
sInverse[{{a_, b_}, {c_, d_}}] := 
      {{Inverse[a], -Inverse[a] . b . Inverse[d]}, {0, Inverse[d]}}

However you can recursively define this function in terms of itself:

sInverse[{{a_, b_}, {c_, d_}}] := 
 {{sInverse[a], -sInverse[a].b.sInverse[d]}, {0, sInverse[d]}}

This however will crap out on you completely if you actually try this. We have to add some stopping condition to this:

sInverse[{{a_, b_}, {c_, d_}}] /;
  (! ListQ[a]) && (! ListQ[b]) && (! ListQ[c]) && (! ListQ[d]) := 
 a d - b c {{d, -b}, {-c, a}}

This itself also has some issues. Notably, when it reaches the base case you're going to get expressions like -(-a b + c d) . {{m1, m2}, {m3, m4}} . (-p q + r s).

As it stands, this will not get you the exact expression you need. To get things to come out nicely you're going to have to include some more logic in there to flatten the lists out and account for multiplication of matrices and non-matrices. I'm not quite sure how it can be done at this moment, but when I think of it I'll add it.