fminsearch not minimizing function

I am working on a problem of fitting a basic 1D model of optical fluence attenuation in a multilayered biological tissue to experimental data. The model consists of a layer of fat and a layer of venous blood, and is written in terms of their optical and physiological properties. Furthermore, the model is fitted to the experimental data at multiple optical wavelengths (19) to account for the wavelength dependence of optical fluence in biological tissue. I wish to minimise the sum of the square of the difference between the model and the experimental data, and solve for the concentrations of oxy- and deoxyhaemoglobin, CHb and CHbO2, in the vein. To achieve this I use the following code:

 %
error = @(F0_Vein,mua_HbO2,mua_Hb,e_HbO2,e_Hb,musr,Cb_Fat,sO2_Fat,CHb,CHbO2) sum(abs(gamma*(2.303*(e_Hb.*CHb + e_HbO2.*CHbO2)).*F0_Vein.*exp(...
        - Z_F.*sqrt(3.*(Cb_Fat.*(sO2_Fat.*mua_HbO2 + (1-sO2_Fat).*mua_Hb)).*(musr + Cb_Fat.*(sO2_Fat.*mua_HbO2 + (1-sO2_Fat).*mua_Hb)))...
        - Z_V.*sqrt(3*(2.303*(e_Hb.*CHb + e_HbO2.*CHbO2)).*(musr + 2.303*(e_Hb.*CHb + e_HbO2.*CHbO2)))) - PA_Vein/trans_Vein).^2);
 %                             
options = optimset('PlotFcns',@optimplotfval,'MaxFunEvals',1000,'MaxIter',1000,'Display','iter');
[fitted_Param,fval,exitflag,output] = fminsearch(@(day) error(F0_Vein,mua_HbO2,mua_Hb,e_HbO2,e_Hb,musr,Cb_Fat,sO2_Fat,day(1),day(2)),[1;1],options);

Which passes the function ‘error’ to fminsearch and solves for CHb and CHbO2 in the vein with an initial guess of [1,1]. However, the code always returns the initial guess values no matter what they are, and fminsearch indicates that the function converged to a solution (exitflag = 1). If anyone could shed some light on what’s wrong with the above code, or if there’s a more suitable minimisation procedure I could use, it would be greatly appreciated.

Also please note since the fluence is modelled at 19 wavelengths, the optical parameters F0_Vein,mua_HbO2,mua_Hb,e_HbO2,e_Hb,musr & PA_Vein are all vectors with 19 elements. I think Matlab may then consider the fluence to be a function of at least 7 variables, and this may be causing an error, but I am not sure. Gamma, Cb_Fat, sO2_Fat, and trans_Vein are known constants.

I have attached the file 'optical_data.mat' to enable the above code to be run.

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[Star Strider]

I usually use the genetic algorithm ga function for complicated problems such as your. I tweaked your ‘error’ funciton to use norm rather than sum-of-squares, and with it:

error = @(F0_Vein,mua_HbO2,mua_Hb,e_HbO2,e_Hb,musr,Cb_Fat,sO2_Fat,CHb,CHbO2) norm(abs(gamma*(2.303*(e_Hb.*CHb + e_HbO2.*CHbO2)).*F0_Vein.*exp(...
        - Z_F.*sqrt(3.*(Cb_Fat.*(sO2_Fat.*mua_HbO2 + (1-sO2_Fat).*mua_Hb)).*(musr + Cb_Fat.*(sO2_Fat.*mua_HbO2 + (1-sO2_Fat).*mua_Hb)))...
        - Z_V.*sqrt(3*(2.303*(e_Hb.*CHb + e_HbO2.*CHbO2)).*(musr + 2.303*(e_Hb.*CHb + e_HbO2.*CHbO2)))) - PA_Vein/trans_Vein));

the best parameter estimates it returned were:

day =
   0.000402374029160   0.001968620181084

fval =
     4.185954018933532e+04

with ‘fval’ being the final value of ‘error’, so it seems initial estimates of [1,1] may have been too far from the global minimum. I am not certain what units you are using for ‘CHb’ and ‘CHbO2’. If you want to constrain them to specific ranges, ga permits that. My ga call is:

vftns = @(day) error(F0_Vein,mua_HbO2,mua_Hb,e_HbO2,e_Hb,musr,Cb_Fat,sO2_Fat,day(1),day(2));
opts = optimoptions('ga', 'PopulationSize',5000, 'InitialPopulationMatrix', randi(1000, 5000, 2)*1E-3, 'PlotFcn','gaplotbestf');
[day,fval] = ga(vftns, 2, [], [], [], [], [1 1]*1E-6, [1 1]*Inf, [], [], opts)

[aidan]

Dear Star Strider,

Thank you for your answer. The initial guess of [1,1] is indeed too far from the global minimum. I lowered the initial guess to [1e-5.1e-5] and found I get almost identical results with fminsearch as you did with ga.

CHb and CHbO2 are in units of mol L-1, and are set in the simulation at 0.000583 and 0.00175, so the fitted parameters are a little off in an absolute sense, but its their relative value that is important.

However, when I run your code for ga I get the following error, maybe because we are using different versions of Matlab?

 Error using optimoptions (line 114)
 Invalid solver specified. Provide a solver name or handle (such as 'fmincon' or @fminunc).
Type DOC OPTIMOPTIONS for a list of solvers.

[Star Strider]

As always, my pleasure.

I am using R2018a. You may need to use gaoptimset (link) instead. I believe most of the options have the same names.