OutputFcn for selfmade function to show iteration steps

Hello there,

I have written my own Levenberg-Marquardt Algorithm as shown below.

Now I want to show every step of the Iteration. Normally I do this with OutputFcn, but this won't work because OutputFun is only for Matlab functions (so I think).

How can I do this otherwise? Is there another option to integrate an OutputFcn? Otherwise I would need to save every x of the iteration, but I don't know how.

(By the way, there is a mistake in the Algorithm, where if eps_mu <=beta 0 x = x-s is wrong. There should be x(k+1) = x(k) which I also don't know how to do, but I am working on this ;) )

Thanks für every help :)

function x=levenberg_marquardt(f,dF,x0,mu0,beta0,beta1,maxit,tol)
    n=length(x0);
    k=0;
    mu=mu0;
    x=x0;
    s=-(dF(x0)'*dF(x0) + mu^2*eye(n))\(dF(x0)'*f(x0));
    while norm(s)>tol && k<maxit
        for k=k+1
             s=-(dF(x)'*dF(x)+mu^2*eye(n))\(dF(x)'*f(x));
             eps_mu=0.5*((f(x)'*f(x)-(f(x+s)'*f(x+s)))/(-s'*dF(x)'*f(x)));
            if eps_mu <= beta0
                x=x-s;
                mu=mu/2;
            else
                x=x+s;
              if eps_mu >= beta1
                  mu=2*mu;
              end
            end
        end
    end
end

[Walter Roberson]

What is it that you want to display at each point of the algorithm, and what form do you want to display it?

[Nico Lange]

Sorry, didn't know this was important.

I create ten random points of a circle and solve this with lsqnonlin. I want to compare the iterations of lsqnonlin with my own Levenberg-Marquardt. Full Code below.

So i want to display every circle of the iteration of Levenberg-Marquardt in the same figure, like i've done with the function outfunlsq at the end of the code.

Sorry for bad english, not my first language.

clear;clc;

phi1=1:((2*pi/10.71)):2*pi;                       
r = -0.5 + (0.5 + 0.5) * rand(1,10);   
xrand = sin(phi1) + 0.25 * r;           
yrand = cos(phi1) + 0.25 * r;           

figure (1)                             
scatter(xrand,yrand,200,'.k','linewidth',2);
axis equal;
grid on;
hold on

f=@(x)((xrand-x(1)).^2+(yrand-x(2)).^2-x(3).^2)';
x0=[0.5,0.5,0.5]';                         
KreisFit = lsqnonlin(f,x0);             
                                       
figure (1)                              
phi2= linspace(0,2*pi,100);    
xFitlsq = KreisFit(3)*cos(phi2) + KreisFit(1);
yFitlsq = KreisFit(3)*sin(phi2) + KreisFit(2); 
plot(xFitlsq,yFitlsq,'k-','linewidth',3)
xlabel ('x-Achse');
ylabel ('y-Achse');
title ('Punktwolke für n=10 und Iterationsschritte durch lsqnonlin');
lgd=legend ('Punktwolke n=10','Ausgleichskreis lsqnonlin','Location','southoutside');
lgd.NumColumns = 3;
hold on

eps=0.1;
maxit = 100;
xlsq = x0;

optionslsq = optimoptions('lsqnonlin','Display','iter',"MaxIterations",maxit, ...
"FunctionTolerance",eps, 'OutputFcn', @outfunlsq);
[xlsq]=lsqnonlin(f, xlsq, [], [], optionslsq);

dF = @(x) [(2*(sqrt((xrand-x(1)).^2+(yrand-x(2)).^2)-x(3)).*(x(1)-xrand))...
           ./(sqrt((xrand-x(1)).^2+(yrand-x(2)).^2));
           (2*(sqrt((xrand-x(1)).^2+(yrand-x(2)).^2)-x(3)).*(x(2)-yrand))...
           ./(sqrt((xrand-x(1)).^2+(yrand-x(2)).^2));
           (2*(x(3)-sqrt((xrand-x(1)).^2+(yrand-x(2)).^2)))]';

mu0=1;
beta0 = 0.3;
beta1 = 0.9;
tol = 1*10^-6;
xlev=levenberg_marquardt(f,dF,x0,mu0,beta0,beta1,maxit,tol);

function x=levenberg_marquardt(f,dF,x0,mu0,beta0,beta1,maxit,tol)
    n=length(x0);
    k=0;
    mu=mu0;
    x=x0;
    s=-(dF(x0)'*dF(x0) + mu^2*eye(n))\(dF(x0)'*f(x0));
    while norm(s)>tol && k<maxit && mu^2>tol
        for k=k+1
             s=-(dF(x)'*dF(x)+mu^2*eye(n))\(dF(x)'*f(x));
             eps_mu=0.5*((f(x)'*f(x)-(f(x+s)'*f(x+s)))/(-s'*dF(x)'*f(x)));
            if eps_mu <= beta0
                x=x-s;
                mu=mu/2;
            else
                x=x+s;
              if eps_mu >= beta1
                  mu=2*mu;
              end
            end
        end
    end
end

function stop = outfunlsq(xlsq,optimvalues,state)
    stop = false;
    switch state
        case 'iter'
            if (optimvalues.iteration) > 0
            fprintf(['Iterationen = ', num2str(optimvalues.iteration)])
            xlsq
            figure (1)
            phi = linspace(0,2*pi,100);
            xlsqplot = xlsq(3) * cos(phi) + xlsq(1);
            ylsqplot = xlsq(3) * sin(phi) + xlsq(2);
            p=plot(xlsqplot, ylsqplot, 'Displayname',['Iteration',...
                num2str(optimvalues.iteration)]);
            if optimvalues.iteration == 1
                p.LineWidth = 3;
            end
            hold on
            legend
            end
    end
end

[Nico Lange]

So, no more comment or answer? Too much code?

[Walter Roberson]

In the time since then, I have had activity on more than 70 Questions. I have been working so much on responding to people the last while that I am very behind in my Christmas preparations..

[Nico Lange]

It was not meant as an insult to you, I am sorry,

It was just a general question, if maybe someone else has an answer.

Answer 1

[Walter Roberson]

You are permitted to call your plotting function manually.

I kept the interface to the plot function compatible with having it called from lsqnonlin . I added an extra optimvalues field whose presense or absence is detected for backwards compatibility. The field indicates which line to update, so that you can have multiple lines on the same update, as you specifically requested that the LM updates go on the same figure as the other updates (presumably for comparison.)

This code assumes that every iteration you are wanting to update the existing circle rather that draw a new one -- since each iteration you would be getting a better-refined position.

If you wanted to add to the plot, then probably the easiest way would be to use animatedline() and add in new coordinates with an extra nan on the end of them to cause a line break from the next set of circles. If you are wanting to add to the plot, then I worry that the legend is going to overflow.

function x=levenberg_marquardt(f,dF,x0,mu0,beta0,beta1,maxit,tol)
    n=length(x0);
    k=0;
    mu=mu0;
    x=x0;
    s=-(dF(x0)'*dF(x0) + mu^2*eye(n))\(dF(x0)'*f(x0));
    while norm(s)>tol && k<maxit && mu^2>tol
        for k=k+1
             s=-(dF(x)'*dF(x)+mu^2*eye(n))\(dF(x)'*f(x));
             eps_mu=0.5*((f(x)'*f(x)-(f(x+s)'*f(x+s)))/(-s'*dF(x)'*f(x)));
            if eps_mu <= beta0
                x=x-s;
                mu=mu/2;
            else
                x=x+s;
              if eps_mu >= beta1
                  mu=2*mu;
              end
            end
            outfunlsq(x, struct('iteration', k, 'plotnum', 2), 'iter');
        end
    end
end

function stop = outfunlsq(xlsq,optimvalues,state)
    stop = false;
    
    persistent p
    
    switch state
        case 'iter'
            if (optimvalues.iteration) > 0
                fprintf(['Iterationen = ', num2str(optimvalues.iteration)])
                xlsq
                phi = linspace(0,2*pi,100);
                xlsqplot = xlsq(3) * cos(phi) + xlsq(1);
                ylsqplot = xlsq(3) * sin(phi) + xlsq(2);
                if isfield(optimvalues, 'plotnum')
                    plotnum = optimvalues.plotnum;
                else
                    plotnum = 1;
                end
                if length(p) < plotnum || ~isvalid(p(plotnum))
                    fig = figure(1);
                    ax = gca(fig);
                    specs = {'r', 'b'};
                    p(plotnum) = line(ax, nan, nan, specs{plotnum}, 'LineWidth', 3);
                    hold(ax, 'on')
                    legend(ax, 'show')
                end
                set(p(plotnum), 'XData', xlsqplot, 'YData', ylsqplot, ...
                    'DisplayName', "plot#" + plotnum + " Iteration " + optimvalues.iteration);
                drawnow()
            end
    end
end

[Nico Lange]

Thank you :)