MATLAB--高斯牛顿法

利用高斯牛顿法求解非线性问题，包含了利用LDLT分解法求解Ax=b的方程:

clear all;
clc
%%演示牛顿高斯法--非线性优化

%真实值
ar = 1.0;
br = 2.0;
cr = 1.0;
x = 1:2:100;
x = x./100;
nosie = 0.9*randn(1,size(x,2));%产生高斯噪音
y_real = exp(ar*x.^2+br*x+cr)+nosie;

%估计初始值
ae = 2.0;
be = 4.0;
ce = 3.0;
itr = 20;%迭代次数
num_itr = 0;%初始迭代次数
esp = 0.01;
%构造函数f(a,b,c)
syms x_ a b c;
arg = [a;b;c];
f(a,b,c,x_) = exp(a*x_.^2+b*x_+c);
Jacobi = jacobian(f,arg);
%高斯牛顿法
last_cost = 0;
for i=1:itr
    H = 0;
    b = 0;
    cost = 0;
    J = [0 0 0];
    for j = 1:size(x,2)
        error = y_real(j)-f(ae,be,ce,x(j));
        
        %雅可比矩阵
         J(1) = -x(j)^2*exp(ae*x(j)^2 + be*x(j) + ce);
         J(2) =  -x(j)*exp(ae*x(j)^2 + be*x(j) + ce);
         J(3) = -exp(ae*x(j)^2 + be*x(j) + ce);
        
         %用MATLAB库函数求出的Jacobian
         %J_ = -Jacobi(ae,be,ce,x(j));%因为y_real(j)-f(ae,be,ce,x(j))，所以要添个负号
         %for k=1:3
         %   J(k) = J_(k);
         %end
       
        H = H + J'*J;
        b = b+ -J'*error;
        cost = cost + error^2;
    end

    %结束条件
    if abs((cost-last_cost)/cost)<esp
        break
    end
    num_itr = i;
    last_cost = cost;
    fprintf('%d,%f\n',i,cost);
    
    %ldlt分解求解dx
    [L,D] = ldl(H);
    bH = sum(b,2);
    dx = L'\(D\(L\bH));
    
    %矩阵直接求逆求解dx
    %dx = H^-1*b;

    ae = vpa(ae+dx(1),2);%放低精度，加快运算速度
    be = vpa(be+dx(2),2);
    ce = vpa(ce+dx(3),2);
end
plot(x,f(ae,be,ce,x));
hold on;
grid on;
scatter(x,y_real,100,'black','.');
fprintf('迭代次数为：%d\n',num_itr);

运行结果：