多項式で行列を代入するときの問題

Jan 12 2021

例:let

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

エラーが発生します:

TypeError:<class'sympy.matrices.immutable.ImmutableDenseMatrix '>と<class'sympy.core.numbers.One'>を追加できません

最初の2つの項は行列になりますが、最後の項(定数項)は行列ではなくスカラーであるため、これは非常に妥当です。この状況を改善するために、多項式を次のように変更しました。

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

同じエラーが続く。xをMに置き換えて、式を手動で入力する必要がありますか?この例では、多項式には3つの項しかありませんが、実際には数十の項に遭遇します(多変量多項式)。

回答

1 wsdookadr Jan 12 2021 at 12:40

したがって、質問は主に行列多項式の概念を中心に展開していると思います。

(ここで、Pは多項式、Aは行列です)

これは、自由項が数値であり、行列である残りの項と加算できないことを意味していると思います。事実上、これら2つのタイプの間で加算演算は定義されていません。

TypeError:<class'sympy.matrices.immutable.ImmutableDenseMatrix '>と<class'sympy.core.numbers.One'>を追加できません

ただし、これは、特定の行列の行列多項式を評価する関数を定義することで回避できます。ここでの違いは、私たちが使用していることである行列の累乗我々は正しく行列多項式の自由項を計算して、ある単位行列に必要な形状のを:a_0 * II=A^0

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)

出力:

この例では、多項式には3つの項しかありませんが、実際には数十の項に遭遇します(多変量多項式)。

上記の関数は、次のように、係数がゼロ以外の単項式を抽出eval_poly_matrixする.monoms()方法を使用して、多変量多項式で機能するように拡張できます。

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))

注:一部の健全性チェック、エッジケースの処理、および最適化が可能です。

  1. 多項式に存在する変数の数は、パラメーターとして渡される行列の数に関連します(前者は後者より大きくなることはありません。それを処理するために存在する必要があるロジックよりも少ない場合は、ケースのみを処理しました。 2つが等しい場合)
  2. 行列多項式の定義に従って、すべての行列は正方形である必要があります
  3. 多変量バージョンについての議論ホーナーのルールがでていますこの質問のコメント。これは、行列の乗算の数を最小限に抑えるのに役立つ場合があります。
  4. 行列の多項式x*yy*x、行列の乗算が非可換であるためとは異なるという事実を処理します。どうやらsympyのポリ関数は非可換変数をサポートしていない ようですが、でシンボルを定義することができ、そこに前進する方法があるようですcommutative=False

上記の4番目のポイントについては、SymPyでMatrix式がサポートされており、これがここで役立ちます。

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)

出力:

この機能の詳細については、#18555をフォローしてください。