共轭梯度法避坑指南：从数学推导到NumPy实现中的5个关键细节

共轭梯度法避坑指南：从数学推导到NumPy实现中的5个关键细节

在数值计算的世界里，共轭梯度法（Conjugate Gradient Method）常常被描述为一种“优雅”的算法——它兼具理论上的简洁性和实践中的高效性，尤其适合处理大规模稀疏对称正定线性系统。然而，这种优雅的表象之下，却隐藏着不少让开发者头疼的“坑”。许多人在初次实现时，往往会被其数学上的完美推导所迷惑，以为照着公式敲出代码就能万事大吉，结果却在真实数据上遭遇收敛失败、数值不稳定甚至结果完全错误的尴尬局面。

我自己在工程实践中就曾多次踩坑。记得有一次，在处理一个有限元分析生成的大型刚度矩阵时，我严格按照教科书上的FR（Fletcher-Reeves）公式实现了算法，迭代了几百步后残差仍然居高不下。调试了整整两天，最后才发现问题出在了一个极其细微的地方：迭代过程中方向向量失去了共轭性，而根源竟然是浮点误差的累积。那次经历让我深刻意识到，从数学公式到可用的代码，中间隔着一条名为“数值稳定性”的鸿沟。

这篇文章正是为那些已经理解共轭梯度法基本数学原理，正准备或正在将其付诸实践的中级开发者而写。我不会重复教科书上的标准推导，而是聚焦于那些在实现过程中真正影响算法表现、却又容易被忽略的五个关键细节。我们将深入探讨矩阵条件数如何悄无声息地破坏收敛、迭代终止条件设定的艺术、不同变体（FR/PRP/HS）在实际场景中的选择策略，以及如何设计有效的重启机制。更重要的是，我会结合具体的NumPy代码示例，展示如何识别和绕过这些陷阱，让你的共轭梯度法实现真正变得健壮可靠。

1. 矩阵条件数：隐形的性能杀手与预处理艺术

几乎所有关于共轭梯度法的教科书都会提到，它要求系数矩阵A是对称正定的。这个条件看似简单，但在实际应用中，满足对称正定只是拿到了入场券，真正决定算法表现的是矩阵的条件数。条件数衡量了矩阵对输入误差的敏感程度，条件数越大，矩阵越“病态”，共轭梯度法的收敛速度就会越慢，甚至停滞不前。

1.1 条件数如何影响收敛速度

从理论上讲，共轭梯度法在精确算术下的收敛速度与矩阵的特征值分布密切相关。一个经典的误差估计表明，经过k次迭代后，误差的A-范数满足：

||e_k||_A ≤ 2 * [ (√κ - 1)/(√κ + 1) ]^k * ||e_0||_A

其中κ是矩阵A的条件数（最大特征值与最小特征值之比）。从这个不等式可以清晰地看到，κ越大，括号内的值越接近1，收敛速度就越慢。当κ达到10^6甚至更高时，算法可能需要数千次迭代才能达到可接受的精度。

在实际代码中，即使矩阵理论上是正定的，浮点运算的舍入误差也可能导致算法实际“看到”的条件数比理论值更差。特别是当矩阵来自物理问题的离散化（如有限元、有限差分）时，网格细化往往会导致条件数随着问题规模的增大而急剧恶化。

下面这个简单的例子展示了条件数对迭代次数的影响：

import numpy as np
import matplotlib.pyplot as plt

def test_condition_number_impact():
    """测试不同条件数对共轭梯度法收敛速度的影响"""
    np.random.seed(42)
    results = []
    
    for cond in [10, 100, 1000, 10000]:
        # 构造一个具有指定条件数的对称正定矩阵
        n = 50
        U, _ = np.linalg.qr(np.random.randn(n, n))
        # 创建特征值从1到cond的矩阵
        eigenvalues = np.linspace(1, cond, n)
        A = U @ np.diag(eigenvalues) @ U.T
        
        # 真实解和右端项
        x_true = np.random.randn(n)
        b = A @ x_true
        
        # 标准共轭梯度法实现
        x = np.zeros_like(b)
        r = b - A @ x
        p = r.copy()
        rsold = r @ r
        
        residuals = []
        for i in range(n * 2):  # 最多迭代2n次
            Ap = A @ p
            alpha = rsold / (p @ Ap)
            x = x + alpha * p
            r = r - alpha * Ap
            rsnew = r @ r
            residuals.append(np.sqrt(rsnew))
            
            if rsnew < 1e-12:
                break
                
            beta = rsnew / rsold
            p = r + beta * p
            rsold = rsnew
        
        results.append({
            'condition_number': cond,
            'iterations': len(residuals),
            'final_residual': residuals[-1] if residuals else None
        })
    
    return results

运行这个测试，你会发现当条件数从10增加到10000时，达到相同精度所需的迭代次数可能增加一个数量级。这不仅仅是理论上的推演，而是每个实现者都会遇到的现实问题。

1.2 预处理技术：从理论到实践的选择

既然条件数问题是根源，那么解决方案就是预处理。预处理的基本思想是找到一个矩阵M，使得M^{-1}A的条件数远小于A本身的条件数，然后求解等价系统M^{-1}Ax = M^{-1}b。好的预处理子M应该满足几个条件：1) M^{-1}A的条件数小；2) M易于求逆或求解；3) 构造M的代价不能太高。

注意：预处理共轭梯度法（PCG）的实现需要特别注意矩阵的对称性。即使A和M都是对称正定的，M^{-1}A通常不再对称。因此在实际算法中，我们采用不同的形式来保持对称性。

常见的预处理子包括：

下面是一个带雅可比预处理的共轭梯度法实现示例：

def preconditioned_conjugate_gradient(A, b, M_inv, max_iter=1000, tol=1e-10):
    """
    带预处理的共轭梯度法实现
    
    参数:
        A: 对称正定矩阵 (n x n)
        b: 右端向量 (n,)
        M_inv: 预处理矩阵的逆或预处理算子
        max_iter: 最大迭代次数
        tol: 收敛容差
    
    返回:
        x: 近似解
        residuals: 残差历史
    """
    n = len(b)
    x = np.zeros(n)
    
    # 初始残差
    r = b - A @ x
    # 预处理残差
    z = M_inv @ r if callable(M_inv) else M_inv @ r
    
    p = z.copy()
    rz_old = r @ z
    
    residuals = [np.linalg.norm(r)]
    
    for k in range(max_iter):
        Ap = A @ p
        alpha = rz_old / (p @ Ap)
        
        x = x + alpha * p
        r = r - alpha * Ap
        
        residual_norm = np.linalg.norm(r)
        residuals.append(residual_norm)
        
        if residual_norm < tol:
            break
        
        # 预处理步骤
        z = M_inv @ r if callable(M_inv) else M_inv @ r
        rz_new = r @ z
        
        beta = rz_new / rz_old
        p = z + beta * p
        rz_old = rz_new
    
    return x, residuals

在实际选择预处理子时，需要权衡效果和成本。对于小规模稠密矩阵，可能不需要预处理；对于大规模稀疏矩阵，不完全Cholesky通常是首选；而对于来自特定物理问题的矩阵，基于问题知识的定制预处理子往往能取得最佳效果。

2. 迭代终止条件：不仅仅是残差范数那么简单

在教科书的示例中，迭代终止条件通常被简化为‖r_k‖ < ε，其中r_k是第k步的残差。但在实际实现中，这个看似简单的判断蕴含着多个陷阱。

2.1 绝对容差与相对容差的权衡

使用绝对残差‖r_k‖ < ε的主要问题是，ε的选择与问题的尺度相关。对于‖b‖很大的问题，即使相对误差很大，绝对残差也可能很小；反之，对于‖b‖很小的问题，即使解已经很精确，绝对残差也可能达不到要求。

更健壮的终止条件应该考虑相对残差：

终止条件1: ‖r_k‖ / ‖b‖ < ε_rel
终止条件2: ‖r_k‖ / ‖r_0‖ < ε_rel

但即使这样还不够。考虑下面这种情况：当解x的范数很大时，即使相对残差很小，绝对误差‖x - x*‖可能仍然很大。因此，更完整的终止条件应该结合多种判断：

def should_stop(r, b, x, A, k, tol_abs=1e-10, tol_rel=1e-8, max_iter=1000):
    """
    综合判断是否应该终止迭代
    
    参数:
        r: 当前残差
        b: 右端项
        x: 当前解
        A: 系数矩阵
        k: 当前迭代次数
        tol_abs: 绝对容差
        tol_rel: 相对容差
        max_iter: 最大迭代次数
    
    返回:
        bool: 是否终止
        str: 终止原因
    """
    r_norm = np.linalg.norm(r)
    b_norm = np.linalg.norm(b)
    
    # 超过最大迭代次数
    if k >= max_iter:
        return True, "达到最大迭代次数"
    
    # 绝对残差足够小
    if r_norm < tol_abs:
        return True, "绝对残差满足要求"
    
    # 相对残差足够小
    if b_norm > 0 and r_norm / b_norm < tol_rel:
        return True, "相对残差满足要求"
    
    # 残差不再显著下降（停滞检测）
    if k > 10:
        # 可以记录历史残差，检查是否停滞
        pass
    
    return False, "继续迭代"

2.2 向后误差与向前误差

在数值线性代数中，我们还需要区分向后误差和向前误差：

• 向后误差：测量的是解满足方程的程度，即‖b - Ax‖
• 向前误差：测量的是解与真实解的差距，即‖x - x*‖

共轭梯度法天然地控制向后误差，但用户通常关心的是向前误差。对于病态问题，即使向后误差很小，向前误差也可能很大。一个实用的技巧是估计条件数，然后根据条件数调整容差要求。

def estimate_condition_number(A, iterations=20):
    """
    使用幂迭代法估计对称正定矩阵的条件数
    
    参数:
        A: 对称正定矩阵
        iterations: 幂迭代次数
    
    返回:
        条件数的估计值
    """
    n = A.shape[0]
    
    # 估计最大特征值
    v = np.random.randn(n)
    v = v / np.linalg.norm(v)
    
    for _ in range(iterations):
        Av = A @ v
        v = Av / np.linalg.norm(Av)
    
    lambda_max = v @ A @ v
    
    # 估计最小特征值（通过A的逆）
    # 注意：这里使用共轭梯度法求解来避免显式求逆
    b = np.random.randn(n)
    x, _ = preconditioned_conjugate_gradient(A, b, np.eye(n), max_iter=100)
    lambda_min_estimate = 1.0  # 简化处理
    
    # 更准确的方法需要更复杂的实现
    return lambda_max / lambda_min_estimate

2.3 实用的终止策略

基于以上分析，我推荐在实际实现中使用如下的综合终止策略：

1. 主要终止条件：相对残差‖r_k‖/‖b‖ < ε_rel
2. 辅助终止条件：绝对残差‖r_k‖ < ε_abs
3. 停滞检测：如果连续多次迭代残差下降小于某个因子（如1e-3），则警告或终止
4. 最大迭代次数限制：通常设为2n，但对于预处理良好的系统可以更少
5. 时间限制：对于实时应用，设置最大运行时间

下面是一个完整的实现示例：

def robust_conjugate_gradient(A, b, x0=None, max_iter=None, tol=1e-8, 
                             stagnation_threshold=1e-3, stagnation_window=10):
    """
    带有健壮终止条件的共轭梯度法
    
    参数:
        A: 对称正定矩阵
        b: 右端向量
        x0: 初始猜测（可选）
        max_iter: 最大迭代次数，默认为2n
        tol: 相对容差
        stagnation_threshold: 停滞检测阈值
        stagnation_window: 停滞检测窗口大小
    """
    n = len(b)
    if max_iter is None:
        max_iter = 2 * n
    
    if x0 is None:
        x = np.zeros(n)
    else:
        x = x0.copy()
    
    r = b - A @ x
    p = r.copy()
    rsold = r @ r
    
    b_norm = np.linalg.norm(b)
    if b_norm == 0:
        b_norm = 1
    
    residuals = [np.sqrt(rsold) / b_norm]
    stagnation_count = 0
    
    for k in range(max_iter):
        Ap = A @ p
        alpha = rsold / (p @ Ap)
        x = x + alpha * p
        r = r - alpha * Ap
        rsnew = r @ r
        
        # 计算相对残差
        rel_residual = np.sqrt(rsnew) / b_norm
        residuals.append(rel_residual)
        
        # 检查停滞
        if k >= stagnation_window:
            improvement_ratio = residuals[-stagnation_window] / residuals[-1]
            if improvement_ratio < 1 + stagnation_threshold:
                stagnation_count += 1
                if stagnation_count > 3:
                    print(f"警告：迭代在{k}步后停滞")
                    break
            else:
                stagnation_count = 0
        
        # 检查收敛
        if rel_residual < tol:
            break
        
        if rsnew < 1e-30:  # 防止除零
            break
        
        beta = rsnew / rsold
        p = r + beta * p
        rsold = rsnew
    
    return x, residuals

3. 重启策略：何时清零方向向量

标准的共轭梯度法理论保证，对于n维二次函数，最多n次迭代就能得到精确解（在精确算术下）。但在实际浮点运算中，由于舍入误差的累积，共轭方向会逐渐失去共轭性，导致算法需要更多迭代，甚至无法收敛到足够精度。这就是为什么需要重启策略。

3.1 为什么需要重启

在理想情况下，共轭梯度法生成的方向向量应该是A-共轭的，即对于i≠j，有p_i^T A p_j = 0。但在浮点运算中，这个性质会逐渐被破坏。主要原因包括：

1. 舍入误差累积：每次迭代中的浮点运算都会引入微小误差
2. 线性依赖：数值上，方向向量可能变得近似线性相关
3. 正交性丢失：残差向量的正交性逐渐丧失

当这些数值问题积累到一定程度时，算法实际上是在低维子空间中“打转”，无法进一步降低残差。重启策略的核心思想是：定期抛弃历史信息，重新从当前点的负梯度方向开始，相当于“重置”算法的状态。

3.2 常见的重启策略

实践中主要有三种重启策略：

1. 固定间隔重启 最简单的策略是每k次迭代重启一次，通常k取n（问题维度）或更小。对于大规模问题，n可能很大，这时可以取一个固定值如50或100。

def conjugate_gradient_with_restart(A, b, restart_interval=50, max_iter=1000, tol=1e-8):
    """带固定间隔重启的共轭梯度法"""
    n = len(b)
    x = np.zeros(n)
    r = b - A @ x
    b_norm = np.linalg.norm(b)
    
    residuals = []
    total_iterations = 0
    
    while total_iterations < max_iter:
        # 开始新一轮（重启）
        p = r.copy()
        rsold = r @ r
        
        for k in range(min(restart_interval, max_iter - total_iterations)):
            Ap = A @ p
            alpha = rsold / (p @ Ap)
            x = x + alpha * p
            r = r - alpha * Ap
            rsnew = r @ r
            
            rel_residual = np.sqrt(rsnew) / b_norm
            residuals.append(rel_residual)
            total_iterations += 1
            
            if rel_residual < tol:
                return x, residuals
            
            if rsnew < 1e-30:
                break
            
            beta = rsnew / rsold
            p = r + beta * p
            rsold = rsnew
        
        # 如果达到重启间隔但未收敛，继续下一轮
        if np.sqrt(rsnew) / b_norm < tol:
            break
    
    return x, residuals

2. 自适应重启 更智能的策略是根据算法表现决定何时重启。一个常见的方法是监控残差下降速度：

def adaptive_restart_conjugate_gradient(A, b, max_iter=1000, tol=1e-8, 
                                        patience=20, improvement_threshold=0.1):
    """带自适应重启的共轭梯度法"""
    n = len(b)
    x = np.zeros(n)
    r = b - A @ x
    b_norm = np.linalg.norm(b)
    
    residuals = []
    no_improvement_count = 0
    best_residual = float('inf')
    
    iteration = 0
    while iteration < max_iter:
        # 开始新一轮
        p = r.copy()
        rsold = r @ r
        
        inner_iter = 0
        while iteration < max_iter and inner_iter < n:
            Ap = A @ p
            alpha = rsold / (p @ Ap)
            x = x + alpha * p
            r = r - alpha * Ap
            rsnew = r @ r
            
            rel_residual = np.sqrt(rsnew) / b_norm
            residuals.append(rel_residual)
            
            iteration += 1
            inner_iter += 1
            
            # 检查收敛
            if rel_residual < tol:
                return x, residuals
            
            # 检查改进情况
            if rel_residual < best_residual * (1 - improvement_threshold):
                best_residual = rel_residual
                no_improvement_count = 0
            else:
                no_improvement_count += 1
            
            # 如果长时间没有改进，考虑重启
            if no_improvement_count >= patience:
                print(f"在{iteration}次迭代后重启")
                break
            
            if rsnew < 1e-30:
                break
            
            beta = rsnew / rsold
            p = r + beta * p
            rsold = rsnew
    
    return x, residuals

3. 基于正交性的重启 最理论化的方法是监控方向向量的共轭性或残差的正交性。当这些性质破坏到一定程度时重启：

def orthogonal_based_restart(A, b, orthogonality_threshold=1e-8, max_iter=1000, tol=1e-8):
    """基于正交性检测的共轭梯度法"""
    n = len(b)
    x = np.zeros(n)
    r = b - A @ x
    b_norm = np.linalg.norm(b)
    
    # 存储历史方向向量用于正交性检查
    P = []  # 方向向量列表
    residuals = []
    
    iteration = 0
    while iteration < max_iter:
        p = r.copy()
        
        # 检查与历史方向的正交性
        if P:
            for prev_p in P[-5:]:  # 只检查最近几个方向
                orthogonality = abs(p @ A @ prev_p) / (np.linalg.norm(p) * np.linalg.norm(prev_p))
                if orthogonality > orthogonality_threshold:
                    print(f"在{iteration}次迭代检测到正交性破坏: {orthogonality:.2e}")
                    # 清空历史，重启
                    P = []
                    break
        
        P.append(p.copy())
        if len(P) > 10:  # 只保留最近10个方向
            P.pop(0)
        
        rsold = r @ r
        
        # 标准CG迭代
        Ap = A @ p
        alpha = rsold / (p @ Ap)
        x = x + alpha * p
        r = r - alpha * Ap
        rsnew = r @ r
        
        rel_residual = np.sqrt(rsnew) / b_norm
        residuals.append(rel_residual)
        iteration += 1
        
        if rel_residual < tol:
            break
        
        if rsnew < 1e-30:
            break
        
        beta = rsnew / rsold
        p = r + beta * p
        rsold = rsnew
    
    return x, residuals

3.3 重启策略的选择建议

选择哪种重启策略取决于具体问题：

• 固定间隔重启：最简单，适合大多数常规问题。对于n维问题，重启间隔设为n通常是个不错的选择。
• 自适应重启：更智能，能适应不同问题的特性，但需要调整参数。
• 基于正交性的重启：理论最优雅，但计算开销较大，适合对精度要求极高的问题。

在实际应用中，我通常从固定间隔重启开始，如果发现收敛速度不理想，再尝试自适应策略。对于特别病态的问题，基于正交性的重启可能更可靠。

4. 变体选择：FR、PRP、HS等算法的实战对比

共轭梯度法有多种变体，主要区别在于β_k的计算公式。最常见的三种是：

• FR（Fletcher-Reeves）: β_k = (r_{k+1}^T r_{k+1}) / (r_k^T r_k)
• PRP（Polak-Ribière-Polyak）: β_k = (r_{k+1}^T (r_{k+1} - r_k)) / (r_k^T r_k)
• HS（Hestenes-Stiefel）: β_k = (r_{k+1}^T (r_{k+1} - r_k)) / (p_k^T (r_{k+1} - r_k))

这些变体在理论上各有优劣，但在实际应用中，选择哪个往往取决于具体问题。

4.1 不同变体的数学特性对比

4.2 实际性能测试

让我们通过一个实际测试来比较这些变体的表现：

def conjugate_gradient_variant(A, b, variant='FR', max_iter=1000, tol=1e-12):
    """不同变体的共轭梯度法实现"""
    n = len(b)
    x = np.zeros(n)
    r = b - A @ x
    p = r.copy()
    
    rsold = r @ r
    residuals = [np.sqrt(rsold)]
    
    for k in range(max_iter):
        Ap = A @ p
        alpha = rsold / (p @ Ap)
        x = x + alpha * p
        r = r - alpha * Ap
        rsnew = r @ r
        
        residuals.append(np.sqrt(rsnew))
        
        if rsnew < tol:
            break
        
        # 根据变体计算beta
        if variant == 'FR':
            beta = rsnew / rsold
        elif variant == 'PRP':
            delta_r = r - (b - A @ (x - alpha * p))  # 计算r_{k+1} - r_k
            beta = (r @ delta_r) / rsold
        elif variant == 'HS':
            delta_r = r - (b - A @ (x - alpha * p))
            beta = (r @ delta_r) / (p @ delta_r)
        else:
            raise ValueError(f"未知变体: {variant}")
        
        # 确保beta非负（对FR和PRP的修正）
        if variant in ['FR', 'PRP']:
            beta = max(beta, 0)
        
        p = r + beta * p
        rsold = rsnew
    
    return x, residuals

def compare_variants():
    """比较不同变体的性能"""
    np.random.seed(42)
    
    # 创建测试问题
    n = 100
    # 构造一个条件数较大的矩阵
    U, _ = np.linalg.qr(np.random.randn(n, n))
    eigenvalues = np.logspace(0, 6, n)  # 条件数10^6
    A = U @ np.diag(eigenvalues) @ U.T
    
    x_true = np.random.randn(n)
    b = A @ x_true
    
    variants = ['FR', 'PRP', 'HS']
    results = {}
    
    for variant in variants:
        x, residuals = conjugate_gradient_variant(A, b, variant=variant, max_iter=2*n)
        results[variant] = {
            'final_residual': residuals[-1],
            'iterations': len(residuals),
            'residuals': residuals
        }
        error = np.linalg.norm(x - x_true) / np.linalg.norm(x_true)
        print(f"{variant}: {len(residuals)}次迭代，相对误差={error:.2e}")
    
    return results

运行这个测试，你可能会观察到：

• FR：最稳定，但收敛速度可能不是最快的
• PRP：通常收敛更快，但偶尔会不稳定
• HS：在条件数适中时表现良好，但病态问题时可能数值不稳定

4.3 混合策略与自适应选择

在实际应用中，我经常使用一种混合策略：开始时使用PRP以获得快速收敛，当检测到可能的不稳定时切换到FR。下面是一个自适应选择变体的实现：

def adaptive_cg(A, b, max_iter=1000, tol=1e-12):
    """自适应选择变体的共轭梯度法"""
    n = len(b)
    x = np.zeros(n)
    r = b - A @ x
    p = r.copy()
    
    rsold = r @ r
    residuals = [np.sqrt(rsold)]
    
    # 初始使用PRP
    current_variant = 'PRP'
    variant_switch_count = 0
    
    for k in range(max_iter):
        Ap = A @ p
        alpha = rsold / (p @ Ap)
        x_new = x + alpha * p
        r_new = r - alpha * Ap
        rsnew = r_new @ r_new
        
        # 检查步长是否合理
        if alpha < 1e-14 or alpha > 1e14:
            print(f"第{k}次迭代：步长异常({alpha:.2e})，切换到FR")
            current_variant = 'FR'
            variant_switch_count += 1
        
        # 检查残差是否下降
        if rsnew > rsold * 0.999:  # 残差几乎没有下降
            print(f"第{k}次迭代：残差下降缓慢，切换到FR")
            current_variant = 'FR'
            variant_switch_count += 1
        
        # 计算beta
        if current_variant == 'FR':
            beta = rsnew / rsold
        elif current_variant == 'PRP':
            delta_r = r_new - r
            beta_prp = (r_new @ delta_r) / rsold
            # PRP+修正：确保beta非负
            beta = max(beta_prp, 0)
        else:
            beta = rsnew / rsold  # 默认回退到FR
        
        # 更新
        x = x_new
        r = r_new
        p = r + beta * p
        rsold = rsnew
        
        residuals.append(np.sqrt(rsnew))
        
        if rsnew < tol or variant_switch_count > 3:
            break
    
    return x, residuals

这种自适应策略结合了不同变体的优点：在算法表现良好时使用PRP加速收敛，在出现数值问题时自动切换到更稳定的FR。

5. 收敛性测试与调试技巧

实现共轭梯度法后，如何验证它的正确性和评估它的性能？这就需要系统的测试方法。我经常看到开发者只用一个简单的例子测试，然后就认为算法正确，这往往会在实际应用中埋下隐患。

5.1 构建全面的测试套件

一个完整的测试套件应该包含以下几种类型的测试问题：

1. 理想条件测试：小规模、良态矩阵，验证基本正确性
2. 病态问题测试：高条件数矩阵，测试数值稳定性
3. 大规模稀疏测试：模拟实际应用场景
4. 边界情况测试：零向量、奇异近似的矩阵等

class CGTestSuite:
    """共轭梯度法测试套件"""
    
    @staticmethod
    def test_well_conditioned():
        """良态矩阵测试"""
        np.random.seed(42)
        n = 50
        # 创建条件数约10的矩阵
        A = np.random.randn(n, n)
        A = A.T @ A + np.eye(n) * 0.1  # 确保正定
        x_true = np.random.randn(n)
        b = A @ x_true
        
        x, residuals = robust_conjugate_gradient(A, b, tol=1e-12)
        error = np.linalg.norm(x - x_true) / np.linalg.norm(x_true)
        
        print(f"良态测试：相对误差={error:.2e}，迭代次数={len(residuals)}")
        return error < 1e-10
    
    @staticmethod
    def test_ill_conditioned():
        """病态矩阵测试"""
        np.random.seed(42)
        n = 100
        # 创建条件数约1e6的矩阵
        U, _ = np.linalg.qr(np.random.randn(n, n))
        eigenvalues = np.logspace(0, 6, n)  # 从1到10^6
        A = U @ np.diag(eigenvalues) @ U.T
        x_true = np.random.randn(n)
        b = A @ x_true
        
        # 不带预处理
        x1, res1 = robust_conjugate_gradient(A, b, tol=1e-6, max_iter=500)
        error1 = np.linalg.norm(x1 - x_true) / np.linalg.norm(x_true)
        
        # 带对角预处理
        M_inv = np.diag(1.0 / np.diag(A))
        x2, res2 = preconditioned_conjugate_gradient(A, b, M_inv, tol=1e-6, max_iter=500)
        error2 = np.linalg.norm(x2 - x_true) / np.linalg.norm(x_true)
        
        print(f"病态测试（无预处理）：相对误差={error1:.2e}，迭代次数={len(res1)}")
        print(f"病态测试（有预处理）：相对误差={error2:.2e}，迭代次数={len(res2)}")
        
        return error2 < 1e-5
    
    @staticmethod
    def test_sparse_large():
        """大规模稀疏矩阵测试"""
        from scipy import sparse
        import time
        
        np.random.seed(42)
        n = 1000
        # 创建稀疏对称正定矩阵（类似有限元刚度矩阵）
        A = sparse.random(n, n, density=0.01, random_state=42)
        A = A @ A.T + sparse.eye(n) * 0.1
        A = A.tocsr()
        
        x_true = np.random.randn(n)
        b = A @ x_true
        
        # 转换为稠密矩阵用于标准CG（仅用于小规模测试）
        if n <= 1000:
            A_dense = A.toarray()
            start = time.time()
            x, residuals = robust_conjugate_gradient(A_dense, b, tol=1e-8)
            time_dense = time.time() - start
            error = np.linalg.norm(x - x_true) / np.linalg.norm(x_true)
            print(f"稠密实现：时间={time_dense:.3f}s，误差={error:.2e}，迭代次数={len(residuals)}")
        
        return True
    
    @staticmethod
    def test_edge_cases():
        """边界情况测试"""
        results = []
        
        # 测试1：b是零向量
        A = np.eye(5)
        b = np.zeros(5)
        x, _ = robust_conjugate_gradient(A, b)
        results.append(np.allclose(x, 0))
        
        # 测试2：A是单位矩阵
        A = np.eye(10)
        b = np.ones(10)
        x, _ = robust_conjugate_gradient(A, b, tol=1e-12)
        results.append(np.allclose(x, b))
        
        # 测试3：初始猜测就是解
        A = np.random.randn(5, 5)
        A = A.T @ A + np.eye(5)
        x_true = np.random.randn(5)
        b = A @ x_true
        x, residuals = robust_conjugate_gradient(A, b, x0=x_true.copy())
        results.append(len(residuals) <= 2)  # 应该立即收敛
        
        print(f"边界测试通过：{all(results)}")
        return all(results)
    
    @staticmethod
    def run_all_tests():
        """运行所有测试"""
        print("=" * 60)
        print("共轭梯度法测试套件")
        print("=" * 60)
        
        tests = [
            ("良态矩阵测试", CGTestSuite.test_well_conditioned),
            ("病态矩阵测试", CGTestSuite.test_ill_conditioned),
            ("大规模稀疏测试", CGTestSuite.test_sparse_large),
            ("边界情况测试", CGTestSuite.test_edge_cases),
        ]
        
        all_passed = True
        for test_name, test_func in tests:
            try:
                print(f"\n执行测试：{test_name}")
                passed = test_func()
                status = "通过" if passed else "失败"
                print(f"结果：{status}")
                all_passed = all_passed and passed
            except Exception as e:
                print(f"测试异常：{e}")
                all_passed = False
        
        print("\n" + "=" * 60)
        print(f"总体结果：{'所有测试通过' if all_passed else '有测试失败'}")
        print("=" * 60)
        
        return all_passed

5.2 调试与性能分析工具

当算法出现问题时，如何调试？以下是我常用的调试技巧：

1. 残差历史分析 绘制残差随迭代次数的变化图，可以揭示很多问题：

def analyze_convergence(A, b, solver_func, **kwargs):
    """分析收敛行为"""
    x, residuals = solver_func(A, b, **kwargs)
    
    import matplotlib.pyplot as plt
    
    plt.figure(figsize=(10, 6))
    plt.semilogy(residuals, 'b-', linewidth=2, label='残差')
    plt.xlabel('迭代次数')
    plt.ylabel('相对残差')
    plt.title('共轭梯度法收敛历史')
    plt.grid(True, alpha=0.3)
    plt.legend()
    
    # 添加理想收敛线（线性收敛）
    if len(residuals) > 1:
        # 估计收敛速率
        rates = []
        for i in range(1, len(residuals)):
            if residuals[i-1] > 0:
                rate = residuals[i] / residuals[i-1]
                rates.append(rate)
        
        avg_rate = np.mean(rates[-10:]) if len(rates) > 10 else np.mean(rates)
        plt.plot([0, len(residuals)], 
                 [residuals[0], residuals[0] * (avg_rate ** len(residuals))],
                 'r--', alpha=0.5, label=f'估计收敛速率: {avg_rate:.3f}')
        plt.legend()
    
    plt.show()
    
    # 打印统计信息
    print(f"最终残差: {residuals[-1]:.2e}")
    print(f"迭代次数: {len(residuals)}")
    print(f"平均收敛速率: {avg_rate:.3f}" if 'avg_rate' in locals() else "")
    
    return x, residuals

2. 正交性检查 定期检查方向向量的共轭性和残差的正交性：

def check_orthogonality(A, x_history, r_history, p_history):
    """检查共轭梯度法迭代过程中的正交性"""
    k = len(p_history)
    if k < 2:
        return
    
    print("\n正交性检查:")
    
    # 检查方向向量的A-共轭性
    print("方向向量的A-共轭性:")
    for i in range(k):
        for j in range(i+1, k):
            val = p_history[i] @ A @ p_history[j]
            if abs(val) > 1e-8:
                print(f"  p_{i}^T A p_{j} = {val:.2e} (应该接近0)")
    
    # 检查残差的正交性
    print("\n残差的正交性:")
    for i in range(k):
        for j in range(i+1, k):
            val = r_history[i] @ r_history[j]
            if abs(val) > 1e-8:
                print(f"  r_{i}^T r_{j} = {val:.2e} (应该接近0)")

3. 数值稳定性监控 在迭代过程中监控关键数值量的变化：

def monitored_cg(A, b, tol=1e-12, max_iter=1000):
    """带监控的共轭梯度法"""
    n = len(b)
    x = np.zeros(n)
    r = b - A @ x
    p = r.copy()
    
    rsold = r @ r
    b_norm = np.linalg.norm(b)
    
    # 监控数据
    monitor_data = {
        'residuals': [np.sqrt(rsold)/b_norm],
        'alphas': [],
        'betas': [],
        'orthogonality': [],
        'direction_norms': [np.linalg.norm(p)]
    }
    
    for k in range(max_iter):
        Ap = A @ p
        alpha = rsold / (p @ Ap)
        
        # 检查alpha的合理性
        if abs(alpha) < 1e-14:
            print(f"警告：第{k}次迭代，alpha={alpha:.2e}过小")
        elif abs(alpha) > 1e14:
            print(f"警告：第{k}次迭代，alpha={alpha:.2e}过大")
        
        x = x + alpha * p
        r = r - alpha * Ap
        rsnew = r @ r
        
        rel_residual = np.sqrt(rsnew) / b_norm
        monitor_data['residuals'].append(rel_residual)
        monitor_data['alphas'].append(alpha)
        
        if rel_residual < tol:
            break
        
        beta = rsnew / rsold
        monitor_data['betas'].append(beta)
        
        # 检查方向向量范数
        p_new = r + beta * p
        dir_norm = np.linalg.norm(p_new)
        monitor_data['direction_norms'].append(dir_norm)
        
        if dir_norm < 1e-10:
            print(f"警告：第{k}次迭代，方向向量范数过小: {dir_norm:.2e}")
        
        p = p_new
        rsold = rsnew
    
    return x, monitor_data

5.3 性能优化技巧

在确保正确性的基础上，我们还可以从性能角度优化实现：

1. 内存访问优化 共轭梯度法的主要计算开销是矩阵-向量乘法。对于稀疏矩阵，使用合适的存储格式（CSR、CSC等）可以大幅提升性能：

def sparse_cg(A_csr, b, tol=1e-8, max_iter=1000):
    """针对稀疏矩阵优化的共轭梯度法"""
    n = len(b)
    x = np.zeros(n)
    r = b - A_csr.dot(x)
    p = r.copy()
    
    rsold = r @ r
    b_norm = np.linalg.norm(b)
    
    for k in range(max_iter):
        # 使用稀疏矩阵乘法
        Ap = A_csr.dot(p)
        
        alpha = rsold / (p @ Ap)
        x = x + alpha * p
        r = r - alpha * Ap
        rsnew = r @ r
        
        if np.sqrt(rsnew) / b_norm < tol:
            break
        
        beta = rsnew / rsold
        p = r + beta * p
        rsold = rsnew
    
    return x

2. 并行计算 对于大规模问题，可以考虑并行化矩阵-向量乘法和向量运算：

def parallel_cg(A, b, tol=1e-8, max_iter=1000):
    """简单的并行化共轭梯度法（概念示例）"""
    from multiprocessing import Pool
    import numpy as np
    
    n = len(b)
    x = np.zeros(n)
    r = b - A @ x  # 这里需要并行化的矩阵乘法
    p = r.copy()
    
    rsold = r @ r
    b_norm = np.linalg.norm(b)
    
    for k in range(max_iter):
        # 在实际实现中，这里应该使用并行矩阵乘法
        Ap = A @ p  # 这行应该并行化
        
        alpha = rsold / (p @ Ap)
        x = x + alpha * p
        r = r - alpha * Ap
        rsnew = r @ r
        
        if np.sqrt(rsnew) / b_norm < tol:
            break
        
        beta = rsnew / rsold
        p = r + beta * p
        rsold = rsnew
    
    return x

3. 混合精度计算 在某些情况下，使用混合精度（如单精度存储矩阵，双精度进行累加）可以在不损失太多精度的情况下减少内存使用和计算时间：

def mixed_precision_cg(A, b, tol=1e-8, max_iter=1000):
    """混合精度共轭梯度法示例"""
    # 将矩阵转换为单精度以减少内存和计算开销
    A_single = A.astype(np.float32)
    
    # 但保持向量为双精度以保证精度
    x = np.zeros_like(b, dtype=np.float64)
    r = b - A_single.astype(np.float64) @ x
    p = r.copy()
    
    rsold = r @ r
    b_norm = np.linalg.norm(b)
    
    for k in range(max_iter):
        # 矩阵乘法用单精度
        Ap_single = A_single @ p.astype(np.float32)
        Ap = Ap_single.astype(np.float64)
        
        alpha = rsold / (p @ Ap)
        x = x + alpha * p
        r = r - alpha * Ap
        rsnew = r @ r
        
        if np.sqrt(rsnew) / b_norm < tol:
            break
        
        beta = rsnew / rsold
        p = r + beta * p
        rsold = rsnew
    
    return x

通过这些测试和优化技巧，你不仅能确保共轭梯度法实现的正确性，还能在实际应用中获得更好的性能。记住，一个健壮的实现需要经过各种边界情况和极端条件的考验，而性能优化则需要在保证数值稳定性的前提下进行。

在实际项目中，我通常会先实现一个基础版本，通过完整的测试套件验证其正确性，然后根据具体应用场景（矩阵特性、规模、精度要求等）选择合适的预处理方法、重启策略和变体，最后在必要时进行性能优化。这样的分层实现策略既能保证可靠性，又能满足性能需求。