Python实现共轭梯度法

共轭梯度法（Conjugate Gradient Method）是一种求解无约束最优化问题的迭代方法，特别适用于对称正定线性系统的求解。它结合了最速下降法和牛顿法的优点，能够在较少的迭代次数内达到较高的精度。

一、公式解释

共轭梯度法的关键在于构造一组共轭方向，使得在每个方向上进行一次线搜索后，都能使目标函数值有显著的下降。假设我们要求解的无约束优化问题为：

min f(x)

其中 x 是一个向量，f(x) 是目标函数。共轭梯度法的迭代公式如下：

x_{k+1} = x_k + α_k * d_k

其中：

• x_k 是第 k 次迭代的解；
• α_k 是第 k 次迭代的步长，通过线搜索确定；
• d_k 是第 k 次迭代的搜索方向，满足共轭性质；
• x_{k+1} 是第 k+1 次迭代的解。

搜索方向 d_k 的构造是共轭梯度法的核心。通常，初始搜索方向 d_0 选择为负梯度方向 -∇f(x_0)。在每次迭代中，新的搜索方向 d_k 通过以下公式更新：

d_k = -∇f(x_k) + β_k * d_{k-1}

其中 β_k 是一个参数，用于调整搜索方向，常见的计算方式有 Fletcher-Reeves 公式和 Polak-Ribiere 公式等。

二、应用场景

共轭梯度法在实际应用中具有广泛的用途，特别是在大规模线性系统、图像处理、机器学习等领域。由于它能够在不计算Hessian矩阵的情况下充分利用目标函数的二次性质，因此在处理大规模问题时具有较高的效率。

三、示例代码

以下是一个使用Python实现共轭梯度法的简单示例，用于求解二次函数的最小值问题：

import numpy as np  
  
def conjugate_gradient(A, b, x0, tol=1e-6, max_iter=100):  
    """  
    共轭梯度法求解 Ax = b  
    :param A: 对称正定矩阵  
    :param b: 右侧向量  
    :param x0: 初始解向量  
    :param tol: 收敛阈值  
    :param max_iter: 最大迭代次数  
    :return: 最优解向量 x  
    """  
    x = x0.copy()  
    r = b - np.dot(A, x)  
    p = r.copy()  
    rsold = np.dot(r, r)  
  
    for k in range(max_iter):  
        Ap = np.dot(A, p)  
        alpha = rsold / np.dot(p, Ap)  
        x += alpha * p  
        r -= alpha * Ap  
        rsnew = np.dot(r, r)  
        if np.sqrt(rsnew) < tol:  
            break  
        p = r + (rsnew / rsold) * p  
        rsold = rsnew  
  
    return x  
  
# 示例  
A = np.array([[4, 1], [1, 3]], dtype=float)  
b = np.array([1, 2], dtype=float)  
x0 = np.zeros_like(b)  
  
x_star = conjugate_gradient(A, b, x0)  
print("最优解向量 x:", x_star)

在这个示例中，我们定义了一个conjugate_gradient函数，用于求解对称正定线性系统 Ax = b。函数接受矩阵 A、右侧向量 b、初始解向量 x0、收敛阈值 tol 和最大迭代次数 max_iter 作为参数，并返回最优解向量 x。在函数内部，我们实现了共轭梯度法的迭代过程，包括计算搜索方向、步长和更新解向量等步骤。最后，我们通过一个简单的示例展示了如何使用该函数求解二次函数的最小值问题。