matlab用最速下降法（梯度法）计算Rosenbrock函数，求程序代码

Rosenbrock函数

实穗锋纤现代码：

clc,clear all

format long g

x0=[00]

fun=@func

gfun=@gfunc

[x,val,k]=grad(fun,gfun,x0)   %最速下降法（梯度法）

目标函数

function f=func(x)

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

end

梯度函数

function g=gfunc(x)

g=[400*x(1)*(x(1)^2-x(2))-x(2)+2*(x(1)-1)-200*(x(1)^2-x(2))]

end

运行结果

如有问题，可以私信于我。

用遗传算法ga()求Rosenbrock函数的结猜仿果，与基埋用上述方法的结果相接近。

设地球物理数据和模型参数之间满足以下非线性关系:

d=f(m) (8.1)

其中:f表示非线性算子d、m都是列向量。

建立如下目标函数:

φ(m)=[d-f(m)]2=min (8.2)

目标函数在模型mi处的梯度为

地球物理反演教程

梯度法的模型修改量是目标函数的负梯度:

mi+1=mi+Δmi=mi-λgi(8.4)

其中:λ为步长因子穗孝，是一个数，用来控制修改量的大小g、m都为列向量。

下面推导λ的计樱族稿算公式。

将式(8.2)目标函数φ(m)按泰勒公式展开，并略去高次项得

地球物理反演教程

将式(8.4)中的Δmi=-λgi带入式(8.5)得

地球物理反演教程

设经过修改模型后，目标函数φ(mi+1)为零，有

地球物理反演教程

由上式可推出步长因子λ的计算公式:

地球物理反演教程

给定初始模型mi后，首先计算出梯度gi，然后按式(8.8)计算步长因子，最后按式(8.4)修改模型。如果:

φ(mi+1)＜φ(mi) (8.9)

则说明修正量合适，采用新模型继续迭代。否则减小λ后再计算，一般λ减小一半。

梯度法的计算过程如下:

(1)给定初始模型m0

(2)进行正演计算

(3)判断是否满足精度要求，是则反演结束，否则进行第(4)步

(4)按照式(8.4)修改模型，转第(2)步。

一般反演精度采用实测数据和理论数据的相对均方差来量度。

因为目标函数的梯度就是φ值下降最快的方向，所以梯度法又称为“最速下降法”。下面用一个简单的例子来说明梯度法的原理。设有如下一维目标函数:

φ(x)=f(x) (8.10)

从图8.1可见，x0为目标函数的极小值点。g1为x1处的梯度，g2为x2处的梯度。如果初始模型为x1，模型修改量应该为正值才能使目标函数向最小值前进。从图上可知g1为负值，负梯度为正，满足修改方向。同理如果初始模型为x2，模型修改量应该为负值。从图上可知g2为正值，负梯度为负脊孝值，满足修改量为负值的要求。

图8.1 一维目标函数示意图

从这个例子容易看出即使初始模型远远偏离极小值点，只要按照负梯度方向修改模型参数，总能使目标函数达到极小值点。但是上图的极小值点只有一个，容易达到全局极小，如果目标函数具有多个极小值点，那么初始模型的选择就很关键了，选的不好容易陷入局部极小。此外在极小值点附近梯度法反演收敛的速度将会很慢。因此一般在反演的开始采用梯度法，在反演的后期采用其他收敛速度快的反演方法，如前面所介绍的最小二乘法(或称为高斯-牛顿法)。

图8.2 最小二乘法和梯度法修正方向示意图

最小二乘法和梯度法在极小值点附近的模型修正方向如图8.2所示[10]。这个图形将形象的说明为何梯度法在极小值点附近收敛速度慢。

图8.2是二维的简单情况，目标函数是个椭圆面。在初始模型m0处梯度法的修正方向是最速下降方向，也就是和等值线的切线垂直的方向，可见它的方向偏离椭圆的中心极小值点。而最小二乘法(高斯-牛顿法)是解椭圆函数最优化问题的精确方法[6]，它的修正方向将会指向椭圆的中心极小值点。因此在接近极小值点附近最小二乘法的收敛速度要快于梯度法。

为了克服最速下降法收敛慢的缺点，1964年Fletcher和Reeves提出了无约束极小化的共轭梯度法，它是直接从Hestenes和Stiefel(1952)解线性方程组的共轭梯度法发展而来。共轭梯度法使最速下降方向具有共轭性，提高了算法的有效性和可靠性[6]。

梯度法的关键是计算目标函数的梯度，最终还是会归结为计算观测数据对模型参数的偏导数。在一维反演时可以用有限差分法进行偏导数的计算，在高维反演时可以采用其他快速计算偏导数的方法，如第9章将要介绍的利用互换定理计算二维直流电测深偏导数矩阵。

---