【Matlab学习手记】非线性拟合的信頼域方法

• 信頼域原理

    常规的线搜索方法，选定一个初始点，产生一个位移方向（搜索方向，如梯度方向），然后确定位移长度（搜索长度），以此来更新初始点位置。 信頼域方法，给定一个“信頼域半径”作为位移长度，并以当前迭代点为中心，以此“上界”为半径确定一个称之为“信頼域”的闭球区域；然后通过求解“信頼域子问题”（目标函数的二次近似模型）的最优解来确定“候选位移”；若候选位移能使目标函数有充分的下降量，则接受该候选位移，并保持或扩大信頼域半径，继续迭代，否则，说明二次模型与目标函数模型的近似度不够理想，则需要缩小信頼域半径，再通过求解新的信頼域内的子问题确定新的候选位移；如此重复，直到满足迭代终止条件。

    信頼域子问题：目标函数的二次Taylor近似，在邻域内满足最优解：

    在此基础上，产生了一系列改进方法，主要分为以下几类：信頼域子问题需要计算Hesse矩阵，计算量大，考虑将其近似的拟牛顿法（最主要的是BFGS）；子问题求解失败时，即候选位移过大，采用线搜索技术给出一个位移；信頼域子问题的求解方法；信頼域半径的更新技巧。 

•  求解信頼域子问题的Dog-Leg法

•  实例分析

    给一个简单的例子，有四个待拟合参数为[5, 1, 10, 2]；

t = (-10:10)';
Et = 5*sin(t) + 10*cos(2*t);

    运行以下程序，经过23次迭代，程序收敛到[5, 1, 10, 2]；

function Opt_Trust_Region()
clear; clc
t = (-10:10)';
Et = 5*sin(t) + 10*cos(2*t);
x = [1.6, 1.4, 6.2, 1.7]';
maxiter = 100;
iter = 0;
Delta = 1;
g = Jfun(x, t)'*fun(x, t, Et);
while iter < maxiter
    iter = iter + 1;
    if iter == maxiter || norm(fun(x, t, Et), Inf) < 1e-6 || norm(g, Inf) < 1e-6 || Delta <= 1e-6 * (norm(x) + 1e-6)
        disp(iter);
        disp(x);
        break;
    end
    J = Jfun(x, t);
    f = fun(x, t, Et);
    g = J'*f;
    alpha = sumsqr(g)/sumsqr(J*g);
    hsd = -g;
    % Dog Leg Method
    if norm(hgn) <= Delta
        hdl = hgn;
    elseif norm(alpha * hsd) >= Delta
        hdl = Delta / norm(hsd) * hsd;
    else
        a = alpha * hsd;
        b = hgn;
        c = a' * (b - a);
        beta = (-c + sqrt(c^2 + sumsqr(b - a) * (Delta^2 - sumsqr(a)))) / sumsqr(b - a);
        hdl = a + beta*(b - a);
    end
    xnew = x + hdl;
    dF = Fun(x, t, Et) - Fun(xnew, t, Et);
    if norm(hgn) <= Delta
        dL = Fun(x, t, Et);
    elseif norm(alpha * hsd) >= Delta
        dL = Delta * (2 * sumsqr(alpha * g) - Delta) / ( 2*alpha);
    else
        dL = 0.5*alpha*(1 - beta)^2*sumsqr(g) + beta*(2 - beta)*Fun(x, t, Et);
    end
    rou = dF/dL;
    if rou > 0
        x = xnew;
    end
    if rou > 0.75
        Delta = max(Delta, 3*norm(hdl));
    elseif rou < 0.25
        Delta = Delta / 2;
    end
end

function x = Left_Divide(A, y)
[u, s, v] = svd(A);
s(s < 1e-10) = 0;
for i = 1 : length(s)
    if s(i, i) > 0
        s(i, i) = 1 / s(i, i);
    end
end
x = v*s*u'*y;

function y = Fun(x, t, Et)
y = 0.5*sumsqr(fun(x, t, Et));

function y = fun(x, t, Et)
y = x(1)*sin(x(2)*t) + x(3)*cos(x(4)*t) - Et;

function y = Jfun(x, t)
y = [sin(x(2)*t), x(1)*t.*cos(x(2)*t), cos(x(4)*t), -x(3)*t.*sin(x(4)*t)];

补充说明：初始值对结果的影响特别大。