larger computation time for every iteration of levenberg marquardt algorithm?

2 vues (au cours des 30 derniers jours)
asim asrar le 18 Oct 2022
Réponse apportée : Alvaro le 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
0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

Réponses (1)

Alvaro le 25 Jan 2023
There is an implementation of the Levenberg–Marquardt algorithm in the Optimization Toolbox that you might find useful.
0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

Catégories

En savoir plus sur Image Arithmetic dans Help Center et File Exchange

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by