基本概念
共轭:
Q 是一个对称实矩阵,对于方向向量d1,d2,…dm,∀i≠j,dTiQdj=0 则他们关于 Q 共轭Q正定
如果对于矩阵Q,Q>0 ,若一组向量 d1,d2,…dm,m≤n−1 ,关于 Q 共轭,则他们线性无关。(直接用定义可以证明)- 内积空间
定义内积<di,dj>=dTiQdj 则这是一个欧氏空间, 即他为我们提供了求线性无关组的共轭方向的算法。并且两组向量关于 Q,Q>0 则他们在如上定义的内积空间中正交
基本共轭方向算法
对于
n
维二次型函数
给定一组初始点 xx0,d0,d1,…,dn−1 , 其中 di 关于 Q 共轭则有如下迭代法则:
可以证明,对于 n 维二次型可以在
下面证明一个非常强的引理
引理1
在共轭方向算法中所有的
k,0≤k≤n−1,0≤i≤k
,都有
由以上引理可以证明任意的
αk
,满足
设 ϕk(αk)=ϕk(xxk+αkddk)
由于 ϕk 是关于 αk 的凸二次函数, ϕk 有唯一极小值点
共轭梯度算法
- 选择 x0 , 计算 d0,g0
- 计算 αk=−ggTkddkddTkQddk
- 计算 xxk+1=xxk+αkddk
- 计算 gk 判断是否停止
- 计算 βk=ggTk+1QddkddTkQddk
- 计算 dk+1=−gk+1+βkdk
可以利用归纳法证明,搜索方向
d0,d1,…,dn−1
是共轭方向
首先证明
dT0Qd1=0
非二次型问题中的共轭梯度
由于二阶导数计算相当费时,而
Q
的计算只在
alpha
αk=f(xxk+αddk) 展开一维搜索
beta
- Hestenes-Stiefel 公式
xk+1两侧同时乘上Q减去b,得gk+1Qdk=xk+αkdk=gk+αkQdk因此=gk+1−gkαk
替换 beta 中的 Qdk 得到
βk=gTk+1[gk+1−gk]dTk[gk+1−gk] - P-R 公式
将H-S中的分母展开由引理 gTk+1dk=0 ,在迭代等式 dk=−gk+βk−1dk−1 左乘 gTk ,可得
gTkdk=−gTkgk+βk−1gTkdk−1=−gTkgk
带入得
βk=gTk+1[gk+1−gk]gTkgk - F-R修正
将P-R公式分子展开由引理的
βk=gTk+1gk+1gTkgk
Powell 证明用F-R公式计算 β 共轭梯度算法性能比较突出.
代码实现
def cg_gradient(fun, grad, x0, args=(), g_args=(), tol=1e-8, max_iter=5000):
alpha = lambda a, x_k, d: fun(*((x_k + a * d,) + args))
g0 = grad(*((x0,) + g_args))
d0 = -g0
for _ in range(max_iter):
a_k = minimize_scalar(alpha, bounds=(0, 100), args=(x0, d0), tol=1e-4)
x0 = x0 + a_k.x * d0
g_k = grad(*((x0,) + g_args))
if is_stop(g_k, np.zeros(g_k.shape), tol):
break
beta = np.sum(g_k ** 2) / np.sum(g0 ** 2) # Fletcher-Reeves 公式
g0 = g_k
d0 = -g_k + beta * d0
if _ % (len(x0) + 5) == 0:
d0 = -g_k
return OptimizeResult({'x': x0, 'fun': fun(*((x0,) + args)), 'jac': grad(*((x0,) + g_args)), 'nit': max_iter - _})
代码实现中有一个trick 即当迭代次数达到一定的次数后会将其重设为梯度的复方向
下面是在rosen函数上实验的结果
cg res
fun: 1.5458689649160323e-22
jac: array([ -1.10458309e-11, -6.90558721e-12])
nit: 4983
x: array([ 1., 1.])
迭代次数只有梯度下降的一半
完整代码https://github.com/DylanFrank/optimize
转载链接http://blog.csdn.net/Dylan_Frank/article/details/78270326