代码之家 › 专栏 › 技术社区 › FrankyBravo

矩阵方程求解器Python

equation sympy matrix python

FrankyBravo · 技术社区 · 7 年前

问题是:我试图求解一个二阶矩阵方程:

其中X(要查找)和C(已知)的大小为[nxn]。(n约为1000)。 C是已知的对称协方差矩阵。(X也应该是对称的)

这是我的代码:

from sympy import solve
from sympy import Indexed, IndexedBase, Tuple
import numpy as np

X = IndexedBase('X',shape=(n,n))
eqs = Tuple(np.dot(X,X)-np.dot(C,X)-np.eye(n))
solve(eqs, X)

这样做对吗?我的代码需要很长时间。我正在寻找任何类型的算法,可以帮助我有效地解决这类方程。

2 回复 | 直到 7 年前

lockcmpxchg8b 7 年前

大多数象征性工作都可以手工完成:

X^2 - CX - I = 0
-> X^2 + 2EX - I = 0          // sub C = -2E
-> X^2 + 2EX + E^2 - I = E^2  //add E^2 to both sides (i.e., complete the square)
-> (X + E)^2 = E^2 + I        //simplify and add I to both sides
-> X+E = +/-(E^2 + I)^(1/2)   //take square root (now we may have more than one answer)
-> X = -E +/- (E^2 + I)^(1/2) //subtract E from both sides

矩阵平方根可能是,也可能不是你想要符号化求解的东西。Symphy当然可以让你用符号来表示它,但在我的尝试中(在MinGW64上的Python3中),它已经证明无法用数字计算它。

矩阵C是对称的,因此我们可以检查平方根下的项(即 1/2

每 Wikipedia (Symmetric Matrix) :

两个对称矩阵的和和与差也是对称的。
给定对称矩阵A和B,则AB是对称的当且仅当A和B交换。
每个实对称矩阵都可以通过实正交相似性对角化。

每 Wikipedia (Square Root Of A Matrix)::Explicit Formulas::ByDiagonalization

对于任何可对角化矩阵 A=VDV^(-1) 然后 A^(1/2) = VD^(1/2)V^(-1)

正在从 C 我们问,是吗 E^2+I 可对角化(以便有一个简单的矩阵平方根显式公式)? C 是对称的,并且 E = -(1/2)C ; 标量乘法不会改变 C 因为它影响每个细胞;因此 E 是对称的。 E^2 = (E * E) 通勤,所以 E^2 是对称的。最后 I 是对称的,所以 (E^2 + I) 是对称的。

因此,可以使用通过上述(4)的对角化得到的平方根。对角线矩阵的平方根是通过取对角线上元素的平方根来计算的。在这里,你可能会遇到另一个问题,如果这些因素是否定的,你的答案将是复杂的。对于每个平方根,也可能有多个答案,可能会给您一些需要考虑的答案。这很可能是SymPy未能给出数字答案的原因。

user6655984 user6655984 7 年前

你的代码不对。NumPy用于数值计算,它不会创建表示方程左侧的SymPy对象。这不会帮助你得到解析解。这是一个用辛方法求解矩阵系统的例子;它是2乘2,而不是1000乘1000。

import sympy as sym
X = sym.Matrix(sym.MatrixSymbol('X', 2, 2))
covar = sym.Matrix([[2, 1], [1, 3]])
sym.solve([X**2 - covar*X - sym.eye(2), X-X.T], X)

注意,辛矩阵的乘法就是 * . 第一个方程是你写的,第二个方程要求X是对称的( X.T 是X的转置)。

然而,3乘3的情况已经存在问题,1000乘1000是完全没有希望的。一个人并不是简单地通过把它扔给Symphy来求解一个由500000个非线性方程组成的系统。

你可以尝试使用SciPy的多变量求解器来获得一些数值解,但它将是众多数值解中的一个。矩阵方程的正确方法如下 X**2 - C*X - I = 0 不是把它们扔向电脑;这是做数学。