Infinite Recursion in own Levenberg-Marquardt Code

Hello there,

I am writing an Levenberg-Marquardt Algorithm to calculate a circle out of ten points.

x1 and y1 are 10 points of a circle randomly generated. I am this far:

f=@(x)((x1-x(1)).^2+(y1-x(2)).^2-x(3).^2);   %Function of the circle with x, y and r as x(1), x(2), x(3)
DF = @(x) [2.*((sqrt((x1-x(1)).^2 + (y1-x(2)).^2) - x(3)).*(x(1)-x1))./(sqrt((x1-x(1)).^2 + (y1-x(2)).^2));
           2.*((sqrt((x1-x(1)).^2 + (y1-x(2)).^2) - x(3)).*(x(2)-y1))./(sqrt((x1-x(1)).^2 + (y1-x(2)).^2));
           2.*(x(3)-sqrt((x1-x(1)).^2 + (y1-x(2)).^2))];    % DF = Derivates of x, y, r of f
x0=[0,0,1]

mu0=1;      %some tolerances
beta0 = 0.3;
beta1 = 0.9;
maxit = 100;
tol = 1*10^-5;
    %algorithm
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
    [s,mu] = korrektur(f,DF,x,mu,beta0,beta1);
    x=x+s;
    k=k+1;
end

function[s,mu] = korrektur(f,DF,x,mu,beta0,beta1)
    n=length(x);
    s=-(DF(x)*DF(x)'+mu^2*eye(n))\(DF(x)*f(x)');     %%%%%%%%% Error appears because of this line %%%%%%%%%%%
    eps_mu=(f(x)*f(x)'-f(x+s)*f(x+s)')\(f(x)*f(x)'-(f(x)'+DF(x)'*s)*(f(x)'+DF(x)'*s)');
    if eps_mu <= beta0
     [s,mu]=korrektur(f,DF,x,mu,beta0,beta1);
    elseif eps_mu >= beta1
     mu=mu/2;
    end
end

Everytime I run this, i get this Error because of the marked line above.

Out of memory. The likely cause is an infinite recursion within
the program.

Can someone explain me how to fix this?

Thanks for the help

Answer 1

[Ameer Hamza]

Inside the function korrektur, you are again calling korrektur with the same value of input parameter

[s,mu]=korrektur(f,DF,x,mu,beta0,beta1);

So the recursion never breaks. Maybe, the correct is to pass eps_mu instead of mu

[s,mu]=korrektur(f,DF,x,eps_mu,beta0,beta1);

[Nico Lange]

Thanks again :)

[Nico Lange]

Could you help me again? I solved the problem before by changing mu to 2*mu, but I won't get the result i want.

The algorithm works in that way, that i get a circle in the end. r should be around 1, but with my Levenberg-Marquardt r ist about 1*10^-20

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,80,'ob','linewidth',1);    
axis ([-1.5 1.5 -1.5 1.5]);
grid on;
hold on

f=@(x)((xrand-x(1)).^2+(yrand-x(2)).^2-x(3).^2)';
x0=[0,0,1];                         
KreisFit = lsqnonlin(f,x0);        
                                     
figure (1)                             
phi = linspace(0,2*pi,100);             
xFit = KreisFit(3)*cos(phi) + KreisFit(1); 
yFit = KreisFit(3)*sin(phi) + KreisFit(2);
plot(xFit,yFit,'b-','linewidth',3)
hold off
xlabel ('x-Achse');
ylabel ('y-Achse');
title ('Punktwolke für n=10');
legend ('n=10','Ausgleichskreis','Location','northeast');

options = optimoptions(@lsqnonlin,'Display','iter-detailed','MaxfunEvals',35);
[x1abw,resnorm,residual,exitflag,output] = lsqnonlin(f,x0,[],[],options);

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;
maxit = 100;
tol = 1*10^-6;

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
    [s,mu] = korrektur(f,DF,x,mu,beta0,beta1);
    x=x+s;
    k=k+1;
end

figure (2)                              
phi = linspace(0,2*pi,100);             
xFit1neu = x(3)*cos(phi) + x(1);    
yFit1neu = x(3)*sin(phi) + x(2);    
plot(xFit1neu,yFit1neu,'b-','linewidth',3)
axis ([-1.5 1.5 -1.5 1.5]);
grid on;
xlabel ('x-Achse');
ylabel ('y-Achse');
title ('Punktwolke für n=10');
legend ('Ausgleichskreis 1','Location','northeast');

function[s,mu] = korrektur(f,DF,x,mu,beta0,beta1)
    n=length(x);
    s=-(DF(x)'*DF(x)+mu^2*eye(n))\(DF(x)'*f(x));
    eps_mu=(f(x)'*f(x)-f(x+s)'*f(x+s))/(f(x)'*f(x)-(f(x)+DF(x)*s)'*(f(x)+DF(x)*s));
    if eps_mu <= beta0
         [s,mu]=korrektur(f,DF,x,2*mu,beta0,beta1); %Changed mu here to 2*mu
    elseif eps_mu >= beta1
     mu=mu/2;
    end
end

My 'mu' in the end is about 1*10^10 and s is about 1*10^-20 and i can't explain why. The result should be similar to result with lsqnonlin. I get a circle in the end, but a very very small circle :D

Answer 2

[Steve Eddins]

I see that korrektur is calling itself recursively. I noticed also that the recursive call appears to be identical to the top-level call:

[s,mu]=korrektur(f,DF,x,mu,beta0,beta1);

It looks like none of the quantities f, DF, x, mu, beta0, and beta1 are altered before the recursive call. That implies that eps_mu is always going to be same, which would explain why the elseif branch never gets executed to stop the recursion.

I think you need to examine this recursion more carefully to make sure that it is doing what you intend and that it will terminate successfully.

[Nico Lange]

Thank you :)