每日 13

共轭梯度法：
目标： 求解线性方程组 Ax = b，其中 A 是对称正定矩阵。
核心思想：
共轭方向： 选择一组关于矩阵 A 共轭的搜索方向，避免重复搜索，加速收敛。
残差最小化： 每次迭代沿当前共轭方向调整解，使残差（r = b - Ax）的范数最小。
算法步骤：
初始化：
初始解：x₀
初始残差：r₀ = b - Ax₀
初始搜索方向：p₀ = r₀
迭代更新（k = 0, 1, 2, ... 直到收敛）：
计算步长：αₖ = (rₖᵀ rₖ) / (pₖᵀ A pₖ)
更新解：xₖ₊₁ = xₖ + αₖ pₖ
更新残差：rₖ₊₁ = rₖ - αₖ A pₖ
计算系数：βₖ = (rₖ₊₁ᵀ rₖ₊₁) / (rₖᵀ rₖ)
更新搜索方向：pₖ₊₁ = rₖ₊₁ + βₖ pₖ
终止条件：
当 ||rₖ|| 小于预设阈值时停止迭代。
数学原理：
共轭性条件： 搜索方向满足 pᵢᵀ A pⱼ = 0 (i ≠ j)，确保每个方向上的误差独立最小化。
Krylov 子空间： 共轭梯度法在 Krylov 子空间 Kₖ(A, r₀) = span{r₀, Ar₀, ..., Aᵏ⁻¹r₀} 中寻找最优解。
收敛性分析：
理论收敛步数： 对于 n × n 的矩阵，最多 n 步可得到精确解（实际中因舍入误差可能需要更多步）。
误差上界： ||xₖ - x||_A ≤ 2((κ - 1) / (κ + 1))ᵏ ||x₀ - x||_A，其中 κ = λ_max / λ_min 是矩阵的条件数。 条件数越小，收敛越快。

Python 代码示例：
import numpy as np

def conjugate_gradient(A, b, x0, max_iter=1000, tol=1e-6):
x = x0
r = b - A @ x
p = r.copy()
rsold = r.dot(r)

for i in range(max_iter):  
    Ap = A @ p  
    alpha = rsold / p.dot(Ap)  
    x += alpha * p  
    r -= alpha * Ap  
    rsnew = r.dot(r)  
    if np.sqrt(rsnew) < tol:  
        break  
    beta = rsnew / rsold  
    p = r + beta * p  
    rsold = rsnew  
return x