前言
博主一头小山猪
目前已开放所有文章:小山猪——经典算法专栏
活动地址:CSDN21天学习挑战赛
G-N法求解非线性最小二乘优化问题
正文
L-M法求解非线性最小二乘优化问题
L-M法 (Levenberg-Marquardt法)原理
当矩阵 ( J k ) T J k \left(J_{k}\right)^{T} J_{k} (Jk)TJk 为病态矩阵时,用G-N算法可能得不到正确的解,甚至当 ( J k ) T J k \left(J_{k}\right)^{T} J_{k} (Jk)TJk 不可逆时, 这时G-N算法就无法计算下去。L-M算法通过采用系数矩阵阻尼的方法改造矩阵 ( J k ) T J k \left(J_{k}\right)^{T} J_{k} (Jk)TJk 的性 态,使算法能够进行下去。
L-M算法中有两个主要的步骤: 一是解方程 [ Q + μ I ] Δ x = − ∇ S ( x ( k ) ) [Q+\mu I] \Delta x=-\nabla S\left(x^{(k)}\right) [Q+μI]Δx=−∇S(x(k)) (其中 μ \mu μ 为阻尼系 数)求出自变量的增量,另一个是阻尼系数 μ \mu μ 的调整算法。
算法步骤
用L-M法求解非线性最小二乘优化问题 min S ( x ) \min S(x) minS(x) 的算法过程如下:
【1】给定初始点 x ( 0 ) x^{(0)} x(0) ,选取参数 β ∈ ( 0 , 1 ) , μ > 1 \beta \in(0,1) , \mu>1 β∈(0,1),μ>1 及精度 ε > 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】令 Q = ( ∇ f ( x ( k ) ) ) T ∇ f ( x ( k ) ) Q=\left(\nabla f\left(x^{(k)}\right)\right)^{T} \nabla f\left(x^{(k)}\right) Q=(∇f(x(k)))T∇f(x(k)) 解方程 [ Q + μ I ] Δ x = − ∇ S ( x ( k ) ) [Q+\mu I] \Delta x=-\nabla S\left(x^{(k)}\right) [Q+μI]Δx=−∇S(x(k)) ;
【6】令 x ( k + 1 ) = x ( k ) + Δ x x^{(k+1)}=x^{(k)}+\Delta x x(k+1)=x(k)+Δx ,计算终止条件是否满足,不满足转【7】;
【7】若 S ( x ( k + 1 ) ) < S ( x ( k ) ) + β ( ∇ S ( x ( k ) ) ) T Δ x S\left(x^{(k+1)}\right)<S\left(x^{(k)}\right)+\beta\left(\nabla S\left(x^{(k)}\right)\right)^{T} \Delta x S(x(k+1))<S(x(k))+β(∇S(x(k)))TΔx ,令 μ = μ / v \mu=\mu / v μ=μ/v ,转【8】,否则令 μ = μ ∗ v \mu=\mu * v μ=μ∗v ,转【5】;
【8】令 k = k + 1 k=k+1 k=k+1 ,转【2】。
代码实现
L-M法函数如下
function [x,minf] = minLM(f,x0,beta,u,v,var,eps)
format long;
if nargin == 6
eps = 1.0e-6;
end
S = transpose(f)*f;
k = length(f);
n = length(x0);
x0 = transpose(x0);
A = jacobian(f,var);
tol = 1;
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;
Q = transpose(Ax)*Ax;
while 1
dx = -(Q+u*eye(size(Q)))\gSx;
x1 = x0 + dx;
for i=1:k
Fx1(i,1) = Funval(f(i),var,x1);
end
Sx1 = Funval(S,var,x1);
tol = norm(dx);
if tol<=eps
break;
end
if Sx1 >= Sx+beta*transpose(gSx)*dx
u = v*u;
continue;
else
u = u/v;
break;
end
end
x0 = x1;
end
x = x0;
minf = Funval(S,var,x);
format short;