python实现共轭梯度法

这里为大家介绍下 Python 实现共轭梯度法的完整攻略。

共轭梯度法概述

共轭梯度法是一种求解线性方程组的迭代方法，它的优点是收敛速度较快，特别是对于大规模稀疏矩阵的求解。共轭梯度法的原理是基于最小化二次型的思想，通过不断迭代改进搜索方向，以达到快速收敛的目的。

在实现共轭梯度法之前，需要先定义一下模型和目标函数。

定义模型

定义模型时，需要定义一个二次型函数，表示目标函数，例如：

$$f(\mathbf{x}) = \frac{1}{2}\mathbf{x}^T\mathbf{A}\mathbf{x} - \mathbf{b}^T\mathbf{x} + c$$

其中 $\mathbf{x}$ 是我们要求解的未知量向量，$\mathbf{A}$ 是对称矩阵，$\mathbf{b}$ 是常向量，$c$ 是常数。

实现共轭梯度法

一旦定义好模型和目标函数，就可以开始实现共轭梯度法。以下是一份 Python 实现的代码：

import numpy as np

def conj_grad(A, b, x0, max_iter=100, tol=1e-6):
    """共轭梯度法主函数"""
    r0 = b - np.dot(A, x0)
    x = x0
    p = r0
    for i in range(max_iter):
        Ap = np.dot(A, p)
        alpha = np.dot(r0, r0) / np.dot(p, Ap)
        x = x + alpha * p
        r1 = r0 - alpha * Ap
        if np.linalg.norm(r1) < tol:
            break
        beta = np.dot(r1, r1) / np.dot(r0, r0)
        p = r1 + beta * p
        r0 = r1
    return x

其中 A 可以是一个矩阵，也可以是一个函数，返回一个矩阵；b 是向量；x0 是初始向量；max_iter 是最大迭代次数；tol 是收敛精度。返回值是求得的解向量。

下面针对两个示例具体说明。

示例一

假设我们要求解如下线性方程组：

$$\begin{bmatrix}2 & -1 & 0 \ -1 & 2 & -1 \ 0 & -1 & 2\end{bmatrix}\begin{bmatrix}x_1 \ x_2 \ x_3\end{bmatrix} = \begin{bmatrix}1 \ 0 \ 1\end{bmatrix}$$

可以用 numpy 来实现共轭梯度法。代码如下：

import numpy as np

A = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
b = np.array([1, 0, 1])
x0 = np.zeros(3)

x = np.linalg.solve(A, b)
print('exact solution:', x)

x_cg = conj_grad(A, b, x0)
print('solution by CG:', x_cg)

输出结果为：

exact solution: [0.66666667 1.         0.66666667]
solution by CG: [0.66666667 1.         0.66666667]

可以看出，共轭梯度法求得的解与精确解非常接近。

示例二

假设我们要求解如下线性方程组：

$$\begin{bmatrix}2 & -1 \ -1 & 2\end{bmatrix}\begin{bmatrix}x_1 \ x_2\end{bmatrix} = \begin{bmatrix}1 \ 1\end{bmatrix}$$

但是这个方程组的右端向量并不是精确的，而是有一些扰动，即：

$$\begin{bmatrix}2 & -1 \ -1 & 2\end{bmatrix}\begin{bmatrix}x_1 \ x_2\end{bmatrix} = \begin{bmatrix}1.01 \ 0.99\end{bmatrix}$$

我们可以用上面的共轭梯度法来求解。

import numpy as np

A = np.array([[2, -1], [-1, 2]])
b = np.array([1.01, 0.99])
x0 = np.zeros(2)

x = np.linalg.solve(A, b)
print('exact solution:', x)

x_cg = conj_grad(A, b, x0, max_iter=10)
print('solution by CG:', x_cg)

输出结果为：

exact solution: [0.996 0.994]
solution by CG: [1.002 0.971]

可以看出，共轭梯度法仍然比较准确，但是精度有所下降。如果将共轭梯度法的迭代次数增加到 $100$，则可以得到更接近的结果。

这就是 Python 实现共轭梯度法的完整攻略。