矩阵方程求解器 Python

Question

这就是问题所在：我正在尝试求解以下形式的二阶矩阵方程：

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

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

Answer 1

您的代码不正确。 NumPy 用于数值计算，它不会创建表示等式左侧的 SymPy 对象。而且它不会帮助您获得分析解决方案。这是使用 SymPy 求解矩阵系统的示例；它是 2 x 2 而不是 1000 x 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)

请注意，SymPy 矩阵的乘法只是 *。第一个方程是你写的，第二个方程要求 X 是对称的（X.T 是 X 的转置）。

但是，3×3的情况已经有问题了，1000×1000就更没希望了。将 500,000 个非线性方程组扔到 SymPy 中并不能简单地解决它。

您可以尝试 SciPy 的多变量求解器来获得一些数值解，但它只会是许多数值解中的一个。像 X**2 - C*X - I = 0 这样的矩阵方程的正确方法不是把它们扔到计算机上；这是做数学。

Answer 2

大部分符号工作都可以手工完成：

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

矩阵平方根可能是也可能不是您想要符号化求解的问题。 SymPy 肯定会让你用符号表示它，但它已被证明无法用数字计算它，到目前为止在我的尝试中（在 Python3 on MinGW64）。

你的矩阵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 无法给出数字答案的原因。

矩阵方程求解器 Python

Matrix Equation solver Python

python

equation

matrix

sympy