3

I have a Sympy matrix A and a polynomial expression P, and I would like to compute P(A).

Here is an example:

x = Symbol('x')
P = x**2 - 3*x + 5
A = Matrix([ [1,3], [-1,2] ])
P.subs(x, A)

I expect Sympy to compute A**2 - 3*A + 5*eye(2) (in the example, the result is the zero matrix).

But this fails with an error message:

AttributeError: ImmutableMatrix has no attribute as_coeff_Mul.

Is there any way to obtain what I want ?

Edit: I've tried to convert P into Sympy's polynomial class, and substitute afterwards, but the result is useless:

Poly(P).subs(A)
Poly(Matrix([ [ 1, 3], [-1, 2]])**2 - 3*Matrix([ [ 1, 3],        
    [-1, 2]]) + 5, Matrix([ [ 1, 3], [-1, 2]]), domain='ZZ')

I can get the correct result with the following function :

def poly_matrix(P, A):
    coeffs = Poly(P).all_coeffs()[::-1]
    res = zeros(A.rows)
    for i in range(len(coeffs)):
        res += coeffs[i]*(A**i)
    return res

But I'm still looking for a more efficient built-in option.

Adrien
  • 469
  • 3
  • 11

3 Answers3

2

Matrix expressions are still a bit buggy, hopefully they'll get fixed in the future.

Anyway, here's an alternative way to perform the substitution:

In [1]: x = MatrixSymbol('x', 2, 2)

In [2]: P = x**2 - 3*x + 5*eye(2)

In [3]: P
Out[3]: 
                 2
⎡5  0⎤ + -3⋅x + x 
⎢    ⎥            
⎣0  5⎦            

In [4]: A = Matrix([ [1,3], [-1,2] ])

In [5]: P.subs(x, A)
Out[5]: 
                               2
⎡5  0⎤ + -3⋅⎡1   3⎤ + ⎛⎡1   3⎤⎞ 
⎢    ⎥      ⎢     ⎥   ⎜⎢     ⎥⎟ 
⎣0  5⎦      ⎣-1  2⎦   ⎝⎣-1  2⎦⎠ 

In [6]: P.subs(x, A).doit()
Out[6]: 
                   2
⎡2  -9⎤ + ⎛⎡1   3⎤⎞ 
⎢     ⎥   ⎜⎢     ⎥⎟ 
⎣3  -1⎦   ⎝⎣-1  2⎦⎠

Here it appears that MatPow isn't able to perform the .doit() operation. It's probably another bug.

In output number 5 there's also a bug with the printer (+ -3 should just be -3).

I really hope that someone will eventually fix all these problems.

Francesco Bonazzi
  • 1,917
  • 9
  • 10
0

Try to evaluate each term of the polynomial.

(x**2).subs(x,A) - (3*x).subs(x,A) + 5*(eye(2))

This will evaluate your expression.

rresol
  • 323
  • 4
  • 20
0

To evaluate P at A, you can substitute x with a Matrix and convert the constant term by multiplying by eye(2):

P_ = Poly(P, x)
(P_ - P_.coeff_monomial(1)).as_expr().subs(x, A) * eye(2) + P_.coeff_monomial(1) * eye(2) # P(A)

The first multiplication by eye(2) ensures the first term in the sum is a matrix, even if P is just a constant, i.e. P_ - P_.coeff_monomial(1) == 0.

llf
  • 610
  • 10
  • 11