前言
博主一头小山猪
目前已开放所有文章:小山猪——经典算法专栏
活动地址:CSDN21天学习挑战赛
G-N法求解非线性最小二乘优化问题
正文
修正G-N法求解非线性最小二乘优化问题
修正G-N法求解非线性最小二乘优化问题
为了克服 G − N G-N G−N 的缺点,有以下两种修正方案:
(1) 在已经得到一个近似极小点 x k x^{k} xk 后,计算 x k + 1 = x k − α k [ ( J k ) T J k ] − 1 ( J k ) T f ( x k ) x^{k+1}=x^{k}-\alpha_{k}\left[\left(J_{k}\right)^{T} J_{k}\right]^{-1}\left(J_{k}\right)^{T} f\left(x^{k}\right) xk+1=xk−αk[(Jk)TJk]−1(Jk)Tf(xk) ,其中 α k \alpha_{k} αk 的大小用一维无约束优化方法求解 min S ( x k + 1 ) \min S\left(x^{k+1}\right) minS(xk+1) 确定,由于每一步都要进行一维搜 索,计算量很大;
(2) 由于 v k = − [ ( J k ) T J k ] − 1 ( J k ) T f ( x k ) v^{k}=-\left[\left(J_{k}\right)^{T} J_{k}\right]^{-1}\left(J_{k}\right)^{T} f\left(x^{k}\right) vk=−[(Jk)TJk]−1(Jk)Tf(xk) 为目标函数的下降方向,选取很小的正数 δ k \delta_{k} δk , 使得 S ( x k + δ k v k ) < S ( x k ) S\left(x^{k}+\delta_{k} v^{k}\right)<S\left(x^{k}\right) S(xk+δkvk)<S(xk) ,此种方案计算量比较小,可以用比较简单的方法来确定 δ k \delta_{k} δk 的 取值。
因此,根据第二种方案可编制下面的修正G-N算法。
算法步骤
用修正 G − N \mathrm{G}-\mathrm{N} G−N 法求解非线性最小二乘优化问题 min S ( x ) \min S(x) minS(x) 的算法过程如下:
【1】给定初始点 x ( 0 ) x^{(0)} x(0) ,及精度 ε > 0 \varepsilon>0 ε>0 ,置 k = 0 k=0 k=0 ;
【2】计算 f ( x ( k ) ) , S ( x ( k ) ) f\left(x^{(k)}\right), S\left(x^{(k)}\right) f(x(k)),S(x(k)) ;
【3】计算 ∇ f ( x ( k ) ) \nabla f\left(x^{(k)}\right) ∇f(x(k));
【4】 计算 ∇ S ( x ( k ) ) = ( ∇ f ( x ( k ) ) ) T ∗ f ( x ( k ) ) \nabla S\left(x^{(k)}\right)=\left(\nabla f\left(x^{(k)}\right)\right)^{T} * f\left(x^{(k)}\right) ∇S(x(k))=(∇f(x(k)))T∗f(x(k)) ;
【5】解方程 [ ( ∇ f ( x ( k ) ) ) T ∇ f ( x ( k ) ) ] Δ x = − ∇ S ( x ( k ) ) \left[\left(\nabla f\left(x^{(k)}\right)\right)^{T} \nabla f\left(x^{(k)}\right)\right] \Delta x=-\nabla S\left(x^{(k)}\right) [(∇f(x(k)))T∇f(x(k))]Δx=−∇S(x(k)) ;
【6】令 α = 1 , β = 1 0 − 5 \alpha=1, \beta=10^{-5} α=1,β=10−5 ;
【7】 S ( x ( k ) + α Δ x ) ⩽ S ( x ( k ) ) + 2 β α ( Δ x ) T ∇ S ( x ( k ) ) S\left(x^{(k)}+\alpha \Delta x\right) \leqslant S\left(x^{(k)}\right)+2 \beta \alpha(\Delta x)^{T} \nabla S\left(x^{(k)}\right) S(x(k)+αΔx)⩽S(x(k))+2βα(Δx)T∇S(x(k)) ,令 x ( k + 1 ) = x ( k ) + α Δ x x^{(k+1)}=x^{(k)}+\alpha \Delta x x(k+1)=x(k)+αΔx 转【8】; 否则令 α = α / 2 \alpha=\alpha / 2 α=α/2 ,转【7】;
【8】终止条件是否满足,不满足则令 k = k + 1 k=k+1 k=k+1 ,转【2】。
代码实现
修正G-N法函数如下
function [x,minf] = minMGN(f,x0,var,eps)
format long;
if nargin == 3
eps = 1.0e-6;
end
S = transpose(f)*f;
k = length(f);
n = length(x0);
x0 = transpose(x0);
tol = 1;
A = jacobian(f,var);
while tol>eps
Fx = zeros(k,1);
for i=1:k
Fx(i,1) = Funval(f(i),var,x0);
end
Sx = Funval(S,var,x0);
Ax = Funval(A,var,x0);
gSx = transpose(Ax)*Fx;
dx = -transpose(Ax)*Ax\gSx;
alpha = 1;
while 1
S1 = Funval(S,var,x0+alpha*dx);
S2 = Sx+2*(1.0e-5)*alpha*transpose(dx)*gSx;
if S1>S2
alpha = alpha/2;
continue;
else
break;
end
end
x0 = x0 + alpha*dx;
tol = norm(dx);
end
x = x0;
minf = Funval(S,var,x);
format short;