CG 和 PCG 算法实现和分析

实验理论

共轭梯度法

共轭梯度法用于求解一类特殊但十分常见的线性方程组:

$$ \boldsymbol{Ax=b} $$

其中 $\boldsymbol{A}\in \mathbb{R}^{N\times N}$,$\boldsymbol{x,b}\in \mathbb{R}^{N}$。$\boldsymbol{A}$ 是给定的对称正定(Symmetric Positive Definite,SPD)的实系数方阵,$\boldsymbol{b}$ 是已知的向量,$\boldsymbol{x}$ 是需要求解的向量。

现构造一个关于向量 $\boldsymbol{x}$ 的二次型:

$$ f(\boldsymbol{x})=\frac{1}{2}\boldsymbol{x}^\mathrm{T}\boldsymbol{Ax}-\boldsymbol{b}^\mathrm{T}\boldsymbol{x}+c \tag{1} \label{eq:quadradic} $$

其中 $c$ 是标量常数。由于:

$$ \frac{\partial f}{\partial\boldsymbol{x}}=\boldsymbol{Ax-b} \tag{2} \label{eq:gradient} $$

因此,使得函数 $f(\boldsymbol{x})$ 取得最小值的 $\boldsymbol{x}^\star$,就是方程 $\boldsymbol{Ax=b}$ 的解[1]。现就将原本的求解线性方程的问题,转化为一个优化问题:

$$ \boldsymbol{x}^\star=\mathop{\arg\min}_{\boldsymbol{x}}\left(\frac{1}{2}\boldsymbol{x}^\mathrm{T}\boldsymbol{Ax}-\boldsymbol{b}^\mathrm{T}\boldsymbol{x}+c\right) \tag{3} \label{eq:optim} $$

下面使用迭代的思想解决公式 \ref{eq:optim} 所描述的优化问题。首先给向量 $\boldsymbol{x}$ 随意取一个值 $\boldsymbol{x}_0$,由公式 \ref{eq:gradient} 可知,$\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}_0$ 是函数 $f(\boldsymbol{x})$ 在 $\boldsymbol{x}_0$ 处的梯度方向。根据最朴素的思想,$\boldsymbol{x}$ 每一次沿着梯度的方向下降,最终会收敛到最优值上。这就是最速下降(Gradient Descent)法的基本思想。

由于最速下降法每次迭代中只考虑了在该点的梯度值,所以可能会出现「一叶障目,不见泰山」般的盲目性。在实际计算中,最速下降法的迭代次数往往很大,导致算法性能不尽如人意。共轭梯度法为了克服这个问题,在每一步迭代中不是简单地选取梯度方向作为下山方向,而是在由该步上的梯度方向及其共轭方向(也就是相互正交的方向)所长成的「二维平面」上找到新的最优下山方向[2]。

可以证明的是,如果每一次根据一对共轭方向构造出来的下山方向都是最优的,那么当前步的梯度方向的共轭方向,正是上一步的下山方向。以下是共轭梯度法的推导过程:

给定初始向量 $\boldsymbol{x}_0$,第一步选取的下山方向 $\boldsymbol{p}_0$ 仍是梯度方向 $\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}_0$,之后第 $k+1$ 步($k\ge 1$),就在第 $k$ 步的下山方向 $\boldsymbol{p}_{k-1}$ 和第 $k+1$ 步的梯度方向 $\boldsymbol{r}_{k}$ 所长成的平面上找最优的下山方向,即第 $k+1$ 步的下山方向 $\boldsymbol{p}_{k}$ 为:

$$ \boldsymbol{p}_{k}=\boldsymbol{x}_{k} + \xi\boldsymbol{r}_{k} + \eta\boldsymbol{p}_{k-1} $$

其中 $\xi, \eta \in \mathbb{R}$.

当然 $\boldsymbol{p}_k$ 也不能任意取,需要考虑到公式 \ref{eq:quadradic} 二次型 $f(\boldsymbol{x})$ 的限制。将 $\boldsymbol{p}_k$ 带入 $f(\boldsymbol{x})$,构造一个新函数 $\psi(\xi, \eta)$:

$$ \begin{aligned} \psi(\xi, \eta) =\, &f(\boldsymbol{x}_{k} + \xi\boldsymbol{r}_{k} + \eta\boldsymbol{p}_{k-1})\\
=\, &\frac{1}{2}(\boldsymbol{x}_{k} + \xi\boldsymbol{r}_{k} + \eta\boldsymbol{p}_{k-1})^\mathrm{T}\boldsymbol{A}(\boldsymbol{x}_{k} + \xi\boldsymbol{r}_{k} + \eta\boldsymbol{p}_{k-1})\\
&-\boldsymbol{b}^\mathrm{T}(\boldsymbol{x}_{k} + \xi\boldsymbol{r}_{k} + \eta\boldsymbol{p}_{k-1})+c \end{aligned} $$

$$\begin{aligned} \frac{\partial\psi}{\partial\xi} =\, &\boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{A}(\boldsymbol{x}_{k} + \xi\boldsymbol{r}_{k} + \eta\boldsymbol{p}_{k-1}) - \boldsymbol{b}^\mathrm{T}\boldsymbol{r}_{k}\\
=\, &\boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{A}\boldsymbol{x}_{k} + \xi\boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{A}\boldsymbol{r}_{k} + \eta\boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{A}\boldsymbol{p}_{k-1} - \boldsymbol{b}^\mathrm{T}\boldsymbol{r}_{k}\\
=\, &\boldsymbol{r}_{k}^\mathrm{T}(\boldsymbol{b} - \boldsymbol{r}_k) + \xi\boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{A}\boldsymbol{r}_{k} + \eta\boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{A}\boldsymbol{p}_{k-1} - \boldsymbol{b}^\mathrm{T}\boldsymbol{r}_{k}\\
=\, &-\boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{r}_k + \xi\boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{A}\boldsymbol{r}_{k} + \eta\boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{A}\boldsymbol{p}_{k-1}\\
\frac{\partial\psi}{\partial\eta} =\, &\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}(\boldsymbol{x}_{k} + \xi\boldsymbol{r}_{k} + \eta\boldsymbol{p}_{k-1}) - \boldsymbol{b}^\mathrm{T}\boldsymbol{p}_{k-1}\\
=\, &\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}\boldsymbol{x}_{k} + \xi\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}\boldsymbol{r}_{k} + \eta\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}\boldsymbol{p}_{k-1} - \boldsymbol{b}^\mathrm{T}\boldsymbol{p}_{k-1}\\
=\, &\boldsymbol{p}_{k-1}^\mathrm{T}(\boldsymbol{b} - \boldsymbol{r}_k) + \xi\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}\boldsymbol{r}_{k} + \eta\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}\boldsymbol{p}_{k-1} - \boldsymbol{b}^\mathrm{T}\boldsymbol{p}_{k-1}\\
=\, &-\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{r}_{k} + \xi\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}\boldsymbol{r}_{k} + \eta\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}\boldsymbol{p}_{k-1}\\
=\, &\xi\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}\boldsymbol{r}_{k} + \eta\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}\boldsymbol{p}_{k-1}\\
\end{aligned} $$

令 $\frac{\partial\psi}{\partial\xi}=\frac{\partial\psi}{\partial\eta}=0$,可解得极小值点 $\tilde{\boldsymbol{x}}=\boldsymbol{x}_{k} + \tilde\xi\boldsymbol{r}_{k} + \tilde\eta\boldsymbol{p}_{k-1}$

即:

$$ \begin{cases} \tilde\xi\boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{A}\boldsymbol{r}_{k} + \tilde\eta\boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{A}\boldsymbol{p}_{k-1} - \boldsymbol{r}_{k}^\mathrm{T}\boldsymbol{r}_k=0\\
\tilde\xi\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}\boldsymbol{r}_{k} + \tilde\eta\boldsymbol{p}_{k-1}^\mathrm{T}\boldsymbol{A}\boldsymbol{p}_{k-1}=0 \end{cases} \tag{4} \label{eq:zero} $$

由于 $\boldsymbol{r}_k\ne 0$ (因为 $\boldsymbol{r}_k$ 是残差,若为 0 则应该停止迭代了),所以 $\tilde\xi\ne 0$,因此可取:

$$ \boldsymbol{p}_k=\frac{1}{\tilde\xi}(\tilde{\boldsymbol{x}}-\boldsymbol{x}_k)=\boldsymbol{r}_k+\frac{\tilde\eta}{\tilde\xi}\boldsymbol{p}_{k-1} $$

作为新的下山方向。令 $\beta_{k-1}=\frac{\tilde\eta}{\tilde\xi}$,则可由公式 \ref{eq:zero} 得:

$$ \beta_{k-1}=-\frac{\boldsymbol{r}_k^{\mathrm{T}}\boldsymbol{Ap}_{k-1}}{\boldsymbol{p}_{k-1}^{\mathrm{T}}\boldsymbol{Ap}_{k-1}} \tag{5} \label{eq:beta_k_1} $$

至此就找到了新的下山方向,向量 $\boldsymbol{x}$ 将每次沿着下山方向,以 $\alpha$ 为步长下降,即 $\boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k$,将其带入 $f(\boldsymbol{x})$ 可得:

$$ \begin{aligned} \phi(\alpha_k) &= f(\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k)\\
&= \frac{1}{2}(\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k)^{\mathrm{T}}\boldsymbol{A}(\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k)-\boldsymbol{b}^{\mathrm{T}}(\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k) + c\\
&= \frac{1}{2}\alpha^2\boldsymbol{p}_k^{\mathrm{T}}\boldsymbol{Ap}_k- \alpha\boldsymbol{r}_k^{\mathrm{T}}\boldsymbol{p}_k+f(\boldsymbol{x}_k) \end{aligned} $$

令 $\phi^\prime(\alpha_k)=0$ 可解得:

$$\begin{aligned} \alpha_k &= \frac{\boldsymbol{r}_k^{\mathrm{T}}\boldsymbol{p}_k}{\boldsymbol{p}_k^{\mathrm{T}}\boldsymbol{Ap}_k}\\
&= \frac{\boldsymbol{r}_k^{\mathrm{T}}(\boldsymbol{r}_k+\beta_{k-1}\boldsymbol{p}_{k-1})}{\boldsymbol{p}_k^{\mathrm{T}}\boldsymbol{Ap}_k}\\
&= \frac{\boldsymbol{r}_k^{\mathrm{T}}\boldsymbol{r}_k}{\boldsymbol{p}_k^{\mathrm{T}}\boldsymbol{Ap}_k} \end{aligned} \tag{6} \label{eq:alpha} $$

根据 $\boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k$ 重新整理 $\boldsymbol{r}_{k+1}$ 可得:

$$\begin{aligned} \boldsymbol{r}_{k+1} &= \boldsymbol{b}-\boldsymbol{Ax}_{k+1}\\
&= \boldsymbol{b}-\boldsymbol{A}(\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k)\\
&= \boldsymbol{r}_k-\alpha_k\boldsymbol{A}\boldsymbol{p}_k\\
\boldsymbol{A}\boldsymbol{p}_k &= \frac{1}{\alpha}(\boldsymbol{r}_k-\boldsymbol{r}_{k+1}) \end{aligned} \tag{7} \label{eq:rk} $$

将公式 \ref{eq:rk} 带入公式 \ref{eq:beta_k_1} 得:

$$ \begin{aligned} \beta_{k} &= -\frac{\boldsymbol{r}_{k+1}^{\mathrm{T}}\boldsymbol{Ap}_{k}}{\boldsymbol{p}_{k}^{\mathrm{T}}\boldsymbol{Ap}_{k}}\\
&= -\frac{\frac{1}{\alpha}\boldsymbol{r}_{k+1}^{\mathrm{T}}(\boldsymbol{r}_k-\boldsymbol{r}_{k+1})}{\frac{1}{\alpha}\boldsymbol{p}_{k}^{\mathrm{T}}(\boldsymbol{r}_k-\boldsymbol{r}_{k+1})}\\
&= \frac{\boldsymbol{r}_{k+1}^{\mathrm{T}}\boldsymbol{r}_{k+1}}{\boldsymbol{p}_{k}^{\mathrm{T}}\boldsymbol{r}_k}\\
&= \frac{\boldsymbol{r}_{k+1}^{\mathrm{T}}\boldsymbol{r}_{k+1}}{\boldsymbol{r}_{k}^{\mathrm{T}}\boldsymbol{p}_k}\\
&= \frac{\boldsymbol{r}_{k+1}^{\mathrm{T}}\boldsymbol{r}_{k+1}}{\boldsymbol{r}_{k}^{\mathrm{T}}(\boldsymbol{r}_k+\beta_{k-1}\boldsymbol{p}_{k-1})}\\
&= \frac{\boldsymbol{r}_{k+1}^{\mathrm{T}}\boldsymbol{r}_{k+1}}{\boldsymbol{r}_{k}^{\mathrm{T}}\boldsymbol{r}_k} \end{aligned} $$

由此可得到共轭梯度法的一般过程:首先向量 $\boldsymbol{x}$ 任意取一初值 $\boldsymbol{x}_0$,由 $\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}_0$ 得到残差 $\boldsymbol{r}_0$,第一步的下山方向 $\boldsymbol{p}_0=\boldsymbol{r}_0$,根据公式 \ref{eq:alpha} 更新向量 $\boldsymbol{x}$:

$$ \begin{aligned} \boldsymbol{x}_1 &= \boldsymbol{x}_0+\alpha_0\boldsymbol{p}_0\\
&= \boldsymbol{x}_0+\frac{\boldsymbol{r}_0^{\mathrm{T}}\boldsymbol{r}_0}{\boldsymbol{p}_0^{\mathrm{T}}\boldsymbol{Ap}_0}\boldsymbol{p}_0 \end{aligned} $$

之后的第 $k+1$ 步($k\ge 1$),其上一步的残差 $\boldsymbol{r}_{k}$、新下山方向 $\boldsymbol{p}_{k}$、新的向量 $\boldsymbol{x}_{k+1}$ 分别为:

$$ \begin{aligned} \boldsymbol{r}_{k} &= \boldsymbol{r}_{k-1}-\frac{\boldsymbol{r}_{k-1}^{\mathrm{T}}\boldsymbol{r}_{k-1}}{\boldsymbol{p}_{k-1}^{\mathrm{T}}\boldsymbol{Ap}_{k-1}}\boldsymbol{A}\boldsymbol{p}_{k-1}\\
\boldsymbol{p}_{k} &= \boldsymbol{r}_{k}-\frac{\boldsymbol{r}_{k}^{\mathrm{T}}\boldsymbol{r}_{k}}{\boldsymbol{r}_{k-1}^{\mathrm{T}}\boldsymbol{r}_{k-1}}\boldsymbol{p}_{k-1}\\
\boldsymbol{x}_{k+1} &= \boldsymbol{x}_k+\frac{\boldsymbol{r}_k^{\mathrm{T}}\boldsymbol{r}_k}{\boldsymbol{p}_k^{\mathrm{T}}\boldsymbol{Ap}_k}\boldsymbol{p}_{k} \end{aligned} $$

预优共轭梯度法

共轭梯度法本身已经具有非常好的性能,然而在一些极端情况下(例如系数矩阵 $\boldsymbol{A}$ 的条件数很大),迭代的次数还是很多。预优共轭梯度法是在共轭梯度法的基础上,先使用一个预优子(Preconditioner)对系数矩阵 $\boldsymbol{A}$ 进行预优,用一个性质更好且很接近 $\boldsymbol{A}$ 的矩阵来代替它,以达到减少迭代次数的目的。

令 $\tilde{\boldsymbol{A}}=\boldsymbol{C}^{-1}\boldsymbol{AC}^{-1}$,$\tilde{\boldsymbol{x}}=\boldsymbol{Cx}$,$\tilde{\boldsymbol{b}}=\boldsymbol{C}^{-1}\boldsymbol{b}$,同时令 $\tilde{\boldsymbol{x}}_k=\boldsymbol{Cx}_k$,$\tilde{\boldsymbol{r}}_k=\boldsymbol{C}^{-1}\boldsymbol{r}_k$,$\tilde{\boldsymbol{p}}_k=\boldsymbol{Cp}_k$,并记 $\boldsymbol{M}=\boldsymbol{C}^2$,根据前述的共轭梯度法的一般过程,以下是预优共轭梯度法的一般过程:

首先向量 $\boldsymbol{x}$ 任意取一初值 $\boldsymbol{x}_0$,由 $\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}_0$ 得到残差 $\boldsymbol{r}_0$,第一步的下山方向为方程 $\boldsymbol{Mz}_0=\boldsymbol{r}_0$的解 $\boldsymbol{z}_0$,即 $\boldsymbol{p}_0=\boldsymbol{z}_0$,根据公式 \ref{eq:alpha} 更新向量 $\boldsymbol{x}$:

$$ \begin{aligned} \boldsymbol{x}_1 &= \boldsymbol{x}_0+\alpha_0\boldsymbol{p}_0\\
&= \boldsymbol{x}_0+\frac{\boldsymbol{r}_0^{\mathrm{T}}\boldsymbol{z}_0}{\boldsymbol{p}_0^{\mathrm{T}}\boldsymbol{Ap}_0}\boldsymbol{p}_0 \end{aligned} $$

之后的第 $k+1$ 步($k\ge 1$),其上一步的残差 $\boldsymbol{r}_{k}$、新下山方向 $\boldsymbol{p}_{k}$、新的向量 $\boldsymbol{x}_{k+1}$ 分别为:

$$ \begin{aligned} \boldsymbol{z}_{k} &= \boldsymbol{M}^{-1}\boldsymbol{r}_{k}\\
\boldsymbol{r}_{k} &= \boldsymbol{r}_{k-1}-\frac{\boldsymbol{r}_{k-1}^{\mathrm{T}}\boldsymbol{z}_{k-1}}{\boldsymbol{p}_{k-1}^{\mathrm{T}}\boldsymbol{Ap}_{k-1}}\boldsymbol{A}\boldsymbol{p}_{k-1}\\
\boldsymbol{p}_{k} &= \boldsymbol{z}_{k}-\frac{\boldsymbol{r}_{k}^{\mathrm{T}}\boldsymbol{z}_{k}}{\boldsymbol{r}_{k-1}^{\mathrm{T}}\boldsymbol{z}_{k-1}}\boldsymbol{p}_{k-1}\\
\boldsymbol{x}_{k+1} &= \boldsymbol{x}_k+\frac{\boldsymbol{r}_k^{\mathrm{T}}\boldsymbol{z}_k}{\boldsymbol{p}_k^{\mathrm{T}}\boldsymbol{Ap}_k}\boldsymbol{p}_{k} \end{aligned} $$

现在的问题就剩下如何寻找到一个条件良好的矩阵 $\boldsymbol{M}$。比较常见的选择有对角预优矩阵和不完全 Cholesky 因子预优矩阵。具体来说,系数矩阵 $\boldsymbol{A}$ 的对角预优矩阵 $\boldsymbol{M}$ 可取:

$$ \boldsymbol{M} = \mathop{\mathrm{diag}}(\boldsymbol{A}_{11},\ldots,\boldsymbol{A}_{NN}) $$

对于不完全 Cholesky 因子预优矩阵,是先求系数矩阵 $\boldsymbol{A}$ 的不完全 Cholesky 分解:

$$ \boldsymbol{A}=\boldsymbol{LL}^{\mathrm{T}}+\boldsymbol{R} $$

其中 $\boldsymbol{L}$ 是下三角矩阵,然后用 $\boldsymbol{M}=\boldsymbol{LL}^{\mathrm{T}}$ 作为预优矩阵。

实验过程

本实验中的算法使用 Julia 语言实现,Julia 语言的一个最显著的特性是其可以使用 Unicode 字符作为变量名以及运算符,因此其代码的形式看起来非常接近数学公式。除此之外,Julia 语言的语法与 Matlab 语言非常相似,也内置了丰富的用于矩阵运算的函数,同时 Julia 语言的性能表现十分优异,非常适合用于数值计算方面的算法实现练习。

共轭梯度法

根据前述结论,共轭梯度法的代码实现如下:

"""
    solve_cg(A, b[, kₘ=1000, ϵ=1e-6])

使用共轭梯度法计算线性方程组 Ax = b.

# 输入参数
- `A`: 系数矩阵 A.
- `b`: 已知向量 b.
- `kₘ=1000`: 可选参数,最大迭代次数.
- `ϵ=1e-6`: 可选参数,迭代终止阈值条件.

"""
function solve_cg(A, b, kₘ=1000, ϵ=1e-6)
    n = size(A)[1]
    x = zeros(n, 1)
    k = 0
    r = b - A * x
    ρ = dot(r', r)
    p = r
    while √ρ > ϵ * norm(b) && k ≤ kₘ
        w = A * p
        α = ρ / dot(p', w)
        x = x + α * p
        r = r - α * w

        ρ̃  = ρ
        ρ = dot(r', r)
        β = ρ / ρ̃
        p = r + β * p
        k = k + 1
    end
    return x, steps
end

正确性验证

算法的正确性使用与正确结果的二范数来衡量。具体来说,对于线性方程 $\boldsymbol{Ax}=\boldsymbol{b}$,令 $\boldsymbol{x}_1=\boldsymbol{A}^{-1}\boldsymbol{b}$,共轭梯度法计算的结果为 $\boldsymbol{x}_2$,算法计算的结果与正确答案的结果的误差定义为:

$$ \mathit{loss} = \left\Vert \boldsymbol{x}_1 - \boldsymbol{x}_2 \right\Vert_2 $$

以下是正确性验证和迭代次数计算的代码:

steps = 100
loss = zeros(steps)
iterations = zeros(steps)
times = zeros(steps)
for i = 1:steps
    A = rand_spd(500)
    b = rand_mat(500, 1)
    x₁ = A \ b
    (x₂, it), t = @timed solve_cg(A, b)
    loss[i] = norm(x₁ - x₂)
    iterations[i] = it
    times[i] = t
end
println("[CG] Average time in $(steps) steps: $(mean(times))")
println("[CG] Average loss in $(steps) steps: $(mean(loss))")
println("[CG] Average iter in $(steps) steps: $(mean(iterations))")

以下是代码的执行结果:

[CG] Average time  in 100 steps: 0.044025289690000004
[CG] Average loss  in 100 steps: 1.7706899123596664e-6
[CG] Average iter  in 100 steps: 173.94

从结果中可以看到,共轭梯度法的计算结果与正确答案误差的数量级为 $10^{-6}$,与程序中设置的迭代终止阈值条件的数量级相一致,显示了算法的正确性。

预优共轭梯度法

根据前述结论,不完全 Cholesky 因子的代码实现[3]如下:

function ichol(A)
    a = copy(A)
    n = size(a)[1]
    for k = 1:n
        a[k, k] = √a[k, k]
        for i = (k + 1):n
            if a[i, k] ≠ 0
                a[i, k] = a[i, k] / a[k, k]
            end
        end
        for j = (k + 1):n
            for i = j:n
                if a[i, j] ≠ 0
                    a[i, j] = a[i, j] - a[i, k] * a[j, k]
                end
            end
        end
    end
    LowerTriangular(a)
end

预优共轭梯度法的代码实现如下:

"""
    solve_pcg(A, b[, preconditioner=ichol_preconditioner, kₘ=1000, ϵ=1e-6])

使用预优共轭梯度法计算线性方程组 Ax = b.

# 输入参数
- `A`: 系数矩阵 A.
- `b`: 已知向量 b.
- `preconditioner=ichol_preconditioner`: 可选参数,预优子产生函数,
                                         默认为不完全 Cholesky 因子预优子.
- `kₘ=1000`: 可选参数,最大迭代次数.
- `ϵ=1e-6`: 可选参数,迭代终止阈值条件.

"""
function solve_pcg(A, b, preconditioner=ichol_preconditioner, kₘ=1000, ϵ=1e-6)
    n = size(A)[1]
    x = zeros(n, 1)
    k = 0
    r = b - A * x
    M = preconditioner(A)
    M′ = M ^ -1
    p, ρ̃  = 0, 0
    while √dot(r', r) > ϵ * norm(b) && k < kₘ
        z = M′ * r
        ρ = dot(r', z)
        k = k + 1
        if k == 1
            p = z
        else
            β = ρ / ρ̃
            p = z + β * p
        end
        w = A * p
        α = ρ / dot(p', w)
        x = x + α * p
        r = r - α * w
        ρ̃  = ρ
    end
    return x, k
end

正确性验证

预优共轭梯度法的正确性验证的方法和过程与共轭梯度法类似,以下是测试代码的输出结果:

[PCG With Diag Preconditioner] Average time in 100 steps: 0.04717120826000001
[PCG With Diag Preconditioner] Average loss in 100 steps: 1.798713732524729e-6
[PCG With Diag Preconditioner] Average iter in 100 steps: 173.53
[PCG With ICho Preconditioner] Average time in 100 steps: 0.04692745026000001
[PCG With ICho Preconditioner] Average loss in 100 steps: 7.356997898159087e-14
[PCG With ICho Preconditioner] Average iter in 100 steps: 1.0

可以看到,无论是使用对角预优因子还是不完全 Cholesky 预优因子的预优共轭梯度法的误差的数量级都在 $10^{-6}$ 以下,其中使用不完全 Cholesky 预优因子的更是达到了 $10^{-14}$,原因是该算法在绝大多数情况下只迭代了一次就满足迭代终止阈值条件了,而不是一个逐渐靠近然后跨越的过程,所以其误差会小很多。上述的实验结果也显示了使用两种不同的预优子的预优共轭梯度算法的正确性。

算法运行时间分析

下表是共轭梯度法以及使用两种不同预优子的预优共轭梯度法计算 100 个相同的线性方程所用的时间统计:

CG PCG With Diag PCG With ICho
Average Time (ms) 46.40 47.92 48.91

表中的「Average Time」代表解 100 个线性方程组的平均耗时。从表中可以看到,共轭梯度法比使用两种不同预优子的预优共轭梯度法都快,原因是预优共轭梯度法每次迭代需要额外求解一个线性方程组得到 $\boldsymbol{z}_k$,这个步骤十分耗时,同时这也不利于并行化,只有在系数矩阵 $\boldsymbol{A}$ 性质很不好的时候更能体现预优共轭梯度法的性能优势。

算法迭代步数分析

在共轭梯度法及其衍生算法中,在固定的终止阈值条件下,算法的迭代次数与实际所解的方程的系数矩阵的条件数有关。对于共轭梯度法,其系数矩阵就是原始的 $\boldsymbol{A}$,对于预优共轭梯度法,其实际所解的系数矩阵为 $\tilde{\boldsymbol{A}}=\boldsymbol{M}^{-1}\boldsymbol{A}$.

下表是共轭梯度法以及使用两种不同预优子的预优共轭梯度法计算 100 个相同的线性方程所用的时间和相应矩阵的条件数统计。

CG PCG With Diag PCG With ICho
Average Iterations 174.06 173.96 1.0
Avg Condition Num 657.18 656.69 1.0000000000003393

表中的「Average Iterations」代表解 100 个线性方程组的平均迭代步数,「Avg Condition Num」代表相应矩阵的平均条件数。从表中可以看出,使用对角预优子的预优共轭梯度法带来的条件数衰减量很少,几乎没有改变,表现在迭代次数上就是使用对角预优子的预优共轭梯度法的迭代次数与原始的共轭梯度法的十分接近。而对于使用不完全 Cholesky 预优因子的预优共轭梯度法,其对应的系数矩阵的条件数有了显著的降低,表现在算法的迭代次数上就是其平均迭代次数为 1,这是一个非常好的结果。

至于为什么使用了不完全 Cholesky 预优因子对系数矩阵 $\boldsymbol{A}$ 进行预优了之后, $\tilde{\boldsymbol{A}}=\boldsymbol{M}^{-1}\boldsymbol{A}$ 的条件数会十分接近于 1,以下给出一个简单(不严谨)的说明:

系数矩阵 $\boldsymbol{A}$ 的不完全 Cholesky 分解的形式为:

$$ \boldsymbol{A}=\boldsymbol{LL}^{\mathrm{T}}+\boldsymbol{R} $$

其中 $\boldsymbol{L}$ 为下三角矩阵,同时 $\boldsymbol{M}=\boldsymbol{LL}^{\mathrm{T}}$.

由于预条件处理后的矩阵 $\tilde{\boldsymbol{A}}=\boldsymbol{M}^{-1}\boldsymbol{A}$ 满足以下关系[4]:

$$ \begin{aligned} \tilde{\boldsymbol{A}} &= \boldsymbol{M}^{-1}\boldsymbol{A}\\
&= \boldsymbol{M}^{-1}(\boldsymbol{M} + \boldsymbol{R})\\
&= \boldsymbol{I} - \boldsymbol{M}^{-1}\boldsymbol{R} \end{aligned} $$

由于系数矩阵 $\boldsymbol{A}$ 的不完全 Cholesky 分解是用 $\boldsymbol{LL}^{\mathrm{T}}$ 来近似矩阵 $\boldsymbol{A}$,矩阵 $\boldsymbol{R}$ 可以看作是「残差」,因此矩阵 $\boldsymbol{R}$ 的元素应该都比较小(这样才能近似)。所以 $\tilde{\boldsymbol{A}}=\boldsymbol{I} - \boldsymbol{M}^{-1}\boldsymbol{R}$ 的结果就会非常接近于单位矩阵,其条件数也就十分接近于 1.

实验结论

本实验对共轭梯度法及预优共轭梯度法的算法理论进行了详细的梳理,使用 Julia 语言实现了共轭梯度法,同时分别使用对角预优子和不完全 Cholesky 预优子对系数矩阵进行预优,最后分析对比了经过预优后的算法和原始的 CG 算法,对实验的结果进行了详细的说明和解释,验证了预优算法的有效性。

参考文献

  1. 徐树方, 高立, 张平文. 数值线性代数, 2000.

  2. SHEWCHUK J R, OTHERS. An introduction to the conjugate gradient method without the agonizing pain, 1994.

  3. ANONYMOUS. Incomplete Cholesky factorization. 2018. https://en.wikipedia.org/wiki/Incomplete_Cholesky_factorization.

  4. 费建中. 一种带参数的不完全 Cholesky 分解共轭梯度法. 计算数学, 1988, 1 : 005.