CODEWITHSUNDEEP
sympy
jonaslb
jonaslb

Reputation: 174

Solving a matrix equation containing MatrixSymbols of symbolic size in Sympy?

As an introduction i want to point out that if one has a matrix A consisting of 4 submatrices in a 2x2 pattern, where the diagonal matrices are square, then if we denote its inverse as X, the submatrix X22 = (A22-A21(A11^-1)A12)^-1, which is quite easy to show by hand.

I was trying to do the same for a matrix of 4x4 submatrices, but its quite tedious by hand. So I thought Sympy would be of some help. But I cannot figure out how (I have started by just trying to reproduce the 2x2 result).

I've tried:

import sympy as s

def blockmatrix(name, sizes, names=None):
    if names is None:
        names = sizes
    ll = []
    for i, (s1, n1) in enumerate(zip(sizes, names)):
        l = []
        for j, (s2, n2) in enumerate(zip(sizes, names)):
            l.append(s.MatrixSymbol(name+str(n1)+str(n2), s1, s2))
        ll.append(l)
    return ll

def eyes(*sizes):
    ll = []
    for i, s1 in enumerate(sizes):
        l = []
        for j, s2 in enumerate(sizes):
            if i==j:
                l.append(s.Identity(s1))
                continue
            l.append(s.ZeroMatrix(s1, s2))
        ll.append(l)
    return ll

n1, n2 = s.symbols("n1, n2", integer=True, positive=True, nonzero=True)
M = s.Matrix(blockmatrix("m", (n1, n2)))
X = s.Matrix(blockmatrix("x", (n1, n2)))
I = s.Matrix(eyes(n1, n2))
s.solve(M*X[:, 1:]-I[:, 1:], X[:, 1:])

but it just returns an empty list instead of the result.

I have also tried:

but all without luck.

Can you show how to derive a result with Sympy similar to the one I gave as an example for X22?

The most similar other questions on solving with MatrixSymbols seem to have been solved by working around doing exactly that, by using an array of the inner symbols or some such instead. But since I am dealing with symbolically sized MatrixSymbols, that is not an option for me.

Upvotes: 0

Views: 473

Answers (1)

smichr
smichr

Reputation: 18939

Is this what you mean by a matrix of 2x2 matrices?

>>> a = [MatrixSymbol(i,2,2) for i in symbols('a1:5')]
>>> A = Matrix(2,2,a)
>>> X = A.inv()
>>> print(X[1,1])  # [1,1] instead of [2,2] because indexing starts at 0
a1*(a1*a3 - a3*a1)**(-1)

[You indicated not and pointed out that the above is not correct -- that appears to be an issue that should be resolved.]

I am not sure why this isn't implemented, but we can do the solving manually as follows:

>>> 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))
>>> for i in range(n**2):
...   s = solve(eqs[i],v[i])[0]
...   eqs[i+1:] = [e.subs(v[i],s) for e in eqs[i+1:]]
...
>>> s  # solution for v[3] which is B22
(-a2*a0**(-1)*a1 + a3)**(-1)

You can change n to 3 and see a modestly complicated expression. Change it to 4 and check the result by hand to give a new definition to the word "tedious" ;-)

The special structure of the equations to be solved can allow for a faster solution, too: the variable of interest is the last factor in each term containing it:

>>> 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))
...   s = 1/d.subs(v[i], 1)*-c
...   eqs[i+1:] = [e.subs(v[i], s) for e in eqs[i+1:]]

Upvotes: 1

Related Questions