fmincon gets very slow after few time

Hi all,

I am using fmincon for finding endmembers in each pixel of an image. I have a double for loop which goes pixel by pixel applying fmincon. When the code starts, each pixel optimization is pretty fast (miliseconds), but as the run advances it starts getting slower and slower. It does not have to do with memory allocation issues since I do it. Also, if I close Matlab and reopen it, if I continue in the last pixel it gets fast again for the first pixels, but again slower and slower until it takes lime one minute per pixel.

Here the example of code I am using:

maps = single(zeros(rows,columns,3)); % Create empty maps.

cont = 0;

for row = 1:rows

for col = 1:columns

cont = cont + 1;

C_MSE_subt = fmincon( @(x)SAMSID_subtractive(double(squeeze(cube(row,col,:))),E,x'), init_point, [-1 0 0; 0 -1 0; 0 0 -1 ], [0; 0; 0], [1 1 1], 1, [0 0 0], [1 1 1],[],options);

maps(row,col,:) = reshape(C_MSE_subt,1,1,3);

end

end

The optimizatiob function is pretty simple, just calculating a spectrum and comparing it against the one passed as input:

function SAMSID = SAMSID_subtractive(spectra1,E,C)

spectra2 = subtractive_model(E,C);

[~,~,SAMSID]=samsid(spectra1,spectra2);

It would seem normal to me if every pixel behaved slow. But the few first pixels very fast and then get slower and slower... Also I checked RAM and CPU and none of them are busy.

If I change the optimization function it does not get slow. For instance, using MSE instead of SAMSID metric, it is fast from the beggining to the end:

function MSE = mean_squared_error_subtractive(spectra1,E,C)

spectra2 = subtractive_model(E,C);

MSE = sum((spectra1-spectra2).^2) / size(spectra1,1);

Any ideas of why this could happen?

Thanks in advance!

[Matt J]

What do you mean by "getting slower"? How are you observing changes in speed?

Also, what is samsid()? It does not seem to be a Mathworks-provided function.

[Matt J]

Incidentally, you should get rid of your inequality constraints,

A=[-1 0 0; 0 -1 0; 0 0 -1 ]; b= [0; 0; 0];

They are redundant with your explicit lower bounds lb=[0 0 0] and are just slowing things down.

Answer 1

[Matt J]

I suspect samsid() is just not as smooth or well-conditioned as mse() for some pixels, and therefore takes more iterations to terminate. You can check the fmincon output to see how many iterations are being performed in each pixel and see if it correlates with the slow-down you are observing.

[Matt J]

Miguel's reply moved here.

Hi Matt,

Thanks for answering so quickly. When I say that it gets slower, I mean that in the first iteration of the double for loop it takes few miliseconds (I see that using a disp() function which tells me the percetage of pixels already processed). However, after more and more pixels are processed, the optimization takes longer and longer. Hence, for instante in pixel number 200, it takes 1 minute. This could lead me to think that some pixels take longer than others in the optimization. But if I stop the run in pixel 99, close Matlab, re-open Matlab, anc continue on pixel number 100, again the first pixels will take miliseconds to run, and will get slower and slower as the program continues.

This is the content of the samsid function:

function [sm,sd,samsid] = samsid(sp1,sp2)

for i=1:size(sp1,2)

sm(i) = sam(sp1(:,i)',sp2(:,i)');

sd(i) = sid(sp1(:,i)',sp2(:,i)');

end

samsid = sin(sm).*sd;

It computes a combination of spectral angle mapper and spectral divergence as an spectral error metric (both implemented in Matlab).

It is true that samsid takes longer to compute than MSE. However, what I see not normal is that it does not take the same time if it runs the first pixel or if it runs the pixel number 100. The higher the order the slower it gets.

I hope I am explaining well. Please let me know if I don't. Thanks a lot again for your help.

[Matt J]

But if I stop the run in pixel 99, close Matlab, re-open Matlab, anc continue on pixel number 100, again the first pixels will take miliseconds to run, and will get slower and slower as the program continues.

Are you using the same initial point init_point when you resume from pixel 100 as when you start the loop from pixel 1?

Also, how many iterations are being run by fmincon in each case?

I see that using a disp() function which tells me the percetage of pixels already processed

It would be better if you used tic/toc. disp() is not synchronized with the rest of the code.

[Miguel Ángel]

Yes, I was using the same initial point for every optimization. Right now I followed your optimization tips and included the change of initial point so that for each pixel, the x0 initial point is the solution found in the previous pixel.

I am also using tic,toc functions (sorry I forgot to mention). If I build a pixelwise optimization times map, given that I start the optimization from top-left pixel and advance rightwards and downwards, I see how the optimization times increase for each pixel. The following is a plot of the optimization times of the first 420 pixels processed:

As you can see, the time is increasing considerably in each new pixel. But there is one more funny fact. If I stop the run, clear the Matlab workspace and run the code again, starting from the beggining with the same pixels. The optimization time for the first pixels starts being as slow as the last optimization time done in the previous run of the code, and keeps increasing. See the following plot with the optimization times of the same fisrt 137 pixels as in the previous plot (I stopped it before 420 cause it was pretty slow, but optimization time keeps increasing and increasing for each pixel reaching values of decens of seconds or even minutes if I let it run for hours):

However, if I close Matlab, reopen it, and run the code again exactly in the same way, it starts fast again and get slower in each pixel.

I am also saving and checking the exitflag parameter for all pixels optimization, and all of them are either 1 or 2.

Finally, if I plot the number of iterations taken for each pixel, I see no increase trend in the number of iterations:

[Matt J]

Are the function called by your main code (samsid(), subtractive_model(),etc...) local functions or nested functions? It might be good to see the contents of subtractive_model as well.

Also, if running in Windows, check the Task Manager to see if there is an increase in RAM consumption as the pixel loop progresses.

Also, have you looked at the results of the optimization? Are they all reasonable, notwithstanding the increase in computation time?

[Miguel Ángel]

Here are the codes for the functions called:

SAMSID_subtractive applies the forward model creating a new spetrum by combining the endmembers E in the given concentrations C (this is done in the call to subtractive_model). Then it computes the samsid metric between the given spectrum (spectra1) and the created spectrum (sectra2):

function SAMSID = SAMSID_subtractive(spectra1,E,C)

spectra2 = subtractive_model(E,C);

[~,~,SAMSID]=samsid(spectra1,spectra2);

function mixed_spectra = subtractive_model(E,C)

q = size(E,2);

N = size(C,1);

w = size(E,1);

mixed_spectra = ones(w,N);

for n = 1:N

for Q = 1:q

mixed_spectra(:,n) = mixed_spectra(:,n) .* E(:,Q).^C(n,Q);

end

end

Then samsid just call to two metrics implemented in Matlab to compared both spectra: spectral divergence (sid) and spectral angle mapper (sam):

function [sm,sd,samsid] = samsid(sp1,sp2)

for i=1:size(sp1,2)

sm(i) = sam(sp1(:,i)',sp2(:,i)');

sd(i) = sid(sp1(:,i)',sp2(:,i)');

end

samsid = sin(sm).*sd

I've been checking at the results of a bigger sample data. I've studied the first 14000 pixels of the image (a whole weekend running the code). Apparently, the maximum number of iterations is being reached several times, which makes the code go slow. Yet this does not explain the "sometimes fast - sometimes slow" behaviour when doing it over the same data different times (geting slow as the number of trials increase and fast again after resetting Matlab). Ram is not consumed since results are overwritten at every iteration and all variables are properly allocated in memory before the loops.

Here three figures of these results. Optimization time for each pixel (x axis in log scale):

Number of iterations for each pixel (x axis in lo scale):

Optimized metric obtained in the optimizatin for each ixel (x axis in log scale):

Seems like there is a periodical pattern even though the data does not have any periodicity or any given order.

[Matt J]

Ram is not consumed since results are overwritten at every iteration and all variables are properly allocated in memory before the loops.

Have you verified this in the Task Manager, as I suggested?

Also, are your functions local or nested?

Also, is the output reasonable?

[Miguel Ángel]

I did verify. The task manager has always been open with the resources tab selected. No ram variations. I have 64 GB RAM and not even 15 GB are consumed during the whole process.

Functions are local. They are not nested in the same Matlab script, but they are all as .m files in the same folder of the main code.

It is difficult to say whether it is reasonable or not, since in order to check this I would need to get the entire map and post-process it a bit to get the final presence map. But I only got barely 10 lines of it after more than three days of run. SAMSID metric is a combination of sam and sid metrics. The first one is an angle, and samsid takes the sine of this angle which is bounded in [0, 1]. SID is the spectral information divergence which has no upper bound (the lower the more similar the spectra are). SAMSID multiplies sine(sam)*sid, so also has no bounds.

[Matt J]

Functions are local. They are not nested in the same Matlab script, but they are all as .m files in the same folder of the main code.

If so, they are neither local nor nested.

https://www.mathworks.com/help/matlab/matlab_prog/types-of-functions.html

Answer 2

[Matt J]

A few more performance optimizations that you should consider:

maps = zeros(rows,columns,3,'single');   %<------ post conversion to single is unnecessary

cube=double(cube); %<----pre-cast the whole array to double

cont = 0;

for row = 1:rows
    for col = 1:columns
        cont = cont + 1;
        
        tmp=squeeze(cube(row,col,:)); %<------ avoid doing this in every function call
        
        C_MSE_subt = fmincon( @(x)SAMSID_subtractive(tmp,E,x'), init_point, [], [], [1 1 1], 1, [0 0 0], [1 1 1],[],options);
        maps(row,col,:) = reshape(C_MSE_subt,1,1,3);
        
        init_point(:)=C_MSE_subt(:);  %<----- the previous solution might be very close to the next one
    end
end

[Miguel Ángel]

Thanks for all valuable tips. I implemented them all:

cube = double(cube);

SAMSID = zeros(rows,columns,'single'); % Create empty SAMSID map.

ti = SAMSID; % Create empty optimization time map.

flag = SAMSID;

iter = SAMSID;

cont = 0;

for row = 1:rows

for col = 1:columns

cont = cont + 1;

clc, disp(['Generating pigments map: ' num2str(cont/num_pixels*100) '%...'])

tmp = squeeze(cube(row,col,:));

tic

[C_SAMSID_subt,SAMSID(row,col),flag(row,col),output] = fmincon( @(x)SAMSID_subtractive(tmp,E,x'), init_point, [], [], [1 1 1], 1, [0 0 0], [1 1 1],[],options);

ti(row,col) = toc;

iter(row,col) = output.iterations;

maps(row,col,:) = reshape(C_SAMSID_subt,1,1,3);

init_pont(:) = C_SAMSID_subt(:);

end

end