求解线性方程组的算法（只调用Numpy）

显允
114 阅读5分钟

求解线性方程组的算法（只调用Numpy）

NumPy 是 Python 语言的一个扩展程序库，支持大量的维度数组与矩阵运算，此外也针对数组运算提供大量的数学函数库。此外求解线性方程组是数值计算的核心问题，其一般形式为 Ax=b ,其中 A 是系数矩阵（维度m×n），x是未知数向量（维度n×1），b是常数项向量（维度m×1）。根据A的性质（是否方阵、可逆、正定等），需选择不同的求解算法。本文仅使用 NumPy 详解常用方法。

直接解法（适用于方阵非奇异矩阵）

当 m=n 且 A 可逆（满秩）时，方程组有唯一解，可通过直接法求解。

逆矩阵法

原理：若A可逆，则解为 x=A-1b ,其适用场景为小型方阵（计算量随维度立方增长，数值稳定性较差）。

# 逆矩阵法
import numpy as np
A = np.array([[2,1],[1,3]],dtype=np.float64)
b = np.array([4,5],dtype=np.float64)
A_inv = np.linalg.inv(A)
x = A_inv @ b
print("解：",x)
print("验证：",A@x)
# 输出 解： [1.4 1.2]
# 输出 验证： [4. 5.]
LU分解法

原理：将A分解为下三角矩阵L（对角线为1）和上三角矩阵U，及 A=LU 。
求解分两步：
1、前代法解 Ly = b (得到y)
2、回代法解 Ux = y (得到x)
适用场景：可逆方阵，比逆矩阵更高效（避免直接求逆）。 Numpy无法直接LU分解函数，但 np.linalg.solve 内部基于LU分解实现（适用于可逆方阵）：

# 用solve直接求解（底层LU分解）
x = np.linalg.solve(A,b)
print("解：",x)
# 解： [1.4 1.2]

手动实现简单LU分解（以2×2矩阵为例）：

def lu_solve(A,b):
    n = A.shape[0]
    L = np.eye(n)
    U = A.copy()
    #分解A为L和U
    for i in range(n-1):
        L[i+1:,i] = U[i+1:,i]/U[i,i]
        U[i+1:] = U[i+1:] - L[i+1:, i:i+1] @ U[i:i+1]
     # 前代法解Ly = b
    y = np.zeros(n)
    for i in range(n):
        y[i] = b[i] - L[i, :i] @ y[:i]
    
    # 回代法解Ux = y
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (y[i] - U[i, i+1:] @ x[i+1:]) / U[i, i]
    return x

x = lu_solve(A, b)
print("LU分解解：", x)  # 输出 [1.4 1.2]
Cholesky分解法

原理：若A是正定对称矩阵（A=AT且所有特征值为正），可分解为 A=LLT（L是下三角矩阵，对角线为正）。
求解分两步：
1、前代法解 Ly=b ;
2、回代法解 LTx=y 。
适用场景：正定对称矩阵（如协方差矩阵），计算量是 LU 分解的一半，数值稳定。

import numpy as np
# 正定对称矩阵示例：4x + 2y = 6；2x + 5y = 7
A = np.array([[4, 2], [2, 5]], dtype=np.float64)  # 对称且顺序主子式>0（正定）
b = np.array([6, 7], dtype=np.float64)

# Cholesky分解
L = np.linalg.cholesky(A)

# 前代法解Ly = b
y = np.linalg.solve(L, b)  # 等价于手动前代

# 回代法解L^T x = y
x = np.linalg.solve(L.T, y)  # 等价于手动回代

print("解：", x)  # 输出 [1. 1.]
print("验证：", A @ x)  # 输出 [6. 7.]（与b一致）
适用于任意矩阵的方法（含长方阵）

当 m!=n (超定:m>n;欠定:m<n)或A奇异(秩亏)时,需用最小二乘或伪逆求解。
1、QR分解法（正交三角分解）
正交矩阵（Orthogonal Matrix）是一类特殊的方阵，其行向量和列向量都是标准正交向量组。具体来说，若一个 n*n 实矩阵 Q 满足以下条件：
                                                              QTQ = QQT = I
其中表示的转置矩阵，表示阶单位矩阵，则称为正交矩阵。
正交矩阵的性质：
1、行列式值：正交矩阵的行列式值为1或-1，即： det(Q)=±1
2、向量长度保持：对于任意向量x，正交变换不改变向量的长度，即： ||Qx||=||x||
3、向量间夹角保持：正交变换不改变向量间的夹角，即对任意向量x和y，有： (Qx)T(Qy)=xTy
4、逆矩阵等于转置：正交矩阵的逆矩阵等于其转置矩阵，即： Q-1=QT

则，QR分解是最常见的矩阵分解方式之一。设有任意矩阵且（若可研究），那么一定可以分解成两个矩阵的乘积：表示正交阵（orthogonal matrix），表示上三角阵（upper triangular matrix）。
其中是一个正交阵，是一个上三角阵。根据、形状的不同，QR分解也分两种形式：
注：如无特殊说明，正交阵一般指正交单位阵。
原理：将 A 分解为正交矩阵 Q（ QTQ=I）和上三角矩阵 R，即 A=QR。
对于超定方程组(m > n)，最小二乘解（使 ||Ax-b||2 最小）满足：Rx = QTb（解上三角方程组即可）。
适用场景：任意矩阵（尤其是超定方程组），数值稳定性好。

# 超定方程组示例：3个方程2个未知数 
A = np.array([[1, 2], [3, 4], [5, 6]], dtype=np.float64) # 3x2矩阵 
b = np.array([7, 11, 15], dtype=np.float64) 
# QR分解 
Q, R = np.linalg.qr(A) # Q:3x3正交矩阵；R:3x2上三角（有效部分为2x2） 
# 求解Rx = Q^T b（取R的前2行和Q^T b的前2个元素） 
x = np.linalg.solve(R[:2, :2], (Q.T @ b)[:2]) 
print("最小二乘解：", x) # 输出 [1. 3.] 
print("残差：", np.linalg.norm(A @ x - b)) # 残差接近0（本例构造为精确解）
SVD 分解

将在下一篇文章详细介绍

avatar
文章
阅读
粉丝