How to combine polynomials in matrix operations in Sympy?

Question

I'm doing some matrix operations, sometimes involving matrices whose entries have constant values.

But for some reason, I cannot get the matrix operation to combine the results into a single polynomial, even when the result is simple. for example, consider the following:

from sympy.matrices import *
import sympy

x= sympy.symbol.symbols('x')

Poly_matrix = sympy.Matrix([[sympy.Poly(x, x)],[sympy.Poly(x, x)]])
constant_matrix = sympy.Matrix([[0,1]])
constant_matrix_poly = constant_matrix.applyfunc(lambda val: sympy.Poly(val, x))

# this doesn't combine them
result1 = constant_matrix * Poly_matrix
print result1
>>>> Matrix([[Poly(0, x, domain='ZZ') + Poly(x, x, domain='ZZ')]])

# even THIS doesn't combine them when I convert the constants to Polynomials 
result = constant_matrix_poly * Poly_matrix
print result
>>>> Matrix([[Poly(0, x, domain='ZZ') + Poly(x, x, domain='ZZ')]])

The problem with this, is that when I try to extract the expression, and turn this result into a different polynomial, I get the following error:

# This is where the trouble starts
sympy.Poly(result[0].as_expr(), x)
sympy.Poly(result1[0].as_expr(), x)

And the resulting error is a long traceback, ending with:

PolynomialError: Poly(x, x, domain='ZZ') contains an element of the set of generators.

To be even more specific, it has trouble with result[0].as_expr() because it cannot convert it to an expression using as_expr(), even though it is still an object of type Poly, so it can still use the method as_expr().

Why is it that these polynomials do not automatically get combined into one?
Or is there another way for me to call sympy.Poly(result[0].as_expr(), x)?

EDIT: Here are some questions with a similar error (although sometimes caused by something different):
sympy: PolynomialError: cos(a) contains an element of the generators set Sympy Error when using POLY with SQRT

Answer 1

I stumbled upon this issue some time ago while running some code from a post. After looking into the source code of the matrix multiplication function _eval_matrix_mul on line 174 of dense.py, it turns out that sympy.Add is used to perform addition during the computation process rather than the + operator. Therefore Poly.add is not invoked and a new expression is created for addition every time.

Further inspection into the source code reveals that there is a PolyMatrix class that rewrites the matrix multiplication functions for polynominals. This class does work as expected as shown in the following. However, it does not show up anywhere in the documentation for unknown reasons so use it with caution. The docstring in the linked source code provides basic documentation for the class.

from sympy import Poly
from sympy.abc import x
from sympy.polys.polymatrix import PolyMatrix

mat = PolyMatrix([[Poly(x, x)], [1]])
const_mat = PolyMatrix([[4, 3]])
print(mat.shape, const_mat.shape)
print(mat.T * mat)
print(mat * mat.T)
print(2 * mat)
print(2 * mat + const_mat.T)

Output:

(2, 1) (1, 2)
Matrix([[Poly(x**2 + 1, x, domain='ZZ')]])
Matrix([[Poly(x**2, x, domain='ZZ'), Poly(x, x, domain='ZZ')], [Poly(x, x, domain='ZZ'), 1]])
Matrix([[Poly(2*x, x, domain='ZZ')], [2]])
Matrix([[Poly(2*x + 4, x, domain='ZZ')], [5]])

Another alternative approach is to use Expr.collect which has the same functionality as sympy.collect, as shown in the following:

from sympy import Poly, Matrix
from sympy.abc import x

mat = Matrix([[Poly(x, x)], [1]])
result = mat.T * mat
simplified = result.applyfunc(lambda p: p.collect(x))
print(simplified)

Output:

Matrix([[Poly(x**2 + 1, x, domain='ZZ')]])

Answer 2

If all you want is the characteristic polynomial, use charpoly. This is more efficient than eigenvals, because sometimes symbolic roots can be expensive to calculate.

Sample

lamda = symbols('lamda') p = M.charpoly(lamda) factor(p) 2 (λ - 5) ⋅(λ - 3)⋅(λ + 2)

Source