Clear Filters
Clear Filters

larger computation time for every iteration of levenberg marquardt algorithm?

3 views (last 30 days)
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

Answers (1)

Alvaro
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

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!