【优化理论】共轭梯度下降算法实现

实验目的

实验步骤

此次实验要求采用共轭梯度下降来解决二次优化问题，在解决此问题前，首先分析下什么是共轭梯度法，共轭梯度法能解决什么问题。
我们来看一个线性方程组Ax=b，求解此方程组的过程可以看成是形如公式(1)的优化问题：
arg⁡min⁡x∥Ax−b∥2…………………………………………(1)\mathop{\arg\min_{x}} \parallel Ax-b\parallel^2…………………………………………(1) argxmin∥Ax−b∥2…………………………………………(1)
其又可转化为公式(2)：
∥Ax−b∥2=(Ax−b)T(Ax−b)=xTATAx−2bTAx+∥b∥2………(2)\parallel Ax-b\parallel ^2 =(Ax-b)^T (Ax-b)=x^T A^T Ax-2b^T Ax+\parallel b\parallel ^2………(2) ∥Ax−b∥2=(Ax−b)T(Ax−b)=xTATAx−2bTAx+∥b∥2………(2)
令H=ATA,g=bTA，H正定，H=A^T A,g=b^T A，H正定，H=ATA,g=bTA，H正定，
则上述优化问题等价于公式(3):
arg⁡min⁡xf(x)=12xTHx−gTx……………………………………(3)\mathop{\arg\min_{x}}f(x)=\frac{1}{2} x^T Hx-g^T x……………………………………(3) argxminf(x)=21xTHx−gTx……………………………………(3)
如此一来，求解线性方程组转化为二次优化问题，那么如何求解公式(3)的优化问题呢？让我们先来看一些定义和性质。
定义1. 给定一个正定矩阵?，我们可以定义?−内积：
⟨x,y⟩H=xTHy,∀x,y∈Rn……………………………………(4)⟨x,y⟩_H=x^T Hy,∀x,y∈R^n……………………………………(4) ⟨x,y⟩H=xTHy,∀x,y∈Rn……………………………………(4)
定义2. 设 H 是一个n×n的正定矩阵，x,y∈R^n，如果⟨x,y⟩_H=0,则称 x 与 y 是H-共轭的。
性质1. 设非零向量d_1,d_2,…,d_n 两两H-共轭，则它们线性无关。
通过定义1、2和性质1，我们可以得出如下结论：
空间任意向量x可以用这组向量基表示:x=∑i=1nαidi………………(5)空间任意向量 x可以用这组向量基表示:x= ∑_{i=1}^nα_i d_i………………(5)空间任意向量x可以用这组向量基表示:x=i=1∑nαidi………………(5)
将公式(5)代入公式(3)可得：
gTx=α1gTd1+α2gTd2+⋯+αngTdn…………………………(6)xTHx=∑i=1nαi2⟨di,di⟩H=∑i=1nαi2diTHdi…………………………(7)f(x)=gTx+12xTHx=∑i=1n(12αi2diTHdi−αigTdi)…………………(8)g^T x=α_1 g^T d_1+α_2 g^T d_2+⋯+α_n g^T d_n…………………………(6) \\ x^T Hx=∑_{i=1}^nα_i^2 ⟨d_i,d_i ⟩_H =∑_{i=1}^nα_i^2 d_i^T Hd_i …………………………(7) \\ f(x)=g^T x+\frac{1}{2}x^T Hx=∑_{i=1}^n(\frac{1}{2}α_i^2 d_i^T Hd_i-α_i g^T d_i ) …………………(8) gTx=α1gTd1+α2gTd2+⋯+αngTdn…………………………(6)xTHx=i=1∑nαi2⟨di,di⟩H=i=1∑nαi2diTHdi…………………………(7)f(x)=gTx+21xTHx=i=1∑n(21αi2diTHdi−αigTdi)…………………(8)
现在的目标就是求公式(8)的最小值，如果参与求和的每一项都取最小值，则f (x)达到最小值，所以对公式(8)中求和的每一项进行求导可得解向量的第i个分量：
αi=gTdidiTHdi,i=1,2,3,…,n…………………………………(9)α_i=\frac {g^T d_i}{d_i^T Hd_i }, i=1,2,3,…,n…………………………………(9) αi=diTHdigTdi,i=1,2,3,…,n…………………………………(9)
注意：上面的解法中求每一个α_i是独立的步骤，只用到一个向量d_i，因此我们并不需要一开始就有n个共轭向量，只需要一个向量就可以启动，而每一次都能求解得到解向量的一个分量，所以理论上当循环次数至多到n次(n指H的维度)时，梯度必降为0，求解完毕。

基本实现原理

对于本题具体的求解算法如图1所示：

图 1 共轭梯度法具体算法算法图中循环终止条件是梯度降为0，符合之前理论上的分析，但在实际计算过程中，由于舍入误差或其他噪声常常导致算法无法有限步收敛，因此停止条件修改为公式(10)：

∥η∥=∥∇f(x)∥<ε.…………………………………(10)\parallelη\parallel=\parallel∇f(x)\parallel<ε.…………………………………(10)∥η∥=∥∇f(x)∥<ε.…………………………………(10)

实现代码

代码如下，也可观看我GitHub。

# -*- coding: utf-8 -*-
"""
Created on Sat Nov 17 13:39:12 2018 @author: YLC
"""
import numpy as np
x = np.array([0,0,0,0]).T #.T表示转置，下同
H = np.array([[158,20,90,101],[20,36,46,61],[90,46,306,156],[101,61,156,245]])
g = np.array([8,-5,1,6]).T
def grad(H,x,g): #梯度计算公式，由原方程求导得到  return np.dot(H,x)-g
eta = grad(H,x,g) #梯度
d = -eta #梯度方向
i = 1 #迭代次数
while(np.linalg.norm(eta,ord=2) > 1e-10):  alpha = -np.dot(eta.T,d)/np.dot(np.dot(d.T,H),d)  x = x + np.dot(alpha,d)  eta = grad(H,x,g)  d = -eta + np.dot(np.dot(np.dot(eta.T,H),d)/np.dot(np.dot(d.T,H),d),d)  #print("========================================")  #print("迭代第"+str(i)+"次||eta||的值为:",np.linalg.norm(eta,ord=2))      #print("迭代第"+str(i)+"次alpha的值为:\n",alpha)  #print("迭代第"+str(i)+"次eta的值为:\n",eta)  #print("迭代第"+str(i)+"次d的值为:\n",d)      print("迭代第"+str(i)+"次x的值为:\n",x)  i = i + 1

实验结果及分析

实验结果如下：
迭代第1次x的值为:
[ 0.03675345 -0.0229709 0.00459418 0.02756508]
迭代第2次x的值为:
[ 0.07945184 -0.09736974 -0.03609274 0.03706019]
迭代第3次x的值为:
[ 0.05512625 -0.28123105 0.00293889 0.06041221]
迭代第4次x的值为:
[ 0.03590246 -0.30049431 -0.00770042 0.08940927]
从实验结果可以看出，整个算法迭代了四次结束，符合前面理论的分析，H为4*4的矩阵，因而n=4，迭代4次结束。在具体实现的过程中要注意阈值ε的取值，不同于前一个实验，这里的ε要非常小，其次用python做矩阵点积预算的时候是dot函数，而不是之间用乘号。

发现与收获

初学时，共轭梯度法较于简单的梯度下降更难以理解，对于数学上的公式、定义、性质要都有了解才好理解。这里的共轭梯度法，共轭指的不是复数的共轭，而是全新的定义H-共轭，有了H-共轭才有了线性无关的向量组d，从而转化为参数的最优化问题。同时，更应该掌握的是共轭梯度法的运用，其根本目的是为了解线性方程组，这说明不同的问题可以采取不同的梯度下降手段，从而降低算法的时空开销，提高算法性能。