reducing the computation time for Levenberg marquardt algo of image reconstruction

tic()
clc;
clear all;
close all;
N =100;

[uin]=incident(N);
G=green_f(N);
H=G;
uin = uin(:);
I = eye(N^2,N^2);
im1 = imread('mono.jpg');
im2 = imresize(im1,[512,512]);
im3=im2(100:199,200:299);
im=im2double(im3);
im = imadjust(im,[0 0.8],[0,0.27]);
y=im(:);
f2=ones(N,N);
f2=f2(:);
stepsize = 0.8;
fo = f2;
po = ones(N,N);
po = po(:);
i=1;
tol =10e-4;
fcontinue=1;
obj_old=0;

while fcontinue

iter=i
ss(i) = stepsize;
maxit=20;
y0=rand(N,N);
y0=y0(:);
AA=I-G.*(fo');
[u]= cgs(AA,uin,[],maxit,[],[],y0);
X11=(u.*H).';
obj_new = 0.5 * ((norm(X11*fo - y.*po,2))^2 ) ;
ff1(i)=obj_new;
Jf=((I+fo.*(G/(I - G.*(fo'))))).*u.';
Hesf=real(Jf'*Jf);
r = (u.*H).'*fo - y.*po;
df = real(Jf'*H'*(r));
fn = fo -(inv(Hesf + stepsize*I))*df;
fn=max(fn,0);
ff(i)=norm(df,2);

Jp = (diag(y'));
Hesp = (Jp'*Jp);
dp=(-1*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

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;
else
    stepsize = 2*stepsize;
end
  
    fo = fn;
    po = pn;
ii1(i)=i;

obj_old=obj_new;
i=i+1;
figure(2)
subplot 121
imagesc((reshape(y,[N,N])))
colormap gray
colorbar
title('diffracted field phase')
subplot 122
imagesc((reshape(abs(fn),[N,N])))
colormap gray
colorbar
title('reconstructed field')
end

figure(2)
% subplot 231
% imagesc(abs(reshape(f,[N,N])))
% colormap gray
% colorbar
% title('phantom magnitude')
subplot 232
imagesc(abs(reshape(y,[N,N])))
colormap gray
colorbar
title('scattered field amplitude')
subplot 233
imagesc(angle(reshape(y,[N,N])))
colormap gray
colorbar
title('scattered field phase')
subplot 234
imagesc(abs(reshape(fn,[N,N])))
colormap gray
colorbar
title('reconstructed field amplitude')
subplot 235
imagesc(abs(reshape(pn,[N,N])))
colormap gray
colorbar
title('phase amplitude')
subplot 236
imagesc(angle(reshape(pn,[N,N])))
colormap gray
colorbar
title('phase angle')
figure(4)
subplot 131
plot(ii1,ff1(1:end-1))
subplot 132
plot(ii1,ff(1:end-1))
subplot 133
plot(ii1,pp(1:end-1))
toc();

%%%%%%%%%% QUESTION DESCRIPTION
% 
% above is the code for image reconstruction, 
% i have tried to implement levenberg marquardt algorithm. but the computational time of the algorithm is very high ,
% can someone help me to speed up the code , or can suggest what chages do i need to ddo to make the code 
% computationally efficient.
% 
% the function Green_f, incident is given there in the attachment.

Answer 1

[Yatharth]

Hi,

I understand that your code's computational time is very high and you would like to optimise it.

For the Levenberg-Marquardt Algorithm you can try using Optimization Toolbox, here are some resources that might be helpful. https://www.mathworks.com/help/optim/

https://www.mathworks.com/help/optim/ug/least-squares-model-fitting-algorithms.html#f204

Here are few suggestions which you can incorporate to make your code more efficient:

1. Pre-allocate Arrays: Pre-allocate arrays for ss, ff1, ff, pp, and ii1 to their final size before the loop. This avoids resizing the arrays in each iteration, resulting in improved performance.
2. Vectorize Operations: Utilize vectorized operations whenever possible to avoid unnecessary loops. For instance, instead of looping over pn to normalize its elements, you can use element-wise division by the norm directly.
3. Use Element-wise Operations: Instead of using matrix multiplication (*) and matrix division (/), use element-wise operations (.* and ./) where appropriate. This can significantly speed up calculations.
4. Profile and Optimize: Use MATLAB's profiling tools to identify the most time-consuming parts of the code. Focus on optimizing those sections by using built-in functions or alternative algorithms.

You can read more about Profiling your code here. https://www.mathworks.com/help/matlab/matlab_prog/profiling-for-improving-performance.html

I hope this helps.