矩阵方程求解器 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

大部分象征性的工作都可以手工完成：

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 肯定会让您以符号方式表示它，但在我的尝试中（在 MinGW64 上的 Python3 中），它已被证明无法以数字方式计算它。

您的矩阵 C 是对称的，因此我们可以检查平方根下的项（即1/2 的幂）是否具有明确的计算公式。一些初步事实：

每Wikipedia (Symmetric Matrix)：

1. 两个对称矩阵的和和差又是对称的。
2. 给定对称矩阵 A 和 B，当且仅当 A 和 B 可交换时，AB 是对称的。
3. 每个实对称矩阵都可以通过实正交相似性对角化。

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

4. 对于任何可对角化矩阵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 无法给出数字答案的原因。

Answer 2

您的代码不正确。 NumPy 用于数值计算，它不会创建一个代表方程左侧的 SymPy 对象。而且它不会帮助您获得分析解决方案。这是一个使用 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)

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

但是，3 x 3 的情况已经存在问题，1000 x 1000 完全没有希望。将一个包含 500,000 个非线性方程的系统扔给 SymPy 并不能简单地求解它。

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