简单的共轭梯度法的Python代码示例代码例子 共轭梯度法优缺点

共轭梯度法

最速下降法以及牛顿法都具有其自身的局限性。本文将要介绍的共轭梯度法是介于最速下降法与牛顿法之间的一种无约束优化算法，它具有超线性的收敛速度，而且算法结构简单，容易编程实现。此外，根最速下降法相类似，共轭梯度法只用到了目标函数及其梯度值，避免了二阶导数的计算，从而降低了计算量和存储量，因此它是求解无约束优化问题的一种比较有效而使用的算法。

一、共轭方向发
共轭方向法的基本思想是在求解$n$维正定二次目标函数极小点时产生一组共轭方向作为搜索方向，在线搜索条件下算法之多迭代$n$步即能求得极小点。京故宫适当修正后共轭方向法可以推广到求解一般非二次目标函数情形。下面先介绍共轭方向的概念。
定义1：设$G$是$n$阶对称正定矩阵，若$n$维向量组$d_1,d_2,\cdots,d_m(m\le n)满足$$d_i^TGd_j=0,i\neq j$，则乘$d_1,d_2,\cdots,d_m$是$G$共轭的。
显然，向量组的共轭是正交的推广，即当$G=I$（单位阵）时，上述定义变成了向量组正交的定义。此外，不难证明，对称正定矩阵$G$的共轭向量组必然是线性无关的。
下面我们考虑求解正定二次目标函数极小点的共轭方向法。设
$\min f(x)=\frac{1}{2}x^TGx+b^Tx+c\tag{1}$
其中$G$是$n$阶对称正定阵，$b$为$n$为常向量，$c$为常数。我们有下面的算法：

---

算法1：共轭方向法
0. 给定迭代精度$0\le\epsilon\le 1$和初始点$x_0$。计算$g_0=\nabla f(x_0)$。选取初始方向$d_0$使得$d_0^Tg_0\lt 0$。令$k=0$.

1. 若$||g_k||\le\epsilon$，停算，输出$x^*\approx x_k$.
2. 利用线搜索方法确定步长$\alpha_k$
3. 令$x_{k+1}=x_k+\alpha_kd_k$，并计算$g_{k+1}=\nabla f(x_{k+1})$.
4. 选取$d_{k+1}$满足下降性和共轭性条件：$d_{k+1}^Tg_{k+1}\lt 0,d_{k+1}^TGd_i=0,\quad,i=0,1,\cdots,k$.
5. 令$k=k+1$，转步1

---

该算法的收敛性证明略过，如果感兴趣，可以去查找相应的数值优化专著。这里直接给出结论，在精确线搜索下，算法1求解正定二次目标函数极小化问题，之多在$n$步内即可求得其唯一极小点。这种能在有限步内求得二次函数极小点的性质通常称为二次终止性。

二、共轭梯度法
共轭梯度法是在每一步迭代利用当前点处的最速下降方向来生成关于凸二次函数$f$的海塞阵$G$的共轭方向，并建立求$f$在$\mathbb{R^n}$上的极小点的方法。这一方法最早是由Hesteness和Stiefel于1952年为求解正定线性方程组而提出来的，后经Fletcher等人研究并应用于无约束优化问题取得了丰富的成果，共轭梯度法也因此成为当前求解无约束优化问题的重要算法类。
设函数如（1）式所定义，则$f$的梯度和海塞矩阵为
$g(x)=\nabla f(x)=Gx+b,\quad G(x)=\nabla^2 f(x)=G\tag{2}$
下面我们讨论算法（1）中共轭方向的构造。我们取初始方向$d_0$为初始点$x_0$处的负梯度方向，即
$d_0=-\nabla f(x_0)=-g_0 \tag{3}$
从$x_0$出发沿$d_0$方向进行线搜索得到步长$\alpha_0$，令
$x_1=x_0+\alpha_0d_0$
其中$\alpha_0$满足条件
$\nabla f(x_1)^Td_0=g_1^Td_0 \tag{4}$
在$x_1$处，用$f$在$x_1$的负梯度方向$-g_1$与$d_0$的组合来生成$d_1$，即
$d_1=-g_1+\beta_0d_0 \tag{5}$
然后选取系数$\beta_0$使$d_1$与$d_0$关于G共轭，即令
$d_1^TGd_0 = 0 \tag{6}$
来确定$\beta_0$，将（6）代入（5）得
$\beta_0 = \frac{g_1^TGd_0}{d_0^TGd_0} \tag{7}$
由（2）得
$g_1-g_0=G(x_1-x_0)=\alpha_0Gd_0 \tag{8}$
故由（3）~（5）可得
$g_2^Tg_0=0,g_2^Tg_1=0,d_0^Tg_0=-g_0^Tg_0,d_1^Tg_1=-g_1^Tg_1$
现假设已得到互相共轭得搜索方向$d_1,d_2,\cdots,d_{k-1}$，精确线搜索得到得步长为$\alpha_0,\alpha_1,\cdots,\alpha_{k-1}$，且满足
$\left\{ \begin{array}{rcl} d_{k-1}^TGd_i=0, &i=0,1,\cdots,k-2,\\ d_i^Tg_i=-g_i^Tg_i,&i=0,1,\cdots,k-1,\\ g_k^Tg_i=0,g_k^Td_i=0,&i=0,1,\cdots,k-1. \end{array} \right. \tag{9}$
现令
$d_k=-g_k+\beta_{k-1}d_{k-1}+\sum_{i=0}^{k-1}\beta_{k}^{(i)}d_i\tag{10}$
其中$\beta_{k-1},\beta_k^{(i)}(i=0,1,\cdots,k-2)$得选择要满足
$d_k^TGd_i=0,i=0,1,\cdots,k-1\tag{11}$
用$d_i^TG(i=0,1,\cdots,k-1)$左乘（10）得
$\beta_{k-1}=\frac{g_k^TGd_{k-1}}{d_{k-1}^TGd_{k-1}},\beta_k^{(1)}=\frac{g_k^TGd_i}{d_i^TGd_i},i=0,1,\cdots,k-2\tag{12}$
类似于（8），我们有
$g_{i+1}-g_i=G(x_{i+1}-x_i)=\alpha_iGd_i,i=0,1,\cdots,k-1$
及
$\alpha_iGd_i=g_{i+1}-g_i,i=0,1,\cdots,k-1\tag{13}$
于是由归纳法假设（9）可得
$\beta_k^{(i)}=\frac{g_k^TGd_i}{d_i^TGd_i}=\frac{g_k^T(g_{i+1}-g_i)}{d_i^T(g_{i+1}-g_i)}=0,i=0,1,\cdots,k-2.$
于是，第k步得搜索方向为
$d_k=-g_k+\beta_{k-1}d_{k-1},\tag{14}$
其中$\beta_{k-1}$由（12）确定，即
$\beta_{k-1}=\frac{g_k^TGd_{k-1}}{d_{k-1}^TGd_{k-1}}\tag{15}$
同时有$d_k^Tg_k=-g_k^Tg_k$。这样确定了一组由负梯度方向形成得共轭方向，而把沿着这组方向进行迭代得方向称为共轭梯度法。其证明过程这里略过。
下面我们给出共轭梯度法求解无约束优化问题（1）极小点得算法步骤

---

算法2：共轭梯度法
0. 给定迭代精度$0\le\epsilon\lt 1$和初始点$x_0$。计算$g_0=\nabla f(x_0)$。令$k=0$

1. 若$||g_k||\le\epsilon$，停算，输出$x^)\approx x_k$
2. 计算搜索方向$d_k$:
   KaTeX parse error: Unknown column alignment: s at position 29: …\begin{array}{ls̲c} -g_k,&k=0,\\…
   其中当$k\ge 1$时，$\beta_{k-1}$由（15）确定
3. 利用线搜索方法确定搜索步长$\alpha_k$
4. 令$x_{k+1}=x_k+\alpha_kd_k$，并计算$g_{k+1}=\nabla f(x_{k+1})$
5. 令$k=k+1$，转步1

---

计算公式（15）是由Fletcher和Reeves给出得，故称之为FR公式，算法2也称之为FR共轭梯度法。除FR工诗外，尚有下列著名公式：
$\begin{aligned} \beta_k&=\frac{g_{k+1}^Tg_{k+1}}{-d_k^Tg_k}, (Dixon公式) \\ \beta_k&=\frac{g_{k+1}^Tg_{k+1}}{d_k^T(g_{k+1}-g_k)}, (Dai-Yuan公式) \\ \beta_k&=\frac{g_{k+1}^T(g_{k+1}-g_k)}{d_k^T(g_{k+1}-g_k)}, (Crowder-Wolfe公式) \\ \beta_k&=\frac{g_{k+1}^T(g_{k+1}-g_k)}{g_k^Tg_k}, (Polak,Ribiere,Ployak,PRP公式) \\ \end{aligned}$

三、共轭梯度法得matlab实现
在共轭梯度法得实际使用中，通常在迭代n步或n+1步之后，重新选取负梯度方向作为搜索方向，我们称之为再开始共轭梯度法。这是因为对于一般非二次函数而言，n步迭代后共轭梯度法产生得搜索方向往往不再具有共轭性。而对于大规模问题，常常每$m (m<n或m\ll n)$步就进行再开始。此外，当搜索方向不是下降方向时，也插入负梯度方向所作为搜索方向。
这里给出基于Armijo-rule非精确线搜索得再开始FR共轭梯度法得matlab程序。

function [fmin, xmin] = frcg(fun, gfun, x0, epsilon)

maxk = 5000;
rho = 0.6;
sigma = 0.4;
k = 0;
n = length(x0);
 
 while k < maxk
     g = feval（gfun, x0);
     itern = k - (n+1) * floor(k / (n + 1));
     itern = itern + 1;
     
     if (itern == 1)
         d = -g;
     else
         beta = (g' * g) / (g0' * g0);
         d = -g + beta * d0; gd = g' * d;
         if (gd >= 0.0)
             d = -g;
         end
     end
     
     if (norm(g) < epsilon), break; end
     
     m = 0; mk = 0;
     
     while (m < 20)  % armijo-rule 
         if (feval（fun, x0 + rho^m*d) < feval（fun, x0) + sigma * rho^m*g'*d)
             mk = m; break;
         end
         m = m + 1;
     end
     
     x0 = x0 + rho^mk*d;
     val = feval（fun, x0);
     fprintf('kIter = %d, fmin = %f\n', k, val);
     
     g0 = g; d0 = d;
     k = k+1;
 end
 x = x0;
 val = feval（fun, x);
 xmin = x;
 fmin = val;


end

利用程序求解无约束优化问题
$\min_{x\in\mathbb{R^2}}\quad f(x)=100(x_1^2-x_2)^2+(x_1-1)^2$
该问题由精确解$x^*=(1,1)^T,f(x^*)=0$

求解main函数

x0 = [0, 0]';
epsilon = 1e-4;

[fmin, xmin] = frcg('func', 'gfunc', x0, epsilon);
fprintf('frcg: fmin = %f, xmin = (%f, %f)\n', fmin, xmin(1), xmin(2));

[x, f] = fminsearch('func', x0);
fprintf('build-in search: fmin = %f, xmin = (%f, %f)\n', f, x(1), x(2));

函数定义以及梯度求解

function f = func(x)

f = 100 * (x(1)^2 - x(2))^2 + (x(1) - 1)^2;

end

function grad = gfunc(x)

grad = [400 * x(1) * (x(1)^2-x(2))+2*(x(1)-1); ...
    -200 * (x(1)^2-x(2))];

end

求解结果为：

frcg: fmin = 0.000000, xmin = (0.999921, 0.999841)
build-in search: fmin = 0.000000, xmin = (1.000004, 1.000011)