MATLAB Answers

Infinite Recursion in own Levenberg-Marquardt Code

2 views (last 30 days)
Nico Lange
Nico Lange on 3 Dec 2020
Commented: Nico Lange on 3 Dec 2020
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

Answers (2)

Ameer Hamza
Ameer Hamza on 3 Dec 2020
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);
  2 Comments
Nico Lange
Nico Lange on 3 Dec 2020
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

Sign in to comment.


Steve Eddins
Steve Eddins on 3 Dec 2020
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by