LU分解的数值精度与复数64.

Question

我正在寻找最准确的类似C的实现，用于求解复浮点类型(又名np.complex64(的线性方程组。要进行LU分解。我想首先在 numpy 中实现它，看到一切顺利，然后将实现转换为 C。

这是我目前得到的：

import scipy.linalg as la
import numpy as np

def lu_factor(A):
L, U = np.empty_like(A), np.empty_like(A)
n = A.shape[0]
for k in range(n):
L[k, k] = 1
U[k, k] = A[k, k] - L[k, :k] @ U[:k, k]
for j in range(k + 1, n):
U[k, j] = A[k, j] - L[k, :k] @ U[:k, j]
for i in range(k + 1, n):
L[i, k] = (A[i, k] - L[i, :k] @ U[:k, k]) / U[k, k]
return L, U

def forward_sub(L, b):
x = np.empty_like(b)
for i in range(b.size):
x[i] = (b[i] - L[i, :i] @ x[:i]) / L[i, i]
return x

def backward_sub(U, b):
x = np.empty_like(b)
for i in reversed(range(b.size)):
x[i] = (b[i] - U[i, i + 1:] @ x[i + 1:]) / U[i, i]
return x

def lu_solve(A, b):
L, U = lu_factor(A)
return backward_sub(U, forward_sub(L, b))

b = np.array([6 + 1j, -4 + 2j, 27 + 3j], dtype=np.complex128)
A = np.array([
[1 + 4j, 1 + 5j,  1 + 6j],
[0 + 7j, 2 + 8j,  5 + 9j],
[2 + 1j, 5 + 2j, -1 + 3j],
], dtype=np.complex128)

x_expected = la.lu_solve(la.lu_factor(A), b)
x = lu_solve(A, b)
np.testing.assert_allclose(x_expected, x)

A和b是随机选择的。注意 - 它们的类型是np.complex128(复杂双精度(，并且"朴素"实现的结果足够接近 scipy 实现。

将类型更改为np.complex64(复浮点数(时，我们得到：

Mismatch: 66.7%
Max absolute difference: 1.0612305e-06
Max relative difference: 6.8692873e-07
x: array([ 1.387071-0.680237j,  3.673277+1.09019j , -3.683192-1.225474j], dtype=complex64)
y: array([ 1.387072-0.680236j,  3.673277+1.090189j, -3.683192-1.225474j], dtype=complex64)

e-06对我来说似乎很高。除了"双倍到浮动精度"之外，还有很好的解释吗？numpy/BLAS 到底做了什么而我没有？这可以使用基本操作进行复制吗？

注意：性能对我来说不是问题，我只关心精度。

Answer 1

编写一个健壮而准确的数值线性代数过程可能非常棘手。 我的第一个建议是使用像Lapack这样的成熟库。

如果要编写自己的 LU，则绝对必须至少使用"部分透视"来获得可靠的过程，否则：

如果矩阵中没有正确的排序或排列，则 因子分解可能无法实现。

注意：wikepedia建议了一个当然可以是一个很好的起点的实现。您只需要修改它以支持复数而不是双精度数。

一旦你解决了你的系统，就可以通过执行一两个迭代细化步骤来提高精度。

Compute the residual: r = b − A.x
Solve the system:     A.d = r
Add the correction:   x = x + d

只需使用 LU 代码应用前面的过程即可求解 A.d=r 系统。如果您正在寻找高精度，请不要忽视最后一步，它在实践中通常效果很好。