用MATLAB实现共轭梯度法求解实例

求解无约束优化问题的共轭梯度法

求解无约束优化问题的共轭梯度法李芳梅，姚瑞哲指导教师：李良摘要：本文主要针对无约束优化问题，利用共轭梯度法（CG方法）求解此类问题，并得出其迭代次数及问题的解。

论文对此种方法给出了具体事例，并对例子进行了matlab软件实现。

1.引言共轭梯度法时介于最速下降法与牛顿法之间的一个方法。

它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点。

共轭梯度法最早是由Hestenes和Stiefel（1952）提出来的，用于解正定系数矩阵的线性方程组。

由于共轭梯度法不需要矩阵存储，且有较快的收敛速度和二次终止性等优点，现已广泛应用与实际问题中。

2.共轭梯度法共轭梯度法是共轭方向法的一种。

所谓共轭方向法就是其所有的搜索方向都是相互共轭的方法。

定义设G是n×n对称正定矩阵，d1，d2是n维非零向量。

如果d1T G d2=0，则称向量d1和d2是G-共轭的，简称共轭的。

设d1，d2，…，d m是R n中任一组非零向量，如果d i T G d j=0 （i≠j），则称d1，d2，…，d m是G-共轭的，简称共轭的。

显然，如果d1，d2，…，d m是共轭的，则它们是线性无关的。

算法（共轭梯度法）步1.初始步：给出x0，ε>0,计算r0=G x0-b，令d0=-r0，k:=0. 步2.如果‖r k‖≤ε，停止.步3.计算αk=r k T r k/d k T Gd k,x k+1=x k+αk d k，r k+1=r k+αk Gd k，βk=r k+1T r k+1/r k T r k，d k+1=-r k+1+βk d k.步4.令k:=k+1，转步2.1.共轭梯度法的matlab实现（数值例子）首先建立如下的m.文件function[x,iter]=cgopt(G,b,x0,max_iter,tol)x=x0;fprintf('\n x0=');fprintf('%10.6f',x0);r=G*x-b;%残差d=-r;for k=1:max_iterif norm(r)<=toliter=k-1;fprintf('\n Algorithm finds a solution!');return endalpha=(r'*r)/(d'*G*d);%收敛速度 xx=x+alpha*d; rr=r+alpha*G*d; beta=(rr'*rr)/(r'*r); d=-rr+beta*d; x=xx; r=rr;fprintf('\n x%d=',k); fprintf('%10.6f',x); enditer=max_iter; return然后建立CG_main 的m.文件，带入数值例子function CG_main()12345123451234512345123451023412923277351432312174351512x x x x x x x x x x x x x x x x x x x x x x x x x ++++=⎧⎪+-+-=-⎪⎪-++-=⎨⎪+++-=-⎪---+=⎪⎩ G=[10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15]; b=[12 -27 14 -17 12]'; x0=[0 0 0 0 0]'; max_iter=1000; tol=1e-6;fprintf('\n');fprintf('Conjugate Gradiant Method:\n'); fprintf('=================\n');[x,iter]=cgopt(G,b,x0,max_iter,tol);fprintf('Iterative number:\n %d\n',iter); fprintf('Solution:\n'); fprintf('%10.6f',x);fprintf('\n================\n');4.结论实际上，共轭梯度法是最速下降法的一种改进。

matlab梯度算法

matlab梯度算法Matlab梯度算法在数学和计算机科学中，梯度是指一个多元函数在某一点上的变化率或斜率。

梯度算法是一种优化算法，用于找到函数的最小值或最大值。

在Matlab中，有多种方法可以使用梯度算法来优化函数，包括梯度下降和共轭梯度法。

本文将详细介绍Matlab中的梯度算法，并逐步讲解其原理和应用。

I. 梯度下降法梯度下降法是一种基于迭代的优化算法，通过计算函数的梯度来更新参数的值，以逐步接近函数的最小值。

在Matlab中，可以使用"gradientDescent"函数来实现梯度下降法。

1. 实现梯度下降法首先，我们需要定义一个优化目标函数，例如：f(x) = x^2 + 2x + 1。

然后，定义其梯度函数为g(x) = 2x + 2。

接下来，我们可以使用以下代码来计算梯度下降：matlab定义优化目标函数f = (x) x^2 + 2*x + 1;定义梯度函数g = (x) 2*x + 2;初始化参数x0 = 0;设置学习率和迭代次数alpha = 0.01;iterations = 100;梯度下降法for i = 1:iterationsx0 = x0 - alpha * g(x0);end打印最优解disp(['Optimal solution: ', num2str(x0)]);在这个例子中，我们使用了学习率(alpha)为0.01，迭代次数(iterations)为100。

通过不断更新参数x0的值，最终得到了最优解。

2. 梯度下降法的原理梯度下降法的核心思想是利用函数在当前点的梯度信息来更新参数的值，以便能够向着函数的最小值前进。

具体来说，算法的步骤如下：a. 初始化参数的值：选择一个初始参数的值作为起始点。

b. 计算梯度：计算函数在当前点的梯度，即求解函数关于参数的偏导数。

c. 更新参数：根据当前点的梯度和学习率，通过减去梯度的乘积来更新参数的值。

基于matlab平台的三种迭代法求解矩阵方程

数值分析第二次作业学院：电子工程学院基于matlab平台的三种迭代法求解矩阵方程组求解系数矩阵由16阶Hilbert方程组构成的线性方程组的解，其中右端项为[2877/851,3491/1431,816/409,2035/1187,2155/1423,538/395,1587/1279,573/502,947 /895,1669/1691,1589/1717,414/475,337/409,905/1158,1272/1711,173/244].要求：1）Gauss_Sedel迭代法；2）最速下降法；3）共轭梯度法；4）将结果进行分析对比。

解：根据题目要求，编写了对应算法的matlab程序，求解结果如下：（求解精度为10e-4，最大迭代次数1000）1、方程的解：如下图1所示图1 三种方法求解的结果对比图2 Gause_Sedel算法收敛特性图3 最速下降法收敛特性图3 共轭梯度法收敛特性从图中可以看到，在相同的最大迭代次数和预设求解精度条件下，共轭梯度算法仅需要4次迭代便可求出方程组的解，耗时0.000454秒，而且求出解的精度最高；Gauss_Sedel方法需要465次迭代，耗时0.006779秒，求解精度最差；最速下降法需要398次迭代，耗时0.007595秒，求解精度与共轭梯度算法差不多，因此两者求出的解也几乎相同。

从中可以得出结论，共轭梯度算法无论从求解精度还是求解速度上都优于其他两种，最速下降法在求解精度上几乎与共轭梯度算法持平，但求解速度更慢。

Gauss_Sedel方法在求解精度和速度两方面都最差。

具体的解为:Gauss_Sedel迭代法:(共需465次迭代，求解精度达到9.97e-5) X=[0.995328360833192 1.01431732497804 1.052861239300110.934006974137998 0.931493373808838 0.9665081384030661.00661848511341 1.03799789809258 1.051806903036541.06215849948572 1.04857676431223 1.028561990411131.01999170162638 0.971831831519515 0.9525261666348130.916996019179182].最速下降法:(共需398次迭代，求解精度达到9.94e-5)X=[0.998835379744322 1.01507463472900 0.9825890937201850.980191460759243 0.991245169713628 1.003780222253291.01350884374478 1.01928337905816 1.020859096651941.01930314197028 1.01444777381651 1.007040589892970.998384452250809 0.987399404644377 0.9757678149709120.963209150871750].共轭梯度法:(共需4次迭代，求解精度达到3.98e-5)X=[0.996472751179456 1.02707840189049 0.9776233734098530.973206695321590 0.986133032967607 1.001289025642341.01322158496914 1.02047386502293 1.023009050605651.02163015083975 1.01678089454399 1.009203108638740.999772406055155 0.988443827498859 0.9760941924969490.962844741655005].Matlab程序主程序:clc;clear;%% 本程序用于计算第二次数值分析作业，关于希尔伯特矩阵方程的解，用三种方法，分析并比较，也可推广至任意n维的矩阵方程%%A=hilb(16); %生成希尔伯特系数矩阵b=[2877/851;3491/1431;816/409;2035/1187;2155/1423;538/395;1587/1279;573/502;947/895;166 9/1691;1589/1717;414/475;337/409;905/1158;1272/1711;173/244]; %右端向量M=1000; %最大迭代次数err=1.0e-4; %求解精度[x,n,xx,cc,jingdu]=yakebi_diedai(A,b,err,M); % 雅克比算法求解tic;[x1,n1,xx1,cc1,jingdu1]=gauss_seidel(A,b,err,M); % gauss_seidel算法求解toc;tic;[x2,n2,xx2,jingdu2]=zuisuxiajiangfa(A,b,err,M); % 最速下降法求解toc;tic;[x3,flag,jingdu3,n3]=bicg(A,b,err); % matlab内置双共轭梯度算法求解toc;tic;[x4,xx4,n4,jingdu4]=con_grad(A,b,err,M); % 教材共轭梯度算法求解toc;%% 计算相应结果，用于作图%%num=[1:16]';jie=[num,x1,x2,x4]; % 三者的解对比% 三者的收敛情况对比num1=[1:n1]';fit1=[num1,jingdu1'];num2=[1:n2]';fit2=[num2,jingdu2'];num4=[1:n4]';fit4=[num4,jingdu4'];子函数1(Gause_Sedel算法):function [x,n,xx,cc,jingdu] = gauss_seidel(A,b,err,M)% 利用迭代方法求解矩阵方程这里是高斯赛尔得迭代方法% A 为系数矩阵b 为右端向量err为精度大小返回求解所得向量x及迭代次数% M 为最大迭代次数cc 迭代矩阵普半径jingdu 求解过程的精度n 所需迭代次数xx 存储求解过程中每次迭代产生的解for ii=1:length(b)if A(ii,ii)==0x='error';break;endendD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=(D-L)\U;cc=vrho(B); %迭代矩阵普半径FG=(D-L)\b;x0=zeros(length(b),1);x=B*x0+FG;k=0;xx(:,1)=x;while norm(A*x-b)>errx0=x;x=B*x0+FG;k=k+1;xx(:,k+1)=x;if k>=Mdisp('迭代次数太多可能不收敛！');break;endjingdu(k)=norm(A*x-b);endend子函数2(最速下降算法):function [x,n,xx,jingdu]=zuisuxiajiangfa(A,b,eps,M)% 利用迭代方法求解矩阵方程这里是最速下降迭代方法% A 为系数矩阵b 为右端向量err为精度大小返回求解所得向量x及迭代次数% % M 为最大迭代次数jingdu 求解过程的精度n 所需迭代次数xx 存储求解过程中每次迭代产生的解x0=zeros(length(b),1);r0=b-A*x0;t0=r0'*r0/(r0'*A*r0);x=x0+t0*r0;r=b-A*x;xx(:,1)=x;k=0;while norm(r)>epsr=r;x=x;t=r'*r/(r'*A*r);x=x+t*r;r=b-A*x;k=k+1;xx(:,k+1)=x;if k>=Mdisp('迭代次数太多可能不收敛！');break;endn=k;jingdu(k)=norm(r);endend子函31(共轭梯度法):function [x,xx,n,jingdu]=con_grad(A,b,eps,M)% 利用迭代方法求解矩阵方程这里是共轭梯度迭代方法% A 为系数矩阵b 为右端向量err为精度大小返回求解所得向量x及迭代次数% M 为最大迭代次数jingdu 求解过程的精度n 所需迭代次数xx 存储求解过程中每次迭代产生的解x0=zeros(length(b),1);r0=b-A*x0;p0=r0;% t0=r0'*r0/(r0'*A*r0);% x=x0+t0*r0;% xx(:,1)=x;k=0;x=x0;r=r0;p=p0;while norm(r)>epsx=x;r=r;p=p;afa=r'*r/(p'*A*p);x1=x+afa*p;r1=r-afa*A*p;beta=r1'*r1/(r'*r);p1=r1+beta*p;x=x1;r=r1;p=p1;k=k+1;xx(:,k)=x;if k>=Mdisp('迭代次数太多可能不收敛！');break;endn=k;jingdu(k)=norm(r);endend。

共轭梯度法matlab最优化问题

共轭梯度法是一种在求解最优化问题时常用的算法。

下面是一个在MATLAB 中实现共轭梯度法的简单示例。

请注意，这个示例是为了教学目的而编写的，可能不适用于所有最优化问题。

首先，假设我们有一个目标函数f(x)，我们需要找到使得f(x) 最小化的x。

假设f(x) 是一个二次函数，形式为f(x) = x^T Ax + b^T x + c，其中A 是对称正定矩阵，b 和c 是常数向量和标量。

以下是一个使用MATLAB 实现共轭梯度法的示例代码：```matlabfunction [x, iter] = conjugate_gradient(A, b, x0, tol, max_iter)% A -目标函数的系数矩阵% b -目标函数的常数向量% x0 -初始解% tol -容忍的误差% max_iter -最大迭代次数x = x0;r = b - A*x;p = r;iter = 0;while (norm(r) > tol) && (iter < max_iter)Ap = A*p;alpha = (p'*r) / (p'*Ap);x = x + alpha*p;r = r - alpha*Ap;beta = (r'*r) / (p'*r);p = r + beta*p;iter = iter + 1;endend```这个函数接受一个对称正定矩阵A，一个常数向量b，一个初始解x0，一个容忍的误差tol，和一个最大迭代次数max_iter 作为输入，并返回最优解x 和迭代次数iter。

注意，这个函数没有包括一些可能的特殊情况处理，例如如果A 是奇异的或者接近奇异的，那么这个函数可能无法正确地收敛。

在使用这个函数之前，你可能需要根据你的具体问题对其进行一些修改和增强。

数值计算实例MATLAB实现(附带详细源码)

数值计算实例MATLAB实现附带详细源码1.在化学反应中，A 的一个分子和 B 的一个分子结合形成物质 C 的分子。

若在时刻t 时，物质 C 的浓度为() y t ，则其是下述初值问题的解()()() ,00y k a y b y y '=--=其中k 为正常数，a 和 b 分别表示 A 和 B 的初始浓度。

假设k = 0.01, a =70毫摩/升， b = 50 毫摩/升. 该方程的真解为0.20.2350(1)()75t te y t e---=- （1）自己编写程序，使用四阶经典Runge-Kutta （龙格-库塔法），以步长为0.5h =，在区间[0, 20]上给出() y t 的近似解；（2）列表给出真解和近似解的比较；（3）讨论当t →∞时，近似解的变化趋势，并分析该数值结果。

解：数学原理：四阶经典Runge-Kutta （龙格-库塔法）112341213243(22)6(,)(,)22(,)22(,)m m m m m m m m m m hu u k k k k k f t u h hk f t u k h hk f t u k k f t h u hk +=++++==++=++=++程序设计见附录结果如下表：（3）近似解变化趋势当t→∞时，由以下极限方程可知：0.20.2350(1)()75lim()tttey tey t--→∞⎧-=-⎪⎨⎪⎩随着t→∞，近似值越来越接近真实值，极限的真实值为50，lim()50ty t→∞=，变化趋势也可由一下曲线图表示：感想：四阶Runge-Kutta法计算的结果精度非常好，其结果与真实解误差不大。

2.考虑定义在闭区间[−5, 5]上的函数()2112()5f x x -=+ ；（1）利用等距节点构造次数分别为 n = 4,8,16, 32 的插值多项式()n p x ，并分别画()()()()481632,,,p x p x p x p x ；（2）利用chebyshev 零点构造次数分别为 n = 4,8,16, 32 的插值多项式()n pp x()()()()481632,,,pp x pp x pp x pp x ；（3）画出当 n = 32 时，两种插值多项式的比较图，误差图，并给出相应的误差估计;（4）在这个问题中能观察到龙格现象吗？解：数学原理：拉格朗日插值多项式：001122()()()()()n n n L x l x y l x y l x y l x y =+++011011()()()()(),0,1,2,()()()()k k n k k k k k k k n x x x x x x x x l x k n x x x x x x x x -+-+----==----0()()()nn n in k k k k k j k jj kx x L x l x y y x x ===≠-==-∑∑∏程序设计见附录(1) 利用等距节点构造次数分别为 n = 4,8,16, 32 的插值多项式如下： ()43240.00160.00.0640.60061400p x x x x x ++=++()876542830.00280.00640.02500.02500.00640.00260.000168.001p x x x x x x x x x ++++++=++()1615141312161110987654320.00210.00280.00410.0064 60.01120.02500.09290.09290.02050 0.01120.00640.00410.002.00160180.021.000p x x x x x x x x x x x x x x x x x ++++++++++++++=++()3231302928272632252423222120191817160001600018000210002400028000340004100050006400083001120016100250004350092902906029p x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .x x .=+++++++++++++++++151413121110987654320600929004350025000161001120008300064000500041000340002800024000210001800016x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .+++++++++++++++（2）利用chebyshev 零点构造次数分别为 n = 4,8,16, 32 的插值多项式如下：()43240.00160.00320.00320.0016x x p x x p x =++++()87654328+0.00190.00320.01080.01080.00320.00196=0.0.0106001pp x x x x x x x x x +++++++()161514131211109168765432=0.0016 0.0017 0.0019 0.00230.00320.00520.01080.0403 1.00000.04030.01080.00520.00320.00230.0019 0.0017 0.0016 pp x x x x x x x x x x x x x x x x x ++++++++++++++++()323130292827263225242322212019181700016000160001700017000190002100023000270003200040000520007100108001860040301428pp x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .x x =+++++++++++++++++16151413121110987654320142800403001860010800071000520004000320002700023000210001900017000170001600016.x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .+++++++++++++++++(3)两种插值多项式的比较误差图如下(a)等距插值误差 (b) chebyshev零点插值误差(4) 等距插值在高次插值中能观察到龙格现象，而chebyshev零点插值观察不到龙格现象。

用共轭梯度法解最小二乘问题

用共轭梯度法求解最小二乘问题摘要本文先讨论了求解对称正定线性方程组的共轭梯度法.然后对系数矩阵列满秩的线性方程组运用正则化方法将其转化为对称正定线性方程组后再运用实用共轭梯度法进行求解，最后举例并通过Matlab 程序实现其结果.关键词共轭梯度法；正则化方法；最小二乘问题；Krylov 子空间1 引言在实际的科学与工程问题中，常常将问题归结为一个线性方程组的求解问题，而求解线性方程组的数值解法大体上可分为直接法和迭代法两大类.直接法是指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法.因此，直接法又称为精确法.迭代法则是采取逐次逼近的方法，亦即从一个初始向量出发，按照一定的计算格式，构造一个向量的无穷序列,其极限才是方程组的精确解，只经过有限次运算得不到精确解.当线性方程组的系数矩阵为对称正定矩阵时，我们常用共轭梯度法（或简称CG 法）求解，目前有关的方法与理论已经相当成熟，并且已成为求解大型稀疏线性方程组最受欢迎的一类方法.2 最小二乘问题定义1[1] 给定矩阵m n A R ⨯∈，A 列满秩及向量m b R ∈，确定nx R ∈使得()()2222min min n ny Ry Rb Ax r x r y Ay b ∈∈-===-. 该为问题称为最小二乘问题，简记为LS （Least Squares ）问题，其中()r x 称为残向量.最小二乘问题的解x 又可称为线性方程组Ax b =，m n A R ⨯∈的最小二乘解，即x 在残向量()r x b Ax =-的2范数最小的意义下满足线性方程组 Ax b =，m n A R ⨯∈.3 共轭梯度法考虑线性方程组Ax b =的求解问题，其中A 是给定的n 阶对称正定矩阵，b 是给定的n 维向量，x 是待求解的n 维向量.为此,定义二次泛函()2T T x x Ax b x ϕ=-.定理1[1]设A 对称正定,求方程组Ax b =的解，等价于求二次泛函()x ϕ的极小点. 定理1表明，求解线性方程组问题就转化为求二次泛函()x ϕ的极小点问题.求解二次函数极小值问题，通常好像盲人下山那样，先给定一个初始向量0x ，确定一个下山方向0p ，沿着经过点0x 而方向为0p 的直线00x x p α=+找一个点1000x x p α=+，使得对所有实数α有()()00000x p x p ϕαϕα+≤+，即在这条直线上1x 使()x ϕ达到极小.然后从1x 出发，再确定一个下山的方向1p ，沿着直线11x x p α=+再跨出一步，即找到1α使得()x ϕ在2111x x p α=+达到极小：()()11111x p x p ϕαϕα+≤+.重复此步骤，得到一串012,,,ααα 和 012,,,p p p ，称k p 为搜索方向，k α为步长.一般情况下，先在k x 点找下山方向k p ，再在直线k k x x p α=+上确定步长k α使()(),k k k k k x p x p ϕαϕα+≤+最后求出1k k k k x x p α+=+.然而对不同的搜索方向和步长，得到各种不同的算法.由此，先考虑如何确定k α.设从k x 出发，已经选定下山方向k p .令()()k k f x p αϕα=+()()()2TT k k k k k k x p A x p b x p ααα=++-+()22TT k k k k k p Ap r p x ααϕ=-+，其中k k r b Ap =-.由一元函数极值存在的必要条件有()220TT k k k k f p Ap r p αα'=-=所确定的α即为所求步长k α，即T k kk Tk kr p p Ap α=. 步长确定后，即可算出1k k k k x x p α+=+.此时，只要0Tk k r p ≠，就有()()()()1k k k k k k x x x p x ϕϕϕαϕ+-=+-()2220T kk TTk k k k k kT k kr p p Ap r p p Ap αα=-=-<即()()1k k x x ϕϕ+<.再考虑如何确定下山方向k p .易知负梯度方向是()x ϕ减小最快的方向，但简单分析就会发现负梯度方向只是局部最佳的下山方向，而从整体来看并非最佳.故采用新的方法寻求更好的下山方向——共轭梯度法.定义2[2]若n 维非零向量,x y 满足0T x Ay =其中A 为n 阶对称正定矩阵，则称x 与y 是相互共轭（A -共轭）的. 下面给出共轭梯度法的具体计算过程:给定初始向量0x ，第一步仍选用负梯度方向为下山方向，即00p r =，于是有00010001000,,T T r r x x p r b Ax p Ap αα==+=-.对以后各步，例如第k+1步（k ≥1）,下山方向不再取k r ，而是在过点由向量k r 和1k p -所张成的二维平面21{|,,}k k k x x x r p R πξηξη-==++∈内找出使函数ϕ下降最快的方向作为新的下山方向k p .考虑ϕ在2π上的限制：()1,()k k k x r p ψξηϕξη-=++11()()Tk k k k k k x r p A x r p ξηξη--=++++12()Tk k k b x r p ξη--++.计算ψ关于,ξη的偏导得：()()11112,2,T T T k k k k k k T Tk k k k r Ar r Ap r r r Ap p Ap ψξηξψξηη----∂=+-∂∂=+∂其中最后一式用到了10Tk k r p -=，这可由k r 的定义直接验证.令0ψψξη∂∂==∂∂，即知ϕ在2π内有唯一的极小值点001k k k x x r p ξη-=++，其中0ξ和0η满足00101011,0.T T T k k k k k k T Tkk k k r Ar r Ap r r r Ap p Ap ξηξη----⎧+=⎨+=⎩ 由于0k r ≠必有00ξ≠，所以可取()0101k k k k p x x r p ηξξ-=-=+作为新的下山方向.显然，这是在平面2π内可得的最佳下山方向.令010k ηβξ-=，则可得 1111.T k k k T k k r Ap p Ap β----=-注：这样确定的k p 满足10Tk k p Ap -=，即k p 与1k p -是相互共轭的.总结上面的讨论，可得如下的计算公式：T k kk Tk kr p p Ap α= ， 1k k k k x x p α+=+， 11k k r b Ax ++=-， 1T k kk Tk kr Ap p Ap β+=-， 11k k k k p r p β++=+. 在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称的计算公式.首先来简化1k r +的计算公式：11()k k k k k k k k r b Ax b A x p r Ap αα++=-=-+=-.因为k Ap 在计算k α是已经求出，所以计算1k r +时可以不必将1k x +代入方程计算，而是从递推关系1k k k r b Ap α+=-得到.再来简化k α和k β的计算公式.此处需要用到关系式1110,T T T k k k k k k r r r p r p +-+=== 1,2,k =.从而可导出1111,T T k k k kr r r α+++=-，()111T TTk k k k k k k kkp Ap p r r p r αα+=-=()1111T T k k k k k k kkr r p r r βαα--=+=.由此可得,T k k k T k k r r p Ap α=， 11.T k k k T k kr r r r β++=.从而有求解对称正定方程组的共轭梯度法算法如下：0x =初值00r b Ax =-；0k =while 0k r ≠1k k =+if 1k = 00p r =else21122T Tk k k k k r r r r β-----= 1122k k k k p r p β----=+ end11111T Tk k k k k r r p Ap α-----= 111k k k k x x p α---=+ 111k k k k r r Ap α---=-endk x x =注：该算法每迭代一次仅需要使用系数矩阵A 做一次矩阵向量积运算. 定理2 [1]由共轭梯度法得到的向量组{}i r 和{}i p 具有如下基本性质：（1）0Ti j p r =， 0;i j k ≤<≤ （2）0Ti j r r =， i j ≠，0,;i j k ≤≤ （3）0Ti j p Ap =， i j ≠，0,;i j k ≤≤ （4）000{,,}{,,}(,,1)k k span r r span p p A r k κ==+，其中0000(,,1){,,,}k A r k span r Ar A r κ+=，通常称之为Krylov 子空间.下面给出共轭梯度法全局最优性定理：定理3[1]用共轭梯度法计算得到的近似解k x 满足()(){}00min :(,,)k x x x x A r k ϕϕκ=∈+或{}**00min :(,,)k AAx x x x x x A r k κ-=-∈+，其中Ax=，*x 是方程组Ax b =的解，0(,,)A r k κ是由所定义的Krylov 子空间.定理2表明，向量组0,,k r r 和0,,k p p 分别是Krylov 子空间0(,,1)A r k κ+的正交基和共轭正交基.由此可知，共轭梯度法最多n 步便可得到方程组的解*x .因此，理论上来讲，共轭梯度法是直接法.然而实际使用时，由于误差的出现，使k r 之间的正交性很快损失，以致于其有限步终止性已不再成立.此外，在实际应用共轭梯度法时，由于一般n 很大，以至于迭代()O n 次所耗费的计算时间就已经使用户无法接受了.因此，实际上将共轭梯度法作为一种迭代法使用，而且通常是k r 是否已经很小及迭代次数是否已经达到最大允许的迭代次数max k 来终止迭代.从而得到解对称正定线性方程组的实用共轭梯度法，其算法如下：x =初值0;k =;r b Ax =-T r r ρ=while)()max2band k kε><1k k =+if 1k =p r = else;p r p βρρβ==+ end;;TAp p x x p ωαρωα===+ ;;Tr r r r αωρρρ=-== end算法中，系数矩阵A 的作用仅仅是用来由已知向量p 产生向量Ap ω=，这不仅可以充分利用A 的稀疏性，而且对某些提供矩阵A 较为困难而由已知向量p 产生向量Ap ω=又十分方便的应用问题有益.4 共轭梯度法求解最小二乘问题的正则化方程组（法方程组）定理4[1] 当A 列满秩时，求最小二乘问题的解等价于解T TA Ax A b =.应用共轭梯度法于对称正定方程组T T A Ax A b =来求方程组Ax b =，m nA R ⨯∈()m n ≥且A 列满秩的最小二乘解，即为Krylov 子空间法中的正则化方法.由A 列满秩有T A A 对称正定，则方程组T T A Ax A b =，m nA R ⨯∈()m n ≥存在唯一解.下面给出其实用共轭梯度法的详细算法且算法中不出现计算TA A 情形:x =初值0;;;;T T k b A b m Ax r b A m ====-while)()max2band k kε><1k k =+if 1k =p r = else;p r p βρρβ==+end;;;T Tn Ap A n p x x p ωαρωα====+ ;;Tr r r r αωρρρ=-==end注：算法中采用了两次矩阵向量积来避免出现计算T A A 情形.算例编写实用共轭梯度法的Matlab 程序求解方程组TTA Ax A b =，其中11112231A -⎡⎤⎢⎥-⎢⎥=⎢⎥-⎢⎥-⎣⎦, 1234b ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦. 解先建立conjgrad.m 文件，内容如下:function [x]=conjgrad(A,b,x) A=[1 -1;-1 1;2 -2;-3 1]; b=[1;2;3;4]; x=rand(2,1); k=0;b=A'*b;m=A*x; r=b-A'*m; T=r'*r;While norm(r)>1e-10*norm(b) k=k+1; if k==1 p=r;else B=T/t; p=r+B*p;endn=A*p;w=A'*n; a=T/(p'*w); x=x+a*p; r=r-a*w; t=T; T=r'*r; end然后运行后得 ans =-2.4167 -3.2500即有方程组的数值解 2.41673.2500x -⎡⎤=⎢⎥-⎣⎦.而其精确解可由如下方法求得: 111123111591121229731T A A -⎡⎤⎢⎥----⎡⎤⎡⎤⎢⎥==⎢⎥⎢⎥⎢⎥----⎣⎦⎣⎦⎢⎥-⎣⎦，11123271121314T A b ⎡⎤⎢⎥---⎡⎤⎡⎤⎢⎥==⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎢⎥⎣⎦，则有121597971x x --⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦，解得122912134x x ⎡⎤-⎢⎥⎡⎤=⎢⎥⎢⎥⎣⎦⎢⎥-⎢⎥⎣⎦，即方程精确解为*2912134x ⎡⎤-⎢⎥=⎢⎥⎢⎥-⎢⎥⎣⎦，故可验证通过Matlab 程序求得的数值解2.41673.2500x -⎡⎤=⎢⎥-⎣⎦满足精度要求.5 总结本文首先给出最小二乘问题的定义，随后从盲人下山法开始讨论了共轭梯度法的具体推导过程及其相关性质与算法.继而重点给出正则化方法的实用共轭梯度算法并举例进行检验.最后，需要说明虽然正则化方法是求一般线性方程组Ax b =，m nA R⨯∈()m n ≥且A 列满秩的最小二乘解的一种方法且简单易行，但是也有许多不足之处，如m n >时一般无解；TA A 形成时运算量大，A 中某些信息会丢失；当A 病态时其收敛性速度由于222()()T A A A κκ=很大变得非常之慢等，故为了避免正则化方法的缺点，还可运用残量极小化方法或残量正交化方法等更好的方法来解决此类问题.参考文献[1] 徐树方，高立，张平文.数值线性代数[M].北京：北京大学出版社，2000.139--151. [2] 施光燕，董加礼.最优化方法[M].北京：高等教育出版社，1999.47--52.。

双共轭梯度法matlab_概述及解释说明

双共轭梯度法matlab 概述及解释说明1. 引言1.1 概述引言部分将介绍“双共轭梯度法(Matlab)”，该方法是一种用于解决优化问题的迭代算法，常用于求解大规模线性方程组、最小二乘问题和非线性最优化等。

本文将全面讲解双共轭梯度法的基础知识、算法流程及其在MATLAB中的应用与实现。

1.2 文章结构本文按照以下方式组织：- 第二节将介绍双共轭梯度法的基础知识，包括梯度下降法、共轭梯度法和双共轭梯度法的简介。

- 第三节将详细阐述双共轭梯度法的算法流程及具体步骤解释，包括初始化步骤、迭代更新步骤以及收敛准则和结束条件设定。

- 第四节将以MATLAB为工具，展示双共轭梯度法在实践中的应用与实现举例。

这一部分将给出MATLAB代码编写指导原则，描述一个示例问题，并说明求解过程和结果分析。

- 最后一节是结论与展望，总结了双共轭梯度法的优点和局限性，并提供对未来可能的研究方向的展望和建议。

1.3 目的本文旨在介绍双共轭梯度法的原理、算法流程及其在MATLAB中的实际应用。

读者将通过本文了解如何使用该方法解决优化问题，并深入理解算法背后的理论基础。

同时，本文还将探讨双共轭梯度法存在的局限性，并展望未来可能的研究方向，为相关领域的研究提供参考。

2. 双共轭梯度法基础知识2.1 梯度下降法简介梯度下降法是一种优化算法，用于求解无约束问题的最小值。

其基本思想是通过沿着目标函数的负梯度方向进行迭代更新，以逐步减小目标函数值。

具体而言，对于一个可微分的目标函数f(x)，初始值$x_0$被选为起点，然后通过以下公式进行迭代更新：$$x_{k+1} = x_k - \alpha_k \nabla f(x_k)$$其中$\alpha_k$是步长或学习率，$\nabla f(x_k)$表示在点$x_k$处的梯度（即函数$f(x)$在$x_k$处的导数）。

该过程将重复执行直到满足预设的终止条件。

2.2 共轭梯度法简介共轭梯度法是一种高效的迭代方法，用于解决对称正定线性系统的问题。

计算方法——共轭梯度法求解线性方程组

k T
n-1 do
k
r r ； k T k d Ad
x
k 1
x k d ；
k k
1
1 共轭梯度法求解线性方程组
③
r
k 1
b Ax
k 1
k 1
；
④ 若 || r k 1 || 或 k 1 n ，则输出近似解 x(k+1)，停止；否则，转⑤； ⑤ ⑥
18 9 27 9 1 18 45 0 45 2 ，b A 9 16 0 126 9 27 45 9 135 8
由于该方程组的系数矩阵以及右端项都比较简单，因此采用从 matlab 命令窗口手动输入的方式来输入数据，取计算精度为 10-6，运行过程及结果如图 2 和图 3（由于迭代的初始值为随机产生，因此每次得到的残量图会有所不同，但最终都趋于 0）：
|| r ||2 k k 22 ； || r ||2
d
k 1
r
k 1
k d ；
k
end do
图 1 共轭梯度法求解线性方程组程序框图
1.2 程序使用说明共轭梯度法求解线性方程组的 matlab 程序见附录 1，该程序可以求解系数矩阵为对称正定矩阵的线性方程组。在使用该程序时，可将程序复制到 matlab 命令窗口中直接运行或者复制到编辑窗口中保存运行，运行时刻根据提示输入，直至得到结果。开始运行程序时，会出现提示“请选择系数矩阵、右端项以及系数矩阵阶数的输
图 4 命令窗口显示的运行结果
1.6 1.4 1.2 1
残量
0.8 0.6 0.4 0.2 0 0 20 40 60 80 100

最优化牛顿法最速下降法共轭梯度法matlab代码

最优化⽜顿法最速下降法共轭梯度法matlab代码⽜顿法迭代公式：(1)2()1()[()]()k k k k x x f x f x +-=-??Matlab 代码：function [x1,k] =newton(x1,eps)hs=inline('(x-1)^4+y^2'); 写⼊函数ezcontour(hs,[-10 10 -10 10]); 建⽴坐标系hold on; 显⽰图像syms x y 定义变量f=(x-1)^4+y^2; 定义函数grad1=jacobian(f,[x,y]); 求f 的⼀阶梯度grad2=jacobian(grad1,[x,y]); 求f 的⼆阶梯度k=0; 迭代初始值while 1 循环grad1z=subs(subs(grad1,x,x1(1)),y,x1(2)); 给f ⼀阶梯度赋初值 grad2z=subs(subs(grad2,x,x1(1)),y,x1(2)); 给f ⼆阶梯度赋初值x2=x1-inv(grad2z)*(grad1z)'; 核⼼迭代公式if norm(x1-x2)break;elseplot([x1(1),x2(1)],[x1(2),x2(2)],'-r*'); 画图k=k+1; 迭代继续x1=x2; 赋值endendend优点：在极⼩点附近收敛快缺点：但是要计算⽬标函数的hesse 矩阵最速下降法1. ：选取初始点xo ，给定误差2. 计算⼀阶梯度。

若⼀阶梯度⼩于误差，停⽌迭代，输出3. 取()()()k k p f x =?4. 10t ()(), 1.min k k k k k k k k k k t f x t p f x tp x x t p k k +≥+=+=+=+进⾏⼀维搜索，求，使得令转第⼆步例题：求min (x-2)^4+(x-2*y)^2.初始值（0,3）误差为0.1(1)编写⼀个⽬标函数，存为f.mfunction z = f( x,y )z=(x-2.0)^4+(x-2.0*y)^2;end(2)分别关于x 和y 求出⼀阶梯度，分别存为fx.m 和fy.m function z = fx( x,y ) z=2.0*x-4.0*y+4.0*(x-2.0)^3;end和function z = fy( x,y )z=8.0*y-4.0*x;end（3）下⾯是脚本⽂件，⼀维搜索⽤的是黄⾦分割法Tic 计算时间eps=10^(-4);误差err=10;dt=0.01;x0=1.0;初始值y0=1.0;mm=0;while err>eps 黄⾦分割法dfx=-fx(x0,y0);dfy=-fy(x0,y0);tl=0;tr=1;确定⼀维搜索的区间h=3;nn=0;gerr=10;geps=10^(-4);while gerr>gepstll=tl+0.382*abs(tr-tl);trr=tl+0.618*abs(tr-tl);iff(x0+tll*h*dfx,y0+tll*h*dfy)>f(x0+trr*h*dfx,y0+trr*h*dfy) tl=tll;elsetr=trr;endgerr=abs(tl-tr); 区间的长度之差tt=0.5*(tl+tr);nn=nn+1;步数增加if nn>200 迭代终⽌条件breakendendx0=x0+tt*h*dfx; 重新迭代y0=y0+tt*h*dfy;err=sqrt(fx(x0,y0)^2+fy(x0,y0)^2);mm=mm+1;步数增加if mm>700 迭代步数超过700，终⽌breakendendres=[x0,y0];输出最后的x，y。

用MATLAB实现共轭梯度法求解实例

用MATLAB实现共轭梯度法求解实例介绍共轭梯度法是一种迭代优化算法，用于求解线性方程组或最小化二次函数的问题。

在本文档中，我们将使用MATLAB来实现共轭梯度法，并通过一个实例来演示它的应用。

共轭梯度法共轭梯度法是一种迭代优化算法，用于求解形如Ax=b的线性方程组。

它的主要思想是利用迭代过程中的残差和共轭梯度的特性，逐步逼近方程的解。

共轭梯度法的优点是收敛速度快，尤其适用于大规模稀疏线性方程组的求解。

共轭梯度法的具体步骤共轭梯度法的具体步骤如下： 1. 初始化解向量x和残差r，令初始残差r0等于b减去Ax的乘积，初始搜索方向p等于r0。

2. 迭代更新解向量： - 计算搜索步长α：α等于r的转置乘以r除以p的转置乘以Ap的乘积。

- 更新解向量和残差：x等于x加上α乘以p，r等于r减去α乘以Ap。

- 计算残差的L2范数：如果残差的L2范数小于预设阈值，停止迭代；否则，继续迭代。

- 计算搜索方向的系数β：β等于r的转置乘以r除以上一次迭代的残差r的转置乘以上一次迭代的残差r的乘积。

- 更新搜索方向：p等于r加上β乘以上一次迭代搜索方向p。

3. 输出解向量x，即为线性方程组的解。

实例在这个实例中，我们将使用共轭梯度法来求解以下线性方程组：A = [4, -1, 0; -1, 4, -1; 0, -1, 4]b = [5; 5; 10]首先，我们将初始化解向量x、残差r和搜索方向p：A = [4, -1, 0; -1, 4, -1; 0, -1, 4];b = [5; 5; 10];x = zeros(size(b)); % 初始化解向量xr = b - A*x; % 初始化残差rp = r; % 初始化搜索方向p然后，我们将进行迭代更新解向量：while norm(r) > 1e-6% 设置迭代终止条件为残差的L2范数小于1e-6Ap = A*p; % 计算Apalpha = (r'*r) / (p'*Ap); % 计算搜索步长alphax = x + alpha*p; % 更新解向量xr_new = r - alpha*Ap; % 计算新的残差rbeta = (r_new'*r_new) / (r'*r); % 计算搜索方向系数betap = r_new + beta*p; % 更新搜索方向pr = r_new; % 更新残差rend最后，输出解向量x的值：x共轭梯度法是一种有效的迭代优化算法，用于求解线性方程组。

matlab共轭梯度法求解方程组

主题：matlab共轭梯度法求解方程组近年来，随着科学技术的不断发展，数学建模和计算机仿真成为科学研究和工程技术领域的重要手段。

在实际应用中，我们常常需要解决线性方程组的求解问题，而共轭梯度法作为一种高效的迭代求解方法，广泛应用于信号处理、图像处理、地球物理勘探和优化问题等领域。

本文将介绍如何利用matlab中的共轭梯度法求解线性方程组的基本原理和实际操作方法。

1. 共轭梯度法的基本原理共轭梯度法是一种迭代法，用于求解对称正定线性方程组Ax=b。

该方法的核心思想是通过一系列的迭代操作，逐步逼近方程组的解，直到满足一定的精度要求。

在每一步迭代中，共轭梯度法利用残差和方向向量的共轭性质，不断寻找最优的步长，从而实现方程组的求解。

2. matlab中共轭梯度法的基本调用方法在matlab中，调用共轭梯度法求解线性方程组非常简单。

需要将方程组的系数矩阵A和右端向量b输入到matlab中，然后利用内置函数conjugateGradient进行求解。

具体的调用方法如下：x = conjugateGradient(A, b, x0, maxIter, tol)其中，A为系数矩阵，b为右端向量，x0为初始解向量，maxIter为最大迭代次数，tol为精度要求。

调用完毕后，matlab将返回方程组的近似解x。

3. 共轭梯度法在实际工程中的应用共轭梯度法作为一种高效的求解方法，在工程技术领域得到了广泛的应用。

以图像处理为例，图像处理中经常需要解决大规模的线性方程组，而共轭梯度法能够高效地求解这类问题，提高了图像处理算法的效率和稳定性。

另外，在地球物理勘探中，共轭梯度法也被广泛应用于三维数据的快速处理和解释。

可以说，共轭梯度法在实际工程中发挥着重要的作用。

4. 共轭梯度法的优缺点分析尽管共轭梯度法具有非常高的效率和稳定性，但是该方法也存在一些缺点。

该方法只适用于对称正定的线性方程组，对于一般的线性方程组并不适用。

共轭梯度法的收敛速度受到方程条件数的影响，对于病态问题，可能收敛速度较慢。

共轭梯度法（Python实现）

共轭梯度法（Python实现）共轭梯度法法（Python实现）使⽤共轭梯度法，分别使⽤Armijo准则和Wolfe准则来求步长求解⽅程f(x1,x2)=(x21−2)4+(x1−2x2)2f(x1,x2)=(x21−2)4+(x1−2x2)2的极⼩值import numpy as np# import tensorflow as tfdef gfun(x): # 梯度# x = tf.Variable(x, dtype=tf.float32)# with tf.GradientTape() as tape:# tape.watch(x)# z = fun(x)# return tape.gradient(z, x).numpy() # 这⾥使⽤TensorFlow来求梯度，直接⼿算梯度返回也⾏return np.array([4 * (x[0] - 2) ** 3 + 2 * (x[0] - 2 * x[1]), -4 * (x[0] - 2 * x[1])]).reshape(len(x), 1)def fun(x): # 函数return (x[0] - 2) ** 4 + (x[0] - 2 * x[1]) ** 2def frcg_armijo(x0):maxk = 5000 # 最⼤迭代次数rho = .6 # Armijo准测参数sigma = .4k = 0epsilon = 1e-4n = len(x0) # 输⼊的维度while k < maxk: # 最⼤迭代次数g = gfun(x0) # 计算梯度itern = k % nif itern == 0: # 迭代n（维度）次后，重新选取负梯度⽅向作为搜索⽅向d = -gelse:beta = (g.T @ g) / (g0.T @ g0) # 计算betad = -g + beta * d0gd = g.T @ dif gd >= 0: # 必要条件，要⼩于0，取负梯度⽅向d = -gif np.linalg.norm(g) < epsilon: # 满⾜精度则结束循环breakm = 0mk = 0while m < 20: # 使⽤Armijo搜索（⾮精确线搜索）if fun(x0 + rho ** m * d) < fun(x0) + sigma * rho ** m * g.T @ d:mk = mbreakm += 1x0 = x0 + rho ** mk * dg0 = gd0 = dk += 1val = fun(x0)return x0, val, kdef frcg_wolfe(x0):maxk = 5000 # 最⼤迭代次数k = 0epsilon = 1e-4n = len(x0) # 输⼊的维度while k < maxk: # 最⼤迭代次数g = gfun(x0) # 计算梯度itern = k % nif itern == 0: # 迭代n（维度）次后，重新选取负梯度⽅向作为搜索⽅向d = -gelse:beta = (g.T @ g) / (g0.T @ g0) # 计算betad = -g + beta * d0gd = g.T @ dif gd >= 0: # 必要条件，要⼩于0，取负梯度⽅向d = -gif np.linalg.norm(g) < epsilon: # 满⾜精度则结束循环breakrho = 0.4sigma = 0.5a = 0b = np.infalpha = 1while True:if not ((fun(x0) - fun(x0 + alpha * d)) >= (-rho * alpha * gfun(x0).T @ d)): b = alphaalpha = (a + alpha) / 2continueif not (gfun(x0 + alpha * d).T @ d >= sigma * gfun(x0).T @ d):a = alphaalpha = np.min([2 * alpha, (alpha + b) / 2])continuebreakx0 = x0 + alpha * dg0 = gd0 = dk += 1x = x0val = fun(x)return x, val, kif __name__ == '__main__':x0 = np.array([[0], [3]])x0, val, k = frcg_armijo(x0) # 使⽤armijo准则print(f'近似最优点：\n{x0}\n迭代次数：{k}\n⽬标函数值：{val.item()}')x0 = np.array([[-1.2], [-1]])x0, val, k = frcg_wolfe(x0) # 使⽤wolfe准则print(f'近似最优点：\n{x0}\n迭代次数：{k}\n⽬标函数值：{val.item()}')。

共轭梯度法编程

P90 第二章 11（2）用大M法求解min w=2x1+x2-x3-x4s.t x1-x2+2x3-x4=22x1+x2-3x3+x4=6x1+x2+x3+x4=7xi≥0, i = 1, 2, 3, 4用matlab求解如下：f=[2,1,-1,-1,200,200,200];a=[1 -1 2 -1 1 0 0;2 1 -3 1 0 1 0;1 1 1 1 0 0 1];b=[2 6 7];lb=zeros(7,1);[x,fval,exitflag,output,lambda]=linprog(f,[],[],a,b,lb)Optimization terminated.运行结果如下：x =3.00000.00001.00003.00000.00000.00000.0000fval =2.0000exitflag =1output =iterations: 7algorithm: 'large-scale: interior point'cgiterations: 0message: 'Optimization terminated.'lambda =ineqlin: [0x1 double]eqlin: [3x1 double]upper: [7x1 double]lower: [7x1 double]从上述运行结果可以得出：最优解为x= ，最小值约为f*=2。

P151 第三章 26用共轭梯度算法求f(x) = (x1-1)^2+5*(x2-x1^2)^2的极小点，取初始点x0= 。

用matlab求解如下：function mg=MG()%%共轭梯度法求解习题三第26题%clc;clear;n=2;x=[2 0]';max_k=100;count_k=1;trace(1,1)=x(1);trace(2,1)=x(2);trace(3,1)=f_fun(x);k=0;g1=f_dfun(x);s=-g1;while count_k<=max_kif k==ng0=f_dfun(x);s=-g0;k=0;elser_min=fminbnd(@(t) f_fun(x+t*s),-100,100); x=x+r_min*s;g0=g1;g1=f_dfun(x);if norm(g1)<10^(-6)break;endm=(norm(g1)^2)/(norm(g0)^2);s=-g1+m*s;count_k=count_k+1;trace(1,count_k)=x(1);trace(2,count_k)=x(2);trace(3,count_k)=f_fun(x);k=k+1;endendcount_kxf=f_fun(x)function g=f_dfun(x)g(1,1)=20*x(1)^3-20*x(1)*x(2)+2*x(1)-2;g(2,1)=10*x(2)-10*x(1)^2;function f=f_fun(x)f=5*x(1)^4-10*x(1)^2*x(2)+x(1)^2-2*x(1)+1+5*x(2)^2;运行结果如下：k =59x0 =1.00001.0000f0 =从上述运行结果可以得出：最优解为x= ，最小值约为f*=0。

用共轭梯度法解最小二乘问题

用共轭梯度法求解最小二乘问题摘要本文先讨论了求解对称正定线性方程组的共轭梯度法.然后对系数矩阵列满秩的线性方程组运用正则化方法将其转化为对称正定线性方程组后再运用实用共轭梯度法进行求解，最后举例并通过Matlab 程序实现其结果.关键词共轭梯度法；正则化方法；最小二乘问题；Krylov 子空间1 引言在实际的科学与工程问题中，常常将问题归结为一个线性方程组的求解问题，而求解线性方程组的数值解法大体上可分为直接法和迭代法两大类.直接法是指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法.因此，直接法又称为精确法.迭代法则是采取逐次逼近的方法，亦即从一个初始向量出发，按照一定的计算格式，构造一个向量的无穷序列,其极限才是方程组的精确解，只经过有限次运算得不到精确解.当线性方程组的系数矩阵为对称正定矩阵时，我们常用共轭梯度法（或简称CG 法）求解，目前有关的方法与理论已经相当成熟，并且已成为求解大型稀疏线性方程组最受欢迎的一类方法.2 最小二乘问题定义1[1] 给定矩阵m n A R ⨯∈，A 列满秩及向量m b R ∈，确定nx R ∈使得()()2222min min n ny Ry Rb Ax r x r y Ay b ∈∈-===-. 该为问题称为最小二乘问题，简记为LS （Least Squares ）问题，其中()r x 称为残向量.最小二乘问题的解x 又可称为线性方程组Ax b =，m n A R ⨯∈的最小二乘解，即x 在残向量()r x b Ax =-的2范数最小的意义下满足线性方程组 Ax b =，m n A R ⨯∈.3 共轭梯度法考虑线性方程组Ax b =的求解问题，其中A 是给定的n 阶对称正定矩阵，b 是给定的n 维向量，x 是待求解的n 维向量.为此,定义二次泛函()2T T x x Ax b x ϕ=-.定理1[1]设A 对称正定,求方程组Ax b =的解，等价于求二次泛函()x ϕ的极小点. 定理1表明，求解线性方程组问题就转化为求二次泛函()x ϕ的极小点问题.求解二次函数极小值问题，通常好像盲人下山那样，先给定一个初始向量0x ，确定一个下山方向0p ，沿着经过点0x 而方向为0p 的直线00x x p α=+找一个点1000x x p α=+，使得对所有实数α有()()00000x p x p ϕαϕα+≤+，即在这条直线上1x 使()x ϕ达到极小.然后从1x 出发，再确定一个下山的方向1p ，沿着直线11x x p α=+再跨出一步，即找到1α使得()x ϕ在2111x x p α=+达到极小：()()11111x p x p ϕαϕα+≤+.重复此步骤，得到一串012,,,ααα 和 012,,,p p p ，称k p 为搜索方向，k α为步长.一般情况下，先在k x 点找下山方向k p ，再在直线k k x x p α=+上确定步长k α使()(),k k k k k x p x p ϕαϕα+≤+最后求出1k k k k x x p α+=+.然而对不同的搜索方向和步长，得到各种不同的算法.由此，先考虑如何确定k α.设从k x 出发，已经选定下山方向k p .令()()k k f x p αϕα=+()()()2TT k k k k k k x p A x p b x p ααα=++-+()22TT k k k k k p Ap r p x ααϕ=-+，其中k k r b Ap =-.由一元函数极值存在的必要条件有()220TT k k k k f p Ap r p αα'=-=所确定的α即为所求步长k α，即T k kk Tk kr p p Ap α=. 步长确定后，即可算出1k k k k x x p α+=+.此时，只要0Tk k r p ≠，就有()()()()1k k k k k k x x x p x ϕϕϕαϕ+-=+-()2220T kk TTk k k k k kT k kr p p Ap r p p Ap αα=-=-<即()()1k k x x ϕϕ+<.再考虑如何确定下山方向k p .易知负梯度方向是()x ϕ减小最快的方向，但简单分析就会发现负梯度方向只是局部最佳的下山方向，而从整体来看并非最佳.故采用新的方法寻求更好的下山方向——共轭梯度法.定义2[2]若n 维非零向量,x y 满足0T x Ay =其中A 为n 阶对称正定矩阵，则称x 与y 是相互共轭（A -共轭）的. 下面给出共轭梯度法的具体计算过程:给定初始向量0x ，第一步仍选用负梯度方向为下山方向，即00p r =，于是有00010001000,,T T r r x x p r b Ax p Ap αα==+=-.对以后各步，例如第k+1步（k ≥1）,下山方向不再取k r ，而是在过点由向量k r 和1k p -所张成的二维平面21{|,,}k k k x x x r p R πξηξη-==++∈内找出使函数ϕ下降最快的方向作为新的下山方向k p .考虑ϕ在2π上的限制：()1,()k k k x r p ψξηϕξη-=++11()()Tk k k k k k x r p A x r p ξηξη--=++++12()Tk k k b x r p ξη--++.计算ψ关于,ξη的偏导得：()()11112,2,T T T k k k k k k T Tk k k k r Ar r Ap r r r Ap p Ap ψξηξψξηη----∂=+-∂∂=+∂其中最后一式用到了10Tk k r p -=，这可由k r 的定义直接验证.令0ψψξη∂∂==∂∂，即知ϕ在2π内有唯一的极小值点001k k k x x r p ξη-=++，其中0ξ和0η满足00101011,0.T T T k k k k k k T Tkk k k r Ar r Ap r r r Ap p Ap ξηξη----⎧+=⎨+=⎩ 由于0k r ≠必有00ξ≠，所以可取()0101k k k k p x x r p ηξξ-=-=+作为新的下山方向.显然，这是在平面2π内可得的最佳下山方向.令010k ηβξ-=，则可得 1111.T k k k T k k r Ap p Ap β----=-注：这样确定的k p 满足10Tk k p Ap -=，即k p 与1k p -是相互共轭的.总结上面的讨论，可得如下的计算公式：T k kk Tk kr p p Ap α= ， 1k k k k x x p α+=+， 11k k r b Ax ++=-， 1T k kk Tk kr Ap p Ap β+=-， 11k k k k p r p β++=+. 在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称的计算公式.首先来简化1k r +的计算公式：11()k k k k k k k k r b Ax b A x p r Ap αα++=-=-+=-.因为k Ap 在计算k α是已经求出，所以计算1k r +时可以不必将1k x +代入方程计算，而是从递推关系1k k k r b Ap α+=-得到.再来简化k α和k β的计算公式.此处需要用到关系式1110,T T T k k k k k k r r r p r p +-+=== 1,2,k =.从而可导出1111,T T k k k kr r r α+++=-，()111T TTk k k k k k k kkp Ap p r r p r αα+=-=()1111T T k k k k k k kkr r p r r βαα--=+=.由此可得,T k k k T k k r r p Ap α=， 11.T k k k T k kr r r r β++=.从而有求解对称正定方程组的共轭梯度法算法如下：0x =初值00r b Ax =-；0k =while 0k r ≠1k k =+if 1k = 00p r =else21122T Tk k k k k r r r r β-----= 1122k k k k p r p β----=+ end11111T Tk k k k k r r p Ap α-----= 111k k k k x x p α---=+ 111k k k k r r Ap α---=-endk x x =注：该算法每迭代一次仅需要使用系数矩阵A 做一次矩阵向量积运算. 定理2 [1]由共轭梯度法得到的向量组{}i r 和{}i p 具有如下基本性质：（1）0Ti j p r =， 0;i j k ≤<≤ （2）0Ti j r r =， i j ≠，0,;i j k ≤≤ （3）0Ti j p Ap =， i j ≠，0,;i j k ≤≤ （4）000{,,}{,,}(,,1)k k span r r span p p A r k κ==+，其中0000(,,1){,,,}k A r k span r Ar A r κ+=，通常称之为Krylov 子空间.下面给出共轭梯度法全局最优性定理：定理3[1]用共轭梯度法计算得到的近似解k x 满足()(){}00min :(,,)k x x x x A r k ϕϕκ=∈+或{}**00min :(,,)k AAx x x x x x A r k κ-=-∈+，其中Ax=，*x 是方程组Ax b =的解，0(,,)A r k κ是由所定义的Krylov 子空间.定理2表明，向量组0,,k r r 和0,,k p p 分别是Krylov 子空间0(,,1)A r k κ+的正交基和共轭正交基.由此可知，共轭梯度法最多n 步便可得到方程组的解*x .因此，理论上来讲，共轭梯度法是直接法.然而实际使用时，由于误差的出现，使k r 之间的正交性很快损失，以致于其有限步终止性已不再成立.此外，在实际应用共轭梯度法时，由于一般n 很大，以至于迭代()O n 次所耗费的计算时间就已经使用户无法接受了.因此，实际上将共轭梯度法作为一种迭代法使用，而且通常是k r 是否已经很小及迭代次数是否已经达到最大允许的迭代次数max k 来终止迭代.从而得到解对称正定线性方程组的实用共轭梯度法，其算法如下：x =初值0;k =;r b Ax =-T r r ρ=while)()max2band k kε><1k k =+if 1k =p r = else;p r p βρρβ==+ end;;TAp p x x p ωαρωα===+ ;;Tr r r r αωρρρ=-== end算法中，系数矩阵A 的作用仅仅是用来由已知向量p 产生向量Ap ω=，这不仅可以充分利用A 的稀疏性，而且对某些提供矩阵A 较为困难而由已知向量p 产生向量Ap ω=又十分方便的应用问题有益.4 共轭梯度法求解最小二乘问题的正则化方程组（法方程组）定理4[1] 当A 列满秩时，求最小二乘问题的解等价于解T TA Ax A b =.应用共轭梯度法于对称正定方程组T T A Ax A b =来求方程组Ax b =，m nA R ⨯∈()m n ≥且A 列满秩的最小二乘解，即为Krylov 子空间法中的正则化方法.由A 列满秩有T A A 对称正定，则方程组T T A Ax A b =，m nA R ⨯∈()m n ≥存在唯一解.下面给出其实用共轭梯度法的详细算法且算法中不出现计算TA A 情形:x =初值0;;;;T T k b A b m Ax r b A m ====-while)()max2band k kε><1k k =+if 1k =p r = else;p r p βρρβ==+end;;;T Tn Ap A n p x x p ωαρωα====+ ;;Tr r r r αωρρρ=-==end注：算法中采用了两次矩阵向量积来避免出现计算T A A 情形.算例编写实用共轭梯度法的Matlab 程序求解方程组TTA Ax A b =，其中11112231A -⎡⎤⎢⎥-⎢⎥=⎢⎥-⎢⎥-⎣⎦, 1234b ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦. 解先建立conjgrad.m 文件，内容如下:function [x]=conjgrad(A,b,x) A=[1 -1;-1 1;2 -2;-3 1]; b=[1;2;3;4]; x=rand(2,1); k=0;b=A'*b;m=A*x; r=b-A'*m; T=r'*r;While norm(r)>1e-10*norm(b) k=k+1; if k==1 p=r;else B=T/t; p=r+B*p;endn=A*p;w=A'*n; a=T/(p'*w); x=x+a*p; r=r-a*w; t=T; T=r'*r; end然后运行后得 ans =-2.4167 -3.2500即有方程组的数值解 2.41673.2500x -⎡⎤=⎢⎥-⎣⎦.而其精确解可由如下方法求得: 111123111591121229731T A A -⎡⎤⎢⎥----⎡⎤⎡⎤⎢⎥==⎢⎥⎢⎥⎢⎥----⎣⎦⎣⎦⎢⎥-⎣⎦，11123271121314T A b ⎡⎤⎢⎥---⎡⎤⎡⎤⎢⎥==⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎢⎥⎣⎦，则有121597971x x --⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦，解得122912134x x ⎡⎤-⎢⎥⎡⎤=⎢⎥⎢⎥⎣⎦⎢⎥-⎢⎥⎣⎦，即方程精确解为*2912134x ⎡⎤-⎢⎥=⎢⎥⎢⎥-⎢⎥⎣⎦，故可验证通过Matlab 程序求得的数值解2.41673.2500x -⎡⎤=⎢⎥-⎣⎦满足精度要求.5 总结本文首先给出最小二乘问题的定义，随后从盲人下山法开始讨论了共轭梯度法的具体推导过程及其相关性质与算法.继而重点给出正则化方法的实用共轭梯度算法并举例进行检验.最后，需要说明虽然正则化方法是求一般线性方程组Ax b =，m nA R⨯∈()m n ≥且A 列满秩的最小二乘解的一种方法且简单易行，但是也有许多不足之处，如m n >时一般无解；TA A 形成时运算量大，A 中某些信息会丢失；当A 病态时其收敛性速度由于222()()T A A A κκ=很大变得非常之慢等，故为了避免正则化方法的缺点，还可运用残量极小化方法或残量正交化方法等更好的方法来解决此类问题.参考文献[1] 徐树方，高立，张平文.数值线性代数[M].北京：北京大学出版社，2000.139--151. [2] 施光燕，董加礼.最优化方法[M].北京：高等教育出版社，1999.47--52.。

用MATLAB计算复数的实部、虚部、模、辐角,共轭复数、简单复数方程根及复数的极限

z = 2.1520 + 0.9505i
二、复数的三角函数运算同实数的三角函数运算
三角函数运算函数为 sin(x),cos(x),tan(x),cot(x),sec(x),csc(x),sinh(x),cosh(x),tanh(x),coth(x),sech( x),csch(x) 反三角函数运算函数为 asin(x),acos(x),atan(x),acot(x),asec(x),acsc(x), 例 5 求复数 3+4i 的三角函数 >> z=3+4i; >> sin(z) ans = 3.8537 -27.0168i
三、用 MATLAB 计算复数方程的根及极限
例 6 求方程 z 8 0 的根. >> solve('z^3+8=0') ans = -2 1+i*3^(1/2) 1-i*3^(1/2) 例 8 对复变函数 f ( z ) z , 取 z 0 1 2i syms z z0 >> f=z^2; >> z0=1+2i z0 = 1.0000 + 2.0000i >> limit(f,z,z0) ans = -3+4*i
例 1：求复数 12 2i 的实部、虚部、模、共轭复数和辐角 >> z=-sqrt(12)-2i; >> x=real(z) x = -3.4641 >> y=imag(z) y = -2 >> abs(z) ans = 4.0000 >> conj(z) ans =
1
-3.4641 + 2.0000i >> angle(z) ans = -2.6180 例 2 imag([5-8j,3+4j]) ans = -8 4

用MATLAB实现共轭梯度法求解实例

求解无约束优化问题的共轭梯度法

matlab梯度算法

基于matlab平台的三种迭代法求解矩阵方程

共轭梯度法matlab最优化问题

数值计算实例MATLAB实现(附带详细源码)

用共轭梯度法解最小二乘问题

双共轭梯度法matlab_概述及解释说明

计算方法——共轭梯度法求解线性方程组

最优化牛顿法最速下降法共轭梯度法matlab代码

用MATLAB实现共轭梯度法求解实例

matlab共轭梯度法求解方程组

共轭梯度法（Python实现）

共轭梯度法编程

最新a03共轭梯度法编程汇总

用共轭梯度法解最小二乘问题

用MATLAB计算复数的实部、虚部、模、辐角,共轭复数、简单复数方程根及复数的极限