python - Python中正向替换的数值稳定性

我正在用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

我正在使用numpy和64位浮点。
简易测试工具
我建立了一个简单的测试套件，生成系数矩阵和向量，计算x，然后使用正向或反向替换来恢复b，并将其与已知的有效性值进行比较。
以下代码执行这些检查：

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()

注意，测试矩阵的最大大小是49×49。即使在这种情况下，系统也不能总是计算出合适的解，并且失败率超过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，这是正常的行为吗？因为问题本质上没有结构。
对于各种实际应用，解100个线性方程组怎么样？这些限制仅仅是一个公认的事实吗？例如，优化算法对这些问题自然是健壮的？还是我遗漏了这个问题的一些重要方面？
[0]：出版社，威廉H.数字食谱第三版：科学计算的艺术。剑桥大学出版社，2007年。

最佳答案

没有限制。这是一个非常有成效的练习，我们都意识到了；编写线性解算器并不是那么容易，这就是为什么几乎所有的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);

请注意，由于您是从均匀分布中提取数字，因此当行有更多的可能出现秩亏时，它越来越可能击中病态矩阵（wrt到inversion）。这就是错误变得越来越大的原因。将某个单位矩阵乘以一个标量，所有错误都会返回到eps*n级别。
但最好的办法是，把这个留给经过几十年考验的专家算法写这些真的没那么简单。您可以阅读Fortran代码，例如，dtrsm求解三角形系统。
在python方面，您可以使用scipy.linalg.solve_triangular，它使用lapack中的?trtrs例程。