8

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

makansij
  • 9,303
  • 37
  • 105
  • 183

3 Answers3

2

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')]])
GZ0
  • 4,055
  • 1
  • 10
  • 21
  • Interesting, I didnk't know about `PolyMatrix`. Can you also critique my solution below? – makansij Sep 26 '19 at 21:28
  • 1
    @Sother I've added another approach. – GZ0 Sep 26 '19 at 23:06
  • this is great thanks. So, is there ever any reason to use `Matrix` instead of `PolyMatrix`? – makansij Sep 28 '19 at 23:16
  • 1
    The only possible reason is `PolyMatrix` is not in the official documentation so its status is unknown. If you want to be sure, you may post a question on `sympy`'s GitHub issues page to ask about it. – GZ0 Sep 28 '19 at 23:22
  • Your example works (thanks!). But for some reason when I try to create a function that converts matrices to PolyMatrix objects, this doesn't work. Can I message you in chat to continue this discussion so as not to disturb things here? – makansij Sep 29 '19 at 02:52
  • 1
    @Sother In case you have not seen it, I had replied in the chat. – GZ0 Sep 30 '19 at 14:27
  • @GZ0 "...it turns out that `sympy.Add` is used to perform addition..." I believe you are wrong. `sympy.Add` is used at least since version `1.1.1`. And in version `1.1.1` computation performs without error: `constant_matrix * Poly_matrix` gives you `Matrix([[Poly(x, x, domain='ZZ')]])`. – Vadim Shkaberda Sep 30 '19 at 19:12
  • @GZ0 The issue is how `sympy.Poly` multiplies by `sympy.core.numbers`. In version `1.1.1` multiplication `Poly(x, x, domain='ZZ')` by `sympy.core.numbers.Integer(N)` gave you `N * Poly(x, x, domain='ZZ')`. And since version `1.2` the output became `Poly(N*x, x, domain='ZZ')`. – Vadim Shkaberda Sep 30 '19 at 19:45
  • @VadimShkaberda The fix in the [`PolyMatrix` class](https://github.com/sympy/sympy/blob/master/sympy/polys/polymatrix.py) avoids the use of `Add` and there is an explicit comment on line 75 saying it should not be used. In the latest version `1.4` the multiplication of polynomials has already been simplified during matrix multiplication. It is the additions that remain as expressions. I am not sure why version `1.1.1` would 'work' (I doubt it returns the same for the examples in my answer) even though your comment seems to suggest that it does not even simplify polynomial multiplications. – GZ0 Sep 30 '19 at 20:02
  • 1
    @GZ0 In version `1.1.1` the same answers are returned except `mat.T * mat` which returns `Matrix([[Poly(x**2, x, domain='ZZ') + 1]])`. Indeed, problem is due to the `sympy.Add`. But it's because of behaviour of `core.Add` was changed. In version `1.1.1` output of `sympy.Add(Poly(2*x, x) + Poly(x, x))` is `Poly(3*x, x)`, and since version `1.2` it never simplifies and remains `Poly(x, x) + Poly(2*x, x)`. – Vadim Shkaberda Sep 30 '19 at 20:29
  • 1
    @Sother May I ask why wasn't the bounty awarded? – GZ0 Oct 06 '19 at 15:52
  • Yikes I didn't log into my account for a while I'm so sorry! this answer totally deserves the bounty. – makansij Oct 15 '19 at 02:14
  • @Sother You could still start another bounty and award it, if you want. – GZ0 Oct 15 '19 at 03:41
1

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

champion-runner
  • 1,489
  • 1
  • 13
  • 26
0

Here is one way to do this: 

# even though this prints as two terms (which is kind of weird)

print result[0].as_expr()

#...this collapses it into one
exec('poly_expr = '+str(result[0].as_expr()))
print poly_expr

#...which allows the expressions to be combined

Poly(poly_expr, x)

However, this seems like a very hackish solution. Is there an easier way to do this?

makansij
  • 9,303
  • 37
  • 105
  • 183
  • 1
    I think you may just do `output = Poly(str(result[0]), x)` to get the same result. – GZ0 Sep 26 '19 at 21:56