今天介绍数值计算和优化方法中非常有效的一种数值解法，共轭梯度法。我们知道，在解大型线性方程组的时候，很少会有一步到位的精确解析解，一般都需要通过迭代来进行逼近，而 PCG 就是这样一种迭代逼近算法。

我们先从一种特殊的线性方程组的定义开始，比如我们需要解如下的线性方程组：

Ax=b A x = b
<script type="math/tex; mode=display" id="MathJax-Element-37"> \mathbf{A} \mathbf{x} = \mathbf{b} </script>

这里的 A(n×n) A ( n × n ) <script type="math/tex" id="MathJax-Element-38"> \mathbf{A} (n \times n)</script> 是对称，正定矩阵， b(n×1) b ( n × 1 ) <script type="math/tex" id="MathJax-Element-39"> \mathbf{b} (n \times 1) </script> 同样也是已知的列向量，我们需要通过 A A <script type="math/tex" id="MathJax-Element-40">\mathbf{A}</script> 和 b b <script type="math/tex" id="MathJax-Element-41"> \mathbf{b}</script> 来求解 x(n×1) x ( n × 1 ) <script type="math/tex" id="MathJax-Element-42"> \mathbf{x} (n \times 1) </script>, 这其实是我们熟知的一些线性系统的表达式。

直接求解

首先，我们来看一种直观的解法，我们定义满足如下关系的向量为关于 矩阵 A A <script type="math/tex" id="MathJax-Element-43">\mathbf{A}</script> 的共轭向量，

uTAv=0 u T A v = 0
<script type="math/tex; mode=display" id="MathJax-Element-44"> \mathbf{u}^\mathsf{T} \mathbf{A} \mathbf{v} = 0 </script>

因为矩阵 A A <script type="math/tex" id="MathJax-Element-45">\mathbf{A}</script> 是对称正定矩阵，所以矩阵 A A <script type="math/tex" id="MathJax-Element-46">\mathbf{A}</script> 定义了一个内积空间：

⟨u,v⟩A:=⟨Au,v⟩=⟨u,ATv⟩=⟨u,Av⟩=uTAv ⟨ u , v ⟩ A := ⟨ A u , v ⟩ = ⟨ u , A T v ⟩ = ⟨ u , A v ⟩ = u T A v
<script type="math/tex; mode=display" id="MathJax-Element-47"> \langle \mathbf{u},\mathbf{v} \rangle_\mathbf{A} := \langle \mathbf{A} \mathbf{u}, \mathbf{v}\rangle = \langle \mathbf{u}, \mathbf{A}^\mathsf{T} \mathbf{v}\rangle = \langle \mathbf{u}, \mathbf{A}\mathbf{v} \rangle = \mathbf{u}^\mathsf{T} \mathbf{A} \mathbf{v} </script>

基于此，我们可以定义一组向量 P P <script type="math/tex" id="MathJax-Element-48"> P </script>

P={p1,…,pn} P = { p 1 , … , p n }
<script type="math/tex; mode=display" id="MathJax-Element-49"> P= \left \{\mathbf{p}_1, \dots, \mathbf{p}_n \right \} </script>

其中的向量 p1 p 1 <script type="math/tex" id="MathJax-Element-50">\mathbf{p}_1</script> , p2 p 2 <script type="math/tex" id="MathJax-Element-51"> \mathbf{p}_2</script>, … , pn p n <script type="math/tex" id="MathJax-Element-52">\mathbf{p}_n</script> 都是互为共轭的，那么 P P <script type="math/tex" id="MathJax-Element-53"> P </script> 构成了 Rn R n <script type="math/tex" id="MathJax-Element-54"> \mathbb{R}^{n} </script> 空间的一个基，上述方程的解 x∗ x ∗ <script type="math/tex" id="MathJax-Element-55"> \mathbf{x}_* </script> 可以表示成 P P <script type="math/tex" id="MathJax-Element-56"> P </script> 中向量的线性组合：

x∗=∑i=1nαipi x ∗ = ∑ i = 1 n α i p i
<script type="math/tex; mode=display" id="MathJax-Element-57"> \mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{p}_i </script>

根据上面的表达式，我们可以得到：

Ax∗=∑i=1nαiApipTkAx∗=∑i=1nαipTkApi(Multiply left by pTk)pTkb=∑i=1nαi⟨pk,pi⟩A(Ax∗=b and ⟨u,v⟩A=uTAv)⟨pk,b⟩=αk⟨pk,pk⟩A(uTv=⟨u,v⟩ and ∀i≠k:⟨pk,pi⟩A=0) A x ∗ = ∑ i = 1 n α i A p i p k T A x ∗ = ∑ i = 1 n α i p k T A p i (Multiply left by  p k T ) p k T b = ∑ i = 1 n α i ⟨ p k , p i ⟩ A ( A x ∗ = b  and  ⟨ u , v ⟩ A = u T A v ) ⟨ p k , b ⟩ = α k ⟨ p k , p k ⟩ A ( u T v = ⟨ u , v ⟩  and  ∀ i ≠ k : ⟨ p k , p i ⟩ A = 0 )
<script type="math/tex; mode=display" id="MathJax-Element-58"> \mathbf{A} \mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i \\ \mathbf{p}_k^\mathsf{T} \mathbf{A} \mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{p}_k^\mathsf{T} \mathbf{A} \mathbf{p}_i \quad \text{(Multiply left by } \mathbf{p}_k^\mathsf{T} \text{)} \\ \mathbf{p}_k^\mathsf{T} \mathbf{b} = \sum^{n}_{i=1} \alpha_i \left \langle \mathbf{p}_k, \mathbf{p}_i \right \rangle_{\mathbf{A}} \qquad (\mathbf{Ax_*} = \mathbf{b} \text{ and } \langle \mathbf{u},\mathbf{v} \rangle_\mathbf{A} = \mathbf{u}^\mathsf{T} \mathbf{A} \mathbf{v}) \\ \left \langle \mathbf{p}_k, \mathbf{b} \right \rangle =\alpha_k \left \langle \mathbf{p}_k, \mathbf{p}_k \right \rangle_{\mathbf{A}} \qquad (\mathbf{u}^\mathsf{T} \mathbf{v} = \left \langle \mathbf{u}, \mathbf{v} \right \rangle \text{ and } \forall i \neq k: \left \langle \mathbf{p}_k, \mathbf{p}_i \right \rangle_{\mathbf{A}} = 0 ) </script>
这意味着：

αk=⟨pk,b⟩⟨pk,pk⟩A α k = ⟨ p k , b ⟩ ⟨ p k , p k ⟩ A
<script type="math/tex; mode=display" id="MathJax-Element-59"> \alpha_k =\frac{\left \langle \mathbf{p}_k, \mathbf{b} \right \rangle}{\left \langle \mathbf{p}_k, \mathbf{p}_k \right \rangle_\mathbf{A}} </script>

所以，如果我们要直接求解的，可以先对矩阵 A A <script type="math/tex" id="MathJax-Element-60"> \mathbf{A} </script> 进行特征值分解，求出一系列的共轭向量，然后求出系数，最后可以得到方程的解 x∗ x ∗ <script type="math/tex" id="MathJax-Element-61"> \mathbf{x_*} </script>

迭代求解

上面的方法已经说明， x∗ x ∗ <script type="math/tex" id="MathJax-Element-62"> \mathbf{x}_* </script> 是一系列共轭向量 p p <script type="math/tex" id="MathJax-Element-63"> \mathbf{p} </script> 的线性组合，学过 PCA 的都知道，可以用前面占比高的向量组合进行逼近，而不需要把所有的向量都组合到一起，PCG 也是用到了这种思想，通过仔细的挑选共轭向量 p p <script type="math/tex" id="MathJax-Element-64"> \mathbf{p} </script> 来重建方程的解 x∗ x ∗ <script type="math/tex" id="MathJax-Element-65"> \mathbf{x_*} </script>。

我们先来看下面的一个方程：

f(x)=12xTAx−xTb,x∈Rn f ( x ) = 1 2 x T A x − x T b , x ∈ R n
<script type="math/tex; mode=display" id="MathJax-Element-66"> f(\mathbf{x}) = \tfrac12 \mathbf{x}^\mathsf{T} \mathbf{A}\mathbf{x} - \mathbf{x}^\mathsf{T} \mathbf{b}, \qquad \mathbf{x}\in\mathbf{R}^n </script>

对上面的方程求导，我们可以得到：

D2f(x)=A D 2 f ( x ) = A
<script type="math/tex; mode=display" id="MathJax-Element-67"> \mathrm{D}^2 f(\mathbf{x}) = \mathbf{A} </script>

Df(x)=Ax−b D f ( x ) = A x − b
<script type="math/tex; mode=display" id="MathJax-Element-68"> \mathrm{D} f(\mathbf{x}) = \mathbf{A} \mathbf{x} - \mathbf{b} </script>

可以看到，方程的一阶导数就是我们需要解的线性方程组，令一阶导数为 0，那么我们需要解的就是这样一个线性方程组了。

假设我们随机定义 x x <script type="math/tex" id="MathJax-Element-69"> \mathbf{x} </script> 的一个初始向量为 x0 x 0 <script type="math/tex" id="MathJax-Element-70"> \mathbf{x_0} </script>，那么我们可以定义第一个共轭向量为 p0=b−Ax0 p 0 = b − A x 0 <script type="math/tex" id="MathJax-Element-71"> \mathbf{p}_0 = \mathbf{b} - \mathbf{A} \mathbf{x}_0 </script>, 后续的基向量都是和梯度共轭的，所以称为共轭梯度法。

下面给出详细的算法流程：

这里写图片描述

而 preconditioned conjugate gradient method 与共轭梯度法的不同之处在于预先定义了一个特殊矩阵 M M <script type="math/tex" id="MathJax-Element-72"> \mathbf{M} </script>：

这里写图片描述

参考来源：wiki 百科

https://en.wikipedia.org/wiki/Conjugate_gradient_method#The_preconditioned_conjugate_gradient_method

Logo
DAMO开发者矩阵

DAMO开发者矩阵，由阿里巴巴达摩院和中国互联网协会联合发起，致力于探讨最前沿的技术趋势与应用成果，搭建高质量的交流与分享平台，推动技术创新与产业应用链接，围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐

cover

风电场智能巡检技术方案：无人机+机器人协同实践

avatar DAMO开发者矩阵

工业AI文献研究_industrial-ai-research

以下为本文档的中文说明工业 AI 文献研究技能，运行一个精简但充分的、源感知的学术研究工作流，专注于工业人工智能领域。它的覆盖范围包括预测性维护、智能调度、工业异常检测、智能制造、信息物理系统、工业自动化边缘 AI、机器人与工业交叉等方向，同时也处理数字孪生、工业物联网、工业 4.0、制造 AI 等相关术语。使用场景包括：研究人员需要系统性地调研某个工业 AI 子课题的文献现状时；工程师需要对比不

avatar DAMO开发者矩阵
cover

[论文学习]Unsafer in Many Turns：工具使用型Agent的多轮安全风险基准测试与防御

avatar DAMO开发者矩阵