Python 中前向替换的数值稳定性

Question

我正在用 Python 实现一些基本的线性方程求解器。

我目前已经实现了三角方程组的前向和后向替代（因此求解起来非常简单！），但即使对于大约 50 个方程组（50x50 系数矩阵），解的精度也会变得非常差。

以下代码执行前向/后向替换：

FORWARD_SUBSTITUTION = 1
BACKWARD_SUBSTITUTION = 2
def solve_triang_subst(A: np.ndarray, b: np.ndarray,
                       substitution=FORWARD_SUBSTITUTION) -> np.ndarray:
    """Solves a triangular system via
    forward or backward substitution.

    A must be triangular. FORWARD_SUBSTITUTION means A should be
    lower-triangular, BACKWARD_SUBSTITUTION means A should be upper-triangular.
    """
    rows = len(A)
    x = np.zeros(rows, dtype=A.dtype)
    row_sequence = reversed(range(rows)) if substitution == BACKWARD_SUBSTITUTION else range(rows)

    for row in row_sequence:
        delta = b[row] - np.dot(A[row], x)
        cur_x = delta / A[row][row]
        x[row] = cur_x

    return x

我正在使用numpy64 位浮点数。

简单的测试工具

我已经建立了一个简单的测试套件，它生成系数矩阵和x向量，计算b，然后使用前向或后向替换来恢复x，将其与已知的有效性值进行比较。

以下代码执行这些检查：

import numpy as np
import scipy.linalg as sp_la

RANDOM_SEED = 1984
np.random.seed(RANDOM_SEED)


def check(sol: np.ndarray, x_gt: np.ndarray, description: str) -> None:
    if not np.allclose(sol, x_gt, rtol=0.1):
        print("Found inaccurate solution:")
        print(sol)
        print("Ground truth (not achieved...):")
        print(x_gt)
        raise ValueError("{} did not work!".format(description))

def fuzz_test_solving():
    N_ITERATIONS = 100
    refine_result = True
    for mode in [FORWARD_SUBSTITUTION, BACKWARD_SUBSTITUTION]:
        print("Starting mode {}".format(mode))
        for iteration in range(N_ITERATIONS):
            N = np.random.randint(3, 50)
            A = np.random.uniform(0.0, 1.0, [N, N]).astype(np.float64)

            if mode == BACKWARD_SUBSTITUTION:
                A = np.triu(A)
            elif mode == FORWARD_SUBSTITUTION:
                A = np.tril(A)
            else:
                raise ValueError()

            x_gt = np.random.uniform(0.0, 1.0, N).astype(np.float64)
            b = np.dot(A, x_gt)

            x_est = solve_triang_subst(A, b, substitution=mode,
                                       refine_result=refine_result)
            # TODO report error and count, don't throw!
            # Keep track of error norm!!
            check(x_est, x_gt,
                  "Mode {} custom triang iteration {}".format(mode, iteration))

if __name__ == '__main__':
    fuzz_test_solving()

请注意，测试矩阵的最大尺寸为 49x49。即使在这种情况下，系统也无法始终计算出合适的解决方案，并且失败率会超过 0.1。下面是这样一个失败的例子（这是在做向后代入，所以最大的误差在第0个系数；所有的测试数据都是从[0, 1[]中统一采样的：

Solution found with Mode 2 custom triang iteration 24:
[ 0.27876067  0.55200497  0.49499509  0.3259397   0.62420183  0.47041149
  0.63557676  0.41155446  0.47191956  0.74385864  0.03002819  0.4700286
  0.37989592  0.56527691  0.15072607  0.05659282  0.52587574  0.82252197
  0.65662833  0.50250729  0.74139748  0.10852731  0.27864265  0.42981232
  0.16327331  0.74097937  0.24411709  0.96934199  0.890266    0.9183985
  0.14842446  0.51806495  0.36966843  0.18227989  0.85399593  0.89615663
  0.39819336  0.90445931  0.21430972  0.61212349  0.85205597  0.66758689
  0.1793689   0.38067267  0.39104614  0.6765885   0.4118123 ]
Ground truth (not achieved...)
[ 0.20881608  0.71009766  0.44735271  0.31169033  0.63982328  0.49075813
  0.59669585  0.43844108  0.47764942  0.72222069  0.03497499  0.4707452
  0.37679884  0.56439738  0.15120397  0.05635977  0.52616387  0.82230625
  0.65670245  0.50251426  0.74139956  0.10845974  0.27864289  0.42981226
  0.1632732   0.74097939  0.24411707  0.96934199  0.89026601  0.91839849
  0.14842446  0.51806495  0.36966843  0.18227989  0.85399593  0.89615663
  0.39819336  0.90445931  0.21430972  0.61212349  0.85205597  0.66758689
  0.1793689   0.38067267  0.39104614  0.6765885   0.4118123 ]

我还实现了[0]第2.5节中描述的迭代细化方法，虽然它确实有一点帮助，但对于较大的矩阵来说，结果仍然很差。

MATLAB 健全性检查

我也在MATLAB中做了这个实验，即使在那里，一旦方程超过100个，估计误差就会呈指数级增长。

以下是我用于此实验的 MATLAB 代码：

err_norms = [];
range = 1:3:120;
for size=range

  A = rand(size, size);
  A = tril(A);
  x_gt = rand(size, 1);

  b = A * x_gt;
  x_sol = A\b;

  err_norms = [err_norms, norm(x_gt - x_sol)];
end

plot(range, err_norms);
set(gca, 'YScale', 'log')

这是结果图：

主要问题

我的问题是：鉴于我随机生成 A 矩阵和 x，问题本质上没有结构，这是正常行为吗？

为各种实际应用求解由数百个方程组成的线性系统怎么样？这些限制是否只是一个公认的事实，例如，优化算法对这些问题自然具有鲁棒性？或者我错过了这个问题的一些重要方面？

---

[0]：出版社，William H.《数值食谱》第 3 版：科学计算的艺术。剑桥大学出版社，2007。

Answer 1

没有任何限制。我们都认识到这是一项卓有成效的做法；编写线性求解器并不那么容易，这就是为什么 LAPACK 或其在其他语言中的同类产品几乎总是被完全放心地使用。

您遇到了几乎奇异的矩阵，并且因为您使用的是 matlab 的反斜杠，所以当遇到接近奇点时，您看不到 matlab 在幕后切换到最小二乘解。如果您只是更改A\b为linsolve(A,b)因此限制求解器求解平方系统，您可能会在控制台上看到很多警告。

我没有测试它，因为我不再有许可证，但如果我盲目地写，这应该会向您显示每一步矩阵的条件数。

err_norms = [];
range = 1:3:120;
for i=1:40
  size = range(i);
  A = rand(size, size);
  A = tril(A);
  x_gt = rand(size, 1);

  b = A * x_gt;
  x_sol = linsolve(A,b);

  err_norms = [err_norms, norm(x_gt - x_sol)];
  zzz(i) = rcond(A);
end

semilogy(range, err_norms);
figure,semilogy(range,zzz);

请注意，因为您是从均匀分布中选取数字，所以越来越有可能遇到病态矩阵（即反转），因为行更有可能出现排名不足。这就是为什么误差会越来越大。撒上一些单位矩阵乘以标量，所有错误都应该回到原来的eps*n水平。

但最好将其留给经过数十年测试的专家算法。编写这些内容确实不是那么简单。您可以阅读 Fortran 代码，例如dtrsm求解三角系统。

在 Python 方面，您可以使用scipy.linalg.solve_triangular来自?trtrsLAPACK 的例程。