recursive block_collapse

[Chris]

Dear all,

I'm trying to use block matrices in the following way:

    n1, n2 = symbols('n1 n2', integer=True)

    U1 = MatrixSymbol('U1', n1, n1+n2) 

    U2 = MatrixSymbol('U2', n2, n1+n2)

    D1 = MatrixSymbol('D1', n1, n1) 

    D2 = MatrixSymbol('D2', n2, n2) 

    D = BlockDiagMatrix(D1, D2) 

    U = BlockMatrix([[U1], [U2]])

    K = U * D * U.T

    Kc = sympy.block_collapse(K)

For the record, K looks like this:

In [32]: K
Out[32]:
              ⎡  T    T⎤
⎡U₁⎤⋅⎡D₁  0 ⎤⋅⎣U₁   U₂ ⎦
⎢  ⎥ ⎢      ⎥          
⎣U₂⎦ ⎣0   D₂⎦          

Kc looks as follows:

⎡                T                  T⎤

⎢⎡U₁⎤⋅⎡D₁  0 ⎤⋅U₁   ⎡U₁⎤⋅⎡D₁  0 ⎤⋅U₂ ⎥

⎢⎢  ⎥ ⎢      ⎥      ⎢  ⎥ ⎢      ⎥    ⎥

⎣⎣U₂⎦ ⎣0   D₂⎦      ⎣U₂⎦ ⎣0   D₂⎦    ⎦

In [30]: Kc.blockshape
Out[30]: (1, 2)

What I'm looking for is a way to recursively break down these blocks further (i.e. a blockshape of (2,2)). I tried the following, to no avail:

In [31]: sympy.block_collapse(Kc.blocks[0])
Out[31]:


                T
⎡U₁⎤⋅⎡D₁  0 ⎤⋅U₁
⎢  ⎥ ⎢      ⎥    
⎣U₂⎦ ⎣0   D₂⎦    

Any ideas on how to achieve this?

Thank you!

Chris

[Matthew Rocklin]

Hi Chris, 

In order to simplify these further you'll need to break down U1 and U2 so that they are of compatible shape with D1 and the 0-matrix.  At the moment there is no way to take these larger matrices and simplify them against smaller ones.  This mismatch would be more evident if we had size-sensitive printing.  For example seeing the U_1^T on the right as a big solid block rather than a single letter would make this mismatch easier to spot.

One way to cut the matrices into smaller blocks is with blockcut.  The inputs are input matrix, desired row sizes, desired block sizes.  Output is a blockmatrix of MatrixSlice objects

U1 = blockcut(U1, (n1,), (n1, n2))

U2 = blockcut(U2, (n2,), (n1, n2))

They look like this

In [13]: U1

Out[13]: [U₁[:n₁, :n₁]  U₁[:n₁, n₁:n₁ + n₂]]

After you do this (and go through the rest of your pipeline) then block_collapse reduces things cleanly to a 2x2 block matrix.  All of the indexing can get in the way of clean printing unfortunately.  You might consider creating MatrixSymbols U11, U12, U21, U22 explicitly.

Best,

-Matt

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/4297d8bf-c9a7-4927-bd25-428bc2ea4954%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Aaron Meurer]

Speaking of printing, we need to fix the baseline in the pretty
printer. The printing of MatMul with block matrices looks pretty bad.

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/CAJ8oX-E9B1fuuuirUNYZSPC9ckPAFGjHTxcQxABK%2B%2B4u-U3EsA%40mail.gmail.com.

[Chris]

Hi Matt,

thanks for your help! You're right, the dimensions didn't make sense, so here is the fixed version for future reference.

    n1, n2 = symbols('n1 n2', integer=True)

    S1, S2 = symbols('S1 S2')

    U1 = MatrixSymbol('U1', n1, n1+n2) 

    U2 = MatrixSymbol('U2', n2, n1+n2)

    D = MatrixSymbol('D', n1+n2, n1+n2)

    U = BlockMatrix([[U1], [U2]])

    K = U * D * U.T

    block = sympy.block_collapse(K)

    K_inv = sympy.block_collapse(block.inverse())

Thank you!

Chris