使用NumPy求解线性系统

计算矩阵逆 $A^{-1}$ 是一种理论方法，用于通过计算 $x = A^{-1}b$ 来求解系统 $Ax = b$。然而，在实际计算中，这种方法通常不是最实用或数值最稳定的。求逆运算所需的计算量大于实际需要，并且可能引入更大的数值误差，特别是对于接近奇异（不可逆）的矩阵。

幸运的是，NumPy 提供了一个专门的函数 numpy.linalg.solve，它专门用于直接高效地求解 $Ax = b$ 形式的线性方程组。此函数使用精密的算法（通常基于矩阵分解，如LU分解），这些算法通常比手动计算逆矩阵然后乘以 $b$ 更快，并提供更准确的结果。

使用 numpy.linalg.solve

numpy.linalg.solve 函数接受两个参数 (parameter)：

1. 系数矩阵 A（一个二维NumPy数组）。
2. 常数向量 (vector) b（一个一维或二维NumPy数组）。

它返回解向量 x。

我们考虑一个简单的线性方程组：

$$\begin{cases} 2x_1 + 3x_2 = 8 \\ x_1 + 4x_2 = 9 \end{cases}$$

我们可以将此系统表示为矩阵形式 $Ax = b$，其中：

$$A = \begin{bmatrix} 2 & 3 \\ 1 & 4 \end{bmatrix}, \quad x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, \quad b = \begin{bmatrix} 8 \\ 9 \end{bmatrix}$$

现在，我们使用 NumPy 定义 A 和 b 并求解 x：

import numpy as np

# 定义系数矩阵 A
A = np.array([[2, 3],
              [1, 4]])

# 定义常数向量 b
b = np.array([8, 9])

print("矩阵 A:\n", A)
print("\n向量 b:\n", b)

# 求解系统 Ax = b 中的 x
try:
    x = np.linalg.solve(A, b)
    print(f"\n解向量 x: {x}") # 输出应为 [1. 2.]
except np.linalg.LinAlgError as e:
    print(f"\n无法求解系统: {e}")

输出结果显示解为 $x_1 = 1$ 和 $x_2 = 2$。

验证解

验证解的正确性始终是一个好的习惯。我们可以通过使用计算出的解向量 (vector) $x$ 来计算 $Ax$，并检查结果是否确实等于（或非常接近，由于潜在的浮点数不精确性）原始向量 $b$。

# 通过计算 A @ x 来验证解
verification = A @ x # 等同于 np.dot(A, x)

print(f"\n验证 (A @ x): {verification}")

# 检查 A @ x 是否接近 b
# np.allclose 对于比较浮点数很有用
are_close = np.allclose(verification, b)
print(f"\n计算出的 Ax 是否接近原始 b？ {are_close}")

A @ x 的输出应为 [8. 9.]，与我们的原始向量 b 匹配，并且 np.allclose 应返回 True。

为何偏好 solve 而非 inv？

我们简要比较一下使用 np.linalg.solve 求解系统与先计算逆矩阵再进行乘法的两种方式：

# 方法 1：使用 np.linalg.solve（推荐）
x_solve = np.linalg.solve(A, b)
print(f"\n使用 solve 求解: {x_solve}")

# 方法 2：使用 np.linalg.inv
try:
    A_inv = np.linalg.inv(A)
    x_inv = A_inv @ b # 计算 x = A⁻¹b
    print(f"使用逆矩阵求解: {x_inv}")

    # 检查解是否接近
    print(f"solve 和 inv 的解是否接近？ {np.allclose(x_solve, x_inv)}")

except np.linalg.LinAlgError:
    print("\n矩阵 A 是奇异的，无法计算逆矩阵来求解。")

虽然对于这个简单、性质良好的系统，两种方法都产生相同的结果，但 np.linalg.solve 通常在几个方面更优：

1. 效率： 直接求解系统通常比求逆矩阵再执行矩阵乘法计算成本更低。
2. 数值稳定性： solve 中使用的直接方法通常旨在最大限度地减少浮点误差的影响，从而获得更准确的结果，特别是对于更大或更复杂的系统。计算逆矩阵有时会放大误差。

因此，除非您因其他目的明确需要逆矩阵 $A^{-1}$，否则应始终优先使用 np.linalg.solve(A, b) 来找到系统 $Ax = b$ 的解 $x$。

处理不可逆矩阵

如果矩阵 $A$ 是奇异的（不可逆）会怎样？在这种情况下，系统 $Ax = b$ 可能没有唯一解（它可能有无限解或根本没有解）。如前所述，奇异矩阵没有逆矩阵。

如果您尝试将 np.linalg.solve 与奇异矩阵 A 一起使用，NumPy 将会引发 LinAlgError，表明该矩阵是奇异的，并且此方法无法唯一地求解该系统。

# 奇异矩阵示例
A_singular = np.array([[1, 2],
                       [2, 4]]) # 第 2 行是第 1 行的 2 倍

b_test = np.array([3, 6]) # 此 b 可能导致无限解

print("\n尝试使用奇异矩阵求解:")
try:
    x_singular = np.linalg.solve(A_singular, b_test)
    print(f"解（不应到达此处）: {x_singular}")
except np.linalg.LinAlgError as e:
    print(f"捕获到预期错误: {e}") # 输出: Caught expected error: Singular matrix

这种行为很有帮助，因为它立即表明系统可解性存在问题，这归因于系数矩阵 $A$ 的性质。

总而言之，NumPy 的 linalg.solve 函数提供了一种高效的方法来求解线性方程组，有效地处理底层计算，并在由于奇异矩阵而无法确定唯一解时提醒您。它是 Python 数值计算中这项常见任务的标准工具。