Python 中正向替换的数值稳定性
Numerical Stability of Forward Substitution in Python
我正在 Python.
中实现一些基本的线性方程求解器
我目前已经实现了三角方程组的前向和后向代换（求解起来非常简单！），但是即使有大约 50 个方程组（50x50 系数矩阵），解的精度也变得非常差。
以下代码执行 forward/backward 替换：
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，然后使用正向或反向替换来恢复 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，这是否是正常行为，因为问题中基本上没有结构？
对于各种实际应用，如何求解包含 100 多个方程的线性系统？这些限制是否只是一个公认的事实，例如，优化算法对这些问题具有天然的鲁棒性？还是我遗漏了这个问题的一些重要方面？
[0]：出版社，William H. Numerical recipes 第 3 版：科学计算的艺术。剑桥大学出版社, 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 to inversion），因为行更有可能出现秩不足。这就是错误越来越大的原因。将一些单位矩阵乘以一个标量，所有错误都应该回到 eps*n 水平。
但最好将其留给经过数十年测试的专家算法。写这些真的不是那么简单。您可以阅读 Fortran 代码，例如 dtrsm 求解三角系统。
在 Python 方面，您可以使用 scipy.linalg.solve_triangular，它使用来自 LAPACK 的 ?trtrs 例程。
我正在 Python.
中实现一些基本的线性方程求解器我目前已经实现了三角方程组的前向和后向代换（求解起来非常简单！），但是即使有大约 50 个方程组（50x50 系数矩阵），解的精度也变得非常差。
以下代码执行 forward/backward 替换：
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，然后使用正向或反向替换来恢复 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，这是否是正常行为，因为问题中基本上没有结构？
对于各种实际应用，如何求解包含 100 多个方程的线性系统？这些限制是否只是一个公认的事实，例如，优化算法对这些问题具有天然的鲁棒性？还是我遗漏了这个问题的一些重要方面？
[0]：出版社，William H. Numerical recipes 第 3 版：科学计算的艺术。剑桥大学出版社, 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 to inversion），因为行更有可能出现秩不足。这就是错误越来越大的原因。将一些单位矩阵乘以一个标量，所有错误都应该回到 eps*n 水平。
但最好将其留给经过数十年测试的专家算法。写这些真的不是那么简单。您可以阅读 Fortran 代码，例如 dtrsm 求解三角系统。
在 Python 方面，您可以使用 scipy.linalg.solve_triangular，它使用来自 LAPACK 的 ?trtrs 例程。