求解矩阵方程
求解矩阵方程的一些方法 :-)
引言
说到求解矩阵方程，经典的方法总是高斯消元法。但还有许多其他相当有趣的算法可以解决此类矩阵方程。有些（至少在我看来 :-)) 非常优雅，有些则相当复杂，最终它们都做着同样的事情，所以，我产生了疑问：哪种方法最好？这促使我比较不同的算法。为此，我实现了高斯算法、LU分解、Givens变换消元法、使用Cramer法则的逆矩阵以及使用行列式。然后我发现了Householder变换。
好的，我先说一下：在我不知道Householder变换之前，Givens变换是我的最爱。但现在我知道了Householder变换，我的看法改变了。Householder变换确实是一个伟大的东西，它快速且准确。但我们先从算法开始。
高斯消元法
高斯消元法求解矩阵方程的经典方法是消去矩阵主对角线下方的所有元素，并将（例如）一个3x3矩阵方程化为以下形式：
化为以下形式：
现在，可以从最后一个方程求解x₃，然后用x₃求解第二个方程中的x₂，最后用x₂和x₃求解第一个方程中的x₁。这就是高斯消元法 Carl Friedrich Gauss 的基本功能。它包含两个部分。第一部分是消元，第二部分是求解未知x的值。
高斯消元法存在一些缺点：如果元素 aᵢⱼ 大于 1，消元后，下方的行元素会越大，越往下越大。如果它们小于 1，则越往下越小。
一个示例文本方程：
被化为以下形式：
对于大型矩阵，这可能导致不准确，最坏情况是溢出。存在所谓的“枢轴策略”来改善这种情况，其中一种策略是先对所有行进行预处理，使得一行中所有元素的总和接近 1。
通过这种枢轴，相同的示例文本方程在消元后看起来像这样：
这看起来好多了。但高斯消元法仍然不够准确，而且在行中有零的情况下相当脆弱。在这种情况下，如果可能，必须交换行。如果不可能，则必须停止算法，并且方程无法求解。枢轴策略并非总能奏效。我曾在一个多项式插值算法上尝试过。那里的 aᵢⱼ 值分布差异很大，枢轴策略的影响反而变得相当负面。更多内容将在后面讨论。
使用 Givens 变换进行消元。
在 Givens 变换中，使用旋转矩阵来消除元素。该变换的基本思想是将两个元素 aₚ<0xE2><0x82><0x9B> 和 a<0xE2><0x82><0x9B><0xE2><0x82><0x9B> 视为一个二维向量。该向量绕角度 φ 旋转，角度 φ 是该向量在二维空间中的角度的负值。因此，旋转会将向量带入水平方向，并消除其垂直分量，即元素 aₚ<0xE2><0x82><0x9B>。向量的长度不会改变。因此，aᵢⱼ 元素不会像高斯算法那样增大或缩小。这是一个相当大的优势。
在一个 5x5 的矩阵中，旋转矩阵看起来像这样：
其中 p 和 q 是旋转发生的二维平面的两个维度。矩阵 A 与此旋转矩阵的乘法可以消除元素 a₂₂。需要计算的旋转角度为：
为了将矩阵
化为上三角矩阵，需要消除元素 a₂₁、a₃₁ 和 a₃₂，因此需要计算 3 个不同的旋转矩阵，并且 A（以及解向量 Y）必须乘以它们。之后，可以使用高斯算法中相同的过程来求解矩阵方程。
使用此 Givens 变换后，我上面的示例文本方程在消元后看起来像这样：
Givens 变换对于行中的零值来说弱点要少得多，并且它没有元素增大或缩小（在下方各行中）的缺点。这使其相当稳定和准确。这是我的个人最爱 :-)
使用 Prescott Durand Crout 的 LU 分解
LU 分解是一种相当复杂的方法。它的思想是将方程 Ax = b 的矩阵 A 分解为一个下三角矩阵和一个上三角矩阵，得到 LUx = b。矩阵 L 和 U（对于一个 4x4 矩阵）看起来像这样：
现在进行向前替换：
Ly = b
然后进行向后替换：
Ux = y
可以得到解向量 x，它可以在两个简单的循环中计算出来，如下所示：
对于 i ≥ j
以及 l 的公式：
对于 i > j
Crout 的 LU 分解是一个非常巧妙的算法，它在特殊情况下显示出其性能，例如自然样条插值（参见 www.mosismath.com/NaturalSplines/NaturalSplines.html），其中矩阵是三对角矩阵。除此之外，它的表现与 Givens 变换算法非常相似。它执行的乘法比高斯算法少，因此 uᵢⱼ 和 lᵢⱼ 元素不会像高斯算法那样增长或缩小太多。
使用 LU 分解后，我示例文本方程的 U 和 L 矩阵如下：
U =
L =
uᵢⱼ 和 lᵢⱼ 元素不像高斯算法中的 aᵢⱼ 元素那样增长或缩小。但它们比 Givens 变换中的元素增长或缩小得稍多一些。这使得 LU 分解对于大型矩阵来说比 Givens 变换算法的精度稍差。
使用 Gabriel Cramer 的逆矩阵
我第一次接触这种方法是在学习电气工程时，当时我不得不用我的袖珍计算器求解矩阵方程（那已经是很久以前的事了 :-))。如果我们求矩阵 A 的逆，矩阵方程 Ax = b 可以写成 x = bA⁻¹。这意味着用 A 的逆乘以解向量就可以得到 x 向量的结果。根据 Cramer 的说法，逆矩阵可以通过计算矩阵 A 中每个元素的伴随矩阵的行列式，然后除以整个矩阵的行列式，并改变每个元素的符号来计算。对于一个 3x3 矩阵，它看起来像这样：
使用 Gabriel Cramer 的行列式
我是在学习数学课程时接触到通过行列式求解矩阵方程的方法的，我仍然清晰地记得当时我实现的“意大利面条式代码”有多糟糕 :-) 幸运的是，我改进了我的编程技能，现在可以实现一个更好的解决方案。
Cramer 说，要计算矩阵方程 Ax = b 的值 xᵢ，我们必须将矩阵 A 中的第 i 列替换为 b 向量，计算其行列式，然后除以原始矩阵的行列式。对于一个 3x3 矩阵方程，它看起来像这样：
分母部分是 A 的行列式：
使用 Cramer 的方法求解矩阵方程对于折磨学生来说可能是一个很好的方式（至少对我们来说是这样 :-)），但对于严肃的应用，它们并不真正适用。任何使用行列式的方法都相当慢，并且很快就会达到其极限。所以我没有进一步研究这些算法。但它们仍然值得提及，因为它们确实是经典的（并且仍然有趣地研究）。
使用 Householder 变换
我认为 Householder 变换是数学上的又一杰作。它使用一个变换矩阵，通过一次矩阵乘法就可以消除一行中所有需要消除的元素。这意味着将一个 n*n 矩阵转换为一个右上三角矩阵只需要 n-1 次矩阵乘法。这使得该算法真正快速且准确。
数学思想是将一行视为一个向量 x，并旋转该向量，使得需要消除的所有坐标变为 0，同时不改变其长度 x'。因此，在原始向量和由此产生的旋转向量之间放置一个虚拟平面（他们称之为超平面），该平面用作镜像轴，原始向量在该轴上进行镜像。
这个平面（蓝色区域）的参考是向量 ω。向量 ω 是通过减去结果向量与原始向量的差值，然后将结果除以差值的长度（因为 ω 的长度为 1）来计算的。变换矩阵 U 构建如下：
U = I – 2 * ω * ωᵀ
向量 ω 的计算方法是减去结果向量与原始向量的差值，并将结果除以这个差值的长度，因为 ω 的长度为 1。
现在听起来相当简单。但到达这一步是相当努力的过程 J
有关数学的详细描述，请参阅 www.mosismath.com/Householder/Householder.html。
有关数学的详细描述。
不同算法的比较
为了比较这些算法，我尝试了几种不同的方法。对于我小的 4x4 矩阵方程，没有任何区别。所以我尝试了一个 10x10 的矩阵方程。在我的 64 位 Intel Core 3 计算机上，不同算法之间只存在计算时间上的差异，主要是使用行列式的方法与其他方法之间的差异。关于准确性，它们都在双精度值的最低有效数字范围内。所以我放弃了行列式，继续进行多项式插值。现在情况变得更有趣了。
在多项式插值中，对于 n 个支撑点 (xᵢ, yᵢ) 并寻找插值参数 aᵢ，我们得到一个矩阵方程，如下所示：
如果 x 值大于 1，左上角的元素将是最小的，右下角的元素将是最大的。如果 x 值小于 1，则相反。我认为这是一个很好的算法基准测试 J
现在，这种情况对高斯算法来说并不太好。枢轴无法再起作用，因为它不改变一行中最大元素和最小元素之间的比例。这就是为什么高斯算法很快就达到极限的原因。我用了一个小技巧，将矩阵方程颠倒填充。这很有帮助。但无论如何，当有 12 个支撑点时，它开始插值不准确，当有 14 个点时，游戏结束，只有 Givens 变换和 LU 分解仍然有用。但为了比较所有三种方法，我首先进行了 10 个支撑点的插值。我使用了一个弯曲度不大的数学函数，因为多项式插值不喜欢这样。所以我的测试函数是 y = 10/x²，它会产生一个像这样的图：
X 从 1.0 开始，以 0.5 的步长增加。
关于准确性，比较高斯和 Givens 在构建上三角矩阵后以及 LU 算法分解后的矩阵中的值可能是一个好主意。
高斯算法得到的值在 5.004888 和 3.3026311E-103 之间，比例为：
1.666E103 : 1
这已经说明了为什么该算法很快就达到极限。
Givens 得到 1345167.95889 到 3.214294，比例为：
418496 : 1
LU 分解在 L 矩阵中得到 1.0 到 126，在 U 矩阵中得到 0.001953 到 8268.75，最大比例为：
4233600 : 1
到目前为止，Givens 产生的矩阵中最大和最小元素之间的比例最小，因此精度也最高。好的，10 个点看不出来，甚至 20 个点也看不出来，但如果我们把支撑点的数量增加到 30 个，差异就变得可见了。现在 LU 分解的图形与原始形状的偏差比 Givens 的大得多。
Householder 变换得到 1345167.96 到 -3.16，比例为：
425379 : 1
这比 Givens 变换的比例略大。这应该会导致更大的不准确性，但由于 Householder 变换的计算量是所有比较算法中最小的，它的精度比其他所有算法都要高。它甚至可以插值 40 个支撑点而没有太大偏差，并且它所需的时间是 Givens 变换的一又三分之二多一点。所以它比 Givens 变换好一些，也比其他所有方法都好 :-)
使用 30 个支撑点的 LU 分解创建的图形
使用 30 个支撑点的 Givens 变换创建的图形
使用 40 个支撑点的 Householder 变换创建的图形
Givens 和 Householder 变换可能不像其他算法那样出名，但它们是求解矩阵方程的非常巧妙的方法。在我发现 Householder 变换之前，Givens 变换是我的最爱。但现在我不得不改变我的看法 :-)
有关高斯消元法的详细描述，请参见
www.mosismath.com/Matrix_Gauss/MatrixGauss.html
有关 Givens 变换的详细描述，请参见
www.mosismath.com/RotationMatrix/RotationMatrix.html
有关 LU 分解的详细描述，请参见
www.mosismath.com/Matrix_LU/Martix_LU.html
有关 Cramer 方法的详细描述，请参见
www.mosismath.com/Determinant/Martix_Determinant.html
有关 Cramer 逆矩阵的详细描述，请参见
www.mosismath.com/Matrix_inv/Martix_invers.html
有关 Householder 变换的详细描述，请参见
www.mosismath.com/Householder/Householder.html
使用代码
我的 C# 演示项目包含一个主窗口。算法就实现在里面。我是在 Visual Studio 10 中实现的。希望您喜欢 :-)。