CS205A Lecture 12 Conjugate gradients I; Gradient descent, setup

课程主页：https://graphics.stanford.edu/courses/cs205a-13-fall/schedule.html

这次回顾第十二讲，这一讲介绍了迭代法求解线性系统的思路，引出了共轭梯度法的思想。

优化算法

序列二次规划

回顾之前介绍的优化问题：

$$\begin{aligned} \text { minimize }& f(\vec{x}) & \\ \text { such that } & g(\vec{x})=\vec{0} \\ &h(\vec{x}) \geq \vec{0} \end{aligned}$$

现在讨论如何使用迭代算法求解，首先对目标函数和条件使用泰勒展开：

$$\begin{aligned} f(\vec x_{k+1})&\approx f(\vec x_k)+\frac{1}{2} \vec{d}^T H_{f}\left(\vec{x}_{k}\right) \vec{d}+\nabla f\left(\vec{x}_{k}\right) \cdot \vec{d}+f\left(\vec{x}_{k}\right) \\ g_{i}\left(\vec{x}_{k}\right)+\nabla g_{i}\left(\vec{x}_{k}\right) \cdot \vec{d}&=0 \\ h_{i}\left(\vec{x}_{k}\right)+\nabla h_{i}\left(\vec{x}_{k}\right) \cdot \vec{d} &\ge 0 \end{aligned}$$

所以可以用如下方式更新$\vec x_{k+1}$：

$$\begin{aligned} {\vec{x}_{k+1} = \vec{x}_{k}+\arg \min _{\vec{d}} \frac{1}{2} \vec{d}^T H_{f}\left(\vec{x}_{k}\right) \vec{d}+\nabla f\left(\vec{x}_{k}\right) \cdot \vec{d}+f\left(\vec{x}_{k}\right)} \\ {\text { such that } g_{i}\left(\vec{x}_{k}\right)+\nabla g_{i}\left(\vec{x}_{k}\right) \cdot \vec{d}=0} \\ {h_{i}\left(\vec{x}_{k}\right)+\nabla h_{i}\left(\vec{x}_{k}\right) \cdot \vec{d} \geq 0} \end{aligned}$$

如果只有等号约束，那么可以构造如下拉格朗日乘子：

$$\begin{aligned} \Lambda(\vec{d}, \vec{\lambda}) & = \frac{1}{2} \vec{d}^T H_{f}\left(\vec{x}_{k}\right) \vec{d}+\nabla f\left(\vec{x}_{k}\right) \cdot \vec{d}+f\left(\vec{x}_{k}\right)+\vec{\lambda}^T\left(g\left(\vec{x}_{k}\right)+D g\left(\vec{x}_{k}\right) .\vec{d}\right) \\ \Longrightarrow \vec{0}=\nabla_{\vec{d}} \Lambda &=H_{f}\left(\vec{x}_{k}\right) \vec{d}+\nabla f\left(\vec{x}_{k}\right)+\left[D g\left(\vec{x}_{k}\right)\right]^{T} \vec{\lambda} \end{aligned}$$

其中

$$D g\left(\vec{x}_{k}\right) =[\nabla g\left(\vec{x}_{1}\right), \ldots ,\nabla g\left(\vec{x}_{k}\right)]$$

结合等号约束，得到如下线性系统：

$$\left( \begin{array}{cc}{H_{f}\left(\vec{x}_{k}\right)} & {\left[D g \left(\vec{x}_{k}\right)\right]^{T}} \\ {D g\left(\vec{x}_{k}\right)} & {0}\end{array}\right) \left( \begin{array}{l}{\vec{d}} \\ {\vec{\lambda}}\end{array}\right)=\left( \begin{array}{c}{-\nabla f\left(\vec{x}_{k}\right)} \\ {-g\left(\vec{x}_{k}\right)}\end{array}\right)$$

Barrier Methods

另一种处理约束条件的方法是将约束项转化为能量项，例如

$$f_{\rho}(\vec{x})=f(\vec{x})+\rho\|g(\vec{x})\|_{2}^{2}$$

如果$\rho \to \infty$，那么$g(\vec x)\to 0$。

Chapter 10 Iterative Linear Solvers

回到线性系统的求解

$$A\vec x =\vec b$$

这一讲将介绍如何利用迭代算法求解线性系统。

讨论之前，这里做如下三个假设：

1. $A\in \mathbb R^{n\times n}$
2. $A$是对称矩阵，$A^T=A$
3. $A$是正定矩阵，即对$\vec x \neq 0$，$\vec x^T A\vec x >0$

梯度下降

在上述假设下，求解$A\vec x =\vec b$等价于最小化

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

特别的，对上式求梯度可得

$$\nabla f(\vec{x})=A \vec{x}-\vec{b}$$

令$\nabla f(\vec x) =\vec 0$即可得到

$$A \vec{x}=\vec{b}$$

回顾梯度下降算法：

1. $\vec{d}_{k} =-\nabla f\left(\vec{x}_{k-1}\right)=\vec{b}-A \vec{x}_{k-1}$
2. 定义$\vec{x}_{k} = \vec{x}_{k-1}+\alpha_{k} \vec{d}_{k}$，其中选择$\alpha_k$使得$f\left(\vec{x}_{k}\right)<f\left(\vec{x}_{k-1}\right)$

决定$\alpha_k ​$的值是一维“线性搜索”问题，特别的，在这里我们有

$$\begin{aligned} g(\alpha) &= f(\vec{x}+\alpha \vec{d}) \\ &=\frac{1}{2}(\vec{x}+\alpha \vec{d})^T A(\vec{x}+\alpha \vec{d})- \vec{b}^T(\vec{x}+\alpha \vec{d})+c \\ &=\frac{1}{2}\left(\vec{x}^T A \vec{x}+2 \alpha \vec{x}^T A \vec{d}+\alpha^{2} \vec{d}^T A \vec{d}\right)-\vec{b}^T \vec{x}-\alpha \vec{b}^T \vec{d}+c \\ &=\frac{1}{2} \alpha^{2} \vec{d}^T A \vec{d}+\alpha\left(\vec{x}^T A \vec{d}-\vec{b}^T\vec{d}\right)+\text{const} \end{aligned}$$

关于$\alpha$求导可得

$$\begin{aligned} \frac{d g(\alpha)}{d \alpha} &=\alpha \vec{d}^T A \vec{d}+\vec{x}^T A \vec{d}-\vec{b}^T\vec{d} \\ &=\alpha \vec{d}^T A \vec{d}+\vec{d}^TA \vec{x}-\vec{d}^T\vec{b} \\ &=\alpha \vec{d}^T A \vec{d}+\vec{d}^T(A \vec{x}-\vec{b}) \end{aligned}$$

令上式为$0$可得

$$\alpha=\frac{\vec{d}^T(\vec{b}-A \vec{x})}{\vec{d}^T A \vec{d}}$$

因为

$$\vec{d}_{k}=\vec{b}-A \vec{x}_{k}$$

所以

$$\alpha_{k}=\frac{\vec{d}_{k}^T \vec{d}_{k}}{\vec{d}_{k}^T A \vec{d}_{k}}$$

将上述内容总结，得到如下算法：

$$\begin{aligned} \vec{d}_{k} &=\vec{b}-A \vec{x}_{k-1} \\ \alpha_{k} &=\frac{\vec{d}_{k}^T \vec{d}_{k}}{\vec{d}_{k}^T A \vec{d}_{k}}\\ \vec{x}_{k} &=\vec{x}_{k-1}+\alpha_{k} \vec{d}_{k} \end{aligned}$$

收敛性

下面讨论梯度下降的收敛问题，特别的，这里考虑

$$R_{k} = \frac{f\left(\vec{x}_{k}\right)-f\left(\vec{x}'\right)}{f\left(\vec{x}_{k-1}\right)-f\left(\vec{x}'\right)}$$

其中

$$A \vec x' =\vec b$$

如果

$$R_k <\beta<1$$

那么该算法一定收敛。

现在利用

$$\vec x_k =\vec x_{k-1}+\alpha_k \vec d_k$$

计算$f(\vec x_k )$：

$$\begin{aligned} f\left(\vec{x}_{k}\right) &=f\left(\vec{x}_{k-1}+\alpha_{k} \vec{d}_{k}\right) \\ &=\frac{1}{2}\left(\vec{x}_{k-1}+\alpha_{k} \vec{d}_{k}\right)^T A\left(\vec{x}_{k-1}+\alpha_{k} \vec{d}_{k}\right)-\vec{b}^T\left(\vec{x}_{k-1}+\alpha_{k} \vec{d}_{k}\right)+c \\ &=f\left(\vec{x}_{k-1}\right)+\alpha_{k} \vec{d}_{k}^T A \vec{x}_{k-1}+\frac{1}{2} \alpha_{k}^{2} \vec{d}^T_{k} A \vec{d}_{k}-\alpha_{k} \vec{b}^T \vec{d}_{k} &由f的定义 \\ &=f\left(\vec{x}_{k-1}\right)+\alpha_{k} \vec{d}_{k}^T\left(\vec{b}-\vec{d}_{k}\right)+\frac{1}{2} \alpha_{k}^{2} \vec{d}_{k}^T A \vec{d}_{k}-\alpha_{k} \vec{b}^T \vec{d}_{k} & \vec{d}_{k}=\vec{b}-A \vec{x}_{k-1} \\ &= f\left(\vec{x}_{k-1}\right)-\alpha_{k} \vec{d}^T_{k} \vec{d}_{k}+\frac{1}{2} \alpha_{k}^{2} \vec{d}_{k}^{T} A \vec{d}_{k}\\ &=f\left(\vec{x}_{k-1}\right)-\frac{\vec{d}_{k}^T \vec{d}_{k}}{\vec{d}_{k}^T A \vec{d}_{k}} \vec{d}_{k}^T \vec{d}_{k}+\frac{1}{2} \left(\frac{\vec{d}_{k}^T \vec{d}_{k}}{\vec{d}_{k}^T A \vec{d}_{k}}\right)^2 \vec{d}_{k}^T A \vec{d}_{k} & \alpha_{k} =\frac{\vec{d}_{k}^T \vec{d}_{k}}{\vec{d}_{k}^T A \vec{d}_{k}}\\ &=f\left(\vec{x}_{k-1}\right)-\frac{(\vec{d}_{k}^T \vec{d}_{k})^2}{\vec{d}_{k}^T A \vec{d}_{k}} +\frac{1}{2} \frac{(\vec{d}_{k}^T \vec{d}_{k})^2}{\vec{d}_{k}^T A \vec{d}_{k}}\\ &= f\left(\vec{x}_{k-1}\right)-\frac{1}{2} \frac{(\vec{d}_{k}^T \vec{d}_{k})^2}{\vec{d}_{k}^T A \vec{d}_{k}} \end{aligned}$$

因此

$$\begin{aligned} R_{k} &= \frac{f\left(\vec{x}_{k}\right)-f\left(\vec{x}'\right)}{f\left(\vec{x}_{k-1}\right)-f\left(\vec{x}'\right)}\\ &=\frac{f\left(\vec{x}_{k-1}\right)-\frac{\left(d_{k}^T d_{k}\right)^{2}}{2 \vec{d}_{k}^T A_{k}}-f\left(\vec{x}'\right)}{f\left(\vec{x}_{k-1}\right)-f\left(\vec{x}'\right)}\\ &= 1-\frac{\left(\vec{d}_{k}^T \vec{d}_{k}\right)^{2}}{2 \vec{d}_{k}^T A \vec{d}_{k}\left(f\left(\vec{x}_{k-1}\right)-f\left(\vec{x}'\right)\right)} \end{aligned}$$

接下来计算

$$f\left(\vec{x}_{k-1}\right)-f\left(\vec{x}'\right)$$

带入$A\vec x’=\vec b​$可得

$$\begin{aligned} f\left(\vec{x}_{k-1}\right)-f\left(\vec{x}'\right) &=\left[\frac{1}{2} \vec{x}_{k-1}^T A \vec{x}_{k-1}-\vec{b}^T \vec{x}_{k-1}+c\right]-\left[\frac{1}{2}\left(\vec{x}'\right)^T A\vec x'-\vec{b}^T \vec{x}'+c\right] \\ &=\left[\frac{1}{2} \vec{x}_{k-1}^T A \vec{x}_{k-1}-\vec{b}^T \vec{x}_{k-1}+c\right]-\left[\frac{1}{2}\left(\vec{x}'\right)^T \vec{b}-\vec{b}^T \vec{x}'+c\right] & A\vec x'=\vec b\\ &=\frac{1}{2} \vec{x}_{k-1}^T A \vec{x}_{k-1}-\vec{b}^T \vec{x}_{k-1}+\frac{1}{2} \vec{b}^T A^{-1} \vec{b} & \vec x'=A^{-1}\vec b\\ &=\frac{1}{2}\left(A \vec{x}_{k-1}-\vec{b}\right)^T A^{-1}\left(A \vec{x}_{k-1}-\vec{b}\right) &A对称\\ &=\frac{1}{2} \vec{d}_{k} A^{-1} \vec{d}_{k} &\vec{d}_{k}的定义 \end{aligned}$$

代回$R_k$的表达式可得

$$\begin{aligned} R_{k} &=1-\frac{\left(\vec{d}_{k}^T \vec{d}_{k}\right)^{2}}{2 \vec{d}_{k}^T A \vec{d}_{k}\left(f\left(\vec{x}_{k-1}\right)-f\left(\vec{x}'\right)\right)} \\ &=1-\frac{\left(\vec{d}_{k}^T \vec{d}_{k}\right)^{2}}{\vec{d}_{k}^T A \vec{d}_{k} \cdot \vec{d}_{k}^T A^{-1} \vec{d}_{k}} \\ &=1-\frac{\vec{d}_{k}^T \vec{d}_{k}}{\vec{d}_{k}^T A \vec{d}_{k}} \cdot \frac{\vec{d}_{k}^T \vec{d}_{k}}{\vec{d}_{k}^T A^{-1} \vec{d}_{k}} \\ &\le 1-\left(\min _{\|\vec{d}\|=1} \frac{1}{\vec{d}^T A \vec{d}}\right)\left(\min _{\|\vec d{\| }=1} \frac{1}{\vec{d}^T A^{-1} \vec{d}}\right)\\ &=1-\left(\max _{ \| \vec{d} \|=1} \vec{d}^T A \vec{d}\right)^{-1}\left(\max _{\|\vec{d}\|=1} \vec{d}^T A^{-1} \vec{d}\right)^{-1}\\ &=1-\frac{\sigma_{\mathrm{min}}}{\sigma_{\mathrm{max}}} & \sigma_{\max},\sigma_{\min}分别为A最大和最小的奇异值 \\ &=1-\frac{1}{\text{cond} A} \end{aligned}$$

所以梯度下降的收敛速度取决于$A$的条件数，下图展示了两种情形：

（注意因为$\text{cond} A\ge 1$，所以收敛性是显然的）

共轭梯度法

前几章介绍的求解$A\vec x =\vec b $的算法都是在$O(n^3)$时间内计算出结果，梯度下降法每一轮需要的时间为$O(n^2)$（矩阵乘法），如果迭代超过$n$轮，那么时间就会超过$O(n^3)$，这一讲介绍的共轭梯度可以保证在$n$步内收敛，是一种非常强大的算法。

动机

对

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

稍作变形

$$\begin{aligned} f(\vec{x}) &=\frac{1}{2} \vec{x}^T A \vec{x}-\vec{b}^T \vec{x}+c \\ &=\frac{1}{2}\left(\vec{x}-\vec{x}'\right)^T A\left(\vec{x}-\vec{x}'\right)+\vec{x}^T A \vec{x}'-\frac{1}{2}\left(\vec{x}'\right)^T A \vec{x}'-\vec{b}^T \vec{x}+c \\ &=\frac{1}{2}\left(\vec{x}-\vec{x}'\right)^T A\left(\vec{x}-\vec{x}'\right)+\vec{x}^T \vec{b}-\frac{1}{2}\left(\vec{x}'\right)^T \vec{b}-\vec{b}^T \vec{x}+c&A \vec{x}'=\vec{b} \\ &=\frac{1}{2}\left(\vec{x}-\vec{x}'\right)^T A\left(\vec{x}-\vec{x}'\right)+\text { const} \end{aligned}$$

所以最小化$f(\vec x)$等价于最小化$\vec x$到$\vec x’$的$A$范数：

$$\left(\vec{x}-\vec{x}'\right)^T A\left(\vec{x}-\vec{x}'\right)$$

为叙述方便，这里给出如下定义：

$$\begin{aligned} \| \vec v \|_A^2 &=\vec v^T A\vec v \\ \langle\vec{u}, \vec{v}\rangle&=\left(L^T \vec{u}\right) \cdot\left(L^T\vec{v}\right)=\vec{u}^{T} A \vec{v} \end{aligned}$$

注意这里$A$是对称正定矩阵，所以可以使用Cholesky分解：

$$A= LL^T$$

因此

$$\begin{aligned} f(\vec{x}) &=\frac{1}{2}\left(\vec{x}-\vec{x}'\right)^T A\left(\vec{x}-\vec{x}'\right)+\text {const}\\ &= \frac{1}{2}\left(\vec{x}-\vec{x}'\right)^T LL^T\left(\vec{x}-\vec{x}'\right)+\text {const}\\ &= \|L^T\left(\vec{x}-\vec{x}'\right) \|_2^2 +\text {const} \end{aligned}$$

如果定义

$$\vec y =L^T\vec x , \vec y ' =L^T\vec x'$$

那么

$$\bar f(\vec y) =f(\vec x) =\|\vec{y}-\vec{y}' \|_2^2 +\text {const}$$

如果我们计算出$L$，那么优化上述问题将非常简单，后面将介绍在不计算出$L$的情形下快速求解上述问题。

在讨论之前，注意我们有如下命题：

命题10.1

$$\{\vec w_1,\ldots ,\vec w_n\}是\mathbb R^n中正交基，那么通过在\vec w_1,\ldots ,\vec w_n上的线性搜索，最多经过n次即可最小化\bar f$$

这个结论是显然的，假设最小值点为

$$\vec x ' =\sum_{i=1}^n \alpha_i \vec w_i$$

那么只要依次计算出$\alpha_i$即可。

运用上述思路求解之前介绍的问题，这里只要考虑$L^T \vec x$即可，假设

$$\vec w_i =L^T \vec x_i$$

为正交基，那么对于$i\neq j$，我们有

$$0=\vec{w}_{i} \cdot \vec{w}_{j}=\left(L^T \vec{v}_{i}\right)^T\left(L^T\vec{v}_{j}\right)=\vec{v}_{i}^T\left(L L^T\right) \vec{v}_{j}=\vec{v}_{i}^T A \vec{v}_{j}$$

于是引出如下定义：

定义10.2

$$如果\vec v^T A\vec w =0，那么称\vec v, \vec w为A共轭$$

那么结合之前讨论，我们有如下命题：

命题10.2

$$\{\vec v_1,\ldots ,\vec v_n\}为A共轭，那么通过在\vec v_1,\ldots ,\vec v_n上的线性搜索，最多经过n次即可最小化 f$$