2x2矩阵的数值稳定逆矩阵

Question

我正在用C语言开发一个数值求解器，需要对一个2x2的矩阵进行求逆，然后再与另一个矩阵相乘：

C = B . inv(A)

我一直在使用以下的2x2矩阵求逆的定义：

a = A[0][0];
b = A[0][1];
c = A[1][0];
d = A[1][1];
invA[0][0] = d/(a*d-b*c);
invA[0][1] = -b/(a*d-b*c);
invA[1][0] = -c/(a*d-b*c);
invA[1][1] = a/(a*d-b*c);

在我求解器的前几次迭代中，这个方法似乎能给出正确的结果，但过了一段时间后，数值开始变得很大，最终导致崩溃。

然后，我对比了一个使用SciPy的实现，发现同样的数学计算没有崩溃。我能找到的唯一不同之处是，SciPy的代码使用了scipy.linalg.inv()，它内部调用了LAPACK来进行矩阵求逆。

当我把inv()的调用替换成我上面的计算时，Python版本就崩溃了，所以我很确定这就是问题所在。计算中的小差异逐渐显现出来，这让我觉得这是一个数值问题——对于求逆操作来说，这并不意外。

我使用的是双精度浮点数（64位），希望数值问题不会成为障碍，但显然并不是这样。

不过，我希望在我的C代码中解决这个问题，而不需要调用像LAPACK这样的库，因为我将其移植到纯C的原因就是为了在目标系统上运行。此外，我还想理解这个问题，而不仅仅是调用一个黑箱。最终，如果可能的话，我也希望能用单精度运行。

所以，我的问题是，对于这样一个小矩阵，有没有更稳定的方式来计算A的逆？

谢谢。

编辑：目前我在尝试找出是否可以通过求解C来< a href="http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/" rel="nofollow">避免求逆。

Answer 1

你的代码没问题；不过，它在进行四次减法时可能会出现精度损失的风险。

可以考虑使用一些更高级的技术，比如matfunc.py中使用的方法。那段代码通过QR分解来进行反演，QR分解是用豪斯霍尔德变换实现的。最终的结果还通过迭代精化进一步改善。

Answer 2

计算行列式并不是很稳定。更好的方法是使用带有部分选主元的高斯-约旦消元法，这样你可以更容易地进行计算。

解决一个2x2的方程组

我们来解决这个方程组（先用 c, f = 1, 0，然后用 c, f = 0, 1 来得到逆矩阵）

a * x + b * y = c
d * x + e * y = f

用伪代码来说就是这样

if a == 0 and d == 0 then "singular"

if abs(a) >= abs(d):
    alpha <- d / a
    beta <- e - b * alpha
    if beta == 0 then "singular"
    gamma <- f - c * alpha
    y <- gamma / beta
    x <- (c - b * y) / a
else
    swap((a, b, c), (d, e, f))
    restart

这种方法比用行列式加余子式更稳定（beta 是行列式乘以某个常数，以稳定的方式计算出来）。你也可以算出完全选主元的等效方法（比如可能需要交换 x 和 y，这样第一次除以 a 时，a 是 a, b, d, e 中绝对值最大的数），在某些情况下这可能更稳定，但我发现上面的方法效果很好。

这相当于进行 LU 分解（如果你想保存这个 LU 分解，可以存储 gamma, beta, a, b, c）。

计算 QR 分解也可以明确进行（只要你正确操作，它也是非常稳定的），但速度较慢（还涉及到开平方）。选择权在你。

提高准确性

如果你需要更好的准确性（上面的方法是稳定的，但会有一些舍入误差，和特征值的比值成正比），你可以“求解修正值”。

实际上，假设你用上面的方法解出了 A * x = b 的 x。然后你计算 A * x，发现它并不完全等于 b，有一点小误差：

A * x - b = db

现在，如果你在 A * dx = db 中求解 dx，你会得到

A * (x - dx) = b + db - db - ddb = b - ddb

其中 ddb 是通过数值解 A * dx = db 引入的误差，通常比 db 小得多（因为 db 比 b 小得多）。

你可以重复上述过程，但通常只需一步就能恢复完整的机器精度。

Answer 3

别去求矩阵的逆。几乎总是，你用逆矩阵来做的事情，其实可以更快、更准确地完成，而不需要去求逆。矩阵求逆本身就不太稳定，再加上浮点数，就更容易出问题了。

说 C = B . inv(A) 就等于说你想解方程 AC = B 来找出 C。你可以把每个 B 和 C 拆成两列来解决这个问题。通过解 A C1 = B1 和 A C2 = B2，就能得到 C。