inversion of block matrix of submatrices

[smichr]

This question was raised on SO and I had some difficulty coming up with a solution. It seems like something that should work with SymPy.

Consider a block matrix of submatrices:

>>> subm = [MatrixSymbol(i,2,2) for i in symbols('a:4')]
>>> A = BlockMatrix(reshape(subm, [2])); A
Matrix([
[a0, a1],
[a2, a3]])

If we look at bottom right submatrix after inversion we find that it is only implicitly defined:

>>> A.inv()[1,1]
(Matrix([
[a0, a1],
[a2, a3]])**(-1))[1, 1]

If we create a matrix of the symbols (which are non-commutative) we get a more explicit (but wrong) expression:

>>> Matrix(2,2,subm) <-- also note that this is indistiguishable in printing from BlockMatrix
Matrix([
[a0, a1],
[a2, a3]])
>>> _.inv()[-1]
a0*(a0*a2 - a2*a0)**(-1)

Finally, I solved the equation by hand:

n = 2
v = symbols('b:%s'%n**2,commutative=False)
A = Matrix(n,n,symbols('a:%s'%n**2,commutative=False))
B = Matrix(n,n,v)
eqs = list(A*B - eye(n))  <-- solve for b0-b3
for i in range(n**2):
  c,d = eqs[i].expand().as_independent(v[i])
  assert all(j.args[-1]==v[i] for j in Add.make_args(d))
  vi = 1/d.subs(v[i],1)*-c
  eqs[i] = Eq(v[i], vi)
  eqs[i+1:] = [e.subs(v[i], vi) for e in eqs[i+1:]]

print(eqs[-1])
Eq(b3, (-a2*a0**(-1)*a1 + a3)**(-1))

(and that is the right, explicit expression)

I think that will apply - result and solution process -- provided that the matrices on the diagonal are square. e.g.

>>> a=[randMatrix(i,j) for i,j in [(1,1),(1,2),(2,1),(2,2)]]
>>> Matrix.irregular(2,*a).inv()
Matrix([
[-534/21013, 1948/63039, -317/63039],
[ 628/21013, -953/63039,   58/63039],
[  25/21013, -406/63039, 1943/63039]])
>>> (-a[2]*a[0]**(-1)*a[1] + a[3])**(-1)
Matrix([
[-953/63039,   58/63039],
[-406/63039, 1943/63039]])

@oscargus , I think this is the first sighting of the irregular method in the wild!

BTW, trying to solve eqs with solve did not work with or without manual=True:

>>> solve(eqs,v)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "sympy\solvers\solvers.py", line 1173, in solve
    solution = _solve_system(f, symbols, **flags)
  File "sympy\solvers\solvers.py", line 1785, in _solve_system
    poly = g.as_poly(*symbols, extension=True)
  File "sympy\core\basic.py", line 788, in as_poly
    poly = Poly(self, *gens, **args)
  File "sympy\polys\polytools.py", line 109, in __new__
    opt = options.build_options(gens, args)
  File "sympy\polys\polyoptions.py", line 731, in build_options
    return Options(gens, args)
  File "sympy\polys\polyoptions.py", line 154, in __init__

    preprocess_options(args)
  File "sympy\polys\polyoptions.py", line 152, in preprocess_options
    self[option] = cls.preprocess(value)
  File "sympy\polys\polyoptions.py", line 293, in preprocess
    raise GeneratorsError("non-commutative generators: %s" % str(gens))
sympy.polys.polyerrors.GeneratorsError: non-commutative generators: (b0, b1, b2, b3)
>>> solve(eqs,v,manual=True)
[]

See also #5858

SO issue is here

Note: one would have to back substitution the solutions to create the full explicit inverse of A; what I have shown here is only the forward solving part. A naive approach is to do

for i in range(-1, -1-len(eqs), -1):
  rep = dict([eqs[i].args])
  for j in range(i-1,-1-len(eqs),-1):
    eqs[j] = eqs[j].xreplace(rep)
  eqs[i] = eqs[i].rhs

invA = Matrix(n,n,eqs)

[smichr]

For simplicity, it is sufficient to consider the 2x2 case: a 3x3 matrix could be imagined to be partitioned into a 2x2 and 1x1 on the diagonal and 2x1 and 1x2 on the off-digonal (and similar for any other sized matrix). This wikipage also indicates that the inverse can be computed so it involves different regions of the original matrix being inverted.

[BoltonBailey]

It would also be useful to be able to simplify the determinant of block matrices.