Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
    
The simplest way may be to define a function which does exactly that. OF course in the general case this is complicated, but if you're doing only 2x2 or 2^n x 2^n matrices it shouldn't be terribly hard. –  Mike Bantegui Jan 17 '12 at 23:26
    
@MikeBantegui Any specific suggestions? A function is fine. –  Eli Lansey Jan 17 '12 at 23:27
    
I'll add an answer that demonstrates this. –  Mike Bantegui Jan 17 '12 at 23:31
add comment

3 Answers

up vote 17 down vote accepted

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.

share|improve this answer
add comment

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

partitioned matrix

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]}};

enter image description here

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]}}]
share|improve this answer
add comment

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.

share|improve this answer
    
Sorry for the late comment, but I thought block wise inversion was very well known for 4-blocks, provided that the blocks on the diagonal were square. –  rcollyer Jan 30 '12 at 15:58
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.