有没有更快的方法来解决我的牛顿算法？

Question

我有这个牛顿算法，它似乎要花很长时间才能找到解决方案。以至于我让它坐了一整夜，什么也没发生。我想知道这是我做错了什么还是我遗漏了什么。

Gobs=[0.61 1.14 2.33 4.76 6.65 4.77 2.38 1.13 0.59];
x=[-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.5];
syms m h
p1=[m;h];

F = ((((6.67*p1(1)*p1(2))/((x(1)^2)+p1(2)^2)^(3/2))-Gobs(1))^2+(((6.67*p1(1)*p1(2))/((x(2)^2)+p1(2)^2)^(3/2))-Gobs(2))^2 + (((6.67*p1(1)*p1(2))/((x(3)^2)+p1(2)^2)^(3/2))-Gobs(3))^2 + (((6.67*p1(1)*p1(2))/((x(4)^2)+p1(2)^2)^(3/2))-Gobs(4))^2 + (((6.67*p1(1)*p1(2))/((x(5)^2)+p1(2)^2)^(3/2))-Gobs(5))^2 + (((6.67*p1(1)*p1(2))/((x(6)^2)+p1(2)^2)^(3/2))-Gobs(6))^2 + (((6.67*p1(1)*p1(2))/((x(7)^2)+p1(2)^2)^(3/2))-Gobs(7))^2 + (((6.67*p1(1)*p1(2))/((x(8)^2)+p1(2)^2)^(3/2))-Gobs(8))^2 + (((6.67*p1(1)*p1(2))/((x(9)^2)+p1(2)^2)^(3/2))-Gobs(9))^2);

dfm = diff(F,'m');
dfh = diff(F,'h');
d2fm = diff(dfm,'m'); 
d2fh = diff(dfh,'h');
dfmh = diff(dfm,'h');
dfhm = diff(dfh,'m');
G = [dfm;dfh];
Hnewton = [d2fm dfmh; dfmh d2fh];    

m =1.6;
h =1.7;

p1 = subs(p1);
n=1;
tm(n)=p1(1);
th(n)=p1(2);
for i=1:2
F1 = subs(F);
p2 = p1-(subs(Hnewton)+ (100/n)*eye(2))\subs(G);
m=p2(1);
p=p2(2);
F2 = subs(F);
if (F2)>0.9;
p1 = p2;
m=p2(1);
h=p2(2);
else break
end
n=n+1;
tm(n)=p2(1);
th(n)=p2(2);
end

plot(th,tm,'r-','LineWidth',1.4)

Answer 1

定义要优化的函数的两个版本：

• F_eval 函数句柄类型用于评估F，

• F syms 函数用于梯度和hessian矩阵

要评估梯度G 和hessian 矩阵Hnewton，使用它们对应的函数句柄

% Gradient
G_eval = matlabFunction(G,'Vars',{[m h]});

% Hessian 
Hnewton_eval = matlabFunction(Hnewtonn,'Vars',{[m h]}); 

还有循环停止条件错误：

使用两个连续函数评估之间的差异

阅读下面的代码

clc
clear
syms m h
p1 = [m h];

Gobs=[0.61 1.14 2.33 4.76 6.65 4.77 2.38 1.13 0.59];
x=[-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.5];
% Define the input known array x first 
% Used only for gradient and hessian, syms function  
F =(...
    (((6.67*p1(1)*p1(2))/((x(1)^2)+p1(2)^2)^(3/2))-Gobs(1))^2 ...
    +(((6.67*p1(1)*p1(2))/((x(2)^2)+p1(2)^2)^(3/2))-Gobs(2))^2 ...
    +(((6.67*p1(1)*p1(2))/((x(3)^2)+p1(2)^2)^(3/2))-Gobs(3))^2 ...
    +(((6.67*p1(1)*p1(2))/((x(4)^2)+p1(2)^2)^(3/2))-Gobs(4))^2 ... 
    +(((6.67*p1(1)*p1(2))/((x(5)^2)+p1(2)^2)^(3/2))-Gobs(5))^2 ... 
    +(((6.67*p1(1)*p1(2))/((x(6)^2)+p1(2)^2)^(3/2))-Gobs(6))^2 ...
    +(((6.67*p1(1)*p1(2))/((x(7)^2)+p1(2)^2)^(3/2))-Gobs(7))^2 ...
    +(((6.67*p1(1)*p1(2))/((x(8)^2)+p1(2)^2)^(3/2))-Gobs(8))^2 ... 
    +(((6.67*p1(1)*p1(2))/((x(9)^2)+p1(2)^2)^(3/2))-Gobs(9))^2 ...
   );
% Used for evaluating the function , function handle 
F_eval = @(p)...
    (...
    (((6.67*p(1)*p(2))/((x(1)^2)+p(2)^2)^(3/2))-Gobs(1))^2 ...
    +(((6.67*p(1)*p(2))/((x(2)^2)+p(2)^2)^(3/2))-Gobs(2))^2 ...
    +(((6.67*p(1)*p(2))/((x(3)^2)+p(2)^2)^(3/2))-Gobs(3))^2 ...
    +(((6.67*p(1)*p(2))/((x(4)^2)+p(2)^2)^(3/2))-Gobs(4))^2 ... 
    +(((6.67*p(1)*p(2))/((x(5)^2)+p(2)^2)^(3/2))-Gobs(5))^2 ... 
    +(((6.67*p(1)*p(2))/((x(6)^2)+p(2)^2)^(3/2))-Gobs(6))^2 ...
    +(((6.67*p(1)*p(2))/((x(7)^2)+p(2)^2)^(3/2))-Gobs(7))^2 ...
    +(((6.67*p(1)*p(2))/((x(8)^2)+p(2)^2)^(3/2))-Gobs(8))^2 ... 
    +(((6.67*p(1)*p(2))/((x(9)^2)+p(2)^2)^(3/2))-Gobs(9))^2 ...
   );


% Gradient and Hessian in syms 
G = gradient(F);
Hnewton = hessian(F);

%% Tranform into function handle


G_eval = matlabFunction(G,'Vars',{[m h]});

Hnewton_eval = matlabFunction(Hnewton,'Vars',{[m h]});

%% Starting guess
m = 1.6;
h =1.7;
n = 1;
pold = [m h];

%% Record 
% First column is m, second is h
solution_history = pold;

%% Infinite loop
while true

     % Use function handle only for evaluation
     F1 = F_eval(pold);

     psol = pold -(Hnewton_eval(pold)+ (100/n)*eye(2))\G_eval(pold);

     [l,~] = size(psol);

     F2_best = F_eval(psol(1, :));
     p_best = psol(1, :);

     % Find the solution which minimize the most F
     for i = 1:l
         if  F_eval(psol(i, :)) < F2_best
             F2_best = F_eval(psol(i, :));
             p_best = psol(1, :);
         end
     end

    F2 = F2_best;
    pnew = p_best;
     solution_history = [solution_history; pnew];

     n = n + 1;

     pold = pnew;
     % Insert stopping condition here  
      if n == 5
          break;
      end

end

tm = solution_history(:, 1);
th = solution_history(:, 2);

plot(th,tm,'r-','LineWidth',1.4)