用NumPy解决带约束线性方程组

在科学计算和工程领域，经常遇到需要求解带约束的线性方程组问题。本文介绍了一种利用 NumPy 库高效解决具有线性等式约束的线性方程组 AX=b 的方法，该方法尤其适用于同时满足原始方程和所有线性约束的场景。通过将原始方程组与线性约束方程合并，构建增广系统，并巧妙地运用 np.linalg.lstsq 函数求解，可以获得精确或最佳的最小二乘解。这种方法避免了将约束作为优化惩罚项的常见误区，提供了一种更直接、更高效的解决方案，旨在帮助读者更好地解决实际问题，提升数据分析和模型构建的效率。同时，文章还探讨了该方法在实际应用中的注意事项和最佳实践，确保读者能够充分理解并正确运用该技术。

本文介绍如何利用 NumPy 库高效解决具有线性等式约束的线性方程组 AX=b。通过将原始方程组与线性约束方程合并，形成一个增广系统，然后使用 np.linalg.lstsq 函数求解，可以同时满足原始方程和所有线性约束，获得精确或最佳的最小二乘解。

1. 引言：带约束的线性系统求解挑战

线性方程组 AX=b 在科学计算、工程、统计学等领域无处不在。然而，在许多实际应用中，解向量 X 往往需要满足额外的条件，即线性等式约束。例如，某些变量的和必须为零，或者某些变量之间存在固定的比例关系。

当 X 向量受到这些线性等式约束时，直接求解 AX=b 变得复杂。常见的误区是尝试将约束作为优化问题的惩罚项或使用通用非线性优化器（如 scipy.optimize.minimize），但这可能导致解在满足约束的同时，无法精确满足原始方程 AX=b。本文将介绍一种更直接、更高效的方法，通过构建增广系统并利用 NumPy 的 np.linalg.lstsq 函数来解决这类问题。

2. 线性等式约束的数学表达

假设我们有一个原始的线性方程组 A X = b，其中 A 是一个 m x n 的系数矩阵，X 是一个 n x 1 的未知向量，b 是一个 m x 1 的常数向量。

同时，我们有 k 个线性等式约束，它们可以被表示为 C X = d，其中 C 是一个 k x n 的约束矩阵，d 是一个 k x 1 的常数向量。

示例： 假设 X = [x1, y1, x2, y2, x3, y3, x4, y4] 是一个 8x1 的向量。 给定以下三个约束：

1. 0.5 * (y1 + y2) = 0
2. 0.5 * (x3 + x4) = 0
3. 0.5 * (y3 + y4) = 0

我们可以将这些约束转化为矩阵 C 和向量 d 的形式。

• 对于约束 1：0.5 * x[1] + 0.5 * x[3] = 0 (注意，这里 x[1] 和 x[3] 指的是 X 向量中的第二个和第四个元素，即 y1 和 y2)。 对应的 C 行将是 [0, 0.5, 0, 0.5, 0, 0, 0, 0]，d 的对应元素是 0。
• 对于约束 2：0.5 * x[4] + 0.5 * x[6] = 0 (对应 x3 和 x4)。 对应的 C 行将是 [0, 0, 0, 0, 0.5, 0, 0.5, 0]，d 的对应元素是 0。
• 对于约束 3：0.5 * x[5] + 0.5 * x[7] = 0 (对应 y3 和 y4)。 对应的 C 行将是 [0, 0, 0, 0, 0, 0.5, 0, 0.5]，d 的对应元素是 0。

因此，我们可以构建约束矩阵 AC (对应 C) 和约束向量 bC (对应 d)：

import numpy as np

# 假设 A 和 b 已定义
A = np.array([
    [-261.60,  11.26,      0.0,    0.0,     0.0,    0.0,     0.0,    0.0],
    [   4.07, -12.75,      0.0,    0.0,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,  -158.63,  -5.65,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,    -2.81, -12.14,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0, -265.99,  19.29,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0,   12.59, -12.34,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0,     0.0,    0.0, -166.25, -12.63],
    [    0.0,    0.0,      0.0,    0.0,     0.0,    0.0,   -8.40, -11.14]
])

b = np.array([
     -6.95,
     16.35,
     -0.96,
     16.35,
     19.19,
    -15.85,
    -12.36,
    -15.63]).reshape(-1, 1)

# 构建约束矩阵 AC 和约束向量 bC
AC = np.zeros([3, A.shape[1]]) # 3个约束，X有8个变量
bC = np.zeros((3, 1))

# 0.5 * (y1 + y2) = 0  => x[1] 和 x[3]
AC[0, [1, 3]] = 0.5 
# 0.5 * (x3 + x4) = 0  => x[4] 和 x[6]
AC[1, [4, 6]] = 0.5 
# 0.5 * (y3 + y4) = 0  => x[5] 和 x[7]
AC[2, [5, 7]] = 0.5 

print("约束矩阵 AC:\n", AC)
print("约束向量 bC:\n", bC)

3. 构建增广系统

为了同时解决原始方程组和所有线性等式约束，我们可以将它们合并成一个更大的、增广的线性系统。这个增广系统可以表示为：

$$ \begin{bmatrix} A \ C \end{bmatrix} X = \begin{bmatrix} b \ d \end{bmatrix} $$

或者简化为 A_aug X = b_aug，其中：

• A_aug 是通过垂直堆叠 A 和 C 得到的矩阵。
• b_aug 是通过垂直堆叠 b 和 d 得到的向量。

任何满足 A_aug X = b_aug 的向量 X 都将同时满足原始方程 A X = b 和所有线性等式约束 C X = d。

在 NumPy 中，可以使用 np.vstack 函数来实现矩阵和向量的垂直堆叠：

# 堆叠原始矩阵 A 和约束矩阵 AC
A_augmented = np.vstack([A, AC])

# 堆叠原始向量 b 和约束向量 bC
b_augmented = np.vstack([b, bC])

print("\n增广矩阵 A_augmented 的形状:", A_augmented.shape)
print("增广向量 b_augmented 的形状:", b_augmented.shape)

此时，A_augmented 的形状将是 (m+k) x n，b_augmented 的形状将是 (m+k) x 1。

4. 使用 NumPy 的 lstsq 求解

np.linalg.lstsq 函数是 NumPy 库中用于求解线性最小二乘问题的核心工具。它寻找一个向量 X，使得 ||A_aug X - b_aug||^2 最小。

• 如果增广系统 A_aug X = b_aug 是精确可解的（即存在唯一解或无穷多解），lstsq 将找到一个精确解。
• 如果系统是超定的（方程数多于未知数）但一致的，它也能找到精确解。
• 如果系统是不一致的（无精确解），lstsq 将找到一个“最佳”近似解，即残差平方和最小的解。

对于我们构建的增广系统，lstsq 将直接找到一个 X，它在满足所有线性等式约束的同时，也尽可能地满足原始方程 A X = b。如果原始系统与约束本身是兼容的，它将找到一个精确解。

代码示例：

# 使用 np.linalg.lstsq 求解增广系统
x_solution, residuals, rank, singular_values = np.linalg.lstsq(A_augmented, b_augmented, rcond=None)

print("\n求解得到的 X 向量:\n", x_solution)

rcond=None 参数是推荐的用法，它使用机器精度来确定奇异值的阈值，而不是默认的固定值，这有助于提高数值稳定性。

5. 结果验证与分析

求解得到 x_solution 后，我们应该验证它是否同时满足原始方程和所有约束。

验证原始方程 A X = b：

# 检查是否满足原始方程 A X = b
b_predicted = np.matmul(A, x_solution)
print("\n原始方程左侧 (A * X_solution):\n", b_predicted)
print("原始方程右侧 (b):\n", b)

# 计算残差
original_equation_residuals = b_predicted - b
print("\n原始方程残差:\n", original_equation_residuals)
print("原始方程残差的L2范数:", np.linalg.norm(original_equation_residuals))

验证线性等式约束 C X = d：

# 检查是否满足约束 C X = d
constraints_satisfied = np.matmul(AC, x_solution)
print("\n约束左侧 (AC * X_solution):\n", constraints_satisfied)
print("约束右侧 (bC):\n", bC)

# 计算约束残差
constraint_residuals = constraints_satisfied - bC
print("\n约束残差:\n", constraint_residuals)
print("约束残差的L2范数:", np.linalg.norm(constraint_residuals))

通过观察残差是否接近于零，我们可以判断解的质量。对于线性等式约束，通常期望约束残差非常接近零（在浮点精度范围内）。

6. 注意事项与最佳实践

• 约束类型： 此方法专门适用于线性等式约束。对于线性不等式约束（例如 CX >= d）或非线性约束，np.linalg.lstsq 不适用。在这种情况下，需要使用更通用的优化库，如 SciPy 的 scipy.optimize.linprog（用于线性规划）或 scipy.optimize.minimize（配合适当的方法，如 SLSQP 或 COBYLA）。
• 系统一致性： np.linalg.lstsq 总是会返回一个解。如果原始方程组与约束之间存在内在矛盾，导致系统不一致，lstsq 将找到一个使所有方程（包括原始方程和约束方程）的残差平方和最小的“妥协”解。用户需要理解这一点，并根据实际问题判断解的合理性。
• 数值稳定性： np.linalg.lstsq 内部通常采用奇异值分解 (SVD) 等数值稳定的方法。rcond 参数对于处理病态或秩亏的系统非常有用，可以帮助避免数值问题。
• 大型稀疏系统： 对于非常大且稀疏的矩阵 A 和 C，np.linalg.lstsq 可能不是最高效的选择，因为它会生成一个稠密的增广矩阵。在这种情况下，可以考虑使用 SciPy 的稀疏线性代数模块 (scipy.sparse.linalg) 中的迭代求解器，例如 lsmr 或 least_squares，它们可以利用矩阵的稀疏性。

7. 总结

通过将原始线性方程组 AX=b 与线性等式约束 CX=d 合并成一个增广系统 A_aug X = b_aug，并利用 NumPy 的 np.linalg.lstsq 函数，我们可以高效且准确地求解带线性等式约束的线性系统。这种方法简洁、直接，能够同时满足原始方程和所有约束，是处理这类问题的强大工具。理解其背后的最小二乘原理和适用范围，将有助于在实际应用中做出正确的选择。

文中关于的知识介绍，希望对你的学习有所帮助！若是受益匪浅，那就动动鼠标收藏这篇《用NumPy解决带约束线性方程组》文章吧，也可关注golang学习网公众号了解相关技术文章。