larger computation time for every iteration of levenberg marquardt algorithm?

Question

asim asrar 2022-10-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1829633-larger-computation-time-for-every-iteration-of-levenberg-marquardt-algorithm

回答： Alvaro 2023-1-25

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
 