共轭梯度法C语言(西安交大)
共轭梯度法C语言(西安交大)
#include<stdio.h>#include<math.h>#define N 10 /*定义矩阵阶数*/void main(){int i,j,m,A[N][N],B[N];double X[N],akv[N],dka[N],rk[N],dk[N],pk,pkk,ak,bk;for(i=0;i<N;i++) /*输入系数矩阵A*/ {for(j=0;j<N;j++){if(i==j)A[i][j]=-2;else if(abs(i-j)==1)A[i][j]=1;elseA[i][j]=0;}}for(i=0;i<N;i++) /*输入矩阵B*/ {if(i==0||i==N-1)B[i]=-1;elseB[i]=0;}printf("\n"); /*输出系数矩阵A*/printf("the Maxtr A is\n");for(i=0;i<N;i++){printf("\n");for(j=0;j<N;j++)printf("%3d",A[i][j]);}printf("\n"); /*输出矩阵B*/printf("the Maxtr b is\n");for(i=0;i<N;i++)printf("%3d",B[i]);printf("\n");printf("input the Maxtr X\n"); /*给X0输入一个初始向量*/ for(i=0;i<N;i++)X[i]=0;printf("X chushi:\n");for(i=0;i<N;i++) /*显示给定的初始向量X0*/ printf("%3.0f",X[i]);/*开始计算*/for(i=0;i<N;i++) /*计算rk*/{pk=0.0;for(j=0;j<N;j++)pk+=A[i][j]*X[j];akv[i]=pk;}for(i=0;i<N;i++)rk[i]=B[i]-akv[i];for(i=0;i<N;i++) /*给rO赋值*/dk[i]=rk[i];for(m=0;m<N;m++) /*开始迭代循环*/{for(i=0;i<N;i++) /*计算ak*/{pk=0.0;for(j=0;j<N;j++)pk+=A[i][j]*dk[j];dka[i]=pk;}pk=0.0;pkk=0.0;for(i=0;i<N;i++){pk+=dka[i]*dk[i];pkk+=rk[i]*rk[i];}if(pkk==0) /*如果残差pkk=0,计算结束*/ break;ak=pkk/pk;for(i=0;i<N;i++) /*计算X(k+1)*/X[i]=X[i]+ak*dk[i];for(i=0;i<N;i++) /*计算r(k+1)*/{pk=0.0;for(j=0;j<N;j++)pk+=A[i][j]*X[j];akv[i]=pk;}for(i=0;i<N;i++)rk[i]=B[i]-akv[i];pk=0.0;for(i=0;i<N;i++) /*计算bk*/pk+=rk[i]*rk[i];bk=pk/pkk;for(i=0;i<N;i++) /*计算d(k+1)*/dk[i]=rk[i]+bk*dk[i];}printf("\n");printf("X cacualtate value is:\n"); /*输出运算结果X*/for(i=0;i<N;i++)printf("%f\n",X[i]);}。
共轭梯度法的研究
共轭梯度法的研究
共轭梯度法是一种常用的优化算法，广泛应用于求解大规模线性方程组、最小二乘问题、非线性方程组等问题。
该算法利用了线性代数中共轭向量的性质，使得每次迭代都能够跨越一定的距离，从而快速收敛到最优解。
本文将介绍共轭梯度法的基本原理、迭代公式以及算法的实现细节。
同时，我们还将探讨共轭梯度法在不同问题中的应用，以及其优点和不足之处。
最后，我们将结合实例深入探讨共轭梯度法的实际应用效果，并提出未来的研究方向。
- 1 -。
数据科学优化方法 第7章 共轭梯度方法
7.1 共轭方向方法 7.2 针对正定二次函数的共轭梯度方法 7.3 非线性共轭梯度方法 7.4 数值实验
引入
• 关于无约束最优化问题，我们已经学习了一阶方法 (负梯 度方法) 和二阶方法 (牛顿方法).
1 在实际使用中，负梯度方法的收敛速度往往 很慢
2 牛顿方法的收敛速度虽然很快，但却需要计 算 Hesse 矩阵的逆.
因此，
dkTQ xk x0 0.
dkTQ x* x0 dkTQ x* xk dkT b Qxk dkTrk
对比(7.2)和(7.3)，有 k k，由此定理得证.
共轭方向方法的计算步骤及性质
• 当 Q为对角阵，每次沿坐标轴方向迭代可以依次确定极小 点 x*的一个分量.
1. A有 n 个实特征根 2. 属于不同特征值的特征向量正交 3. A可正交对角化
根据该定理，从任意初始点出发，依次沿 Q 的 n 个正交特征 向量方向做精确线搜索，便可得到正定二次函数问题的极小值.
方法的引入
min
f
(x)
1 2
xT
7/ 3 3
2 /2
33 13 /
/2 2
x
-
733 2来自,13 23
3
x
+
41 4
3
3
图7.2 沿负梯度方向的迭代路径 (实线)和沿正交特征向量方向的迭代路径 (虚线)
方法的引入
• 尽管根据谱定理可以解决正定二次函数的优化问题，但求矩阵 的特征值和特征向量需要付出较高的计算成本.
• 为了能处理一般的最优化问题，尤其是大规模最优化问题，我 们必须另外寻找合适的迭代方向.
对 n 维问题，从 x0 出发，依次沿 n 个坐标轴方向作精确线 搜索，便得到该问题最小值
共轭梯度法的C实现
西安交通大学实验报告课程名称：数值分析上机实验实验名称：共轭梯度法学院：___数学学院______________________ 班级姓名：学号：实验日期 2015 年 05 月 26 日自评成绩：97一、实验目的(1)熟练掌握改进平方根法和共轭梯度法的迭代过程(2)尝试使用自己熟悉的计算机语言解决数学中的问题(3)通过上机实验来巩固课本中所学的知识二、实验内容与结果题目2：共轭梯度法源程序2#include <iostream>using namespace std;double f1(double a[10],double n)//构造第一个求和函数，简化主函数{double s=0;int i;for(i=0;i<n;i++){s=s+a[i]*a[i];}return s;}double f2(double r[10],double a[10][10],int n)//构造第二个求和函数，简化主函数{double b[10],m=0;int i,j;for(i=0;i<n;i++){double s=0;for(j=0;j<n;j++){s=s+r[j]*a[j][i];}b[i]=s;}for(i=0;i<n;i++){m=m+b[i]*r[i];}return m;}double *f3(double a[10][10],double b[10],int n)//构造输出列向量的函数{double *r;r=new double[10];int i,j;for(i=0;i<n;i++){double s=0;for(j=0;j<n;j++){s=s+a[i][j]*b[j];}r[i]=s;}return r;}int main(){int i,j,k,n;double a,b,e,A[10][10],B[10],r[10],r1[10],x[10],d[10];double *ax;cout<<"请输入未知数的个数和计算精度："<<endl; cin>>n>>e;cout<<"请输入起始向量"<<endl;for(i=0;i<n;i++){cin>>x[i];}cout<<"请输入右端项： "<<endl;for(i=0;i<n;i++){cin>>B[i];}cout<<"请输入矩阵： "<<endl;for(i=0;i<n;i++){for(j=0;j<n;j++){cin>>A[i][j];}}ax=f3(A,x,n);for(i=0;i<n;i++) {r[i]=B[i]-ax[i];d[i]=r[i];}if(sqrt(f1(r,n))<=e) {for(i=0;i<n;i++) {cout<<x[i]<<'\t';}}else{for(k=0;k<n;k++){a=f1(r,n)/f2(d,A,n);for(i=0;i<n;i++){x[i]=x[i]+a*d[i];}ax=f3(A,x,n);for(i=0;i<n;i++){r1[i]=B[i]-ax[i];}if(sqrt(f1(r1,n))<=e||k+1==n) {for(i=0;i<n;i++){cout<<"x["<<i<<"]="<<x[i]<<'\n'; }break;}else{b=f1(r1,n)/f1(r,n);for(i=0;i<n;i++){d[i]=r1[i]+b*d[i];r[i]=r1[i];}}}}return 0;}运行结果2三、实验总结（列出实验的结论、收获和存在的问题，必写）通过这次实验，我对简化平方根法和共轭梯度法的计算过程和迭代格式有了熟练地掌握，让我学习到了很多，同时也是一个复习的过程。
共轭梯度法
共轭梯度法
数学上，共轭梯度法是求解特定线性系统的数值解的方法，其中那些矩阵为对称和正定。
共轭梯度法是一个迭代方法，所以它适用于稀疏矩阵系统，因为这些系统对于象乔莱斯基分解这样的直接方法太大了。
这种系统在数值求解偏微分方程时相当常见。
共轭梯度法也可以用于求解无约束的最优化问题。
双共轭梯度法提供了一种处理非对称矩阵情况的推广。
方法的表述
设我们要求解下列线性系统
Ax = b,，
其中n-×-n矩阵A是对称的（也即，A T = A），正定的（也即，x T Ax > 0对于所有非0向量x属于R n），并且是实系数的。
将系统的唯一解记作x*。
最后算法
经过一些简化，可以得到下列求解Ax = b的算法，其中A是实对称正定矩阵。
x
:= 0
k := 0
r
:= b
repeat until r k is "sufficiently small":
k := k + 1
if k = 1
p
:= r0
1
else
end if
x
:= x k-1 + αk p k k
r
:= r k-1 - αk A p k k
end repeat
结果为x k。
共轭梯度方法
共轭梯度方法（Conjugate Gradient Method）是求解线性方程组的一种迭代算法。
该方法适用于求解大型稀疏的对称正定线性方程组，可以显著减少计算量和存储空间。
该方法的主要思想是利用共轭方向（Conjugate Directions）的性质，在有限次迭代中求解方程组的解。
共轭梯度方法的基本步骤如下：
选取一个初值$x_0$，并令$r_0=b-Ax_0$，其中$b$ 为方程组的右端向量，$A$ 为系数矩阵。
计算一个共轭方向$p_0=r_0$，即$p_0$ 与$r_0$ 正交，并满足$Ap_0 \neq 0$。
对于$k=0,1,2,\ldots$，执行以下操作：
a. 计算$\alpha_k=\frac{r_k^Tr_k}{p_k^TAp_k}$。
b. 更新解向量$x_{k+1}=x_k+\alpha_kp_k$。
c. 计算残差向量$r_{k+1}=r_k-\alpha_kAp_k$。
d. 计算$\beta_k=\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k}$。
e. 更新共轭方向$p_{k+1}=r_{k+1}+\beta_kp_k$，即$p_{k+1}$ 与$p_k$ 具有共轭性。
如果残差向量$r_k$ 较小，则停止迭代，输出解向量$x_k$。
共轭梯度方法具有收敛速度快、存储空间小等优点，但对于非对称和非正定的线性方程组，该方法可能不收敛。
同时，该方法也有一些变体，如预处理共轭梯度法、共轭残差法等，可以更好地解决不同类型的线性方程组求解问题。
共轭梯度法 计算化学
共轭梯度法计算化学
共轭梯度法是一种常用于求解线性方程组的迭代方法，也可以用于计算化学中的一些问题。
在计算化学中，共轭梯度法常用于求解分子结构优化、能量最小化、电子结构计算等问题。
其中最常见的应用是求解非常大的线性方程组，如Hartree-Fock方程，即通过求解一个自洽的线性方程组来得到分子的基态电子结构。
共轭梯度法利用了线性方程组的特殊性质，通过迭代计算逼近线性方程组的解。
它不需要事先知道线性方程组的解的精确形式，只需要知道线性方程组的系数矩阵和右端向量。
在每一次迭代中，共轭梯度法利用之前的迭代结果和当前的残差信息来计算一个新的搜索方向，从而逼近解的位置。
通过多次迭代，共轭梯度法可以逐渐接近线性方程组的解。
在化学计算中，尤其是量子化学计算中，共轭梯度法经常用于求解大型稠密矩阵的特征值问题。
这些问题通常涉及求解特征方程，得到能量的本征值和本征函数。
共轭梯度法的迭代性质使得它非常适合处理大型稠密矩阵的特征值计算问题，因为它只需要存储和计算与当前解和残差相关的信息，而不需要存储和计算整个矩阵。
总之，共轭梯度法是一种非常有效的迭代方法，可以用于求解线性方程组以及一些涉及大型矩阵的计算化学问题。
在实际应用中，共轭梯度法经常与其他优化算法结合使用，以获得更高效、更准确的结果。
第二节共轭梯度法
Page 1§3.4 共轭梯度法基本思想Page 2利用目标函数在当前迭代点处的负梯度方向与上一步的搜索方向的适当线性组合，构造下一步的搜索方向，要求这些搜索方向是一系列共轭方向。
由Taylor公式知，一个函数在一点附近的性态与二次函数是很接近的，因此，为了建立有效算法,往往选用二次模型，即先针对正定二次函数建立有效的算法，然后再推广到一般函数上去。
Page 3一、共轭方向及其性质注：若,Q I =则是正交的，因此共轭是正交的推广．定义1:如果：()0,T ij d Q d i j =≠则称m d d d ,,,21"是关于Q 共轭的．设m d d d ,,,21"是nR 中任一组非零向量，n nQ R×∈是n 阶实正定矩阵，Page 4定理1:设Q 为n 阶正定阵，非零向量组m d d d ,,,21"关于Q 共轭，则必线性无关．推论1:设Q 为n 阶正定阵，非零向量组n d d d ,,,21"关于Q 共轭，则它们构成nR 的一组基．n 则推论2:设Q 为阶正定阵，非零向量组n d d d ,,,21"关于Q 共轭.若向量v 与n d d d ,,,21"均正交,.0=vPage 5定理2:设为阶正定阵，向量组Q n 12,,kd d d "关于Q 共轭，对正定二次函数()12T Tf x x Qx b x c =++，由任意1x 开始，1,1,2,,i ii i xx d i k α+=+="则：(１)1()0,1,2,k T if x d i k+∇=="(２)当k n =时，1n x +为f(x)在n R 上的极小点．后得沿k 次精确线搜索依次进行i dPage 6启发：用迭代法求解正定二次函数极小：()12T TMin f x x Qx b x c=++转化为构造关于Q共轭的n个方向问题.如何构造n个关于Q共轭的方向?一种做法是借助迭代点处的负梯度方向来构造共轭梯度法Page 7二、求解正定二次函数的共轭梯度法令11(),.k k kk k df xd λλ++=−∇+为待定组合系数1()()()kTk k k T kd Q f x d Q dλ+∇=左乘(),kTd Q 并使1()0,kTk d Q d+=得：1. 构造共轭方向—用待定系数法1,:,k k kk k x x d αα+=+精确线搜得步长进而得索取：0()d f x =−∇令0;k =Page 8不仅如此, 还可以证明:1()0,1,2,,j Tk d Qdj k+=="12,,,nd d d "可见, 若中途不停机的话, 按上述方法得到的关于Q 是两两共轭的.n 个方向Page 92. 正定二次函数的共轭梯度法算法Step2:如果(),kf xε∇<停止迭代；令1.k k kk xx d α+=+Step3:令1,k k =+转Step2.0()()argmin ()()k T k k kk k T kd f x f x d d Qd ααα>∇=+=−否则沿精确线搜索求步长:kd 111()()(),()k T k k k kk k k T kd Q f xd f x d d Qdλλ+++∇=−∇+=其中并计算:Step1:给出0,(),01,0.nx R d f x k ε∈=−∇≤<<=Page 1012,,,nd d d "可以证明, 若中途不停机的话, 这样得到的关于Q 是两两共轭的.n 个方向再结合定理2, 因此1n x+一定是所求的最优解.3. 收敛性1()()()kTk k k T kd Q f x d Q dλ+∇=（Hestenes-Stiefel 公式）组合系数的选取:可见,该公式中里面含有Hessen矩阵Q, 不能直接应用于求解一般非正定二次函数情形.Page 114. 组合系数的其他形式（１）FR 公式212()()k k kf xf x λ+∇=∇（1964）（2）DM 公式21()-()()k k k Tkf xd f x λ+∇=∇（Fletcher-Reeves 公式）（Dixon-Myers 公式）Page 12组合系数的其他形式（3）PRP 公式()12()()()()k Tk kk kf x f xf x f x λ+∇∇−∇=∇（1969）（Polak-Ribiere-Polyak 公式）Page 13三、求解一般函数的FR 共轭梯度法Step1:给出0,(),01,0.nx R d f x k ε∈=−∇≤<<=否则，由精确线搜索求步长:Step3:0argmin ().k kk f x d ααα>=+进而得1.k k kk x x d α+=+返回Step2.Step2:如果(),kf x ε∇<*;kx x =取停,()21112()()k k k kkf xdf x df x +++∇=−∇+∇以及1,k k =+令Page 14FR 共轭梯度法收敛定理定理4:假定()x f 在有界水平集()(){}nL x R f x f x=∈≤上连续可微，且有下界，那么采用精确线搜索下的FR 共轭梯度法产生的点列{}kx 至少有一个聚点是驻点，即：(1)当{}kx 是有穷点列时，其最后一个点是()x f 的驻点．(2)当{}kx 是无穷点列时，它必有聚点，且任一聚点都是()x f 的驻点．Page 15例2：用共轭梯度法求解：()()22012min 4,1,1Tf x x x x =+=取又由()f x 可看出()()1122201,208x f x x x x ⎛⎞⎛⎞=⎜⎟⎜⎟⎝⎠⎝⎠2008Q ⎛⎞=⎜⎟⎝⎠是一个正定二次函数, 其中解:因共轭梯度法的第一步迭代与最速下降法相同,故由上一节例1知分析：由例1知该函数的极小点是(0,0).Tx ∗=因此在求最优步长和组合系数时,有两种方法.Page 16()01,1:Tx =得()()2,8Td f x=−∇=−−(1) 由令00()0:0.13077df x d d ααα+==得由于较大, 因此还需迭代下去.()1f x ∇10120.738460.13077180.04616x x d α⎛⎞⎛⎞⎛⎞=+=−=⎜⎟⎜⎟⎜⎟−⎝⎠⎝⎠⎝⎠进而有:()()()111.47692,0.36923,1.52237Tf xf x ∇=−∇=Page 17()110d f x dλ=−∇+(2) 先构造下一个搜索方向110()1.476922 1.545080.034080.3692380.09659d f x dλ=−∇+−−−⎛⎞⎛⎞⎛⎞=+=⎜⎟⎜⎟⎜⎟−⎝⎠⎝⎠⎝⎠于是有21020()0.03408()f x f x λ∇===∇"01000()()17.723040.03408()520TT d Q f x d Qd λ∇===="或其中, 组合系数Page 18沿搜索方向用精确线搜索求步长:21110.73846 1.5450800.477940.046160.096590x x d α−⎛⎞⎛⎞⎛⎞=+=+=⎜⎟⎜⎟⎜⎟−⎝⎠⎝⎠⎝⎠于是有令111()0:0.47794df x d d ααα+==得为此, 先计算出表达式11()f x d α+所以迭代终止,因()20,f x ∇=最优点为:()20,0.Tx x ∗==可见: 对于例1, 用共轭梯度法经两次迭代求得最优解.优点:Page 19（1）克服了最速下降法的慢收敛性.（2）建立在二次模型上，对正定二次模型具有二次终止性．（3）算法简单，易于编程，需存储空间小等优点,是求解大规模问题的重要方法．Page 20。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
#include<stdio.h>
#include<math.h>
#define N 10 /*定义矩阵阶数*/
void main()
{
int i,j,m,A[N][N],B[N];
double X[N],akv[N],dka[N],rk[N],dk[N],pk,pkk,ak,bk;
for(i=0;i<N;i++) /*输入系数矩阵A*/ {
for(j=0;j<N;j++)
{
if(i==j)
A[i][j]=-2;
else if(abs(i-j)==1)
A[i][j]=1;
else
A[i][j]=0;
}
}
for(i=0;i<N;i++) /*输入矩阵B*/ {
if(i==0||i==N-1)
B[i]=-1;
else
B[i]=0;
}
printf("\n"); /*输出系数矩阵A*/
printf("the Maxtr A is\n");
for(i=0;i<N;i++)
{
printf("\n");
for(j=0;j<N;j++)
printf("%3d",A[i][j]);
}
printf("\n"); /*输出矩阵B*/
printf("the Maxtr b is\n");
for(i=0;i<N;i++)
printf("%3d",B[i]);
printf("\n");
printf("input the Maxtr X\n"); /*给X0输入一个初始向量*/ for(i=0;i<N;i++)
X[i]=0;
printf("X chushi:\n");
for(i=0;i<N;i++) /*显示给定的初始向量X0*/ printf("%3.0f",X[i]);
/*开始计算*/
for(i=0;i<N;i++) /*计算rk*/
{
pk=0.0;
for(j=0;j<N;j++)
pk+=A[i][j]*X[j];
akv[i]=pk;
}
for(i=0;i<N;i++)
rk[i]=B[i]-akv[i];
for(i=0;i<N;i++) /*给rO赋值*/
dk[i]=rk[i];
for(m=0;m<N;m++) /*开始迭代循环*/
{
for(i=0;i<N;i++) /*计算ak*/
{
pk=0.0;
for(j=0;j<N;j++)
pk+=A[i][j]*dk[j];
dka[i]=pk;
}
pk=0.0;pkk=0.0;
for(i=0;i<N;i++)
{
pk+=dka[i]*dk[i];
pkk+=rk[i]*rk[i];
}
if(pkk==0) /*如果残差pkk=0,计算结束*/ break;
ak=pkk/pk;
for(i=0;i<N;i++) /*计算X(k+1)*/
X[i]=X[i]+ak*dk[i];
for(i=0;i<N;i++) /*计算r(k+1)*/
{
pk=0.0;
for(j=0;j<N;j++)
pk+=A[i][j]*X[j];
akv[i]=pk;
}
for(i=0;i<N;i++)
rk[i]=B[i]-akv[i];
pk=0.0;
for(i=0;i<N;i++) /*计算bk*/
pk+=rk[i]*rk[i];
bk=pk/pkk;
for(i=0;i<N;i++) /*计算d(k+1)*/
dk[i]=rk[i]+bk*dk[i];
}
printf("\n");
printf("X cacualtate value is:\n"); /*输出运算结果X*/
for(i=0;i<N;i++)
printf("%f\n",X[i]);
}。