Is there a faster way to solve my newton algorithm?

Question

I have this Newton algorithm which seems to be taking forever to find a solution. To the point where i let it sit for a whole night and nothing happend. I was wondering if it is something I put wrong or maybe something I am missing.

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)

It's very difficult to follow your code, but keeping everything as symbolic variables will always make code very slow. Once you have found all the functions you need using `diff`, convert them to fast anonymous functions using `matlabFunction`. — David, Dec 12 '19 at 23:44
Oh ok, I will try not to use them so much. I tried converting and it says that I can not use nonscalar arrays, which i tried to convert into cells, but that also did not work. — Miro Feliciano, Dec 13 '19 at 00:04
I tried converting them to function handle, but then I am having problems to adapt the rest of my code. — Miro Feliciano, Dec 13 '19 at 09:55
@ChelseaG.There is no way in matlab to have the derivative of a function handle.`sym` variable is used only for the purpose of calculating the `gradient` and `hessian` of the given function. In function evaluation those two are converted into function handle and making sure their inputs are vector — Adam, Dec 13 '19 at 12:57
@MiroFeliciano use the built-in function `gradient()` and `hessian`. For evaluation they are all converted into function handle, don't use `subs`. Provide the function itself and the known value for the optimization, I will post an appropriate answer. — Adam, Dec 13 '19 at 13:04

Adam · Accepted Answer · 2019-12-14T17:29:36.433

Define two versions of the function to be optimized:

F_eval function handle type used for evaluating F,
F syms function used for the gradient and the hessian matrix

To evaluate the gradient G and the hessian matrix Hnewton, used their corresponding function handle

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

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

Also the loop stopping condition is wrong:

use the difference between two successive function evaluation

Read through the code below


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)

Oh, I did not think about using two versions of the function, thank you! Although I am getting an error in the `p2 = p1 -(Hnewton_eval(p1)+ (100/n)*eye(2))\G_eval(p1);` line, saying that I do not have enough arguments, and I am not understanding why. — Miro Feliciano, Dec 14 '19 at 13:30
Also I edited the main question to include `Gobs` and `x`, which was also missing. P.S. my MatLab version is the R2015a. — Miro Feliciano, Dec 14 '19 at 13:31
@MiroFeliciano answer updated, just missed the curly bracket while using `matlabFunction()`. Your stopping condition is weird. — Adam, Dec 14 '19 at 14:09
My stopping criterion is based on contour levels, so that it stops when it reaches below a certain level. I had to transpose some of the variables you posted because the matrix did not have the same size. My question is though how can I do this iteratively? Because with your while loop i guess it just does it once, and I wanted more times — Miro Feliciano, Dec 14 '19 at 16:57
@MiroFeliciano just updated. The code will iterate 5 times, but you can change the **stopping condition** from `n == 5` to what you want. And in each iteration there might have multiple solution, choose the one with minimum value of `F` — Adam, Dec 14 '19 at 17:33
Thank you very much, it works fine! My only problem with this is that it does not what I want, but perhaps that is my bad formulation — Miro Feliciano, Dec 14 '19 at 19:02
You can reformulate the problem in another post I’ll check it for you. — Adam, Dec 14 '19 at 19:53

Is there a faster way to solve my newton algorithm?

1 Answers1