精华内容
下载资源
问答
  • 共轭方向共轭梯度

    千次阅读 多人点赞 2017-11-16 15:18:46
    共轭方向是介于最速下降和Newton之间的一种方法。克服了最速下降的锯齿...同时,共轭方向的迭代公式比较简单,不必目标函数的Hesse矩阵,比Newton减少了计算量和存储量。是一种比较实用而且有效的方法。

    共轭方向法是介于最速下降法和Newton法之间的一种方法。克服了最速下降法的锯齿现象,从而提高了收敛速度;同时,共轭方向法的迭代公式比较简单,不必求目标函数的Hesse矩阵,比Newton法减少了计算量和存储量。是一种比较实用而且有效的方法。
    在讲共轭方向法和共轭梯度法之前,先对共轭向量进行说明。

    共轭向量及其性质

    定义1:(共轭方向) Qn×n 对称正定矩阵,若 n 维向量空间中的非零向量 p0p1pm1 满足

    pTiQpj=0     ij=01m1(ij)(1)
    则称p0p1pm1Q共轭向量或称向量p0p1pm1Q共轭的(简称共轭),称p0p1pm1的方向是Q共轭方向
    Q=I()时,公式(1)变为
    pTipj=0     ij=01m1(ij)(2)
    即向量p0p1pm1互相正交。由此可见,正交是共轭的一种特殊情况,共轭是正交的推广。
    定理1: 若非零向量p0p1pm1Q共轭的,则他们是线性无关的。
    推论1: 在n维的向量空间中,非零的共轭向量个数不超过n。
    定义2:p0p1pm1Rn中线性无关向量,x0Rn,那么由形式为
    z=x0+i=0m1αipi    α1α2αm1
    的向量构成的集合为由点x0和向量p0p1pm1所生成的线性流行。记为L[x0;p0,p1,,pm1]

    基本思想

    在考虑普通函数之前,我们首先用2元正定二次函数进行讲解。首先考虑如下的正定二次函数

    f(x)=12xTQx+bTx+c(3)
    要求的目标函数f(x)的最优值,根据最速下降法的思想,我们首先选定一个初始点x0,然后沿着该点的最速下降方向p0=f(x0)做直线搜索,得到点x1,由最速下降法的性质可知
    f(x1)Tp0=0(4)
    p0与点x1出的等值线相切。
    在第二次迭代过程中我们不适用f(x1)作为这次迭代的搜索方向,我们想直在第二次迭代之后能直接到达最优点x,那么这次的迭代方向p1应该满足什么条件呢?
    首先根据迭代公式我们有
    x=x1+t1p1(5)
    其中t_1是最优步长因子,显然在未到达最优点x之前,t1是不等于0的。对目标函数求梯度,有
    f(x)=Qx+b(6)
    对于极小点x,我们有
    f(x)=Qx+b=0=Q(x1+t1p1)+b=Qx1+b+Qt1p1=f(x1)+Qt1p1=0
    f(x1)+Qt1p1=0
    在上式两边同乘以pT0,由于公式(4),并且t10我们可以得到
    pT0f(x1)+t1pT0Qp1=>t1p0Qp1==pT000(7)
    由公式(7)我们知道,p0p1Q的共轭向量。
    现在我们假设
    p1=f(x1)+α0p0(8)
    在上式两侧同时乘以pT0Q,得到
    pT0Qf(x1)+α0pT0Qp0=pT0Qp1=0
    解得
    α0=pT0Qf(x1)pT0Qp0
    带入到公式(8)得到
    p1=f(x1)+pT0Qf(x1)pT0Qp0p0(9)

    综上所述,对于n元正定二次函数,我们可以从任意点出发,然后沿着这n个共轭方向最多做n次直线搜索(原因将在下一小节进行讲解),就可以求的目标函数的最优点。

    共轭方向法

    定理2: 假设:
    (1)Qn×n对称正定矩阵;
    (2)非零向量p0,p1,,pm1Q共轭向量;
    (3)对二次目标函数(3)顺次进行m次直线搜索

    xi+1=1s(xipi)   i=01m1
    其中x0Rn是任意选定的初始点。则
    i)pTif(xm)=0   0j<m
    ii)xm是二次函数(3)在线性流行L[x0;p0,p1,,pm1]上的极小点。

    (个人理解,如有错误,请指正)前面已经说过,共轭是正交的推广,对于n维空间,我们可以把n(不超过n)个共轭向量作为n为空间的基,只不过这个基不再是正交的,那么共轭向量法就是对于定义在n为空间中的函数,沿着每一个基的方向做直线搜索,那么最多做n次搜索就能得到最优值。对于线性流行L[x0;p0,p1,,pn1],其实可以认为是经过点x0由基$$p_0,p_1,···,p_{m-1}]$构成的超平面。对于定理2,他所指出的就是最后一个迭代点$x_m$处的梯度与之前的搜索方向垂直。根据定理2我们可以得到如下的推论。

    推论2: 在定理2中,当m=n时,xn是正定二次函数在Rn中的极小点。
    推论3: 在定理2中,p0,p1,,pm1的任意线性组合i=0m1βipi都与f(xm)正交。

    对于正定二次函数(3),从任意点出发,顺次沿nQ共轭方向做直线搜索,最多经过n次就可以到达极小点。
    下面是对于正定二次函数的共轭方向法的算法描述。

    已知:正定二次目标函数f(x)=12xTQx+bTx+c和终止限ϵ
    (1)选定初始点x0和下降方向p0;置k=0
    (2)做直线搜索xk+1=1s(xk,pk)
    (3)判别||f(xk+1)||<ϵ是否满足:若满足,输出xk+1;否则转(4)
    (4)提供共轭方向pk+1,使得

    pTjQpk+1=0   j=01k

    (5)置k=k+1,转(2)

    对于正定二次函数,第(2)步可以采取如下显示计算公式

    xk+1=xkpTkf(xk)pTkQpkpk
    公式推导如下。
    xk+1QXk+1+bf(xk+1)pTkf(xk+1)0tk=xk+tkpk=Qxk+tkQpk+b=f(xk)+tkQpk=pTkf(xk)+tkpTkQpk=pTkf(xk)+tkpTkQpk=pTkf(xk)pTkQpk(10)

    在第4步中,只是说要求得共轭方向,并没有说明如何求共轭方向,不同求共轭方向的方法对应了不同的共轭方向法。

    共轭梯度法

    共轭梯度法是一种共轭方向法,在求每一个迭代点的搜索方向时,与改点的梯度有关,故叫做共轭梯度法。共轭梯度法:在初始点x0的搜索方向p0为初始点x0的负梯度方向f(x0);之后迭代点xk的搜索方向pk为该点的负梯度方向f(xk)与已经得到的搜索方向pk1的线性组合(即pk=gk+αk1pk1)。
    首先对正定二次函数做说明。

    用于正定二次函数的共轭梯度法

    对于目标函数为(3)的正定二次函数来说。我们规定gk=f(xk)
    * 第一个迭代点
    我们可以任意指定一个初始点x0,那么初始点处的搜索方向p0=g0(11),从x0出发沿着p0方向做直线搜索

    x1=x0+t0p0
    ,可以求得(可以参照公式(10)的推导)时的推导)
    t0=pT0g0pT0Qp0=gT0g0pT0Qp0
    x1=x0+gT0g0pT0Qp0p0
    由此我们便得到了第一个迭代点。
    * 第二个迭代点
    由直线搜索的特性我们可以知道p0g1=0(可以想一下最速下降法中的锯齿现象),故p0g1是线性无关的。根据共轭梯度法的特性我们可以得到在点x1处的搜索方向p1
    p1=g1+α0p0(12)
    。再由共轭方向法的特性我们知道p1p1Q共轭方向。故
    pT1Qp0=gT1Qp0+α0pT0Qp0=0
    便可以得到
    α0=gT1Qp0pT0Qp0
    由此我们便确定了迭代点x1处的搜索方向p1,之后便是从x1开始沿搜索方向p1做直线搜索x2=x1+t1p1,根据公式(10)和(12),可以得到
    t1=pT1g1pT1Qp1=gT1g1+α0p0g1pT1Qp1=gT1g1pT1Qp1
    由此我们便得到了第二个迭代点x2
    * 第三个迭代点
    同理,我们可以得到p1g2=0,可以得到点x2处的搜索方向为
    p2=g2+α1p1(13)
    同理,p2p1Q共轭方向,故能够得到
    pt2Qp1=gT2Qp1+α1pT1Qp1=0
    便可得到
    α0=gT2Qp1pT1Qp1
    至此,我们便得到了点x2处的搜索方向p2,然后从x2开始沿p2方向做直线搜索,与上相同便可得到x3=x2+t2p2,求得
    t2=gT2g2pT2Qp2
    至此我们便得到了第三个迭代点x3
    但是上面在求第三个迭代点时有一个小小的问题便是,以上求得p3的方法能保证p2p0共轭吗?即p0p1p2Q共轭向量(根据共轭向量定义我们知道共轭向量之间需要亮亮共轭)吗?答案是肯定的,接下来推导pT2Qp0=0,如下。
    pT2Qp0=gT2Qp0+α1pT1Qp0=gT2Qp0=gT2(g1g0)t0=gT2g1gT2g0t0
    由于g0g1都可以表示成p0p1的线性组合,故有gT2g1=gT2g0=0(因为可以把g0g1看作是由p0p1确定的超平面上线,由于gT2p1=0,故有gT2g1=gT2g0=0)。
    关于为什么Qp0=(g1g0)t0,可以参考公式(10)的推导过程,有f(xk+1)=f(xk)+tkQpk,由此得到。
    * 第k+1个迭代点
    与前面相似,我们可以得到
    pk=gk+αk1pk1(14)
    可以得到
    αk1=gkQpk1ptk1Qpk1
    ,我们便能得到xk处的搜索方向pk,在xk处沿pk方向做直线搜索,同理我们能得到
    tk=gTkgkptkQpk
    同理我们能够证出p0p1pkQ共轭向量。

    所以可以按照上述方法,依次构造出共轭向量,最多经过n次迭代就能找到最优点x
    共轭梯度法的算法描述如下。

    已知:二次函数公式(3),终止限ϵ
    (1)选定初始点x0;计算p0=f(x0);置k=0
    (2)做直线搜索xk+1=1s(xk,pk),或者采用如下公式计算

    tk=f(xk)Tf(xk)pTkQpk(15)
    xk+1=xk+tkpk(16)

    (3)判断||f(xk+1||<ϵ)是否满足要求。满足则输出xk+1停止;否则转(4)
    (4)计算
    αk=f(xk+1)TQpkpTkQpk(17)
    pk+1=f(xk+1)+αkpk(18)

    (5)置k=k+1,转(2)

    用于非二次函数的共轭梯度法

    上面所讲的迭代公式要想适用于非二次函数,就要将迭代公式(17)中的Q去掉,那么根据的推到我们可以得到

    Qpk=f(xk+1)f(xk)tk
    ,由此,我们能得到迭代公式
    αk=gTk+1(gk+1gk)pTk(gk+1gk)(19)
    在公式(14)两边同时乘以gk+1,得
    pTkgk+1=gTkgk+1+αk1pTk1gk1(20)
    由直线搜索的性质,我们知道
    pTkgk+1=0(21)
    于是,公式(19)变为
    gTk+1gk=0(21)
    另外
    pTkgk=gTkgk+αk1pk1gk=gTkgk(22)

    * 将公式(20)(21)(22)带入到公式(19)之后,我们得到
    αk=gTk+1gk+1gTkgk=||gk+1||2||gk||2(23)
    这个公式称为Fletcher-Reeves公式
    * 将公式(21)(22)带入到公式(19)后得到
    αk=gTk+1gk+1pTkpk(24)
    这个公式称为Dixon-Myers公式
    * 将公式(21)(23)带入到公式(19)后得到
    αk=gTk+1(gk+1gk)gTkgk(25)
    这个公式称为Polak-Ribiere公式
    将公式(23)(24)(25)替换公式(17)对应了不同的共轭梯度法。共轭梯度法不要求精确的直线搜索,但是不精确的直线搜索可能会造成之后迭代出来的向量不再是共轭的,这将会降低共轭梯度法的效能,解决方法就是重设初始点,即经过n+1次迭代后得到的xn+1作为初始点,开始新一轮的迭代。
    用于一般函数的Fletcher-Reeves共轭梯度法描述如下。

    已知:目标函数f(x)以及梯度函数g(x),问题的维数N以及H终止准则的终止限ϵ1ϵ2ϵ3以及终止限ϵ
    (1)选定初始点x0;计算p0=f(x0);置k=0(2)线x_{k+1}=1s(x_k,p_k)f_{k+1}=f(x_{k+1}),g_{k+1}=g(x_{k+1})$

    (2)进行直线搜索xk+1=1s(xkpk),计算fk+1=f(xk+1)gk+1=g(xk+1)
    (3)判断H终止准则是否满足。满足:输出xk+1,停止;否则:转(4)
    (4)判断k=n是否成立,即是否已经迭代了n+1次。是:重置x0=xk+1f0=fk+1g0=gk+1p0=g0k=0,然后转(2);不是:转(5)
    (5)按照Fletcher-Reeve公式计算

    αk=||gk+1||2||gk||2
    pk+1=gk+αkpk

    (6)判断pTk+1gk+10是否成立(检查pk+1是不是下降方向,由于实际计算中无法精确求解,可能会出现不是下降方向的情况;在精确求解中,pk+1一定是下降法向)。做三种情况处理:i)若pTk+1gk+1>ϵ,这时则改取pk+1作为搜索方向,并置k=k+1,转(2);ii)若|pTk+1gk+1|ϵ|,则重置x0=xk+1f0=fk+1g0=gk+1p0=g0k=0,转(2);iii)若pTk+1gk+1<ϵ,则pk+1就作为搜索方向,并置k=k+1,转(2)

    展开全文
  • 典型的最优化问题,用最速下降、牛顿共轭梯度法求最小值。
  • matlab共轭梯度法求目标函数的最小极值-共轭梯度-王.rar 我是地球物理专业的一名学生,把自己实习的作业发上来大家分享下
  • 预处理共轭梯度法求线性方程组Ax=b的解,数值计算,求解方程
  • 共轭梯度法求函数极小值,其中用进退法求步长区间,用黄金分割法求最佳步长。
  • 共轭梯度法求最优解

    千次阅读 2011-11-19 15:07:05
     已知 f(x) = (x1-1)2+5(x2-5)2+(x3-1)2+5(x4-5)2 ,用快速下降、牛顿共轭梯度法求 minf(x) 。   共轭梯度代码:     //共轭梯度 //请根据具体题目,修改本程序“//@”所在行的下一行代码。...

     

    题目:

          已知 f(x) = (x1-1)2+5(x2-5)2+(x3-1)2+5(x4-5)  用快速下降法、牛顿法或共轭梯度法求 minf(x) 

     

    共轭梯度法代码:

     

     

    //共轭梯度法
    //请根据具体题目,修改本程序“//@”所在行的下一行代码。
    #include<math.h>
    #include<stdio.h>
    #include<stdlib.h>
    //@ 题目中方程是几元的,此处将LEN设为几
    #define LEN 4
    #define TYPE float
    
    TYPE fAnswer(TYPE *x) {    //求f(x)的值,并返回
    	TYPE f;
    	//@ 题目中的方程写与此
    	f = (x[0] - 1) * (x[0] - 1) + 5 * (x[1] - 5) * (x[1] - 5) + (x[2] - 1) * (x[2] - 1) + 5 * (x[3] - 5) * (x[3] - 5);
    	return f;
    }
    
    void vectorMultiply(TYPE *x, TYPE e) {    //常数e乘x向量,赋值给x向量 【若求负梯度,可用梯度乘-1】
    	int i;
    	for(i = 0; i < LEN; i++) {
    		x[i] = e * x[i];
    	}
    }
    
    void vectorAdd(TYPE *x, TYPE *y, TYPE *z) {    //向量加法操作
    	int i;
    	for(i = 0; i < LEN; i++) {
    		z[i] = x[i] + y[i];
    	}
    }
    
    void vectorSub(TYPE *x, TYPE *y, TYPE *z) {    //向量减法操作
    	int i;
    	for(i = 0; i < LEN; i++) {
    		z[i] = x[i] - y[i];
    	}
    }
    
    void vectoreEqual(TYPE *x, TYPE *y) {    //向量赋值操作
    	int i;
    	for(i = 0; i < LEN; i++) {
    		y[i] = x[i];
    	}
    }
    
    TYPE norm2_graded(TYPE *x) {    //负梯度模的平方
    	int i;
    	TYPE norm2 = 0.0;
    
    	for(i = 0; i < LEN; i++) {
    		norm2 += x[i] * x[i];
    	}
    
    	return norm2;
    }
    
    //@ 对题目的xi分别求偏倒,赋值给stepLength[i]
    void setStepLength(TYPE *stepLength, TYPE *x0) {
    	stepLength[0] = 2 * (1 - x0[0]);
        stepLength[1] = 10 * (5 - x0[1]);
    	stepLength[2] = 2 * (1 - x0[2]);
        stepLength[3] = 10 * (5 - x0[3]);
    }
    
    void DSC(TYPE *x0, TYPE *stepLength) {    //用D.S.C.法求下一个x落点
    	TYPE x1[LEN];
    	TYPE x2[LEN];
    	TYPE x3[LEN];
    	TYPE x4[LEN];
    	TYPE x5[LEN];
    	TYPE xa[LEN];
    	TYPE xb[LEN];
    	TYPE xc[LEN];
    	vectorAdd(x0, stepLength, x1);
    	if(fAnswer(x1) > fAnswer(x0))
    		vectorMultiply(stepLength, -1);
    
    	vectorAdd(x0, stepLength, x1);
    	while(fAnswer(x1) <= fAnswer(x0)) {
    		vectoreEqual(x1, x0);
    		vectorAdd(stepLength, stepLength, stepLength);
    		vectorAdd(x1, stepLength, x1);
    	}
    
    	vectorMultiply(stepLength, 0.5);
    	vectorSub(x0, stepLength, x2);
    	vectoreEqual(x1, x4);
    	vectoreEqual(x0, x3);
    	vectorSub(x1, stepLength, x5);
    
    	if(fAnswer(x5) > fAnswer(x3))
    		vectoreEqual(x3, xb);
    	else
    		vectoreEqual(x5, xb);
    
    	vectorSub(xb, stepLength, xa);
    	vectorAdd(xb, stepLength, xc);
    
    	vectorMultiply(stepLength, (fAnswer(xa) - fAnswer(xc)) / (2 * (fAnswer(xa) - 2 * fAnswer(xb) + fAnswer(xc))));
    	vectorAdd(xb, stepLength, x0);
    	setStepLength(stepLength, xb);
    }
    
    void main() {    //方法主体
    	TYPE x0[LEN];    //初始点
    	TYPE stepLength[LEN];    //步长
    	TYPE stepLength2[LEN];
    	TYPE e = 0.001;    //误差
    	TYPE a = 0.00;
    
    	int i;    //用于循环计数
    	int n = LEN;
    	int count = 0;
    
        printf("请输入x的初始值:\n");
    	for(i = 0; i < LEN; i++) {    //初始化x0数组
    		printf("x%d = ", i+1);
    		scanf("%f", &x0[i]);
    	}
    
    	setStepLength(stepLength, x0);
    
       	for(i = 1; norm2_graded(stepLength) > e * e; i++) {    //共轭梯度法的实现主体
    		DSC(x0, stepLength);
    		if(i != n) {
    			if(norm2_graded(stepLength) <= e * e) break;
    			setStepLength(stepLength2, x0);
    			a = norm2_graded(stepLength2) / norm2_graded(stepLength);
    			vectorMultiply(stepLength, a);
    			vectorAdd(stepLength2, stepLength, stepLength); 
    			i = i + 1;
    		} else {
    			i = 1;
    		}
    		count = count + 1;
    	}
    
    	printf("共轭梯度法循环结束,共循环%d次\n", count);
    
    	printf("使用共轭梯度法获得的最优点为:\n");
    	for(i = 0; i < LEN; i++) {
    		printf("x%d = %f\n", i+1, x0[i]);
    	}
    	printf("minf(x) = %f\n", fAnswer(x0));
    
    	printf("共轭梯度法程序结束!!!\n");
    }
    

     

    总结

            本解法所使用的一维搜索算法仍是改进的D.S.C.法,将步长 的初始值改为负梯度向量,D.S.C.算法中的其它一维参量也采用向量形式替换,使此D.S.C.法可适用于多维搜索。

            本解法仅通过将我发表的最速下降法中的main()函数中最速下降法的实现主体替换为共轭梯度法的实现主体而得到。

            经测试验证,此方法无误。

    展开全文
  • 共轭梯度

    2019-11-28 09:35:21
    方向是在出梯度方向的前提下,添加正则项,使得前后两次方向互为共轭所得出的方向向量。 步长是使得该次下降达到最深位置的步长 代码如下,贴出最速下降和共轭梯度的 matlab 代码,可以比较...

    共轭梯度法

    我们知道最速下降法和共轭梯度法的目的都是一致的,那就是通过变分法求解线性方程组。

    大前提是该线性方程组的系数矩阵必须是对称正定矩阵

    我们做用迭代法解线性方程组的时候只关注两件事情
    1.方向
    2.步长

    方向是在求出梯度方向的前提下,添加正则项,使得前后两次方向互为共轭所得出的方向向量。

    步长是使得该次下降达到最深位置的步长

    代码如下

    function [ X ] = conjugateGradient( A, b, X0, e )
    %conjugateGradient: conjugate gradient method
    %Input   -     A: coefficient matrix
    %        -     b: constant term
    %        -    X0: X0 is the initial value of item
    %        -     e: e is the termination condition of iteration
    %Output  -     X: the solution of equations
    X = X0;
    % r 为剩余向量
    r = b-A*X;
    % d 为下降方向
    %0 步的下降方向为负的梯度方向
    d = r;
    %更新第0步的步长
    alpha = (r.*d)/(d'*A*d);
    %完成第一次迭代
    X = X + alpha.*d;
    %以后的搜索方向都是共轭的方向
    while(norm(r)>e)
        %更新新的梯度方向
        r = b-A*X;
        %更新方向
        %更新上一次的方向在这一次方向的表达式上的系数
        beta = -(r'*A*d)/(d'*A*d);%这里的d是上一次的下降方向
        d = r + beta.*d;
        %方向上的更新完成
        %更新步长
        alpha = (r.*d)/(d'*A*d);
        X = X + alpha.*d;
    end
    
    end

    测试

    输入参数:
    X0 = [0;0;0]
    b = [1;1;1]
    e = eps
    挑一个对角占优的矩阵
    在这里插入图片描述
    结果如下:
    在这里插入图片描述

    验证一下:
    在这里插入图片描述
    测试成功!

    展开全文
  • python实现共轭梯度

    2021-01-20 06:05:02
    共轭梯度是介于最速下降与牛顿之间的一个方法,它仅需利用一阶导数信息,但克服了最速下降收敛慢的缺点,又避免了牛顿需要存储和计算Hesse矩阵并逆的缺点,共轭梯度不仅是解决大型线性方程组最有用的...
  • 共轭梯度(1)--的基本原理之前已经搞明白了,梯度下降的基本原理,当然解释的调度是从函数极值的角度出发的,事实上从这个角度来理解,个人感觉是一个最为直接的理解角度,其完完全全是建立在多变量函数的微分...

    共轭梯度法(1)--的基本原理

    之前已经搞明白了,梯度下降法的基本原理,当然解释的调度是从求函数极值的角度出发的,事实上从这个角度来理解,个人感觉是一个最为直接的理解角度,其完完全全是建立在多变量函数的微分系统中的。事实上这个方法(思想)在实际中是被用于求解线性方程的,当然单纯的梯度下降形式并没有被直接采用,被广泛用于求解对称、正系数矩阵方程组的方法是,基于梯度下降原理实现的共轭梯度法,这一方法在大型稀疏线性方程组的求解中,有极高的效率。但是,之前对于梯度下降思想的描述,好像认为这个方法是用于求函数极值的方法,似乎与线性方程组求解问题无关。因此,对于对称、正定系数矩阵的二次型的理解是很必要的,它能让我们明白矩阵系统中线性方程组与函数极值的关系,能够解释为什么一个求函数极值的思想可以被用于求解线性方程组系统。

    $Ax = b$对应的二次型

    这里就不严格去给出严格的数学表示了,因为对于非数学专业的我们来说,好像没啥必要,知晓核心的思想更为重要。注意本次针对的系数矩阵$A$都是对称、正定的,线性方程组系统$Ax = b$对应的是二次型对于自变量导数为0时的表达式,不严谨地可将其看作是函数吧,即二次型为

    $$

    f(x) = \frac{1}{2} x^{T} A x - b^{T} x + c \tag{1}

    $$

    对于$Ax = b$来说,$c$是一个0向量。我们知道求$f(x)$极小值的方法最直接的一种就是利用其导数为0,获得极值点从而获得极小值,因此,求导得到

    $$

    f'(x) = \frac{1}{2}A^T x + \frac{1}{2}A x - b = 0 \tag{2}

    $$

    因此,求解线性方程组$Ax = b$就等价于求解$f'(x) = 0$这个表达式,那么一阶导数为0这个表达换一种说法,不就是去求二次型$f(x)$的函数的极值嘛!这也就解释了为什么一个求解函数极值的方法却可以应用到线性方程组的求解当中来。

    为了更好的理解线性方程组$Ax = b$的求解与函数极值之间的关系,这里举一个简单的,可以被可视化的例子,若方程组为:

    $$

    \begin{bmatrix}

    2 & 2\\

    2 & 5

    \end{bmatrix}

    \begin{bmatrix}

    x\\

    y

    \end{bmatrix}

    =

    \begin{bmatrix}

    6\\

    3

    \end{bmatrix} \tag{3}

    $$

    这样一个对称正定的方程组,为了待求的未知数为了更符合函数表达,这里用了$x, y$来表示,依据二次型公式,得到

    $$

    f(x, y) = \frac{1}{2}

    \begin{bmatrix}

    x\\

    y

    \end{bmatrix}

    \begin{bmatrix}

    2 & 2\\

    2 & 5

    \end{bmatrix}

    \begin{bmatrix}

    x & y

    \end{bmatrix} -

    \begin{bmatrix}

    6 & 3

    \end{bmatrix}

    \begin{bmatrix}

    x\\

    y

    \end{bmatrix}

    $$

    这样一个$2 \times 2$的线性方程组对应的二次型就是一个二元函数,为了得到表达式,这里利用sympy库进行符号运算。import numpy as np

    from sympy import *

    # 先定义符号变量x, y

    x, y = symbols('x y')

    A = np.array([[2, 2], [2, 5]])

    b = np.array([6, 3])

    # 自变量向量这里用x1表示

    x1 = np.array([x, y])

    f_xy = 1/2 * x1.T @ A @ x1 - b.T @ x1

    init_printing(use_unicode = True) # 更美观的显示

    pprint(f_xy.expand())2 2

    1.0⋅x + 2.0⋅x⋅y - 6⋅x + 2.5⋅y - 3⋅y

    利用函数图像来直观地可视化函数极值情况。import matplotlib.pyplot as plt

    import numpy as np

    from mpl_toolkits.mplot3d import Axes3D

    x = np.linspace(-10, 10, 100)

    y = np.linspace(-10, 10, 100)

    xx, yy = np.meshgrid(x, y)

    f_xy = xx * (xx + yy) - 6 * xx + yy * (xx + 2.5 * yy) - 3 * yy

    # 这里用等值线云图与伪色彩图来可视化,也顺便对比下两者的区别

    fig, ax = plt.subplots(1, 2)

    cf1 = ax[0].pcolormesh(xx, yy, f_xy, cmap='RdBu_r')

    cs1 = ax[0].contour(x, y, f_xy, colors='gray', levels=15)

    ax[0].clabel(cs1, inline=True, colors='k')

    fig.colorbar(cf1, orientation='vertical')

    ax[0].set_title('function by pcolormesh')

    cf2 = ax[1].contourf(x, y, f_xy, cmap='RdBu_r')

    cs2 = ax[1].contour(x, y, f_xy, colors='gray', levels=15)

    ax[1].clabel(cs2, inline=True, colors='k')

    ax[1].set_title('function by contourf')

    fig.savefig('function2.jpg', dpi=600)

    plt.show()

    # 三维的函数图像

    fig1 = plt.figure()

    ax = Axes3D(fig1)

    cd = ax.plot_surface(xx, yy, f_xy, cmap='RdBu_r')

    ax.set_title('the function of 3D image')

    fig1.colorbar(cd)

    fig1.savefig('image3d.jpg', dpi=600)

    plt.show()

    通过函数图像的可视化,可以看到线性方程组本身是对称正定的,那么对应的二次型的函数图像则必定都位于在$z= 0 $平面以上,函数均是大于0的,且从函数的等值线及三维图像中可以看到是存在极小值的。明白了这个关系之后,对于梯度下降法或者共轭梯度法应用于方程组(系数矩阵为对称正定)的求解就不足为奇了。

    Gradient descent method求解方程组

    共轭梯度法的实现全来自于或者说是基于梯度下降法,因此,理解梯度下降法求解线性方程系统很有必要。在之前的梯度下降法的原理中,已经明白了这个方法对于函数极值的求解,是通过沿着梯度方向负向进行的,每一次的计算实质上获得是自变量的增量向量,而现在面对方程组时,由于之前的二次型的内容已经得出了线性方程组与多元函数极值之间的关系,即$f'(x_{i}) = Ax_{i} - b = 0$。那么梯度下降法可以被这样描述与构建:每一次产生的近似值$x_{i}$向量,近似解自然会有一个残差向量(residual vector)为$r_{i}$, 下表$i$表示第几步的迭代。

    $$

    r_{i} = b - A x_{i}\\

    r_{i} = - f'(x_{i}) \tag{4}

    $$

    式(4)中的$r_{i}$不就是梯度向量嘛,表示了点下一步搜索或者移动的方向,当然要想确定自变量的增量大小,还需要知道步长,这一点在梯度下降原理中已经给出了解释;在这里步长设定为$\alpha_{i}$,所以,产生的新的点$x_{i + 1}$为:

    $$

    x_{i + 1} = x_{i} + \alpha_{i} r_{i} \tag{5}

    $$

    到这里好像似乎与之前的原理中阐述并无大的区别,实际上现在问题在于如何确定步长$\alpha_i$,而在之前的介绍中步长是一个固定的,对于迭代计算来说,固定的步长的计算效率显然是不能够令人满意的。对于函数的极小值的确定来说,在这里每一步获得的新的函数值为$f(x_{i+1})$,对于整个求解过来说,只要这个更新的函数值达到最小值那么任务就算是完成了,依据式(5),函数值可以看做是关于步长的函数,为了获得极小值,通过导数为0来完成,即

    $$

    \frac{\mathrm{d}f'(x_{i + 1})}{\mathrm{d} \alpha_{i}} = 0 \tag{6}

    $$

    结合式(5)来看,这是一个复合函数,利用链式求导法则,得到

    $$

    式(6) = f'(x_{i+1}) ^{\mathrm{T}} \cdot \frac{\mathrm{d}x_{i + 1}}{\mathrm{d}\alpha_{i}} = 0 \tag{7}

    $$

    这个式子中,$f'(x)$不就是梯度向量嘛,依据式(4)给出的定义,在这里它也称作残差向量$-r_{i + 1}$,$x_{i + 1}$求导的结果可以根据式(5)为$r_{i}$,所以,要保证函数取得极值,就要使得前后两步的梯度向量(也就是残差向量$r$)保证正交!这里写为内积形式:

    $$

    r_{i + 1} ^ {\mathrm{T}} \cdot r_{i} = 0 \tag{8}

    $$

    那么依据这个关系,带入残差向量的表达,得到:

    $$

    (b- A x_{i + 1}) ^ {\mathrm{T}} \cdot r_{i} = 0\\

    (b - A (x + \alpha_{i} r_{i}))^{\mathrm{T}} \cdot r_{i} = 0

    \tag{9}

    $$

    通过这个式子就可以得到步长$\alpha_{i}$的计算公式:

    $$

    \alpha_{i} = \frac{r_{i} ^{\mathrm{T}} r_{i}}{r_{i} ^{\mathrm{T}} A r_{i}} \tag{10}

    $$

    到这里,所有量的表达都得到了解决,那么对于线性方程组系统来说,梯度下降法的迭代过程为:

    $$

    repeat:\\

    r_{i} = b - A x_{i}\\

    \alpha_{i} = \frac{r_{i}^{\mathrm{T}} r_{i}}{r_{i} ^{\mathrm{T}} A r_{i}}\\

    x_{i + 1} = x_{i} + \alpha_{i} r_{i}\\

    until: x_{i}的增量幅度满足设定精度\\

    end

    \tag{11}

    $$

    迭代计算的代码如下:import numpy as np

    # 每一步的残差向量r0都会被以行向量的形式存放在r矩阵中

    def gd_method(A, b, x0, tol, max_time):

    r0 = b - np.dot(A, x0)

    number = 0

    r = np.zeros([max_time, b.shape[0]])

    x = np.zeros([max_time, b.shape[0]])

    while np.linalg.norm(r0, np.inf) > tol:

    r0 = b - np.dot(A, x0)

    r[number, :] = r0

    alpha = np.dot(r0.T, r0) / np.dot(np.dot(r0,A), r0)

    x0 = x0 + alpha * r0

    x[number, :] = x0

    number += 1

    if number > max_time:

    print('warning:the result is not stable!')

    return

    return x, number, r

    这里对一个$3\times3$大小的对称正定矩阵进行验算,梯度下降算法迭代了67次得到结果,$[0.99983945,0.99976565, 1.99978575]$,而这个方程组的精确解为$[1,1,2]$。可以迭代结果还是比较可信与可靠的。A = np.array([[4, -2, -1], [-2, 4, -2], [-1, -2, 3]])

    b = np.array([0, -2, 3])

    x0 = np.array([1, 1, 1])

    x_gd = gd_method(A,b, x0, 0.0001,500)

    print(x_gd[0][x_gd[1]-1,:])

    print(x_gd[1])

    print(x_gd[2])[0.99983945 0.99976565 1.99978575]

    67

    [[-1. -2. 3. ]

    [-0.39130435 0.43478261 0.15942029]

    [ 0.09201888 0.02436289 0.15942029]

    ...

    [ 0. 0. 0. ]

    [ 0. 0. 0. ]

    [ 0. 0. 0. ]]x = np.linalg.solve(A,b)

    print(x)[1. 1. 2.]

    依据上述的分析以及这个算例的计算,实际上对于线性方程组迭代的迭代计算会有一个新的认识。首先,从最为原始的迭代的角度来看,通过每一迭代步中都会产生一个残差向量,似乎在方程组的求解角度来看,残差向量仅仅也就表示一个误差吧,好像也并没有很直观的一个认识。但是从梯度下降的算法角度来看,线性方程组迭代步中的残差向量实际上就是二次型(函数)的负梯度向量,那么也就是说,这个向量$r$不断地在修正点移动的方向,使得移动方向不断的向函数的极小值点处靠近!由于gd_method函数中会记录每一步的参差向量,因为这个是一个3个未知数的方程组,因此残差向量$r$本身就是一个$R^{3}$空间中的向量,同时返回得到的$x$矩阵中存放着每一步的近似解,这样进行可视化,可以看出梯度下降中前后的参差向量确实是正交的。因为三元函数无法进行可视化,不是很容易看出梯度下降的搜寻过程。可能二元函数的化对于理解更方便吧。xs = x_gd[0][0:x_gd[1]-1,:]

    plt.style.use('ggplot')

    fig = plt.figure()

    ax = Axes3D(fig)

    ax.plot(xs[:,0], xs[:,1], xs[:,2], marker='o')

    ax.set_title('point in every step')

    ax.set_xlabel('x')

    ax.set_ylabel('y')

    ax.set_zlabel('z')

    fig.savefig('gdimage.jpg', dpi=600)

    plt.show()

    在上述的分析中,我们会发现实际好像似乎梯度下降这个迭代算法的效率并不是很高,很简单的方程组都要迭代67次,显然这样的效率去求解大型的稀疏方程组是不能令人满意的。这个方法之所以没有被广泛地应用,主要原因是因为每一步的残差向量也就是梯度向量,都必须与相邻步的残差向量保持正交,如果为了便于理解,以二维空间为例,就是每一次的搜索(移动)方向都是相互垂直的,那么间隔的残差向量必定平行,也就说从初始值向精确解移动是呈折线的方式,一个搜索方向会被重复了好多次,这样使得迭代效率并不高。但是梯度下降法实现的思想确是很有意义的,为函数极值的确定与线性方程组的求解之间建立了联系。真正有意义的方法是在这个方法基础上实现的共轭梯度法,它对搜索方向进行了共轭处理,在理解了梯度下降的原理后,再对共轭梯度法进行学习可能更好一点。

    本站文章如非特别说明,均为原创。未经本人同意,请勿转载!转载请务必注明出处!

    展开全文
  • %功能: 用FR共轭梯度求解无约束问题: min f(x) %x0是初始点, fun, gfun分别是目标函数和梯度 %输出: x, val分别是近似最优点和最优值, k是迭代次数. x0 =[2 2]'; maxk=5000; %最大迭代次数 rho=0.6;sigma=0.4; k=0...
  • 共轭梯度学习笔记

    千次阅读 2019-06-27 20:07:28
    共轭梯度是介于最速下降与牛顿之间的一个方法,它仅需一阶导数信息,但克服了最速下降收敛慢的缺点,又避免了牛顿需要存储和计算Hesse矩阵并逆的缺点。共轭梯度的基本思想是把共轭性与最速下降相...
  • main为主函数 fun gfun ggfun分别为输入的...分别为最速下降 牛顿(阻尼)共轭梯度 以及 拟牛顿 F1-4为下降的图示 可以看到牛顿和拟牛顿收敛速度最快 但是牛顿需要矩阵的逆 在实际中 运算量可能较大
  • 优化对象:凸函数 梯度下降 顾名思义,就是沿着与梯度相反...在某一个特定点处将函数泰勒展开(仅保留至二次项),用泰勒展开的二次函数代替原函数最小点。求解最小点的方法就是:令二次函数对变量x的一阶导数等
  • 最优化方法-共轭梯度

    千次阅读 2019-03-31 17:16:23
    其基本思想是把共轭性与最速下降方法相结合,利用已知点处的梯度构造一组共轭方向,并沿这组方向进行搜素,出目标函数的极小点。根据共轭方向基本性质,这种方法具有二次终止性。 对于二次凸函数的共轭梯度: ...
  • %共轭梯度法求二次凸函数的最小值(FR)%此代码不适用于一般函数的最小值function [minf]=myCGMsyms x1 x2;f=fun([x1;x2]);%函数表达式A=[4 2;2 2];%矩阵A即为上面函数的二次型中的参数Ax0=[0 0]';%初始点mu=10^(-...
  • 共轭梯度(Conjugate Gradient)

    千次阅读 2019-03-03 16:34:18
    共轭梯度(Conjugate Gradient)是介于最速下降与牛顿之间的一个方法,它仅需利用一阶导数信息,但克服了最速下降收敛慢的缺点,又避免了牛顿需要存储和计算Hesse矩阵并逆的缺点,共轭梯度不仅是解决...
  • 共轭梯度目录1 共轭梯度原理2 共轭梯度代码实现参考文献和网址 目录 1 共轭梯度原理 每次将一个方向优化到了极小,后面的优化过程将不再影响之前优化方向上的极小值,所以理论上对N维问题极小只用对N个方向...
  • 共轭梯度是介于梯度下降与牛顿之间的一个方法,它仅需利用一阶导数信息,既克服了梯度下降收敛慢的缺点,又避免了牛顿需要存储和计算海塞矩阵并逆的缺点,共轭梯度不仅是解决大型线性方程组最有用的...
  • d0=-g0;%搜索方向alpha=-(g0'*d0)/(d0'*h0*d0);%步长 xk=x0+alpha*d0;...%搜索方向时的系数 dk=-gk+beta*d0;%下一个方向 x0=xk;%更新点g0=gk;%更新所在点的梯度 d0=dk;%更新方向 while g0'*g0>...
  • %精度 %共轭梯度 for i=1:2 if g==0 disp(‘x‘) x break else b=(g‘*a*d)/(d‘*a*d)%β d=-g+b*d; d11=subs(d(1),{x1,x2},{x(1) x(2)}); d22=subs(d(2),{x1,x2},{x(1) x(2)}); d=[d11;d22] af=(-g‘*d)/(d‘*a*d)%...
  • 简单回顾一下,线性共轭梯度是一种不需要矩阵逆或矩阵分解,以迭代的方式求解线性方程的方法,而且可以保证在N次迭代内收敛,其中N是变量的维度。对于维度较大的问题,辅以良好的预处理步骤,算法可以在很少的步...
  • 大规模多输入多输出系统中,最小均方误差信号检测算法是近似最优的,但由于其涉及矩阵逆,计算复杂度随着天线数量增加呈指数增长。提出了低复杂度的预处理共轭梯度...相比直接用共轭梯度,能够更快收敛到最佳值。
  • 牛顿主要解决非线性优化问题,其收敛速度比梯度下降速度更快,其需要解决的问题可以描述为:对于目标函数f(x),在无约束条件下他的最小值minf(x)。 牛顿的主要思想:在现有的极小值估计值的附近对f(x...
  • Optimization_Algorithm梯度下降、牛顿共轭梯度等matlab和python程序:一个空间曲面(3维)的极值点。梯度下降算法速度较慢、迭代次数较大,并且最后的结果是近似值;牛顿利用函数的二阶泰勒展开近似,可以...

空空如也

空空如也

1 2 3 4 5 ... 7
收藏数 136
精华内容 54
关键字:

共轭求法