# larger computation time for every iteration of levenberg marquardt algorithm?

3 views (last 30 days)
asim asrar on 18 Oct 2022
Answered: Alvaro on 25 Jan 2023
the code written below is taking lot of time for each iteration , can someone suggest some way to solve the issue
tic()
clc;
clear all;
close all;
N=100;
[uin]=incident(N);
[G]=green_f(N);
P=50;
f=sc_potc(N,P);
H=G;
uin = uin(:);
f=f(:);
I = eye(N^2,N^2);
A=I-G*diag(f);
maxit=50;
x0 = rand(N,N);
x0=x0(:);
[u_tilda,FLAG,RELRES]=cgs(A,uin,[],maxit,[],[],x0);
y1=H*diag(u_tilda)*f;
y=abs(y1);
f2 =zeros(N,N);
% f2=ones(N,N) + 1i*ones(N,N);
f2=f2(:);
stepsize=0.5;
fo=f2;
po=ones(N,N);
po=po(:);
i=1;
fcontinue=1;
tol = 1e-4;
while fcontinue
iter=i
ss(i)=stepsize;
maxit=20;
y0=rand(N,N);
y0=y0(:);
AA=I-G*diag(fo);
[u]=cgs(AA,uin,[],maxit,[],[],y0);
X11=H*diag(u);
obj_old = 0.5*((norm(X11*fo - (diag(y)*po),2))^2);
ff1(i)=obj_old;
Jf=(I+diag(fo)*(inv(I-G*diag(fo)))*G)*diag(u);
Hesf=(Jf'*Jf);
r=H*diag(u)*fo -diag(y)*po;
df=(Jf'*H'*(r));
fn=fo-(inv(Hesf + stepsize*I))*df;
fn = reshape(fn,[N,N]);
param.tol = 10e-4;
param.verbose = 1;
param.maxit=200;
param.weights=[1,1];
gamma=0.004;
fre=real(fn);
fim=imag(fn);
fre=max(fre,0);
fim=max(fim,0);
fn=fre +1i*fim ;
fn=fn(:);
Jp=(diag(y'));
Hesp=(Jp'*Jp);
dp=(-1*(diag(y'))*r);
pp(i)=norm(dp,2);
pn =po -(inv(Hesp + stepsize*I))*dp;
for j = 1:N^2;
pn(j)=pn(j) / (norm(pn(j)));
end
AA_new=I-G*diag(fn);
z0=rand(N,N);
z0 = z0(:);
[u_new]=cgs(AA_new,uin,[],maxit,[],[],z0);
X11_new = H*diag(u_new);
obj_new =0.5*((norm(X11_new*fn - diag(y)*pn,2))^2);
rel_obj= abs(obj_new - obj_old)/obj_new;
if rel_obj < tol;
fcontinue =0;
break;
end
if obj_new < obj_old
stepsize = stepsize/2;
fo=fn;
po=pn;
else
stepsize = 2*stepsize;
end
ii1(i)=i;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(5)
subplot 121
imagesc((reshape(real(fn),[N,N])))
colormap gray
colorbar
title('recontructed amplitude real')
subplot 122
imagesc((reshape(imag(fn),[N,N])))
colormap gray
colorbar
title('recontructed amplitude imaginary')
pause(0.5)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i=i+1;
end
%%%%%%%%%%%%%%%%%%%%% function sc_pot and green_f and incident are given as
%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%
.
% %save('G.mat','H')
% imagesc(abs(H))
function G=green_f(N,L)
%N=sqrt(N);
%N=100;
um = 1;
%um = 1e-6;
c = 3e8; % m/s
dx = um;
dy = um;
dz = um;
%
% dx = um;
% dy = um;
% dz = um;
dt = 1/4 * dx /c;
lamb =5.32*um;
k = 2*pi / lamb;
% omega region configuration
x = linspace(0,N-1,N);
y=x;
x=x*dx;
y=y*dy;
[X,Y] = meshgrid(x,y);
X=(X(:));
Y=(Y(:));
g=zeros(length(X),length(Y));
%%
for ii=1:length(X)
a=X-X(ii);
g(ii,:)=a;
end
gg=zeros(length(X),length(Y));
for ii=1:length(X)
aa=Y-Y(ii);
gg(ii,:)=aa;
end
r = sqrt(g.^2 + gg.^2 +5^2*dz);
nb=1.33; % refractive index background
siz=[N,N];
%siz=[N^2,N^2]; % size of the region of interest (containing the support of f)
dz=16*lamb/siz(1); % axial discretization step (o have a ROI of 16*lamb)
kdz=2*pi*nb/lamb; % wavenumber
%kg=9.9260
kg=kdz;
p=parpool(20);
nn=size(r,2);
parfor i=1:size(r,1)
%for i=1:size(r,1)
for j=1:nn
G(i,j)= (exp(1i*(kg*r(i,j)*nb)))/(4*pi*r(i,j));
end
end
delete(p);
end
%%%%%%%%%%%%%% incident function %%%%%%%%%%%%%%%%%%%%%
function E1=incident(N)
%%
%N=100;
um = 1;
%um = 1e-6;
c = 3e8; % m/s
dx = um;
dy = um;
dz = um;
x1 = linspace(0,N-1,N);
y1=x1;
z1=x1;
x1=x1*dx;
y1=y1*dy;
z1=z1*dz;
[X1,Y1,Z1] = meshgrid(x1,y1,z1);
lamb=5.32*um; % wavelength
nb=1.33; % refractive index background
kdz=2*pi*nb/lamb; % wavenumber
kg=kdz;
theta = 0;
E = exp(1i*kg*Z1*cos(0));
E1=E(:,:,5);
end
%%%%%%%%%%% sc_potc function %%%%%%%%%%%%%%%%%%
function [f]=sc_potc(N,P)
% HERE WE ARE CONSIDERING 1 UNIT = 10 MICROMETER
% N=75;
% P=38;
shape = [N,N,N];
N1=shape(1);
N2=shape(2);
N3=shape(3);
nb=1.33;
%nb=1.329-0.0516i;
um=1;
%um = 1e-6;
c = 3e8; % m/s
dx = um;
dy = um;
dz = um;
R=(N1/3)*dx;
R1=(N1/2.6)*dx;
lamb=5.32*um; % wavelength
xx = linspace(0,N1-1,N1);
yy = linspace(0,N2-1,N2);
zz = linspace(0,N3-1,N3);
xx=xx*dx;
yy=yy*dy;
zz=zz*dz;
[XX,YY,ZZ] = meshgrid(xx,yy,zz);
center = [(N/2)*dx,(N/2)*dy,(N/2)*dz];
rr = sqrt((center(2)-XX).^2 + (center(1)-YY).^2 + (center(3)-ZZ).^2);
center1 = [(N/2.5)*dx,(N/2.5)*dy,(N/2)*dz];
rr1 = sqrt((center1(2)-XX).^2 + (center1(1)-YY).^2 + (center1(3)-ZZ).^2);
center2 = [(N/1.6)*dx,(N/1.6)*dy,(N/2)*dz];
rr2 = sqrt((center2(2)-XX).^2 + (center2(1)-YY).^2 + (center2(3)-ZZ).^2);
epsr = ones(shape)*(1.329+0.0516i);
mur = ones(shape);
k0=(2*pi*nb)/lamb;
% 2 region
for ii=1:1:size(rr,1)
for jj = 1:1:size(rr,2)
for kk=1:1:size(rr,3)
if rr(ii,jj,kk) <= R
epsr(ii,jj,kk) = epsr(ii,jj,kk)+(0.059+0.0009i);
else if rr(ii,jj,kk) >= R && rr(ii,jj,kk) <= R1
epsr(ii,jj,kk) = epsr(ii,jj,kk)+(0.126 + 0.0159i);
else if rr(ii,jj,kk) >= R1
epsr(ii,jj,kk) = 1.33;
end
end
end
end
end
end
RR1=(N1/8)*dx;
for ii=1:1:size(rr1,1)
for jj = 1:1:size(rr1,2)
for kk=1:1:size(rr1,3)
if rr1(ii,jj,kk) <= RR1
epsr(ii,jj,kk) = epsr(ii,jj,kk)+(0.0465+0.0161i);
end
end
end
end
RR2=(N1/12)*dx;
for ii=1:1:size(rr2,1)
for jj = 1:1:size(rr2,2)
for kk=1:1:size(rr2,3)
if rr2(ii,jj,kk) <= RR2
epsr(ii,jj,kk) = epsr(ii,jj,kk)+(0.0465+0.0161i);
end
end
end
end
eps1 = epsr(:,:,P);
for i1=1:size(eps1)
for j1=1:size(eps1)
%f(i1,j1)=k0^2*((eps1(i1,j1).^2/nb.^2)-1);
f(i1,j1)=k0^2*(eps1(i1,j1).^2 - nb.^2);
end
end
end

Alvaro on 25 Jan 2023
There is an implementation of the Levenberg–Marquardt algorithm in the Optimization Toolbox that you might find useful.

### Categories

Find more on Historical Contests in Help Center and File Exchange

R2022b

### Community Treasure Hunt

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

Start Hunting!