Ondrej and I implemented the following method to the Sympy Matrix class last night:
def get_diag_blocks(self):
sub_blocks = []
def recurse_sub_blocks(M):
i = 1
while i <= M.shape[0]:
if i == 1:
to_the_right = M[0, i:]
to_the_bottom = M[i:, 0]
else:
to_the_right = M[0:i, i:]
to_the_bottom = M[i:, 0:i]
if any(to_the_right) or any(to_the_bottom):
i += 1
continue
else:
sub_blocks.append(M[0:i, 0:i])
if M.shape == M[0:i, 0:i].shape:
return
else:
recurse_sub_blocks(M[i:, i:])
return
recurse_sub_blocks(self)
return sub_blocks
So it essentially finds the upper left sub-block, then recursively calls its self on the remainder of the lower left of the matrix. Our search for subblocks just increments by 1 at each iteration, and this should probably be improved to be a bisect method, but we didn't get that done yet.
This would also be nice to put into numpy, along with a standalone function which creates a block diagonal matrix given a list of square matrices, like this one we also implemented last night:
def block_diag(matrices):
"""
Constructs a block diagonal matrix.
Example:
>>> from sympy import block_diag, symbols
>>> x, y, z = symbols("x y z")
>>> a = Matrix([[1, 2], [2, 3]])
>>> b = Matrix([[3, x], [y, 3]])
>>> block_diag([a, b, b])
[1, 2, 0, 0, 0, 0]
[2, 3, 0, 0, 0, 0]
[0, 0, 3, x, 0, 0]
[0, 0, y, 3, 0, 0]
[0, 0, 0, 0, 3, x]
[0, 0, 0, 0, y, 3]
"""
rows = 0
for m in matrices:
rows += m.rows
A = zeros((rows, rows))
i = 0
for m in matrices:
A[i+0:i+m.rows, i+0:i+m.cols] = m
i += m.rows
return A
~Luke
