共轭梯度法程序
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
一、共轭梯度法
共轭梯度法（Conjugate Gradient）是共轭方向法的一种，因为在该方向法中每一个共轭向量都是依靠赖于迭代点处的负梯度而构造出来的，所以称为共轭梯度法。由于此法最先由Fletcher和Reeves （1964）提出了解非线性最优化问题的，因而又称为FR 共轭梯度法。由于共轭梯度法不需要矩阵存储，且有较快的收敛速度和二次终止性等优点，现在共轭梯度法已经广泛地应用于实际问题中。共轭梯度法是一个典型的共轭方向法，它的每一个搜索方向是互相共轭的，而这些搜索方向d仅仅是负梯度方向与上一次迭代的搜索方向的组合，因此，存储量少，计算方便，效果好。
二、共轭梯度法的原理
设有目标函数
f（X）=1/2X T HX+b T X+c 式1 式中，H作为f（X）的二阶导数矩阵，b为常数矢量，b=[b1,b2,b3,...b n]T 在第k次迭代计算中，从点X（k）出发，沿负梯度方向作一维搜索，得
S(K)=-∆f（X（k））式2 X(k+1)=X(k)+ɑ(k)S(k) 式3
在式中，ɑ（k）为最优步长。
设与S(k）共轭的下一个方向S(k+1)由点S(k）和点X(k+1)负梯度的线性组
合构，即
S (k+1)=-∆f （X （k+1））+β(k)S (k) 式4 根据共轭的条件有
[S (k)]T ∆2f （X (k)）S (k+1)=0 式5 把式2和式4带入式5，得
-[∆f(X (k))]T ∆2f （X (k)）[-∆f （X （k+1））+β(k)S (k) ]=0 式6 对于式1，则在点X (k)和点X (k+1)的梯度可写为
∆f(X (k))=HX (k)+b 式7 ∆f （X （k+1）)=HX (k+1)+b 式8 把上面两式相减并将式3代入得
ɑ(k)H S (k)=∆f （X （k+1）)-∆f(X (k)) 式9 将式4和式9两边分别相乘，并代入式5得
-[∆f （X （k+1）)+β(k)∆f(X (k))]T [∆f （X （k+1）)-∆f(X (k)]=0 式10 将式10展开，并注意到相邻两点梯度间的正交关系，整理后得 β
（k ）
=2
2
||))((||||))1((||k X f k X f ∆+∆ 式11
把式11代入式4和式3，得 S （k+1）=-∆f （X (k)）+β
（k ）
S （k ）
X (k+1)=X （k ）+ɑ（k ）S （k ）
由上可见，只要利用相邻两点的梯度就可以构造一个共轭方向。以这种方式产生共轭方向并进行迭代运算的方法，即共轭梯度法。
三、共轭梯度法的迭代方法步骤：
1、给定的起始点X （0）和迭代计算精度h ；
2、取X （0）的负梯度作为搜索方向
S (0)=-∆f （X (0)） 3、沿方向S(k)进行一维搜索，得 X (k+1)=X （k ）+ɑ（k ）S （k ）
4、进行收敛检验。
若满足||∆f （X (k+1))||≤h,则令
X *=X (k+1），f （X *）=f （X (k+1)） 输出最优解，结束迭代计算；否则，转入（5）
5、若k=n ，则令X (0)=X (k+1)，转入（2）开始新一轮的迭代；否则，转入（6）。
6、构造新的共轭方向 β=2
2
||))((||||))1((||k X f k X f ∆+∆
S （k+1）=-∆f （X (k)）+β（k ）
S （k ）
令k=k+1，转入（3）.
四、共轭梯度法的计算框图如下：
给定X (0), h
S
（0）
=-∆f （X (0)）,k=0
minf(X (0))+ɑS (k)） X （k-1）=X （k ）+ɑ（k ）S （k ）
||∆f （X
（k+1）
）||≤h
k=n ？
Β=2
2||))((||||))1((||k X f k X f ∆+∆
S （k+1）
=-∆f （X (k)）+β
（k ）
S
（k ）
K=k+1
X （0）=X (k+1) Y
N
N
X *
=X (k+1)
f （X *）=f （X (k+1)）
转出
Y
五、共轭梯度法的C++实现程序清单：
#include
#include
#define N 10
#define eps pow(10,-6)
double f(double x[],double p[])
{
double s;
s=x[]*x[]+4*p[]*p[]-1;
return s;
}
/*以下是进退法搜索区间源程序*/
void sb(double *a,double *b,double x[],double p[]) {
double t0,t1,t,h,alpha,f0,f1;
int k=0;
t0=1; /*初始值*/
h=0.01; /*初始步长*/
alpha=1;
f0=f(x,p,t0);
t1=t0+h;
f1=f(x,p,t1);
while(1)
{
if(f1 { h=alpha*h; t=t0; t0=t1; f0=f1; k++; } else { if(k==0) {h=-h;t=t1;} else { *a=t } } t1=t0+h;