共轭梯度法的推导与完整算法

共轭梯度法

学习自知乎：https://www.zhihu.com/question/27157047 and wikipedia and 非线性规划课

简介

在数值线性代数中，共轭梯度法是一种求解对称正定线性方程组Ax=b的迭代方法。

事实上，求解Ax=b等价于求解： m i n ∣ ∣ A x − b ∣ ∣ 2 2 min||Ax-b||_2^2 min∣∣Ax−b∣∣22​ ,将其展开后可以得到： m i n x T A T A x − b T A x + b T b min \quad x^TA^TAx-b^TAx+b^Tb minxTATAx−bTAx+bTb ,也就是等价于求解 m i n 1 2 x T A T A x − b T A x min\quad \frac{1}{2}x^TA^TAx-b^TAx min21​xTATAx−bTAx 。于是解方程问题就转化为了求解二次规划问题(QP)。

共轭梯度法是介于梯度下降法与牛顿法之间的一个方法，是一个一阶方法。它克服了梯度下降法收敛慢的缺点，又避免了存储和计算牛顿法所需要的二阶导数信息。

在n维的优化问题中，共轭梯度法最多n次迭代就能找到最优解（是找到，不是接近），但是只针对二次规划问题。

共轭梯度法的思想就是找到n个两两共轭的共轭方向，每次沿着一个方向优化得到该方向上的极小值，后面再沿其它方向求极小值的时候，不会影响前面已经得到的沿哪些方向上的极小值，所以理论上对n个方向都求出极小值就得到了n维问题的极小值。

算法推导过程

目标函数的标准形式：

min ⁡ x ∈ R n 1 2 x T Q x − b T x \min_{x\in R^n}\frac{1}{2}x^TQx-b^Tx minx∈Rn​21​xTQx−bTx

Q-conjugate: 对于正定矩阵Q，如果非零向量x,y是Q-conjugate的，那么

x T Q y = 0 x^TQy=0 xTQy=0

我们需要找到n个相互Q-conjugate的基向量 d 1 , d 2 , … , d n {d_1,d_2,…,d_n} d1​,d2​,…,dn​，它们相互共轭且线性无关。

因此空间中任意向量x都可以用这组基向量表示：

x = ∑ i = 1 n a i d i x=\sum_{i=1}^na_id_i x=∑i=1n​ai​di​

因此我们的目标函数可以改写为：
min ⁡ a 1 , . . . , a n ∈ R n 1 2 ( ∑ i = 1 n a i d i ) T Q ( ∑ j = 1 n a j d j ) − b T ( ∑ i = 1 n a i d i ) = min ⁡ a 1 , . . . , a n ∈ R n 1 2 ∑ i = 1 n ∑ j = 1 n a i a j d i T Q d j − ∑ i = 1 n a i b T d i \min_{a_1,...,a_n\in R^n}\frac{1}{2}(\sum_{i=1}^na_id_i)^TQ(\sum_{j=1}^na_jd_j)-b^T(\sum_{i=1}^na_id_i) \\ =\min_{a_1,...,a_n\in R^n}\frac{1}{2}\sum_{i=1}^n\sum_{j=1}^na_ia_jd_i^TQd_j-\sum_{i=1}^na_ib^Td_i a1​,...,an​∈Rnmin​21​(i=1∑n​ai​di​)TQ(j=1∑n​aj​dj​)−bT(i=1∑n​ai​di​)=a1​,...,an​∈Rnmin​21​i=1∑n​j=1∑n​ai​aj​diT​Qdj​−i=1∑n​ai​bTdi​
因为 d 1 , d 2 , … d n {d_1,d_2,…d_n} d1​,d2​,…dn​是Q-conjugate，所以由于 d i Q d j = 0 , i f i ≠ j d_iQd_j=0\quad,if\quad i\ne j di​Qdj​=0,ifi​=j，上式可以继续简化为
min ⁡ a 1 , . . . , a n ∈ R n 1 2 ∑ i = 1 n a i 2 d i T Q d i − ∑ i = 1 n a i b T d i \min_{a_1,...,a_n\in R^n}\frac{1}{2}\sum_{i=1}^na_i^2d_i^TQd_i-\sum_{i=1}^na_ib^Td_i a1​,...,an​∈Rnmin​21​i=1∑n​ai2​diT​Qdi​−i=1∑n​ai​bTdi​
将求和符号提出来后，可以得到
min ⁡ a 1 , . . . , a n ∈ R n 1 2 ∑ i = 1 n ( a i 2 d i T Q d i − a i b T d i ) \min_{a_1,...,a_n\in R^n}\frac{1}{2}\sum_{i=1}^n(a_i^2d_i^TQd_i-a_ib^Td_i) a1​,...,an​∈Rnmin​21​i=1∑n​(ai2​diT​Qdi​−ai​bTdi​)
原本我们的问题是寻找一个向量x，使得目标函数值最优。现在我们用 a 1 , a 2 , … a n a_1,a_2,…a_n a1​,a2​,…an​来表示 x = ( x 1 , x 2 , … . , x n ) x=(x_1,x_2,….,x_n) x=(x1​,x2​,….,xn​)，也就是我们的任务变为优化 a 1 , a 2 , … , a n a_1,a_2,…,a_n a1​,a2​,…,an​了。

为了看的清楚点，我们把求和分开来写：
min ⁡ a 1 , . . . , a n ∈ R n 1 2 ( a 1 2 d 1 T Q d 1 + a 1 b T d 1 ) + 1 2 ( a 2 2 d 2 T Q d 2 + a 2 b T d 2 ) + . . . + 1 2 ( a n 2 d n T Q d n + a n b T d n ) \min_{a_1,...,a_n\in R^n}\frac{1}{2}(a_1^2d_1^TQd_1+a_1b^Td_1)+\frac{1}{2}(a_2^2d_2^TQd_2+a_2b^Td_2)+...+\frac{1}{2}(a_n^2d_n^TQd_n+a_nb^Td_n) a1​,...,an​∈Rnmin​21​(a12​d1T​Qd1​+a1​bTd1​)+21​(a22​d2T​Qd2​+a2​bTd2​)+...+21​(an2​dnT​Qdn​+an​bTdn​)
这样，我们分别求第一项 1 2 ( a 1 2 d 1 T Q d 1 + a 1 b T d 1 ) \frac{1}{2}(a_1^2d_1^TQd_1+a_1b^Td_1) 21​(a12​d1T​Qd1​+a1​bTd1​)的最小值，然后第二项，即可。

第一项是一个二次函数，求最小值我们直接求导即可：

a 1 = b T d 1 d 1 T Q d 1 a_1=\frac{b^Td_1}{d_1^TQd_1} a1​=d1T​Qd1​bTd1​​ ，同样可得：

a i = b T d i d i T Q d i a_i=\frac{b^Td_i}{d_i^TQd_i} ai​=diT​Qdi​bTdi​​

所以最优解 x ∗ = ∑ i = 1 n a i d i x^*=\sum_{i=1}^na_id_i x∗=∑i=1n​ai​di​，将上式代入可得

x ∗ = ∑ i = 1 n b T d i d i T Q d i d i x^*=\sum_{i=1}^n\frac{b^Td_i}{d_i^TQd_i}d_i x∗=∑i=1n​diT​Qdi​bTdi​​di​

现在，我们只需要构建Q-conjugate向量组 d 1 , d 2 , … , d n {d_1,d_2,…,d_n} d1​,d2​,…,dn​。有一种称为Gram-Schmidt过程的算法可以用来找到这组向量。在实际算法实现中，我们往往是边求 a i a_i ai​,边求 d i d_i di​。下面给出完整的共轭梯度算法。

完整算法

输入：正定矩阵Q（A），初始向量 x 0 x_0 x0​，向量b

输出：最优解对应的向量 x k + 1 x_{k+1} xk+1​

$r_0:=b-Qx_0 $ (r_i为第i次迭代的误差)

d 0 : = r 0 d_0:=r_0 d0​:=r0​ (d_i是我们要求的共轭向量)

k : = 0 k:=0 k:=0 (k表示第几次迭代)

repeat

   α k : = r k T r k d k T Q d k \alpha_k:=\frac{r_k^Tr_k}{d_k^TQd_k} αk​:=dkT​Qdk​rkT​rk​​ (该项为学习率，是求出来的，对应的是之前说的a_i)

   x k + 1 : = x k + α k d k x_{k+1}:=x_k+\alpha_kd_k xk+1​:=xk​+αk​dk​

   r k + 1 : = r k − α k Q d k r_{k+1}:=r_k-\alpha_kQd_k rk+1​:=rk​−αk​Qdk​

  如果 r k + 1 r_{k+1} rk+1​足够小，则提前退出循环 （也就是认为已经找到最优解了）

   β k : = r k + 1 T r k + 1 r k T r k \beta_k:=\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k} βk​:=rkT​rk​rk+1T​rk+1​​

   d k + 1 : = r k + 1 + β k d k d_{k+1}:=r_{k+1}+\beta_kd_k dk+1​:=rk+1​+βk​dk​ （Gram-Schmidt过程求d_k）

  k:=k+1

end repeat

The result is x k + 1 x_{k+1} xk+1​