Problemas al sustituir una matriz en un polinomio

Jan 12 2021

Ejemplo: dejar

M = Matrix([[1,2],[3,4]]) # and 
p = Poly(x**3 + x + 1)    # then
p.subs(x,M).expand()

da el error:

TypeError: no se puede agregar <class'sympy.matrices.immutable.ImmutableDenseMatrix '> y <class' sympy.core.numbers.One '>

lo cual es muy plausible ya que los dos primeros términos se convierten en matrices pero el último término (el término constante) no es una matriz sino un escalar. Para remediar esta situación cambié el polinomio a

p = Poly(x**3 + x + x**0)    # then

persiste el mismo error. ¿Estoy obligado a escribir la expresión a mano, reemplazando x por M? En este ejemplo, el polinomio tiene solo tres términos, pero en realidad encuentro (polinomios multivariados con) docenas de términos.

Respuestas

1 wsdookadr Jan 12 2021 at 12:40

Entonces creo que la pregunta gira principalmente en torno al concepto de polinomio Matrix :

(donde P es un polinomio y A es una matriz)

Creo que esto quiere decir que el término libre es un número, y no se puede sumar con el resto que es una matriz, efectivamente la operación de suma no está definida entre esos dos tipos.

TypeError: no se puede agregar <class'sympy.matrices.immutable.ImmutableDenseMatrix '> y <class' sympy.core.numbers.One '>

Sin embargo, esto puede evitarse definiendo una función que evalúe el polinomio de la matriz para una matriz específica. La diferencia aquí es que estamos usando la exponenciación de la matriz , por lo que calculamos correctamente el término libre del polinomio de la matriz, a_0 * Idonde I=A^0está la matriz de identidad de la forma requerida:

from sympy import *
x = symbols('x')
M = Matrix([[1,2],[3,4]])
p = Poly(x**3 + x + 1)

def eval_poly_matrix(P,A):
    res = zeros(*A.shape)
    for t in enumerate(P.all_coeffs()[::-1]):
        i, a_i = t
        res += a_i * (A**i)
    return res

eval_poly_matrix(p,M)

Producción:

En este ejemplo, el polinomio tiene solo tres términos, pero en realidad encuentro (polinomios multivariados con) docenas de términos.

La función eval_poly_matrixanterior se puede ampliar para que funcione con polinomios multivariados utilizando el .monoms()método para extraer monomios con coeficientes distintos de cero , así:

from sympy import *
x,y = symbols('x y')
M = Matrix([[1,2],[3,4]])
p = poly( x**3 * y + x * y**2 + y )

def eval_poly_matrix(P,*M):
    res = zeros(*M[0].shape)
    for m in P.monoms():
        term = eye(*M[0].shape)
        for j in enumerate(m):
            i,e = j
            term *= M[i]**e
        res += term
    return res

eval_poly_matrix(p,M,eye(M.rows))

Nota: Es posible realizar algunas comprobaciones de cordura, manejo de casos extremos y optimizaciones:

  1. El número de variables presentes en el polinomio se relaciona con el número de matrices pasadas como parámetros (el primero nunca debe ser mayor que el segundo, y si es menor que lo anterior, es necesario que exista alguna lógica para manejar eso, solo he manejado el caso cuando los dos son iguales)
  2. Todas las matrices deben ser cuadradas según la definición del polinomio de la matriz
  3. Una discusión sobre una versión multivariante de la regla de Horner se presenta en los comentarios de esta pregunta . Esto puede resultar útil para minimizar el número de multiplicaciones de matrices.
  4. Maneje el hecho de que en una matriz el polinomio x*yes diferente y*xporque la multiplicación de matrices no es conmutativa . Aparentemente, las funciones poli en sympy no admiten variables no conmutativas , pero puede definir símbolos con commutative=Falsey parece que hay una forma de avanzar allí.

Acerca del cuarto punto anterior, hay soporte para expresiones Matrix en SymPy, y eso puede ayudar aquí:

from sympy import *
from sympy.matrices import MatrixSymbol

A = Matrix([[1,2],[3,4]])
B = Matrix([[2,3],[3,4]])

X = MatrixSymbol('X',2,2)
Y = MatrixSymbol('Y',2,2)
I = eye(X.rows)

p = X**2 * Y + Y * X ** 2 + X ** 3 - I
display(p)

p = p.subs({X: A, Y: B}).doit()
display(p)

Producción:

Para conocer más novedades sobre esta función, siga el número 18555