MATLAB: Fminsearch not minimizing function

errorfminsearchMATLABminimization

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.

Best Answer

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)
Related Question