Speeding up computation of symbolic determinant in SymPy

Question

I have a 4x4 matrix A with rather long but simple symbolic expressions in each of its entries. About 30 different symbols are involved. By "simple" I mean that these symbols are combined using only addition/subtraction, multiplication/division, and integer powers. By "long" I mean that if I print out the matrix, it covers three or four screens worth.

I need the determinant of this matrix. Or, to be more specific, I know that the determinant is a fourth-order polynomial in one particular symbol, and I need the coefficients of this polynomial. A.det() does not terminate after hours and hours of running, so I need a different approach. Any ideas? So far I've tried to throw various simplify functions at each element of A without any success.

Is there some strategy I can employ to let SymPy be aware of the simple structure of my expressions, or that I know that the result is a polynomial in one of the symbols?

Chaoan

@Meta: No, the entries can be arbitrary sums, products and integer powers of that symbol. – Chaoan

@Meta: No, I allow negative powers of the symbol too, so it's not a polynomial. Maybe SymPy's polynomial module can still handle that though? In general, the entries are ratios of polynomials in the specific symbol. – Chaoan

@Meta: Isn't it kinda confusing to delete comments? Now mine make no sense. Anyway: A typical matrix is 4400 characters long when printed out… I don't think posting it here makes much sense. – Chaoan

Answer 1

Maybe it would work to create the general expression for a 4x4 determinant

In [30]: A = Matrix(4, 4, symbols('A:4:4'))

In [31]: A
Out[31]:
⎡A₀₀  A₀₁  A₀₂  A₀₃⎤
⎢                  ⎥
⎢A₁₀  A₁₁  A₁₂  A₁₃⎥
⎢                  ⎥
⎢A₂₀  A₂₁  A₂₂  A₂₃⎥
⎢                  ⎥
⎣A₃₀  A₃₁  A₃₂  A₃₃⎦

In [32]: A.det()
Out[32]:
A₀₀⋅A₁₁⋅A₂₂⋅A₃₃ - A₀₀⋅A₁₁⋅A₂₃⋅A₃₂ - A₀₀⋅A₁₂⋅A₂₁⋅A₃₃ + A₀₀⋅A₁₂⋅A₂₃⋅A₃₁ + A₀₀⋅A₁₃⋅A₂₁⋅A₃₂ - A₀₀⋅A₁₃⋅A₂₂⋅A₃₁ - A₀₁⋅A₁₀⋅A₂₂⋅A₃₃ + A₀₁⋅A₁₀⋅A₂₃⋅A₃₂ + A₀₁⋅A₁₂⋅A₂₀⋅
A₃₃ - A₀₁⋅A₁₂⋅A₂₃⋅A₃₀ - A₀₁⋅A₁₃⋅A₂₀⋅A₃₂ + A₀₁⋅A₁₃⋅A₂₂⋅A₃₀ + A₀₂⋅A₁₀⋅A₂₁⋅A₃₃ - A₀₂⋅A₁₀⋅A₂₃⋅A₃₁ - A₀₂⋅A₁₁⋅A₂₀⋅A₃₃ + A₀₂⋅A₁₁⋅A₂₃⋅A₃₀ + A₀₂⋅A₁₃⋅A₂₀⋅A₃₁ - A₀₂⋅A₁
₃⋅A₂₁⋅A₃₀ - A₀₃⋅A₁₀⋅A₂₁⋅A₃₂ + A₀₃⋅A₁₀⋅A₂₂⋅A₃₁ + A₀₃⋅A₁₁⋅A₂₀⋅A₃₂ - A₀₃⋅A₁₁⋅A₂₂⋅A₃₀ - A₀₃⋅A₁₂⋅A₂₀⋅A₃₁ + A₀₃⋅A₁₂⋅A₂₁⋅A₃₀

and then substitute in the entries with something like

A.det().subs(zip(list(A), list(your_matrix)))

SymPy being slow to generate a 4x4 determinant is a bug, though. You should report it at https://github.com/sympy/sympy/issues/new.

EDIT (this wouldn't fit in a comment)

It looks like Matrix.det is calling a simplification function. For matrices 3x3 and smaller, the determinant formula is written out explicitly, but for larger matrices, it is computed using the Bareis algorithm. You can see where the simplification function (cancel) is called here, which is necesssary as part of the computation, but end up doing a lot of work because it tries to simplify your very large expressions. It would probably be smarter to only do the simplifications that are needed to cancel terms of the determinant itself. I opened an issue for this.

Another possibility to speed this up, which I'm not sure will work or not, would be to select a different determinant algorithm. The options are Matrix.det(method=alg) where alg is one of "bareis" (the default), "berkowitz", or "det_LU".

Furbelow

Wow! The det method must be doing some implicit simplification, because with your suggestion the whole process is done in seconds. Thanks! – Chaoan

Answer 2

What you describe as a ratio of polynomials is what is known as a rational function: https://en.wikipedia.org/wiki/Rational_function

SymPy's polys module does have ways of representing rational functions although they can be slow especially with lots of variables.

There is a new matrix implementation in sympy 1.7 which is still somewhat experimental but is based on the polys module and can handle rational functions. We can test it here by quickly creating a random matrix:

In [35]: import random                                                                                                            

In [36]: from sympy import random_poly, symbols, Matrix                                                                           

In [37]: randpoly = lambda : random_poly(random.choice(symbols('x:z')), 2, 0, 2)                                                  

In [38]: randfunc = lambda : randpoly() / randpoly()                                                                              

In [39]: M = Matrix([randfunc() for _ in range(16)]).reshape(4, 4)                                                                

In [40]: M                                                                                                                        
Out[40]: 
⎡     2              2            2              2       ⎤
⎢  2⋅z  + 1       2⋅z  + z     2⋅z  + z + 2     x  + 2   ⎥
⎢  ────────     ────────────   ────────────   ────────── ⎥
⎢   2            2                 2             2       ⎥
⎢  y  + 2⋅y     y  + 2⋅y + 1      x  + 1      2⋅z  + 2⋅z ⎥
⎢                                                        ⎥
⎢  2              2                  2        2          ⎥
⎢ y  + y + 1   2⋅x  + 2⋅x + 1       z        z  + 2⋅z + 1⎥
⎢ ──────────   ──────────────     ──────     ────────────⎥
⎢     2            2               2           2         ⎥
⎢  2⋅y  + 2       y  + 2⋅y        y  + 1      x  + x + 2 ⎥
⎢                                                        ⎥
⎢     2           2                2            2        ⎥
⎢  2⋅z  + 2    2⋅z  + 2⋅z + 2     y  + 1     2⋅y  + y + 2⎥
⎢────────────  ──────────────   ──────────   ────────────⎥
⎢   2             2                2             2       ⎥
⎢2⋅z  + z + 1  2⋅x  + 2⋅x + 2   2⋅y  + 2⋅y      x  + 2   ⎥
⎢                                                        ⎥
⎢    2               2            2             2        ⎥
⎢ 2⋅y  + 2⋅y      2⋅y  + y     2⋅x  + x + 1  2⋅x  + x + 1⎥
⎢ ──────────      ────────     ────────────  ────────────⎥
⎢    2              2                 2          2       ⎥
⎣   z  + 2         x  + 2          2⋅y          x  + 2   ⎦

If we convert that to the new matrix implementation then we can compute the determinant using the charpoly method:

In [41]: from sympy.polys.domainmatrix import DomainMatrix                                                                        

In [42]: dM = DomainMatrix.from_list_sympy(*M.shape, M.tolist())                                                                  

In [43]: dM.domain                                                                                                                
Out[43]: ZZ(x,y,z)

In [44]: dM.domain.field                                                                                                          
Out[44]: Rational function field in x, y, z over ZZ with lex order

In [45]: %time det = dM.charpoly()[-1] * (-1)**M.shape[0]                                                                         
CPU times: user 22 s, sys: 231 ms, total: 22.3 s
Wall time: 23 s

This is slower than the approach suggested by @asmeurer above but it produces output in a canonical form as a ratio of expanded polynomials. In particular this means that you can immediately tell if the determinant is zero (for all x, y, z) or not. The time is also taken by the equivalent of cancel but the implementation is more efficient than Matrix.det.

How long this takes is largely a function of how complicated the final output is and you can get some sense of that from the length of its string representation (I won't show the whole thing!):

In [46]: len(str(det))                                                                                                            
Out[46]: 54458

In [47]: str(det)[:80]                                                                                                            
Out[47]: '(16*x**16*y**7*z**4 + 48*x**16*y**7*z**2 + 32*x**16*y**7 + 80*x**16*y**6*z**4 + '

At some point it should be possible to integrate this into the main Matrix class or otherwise to make the DomainMatrix class more publicly accessible.

Roeser