Error in using lsqcurvefit with constrained parameters.

Hi all,

I'm attempting to use the function lsqcurvefit to fit a single parameter, however I am providing the function with two other parameters that have been fixed using this previous post as a guide. Here is my code:

opts2 = optimset('Display', 'off', 'FinDiffType', 'central');
free_water = zeros(X, Y, Z); % The fraction of non-free water (unitless).
for z = 1:length(slices) % For every particular voxel...
    for y = 1:Y
        for x = 1:X
            if mask(x, y, slices(z)) == 0 % Set all values outside of the mask to be 0.
                free_water(x, y, z) = 0;
            else
                S_low_avg_formatted = zeros(low_num_unique, 1);  % Reformat the average signal within a voxel.
                for b = 1:low_num_unique % Vectorize eventually!!!
                    S_low_avg_formatted(b, 1) = S_avg_low(x, y, z, b);
                end
                
                params0 = [0; lambda_perp(x, y, z); deltaLambda(x, y, z)];
                fixedset = [0 1 1];
                fixedvals = [lambda_perp(x, y, z), deltaLambda(x, y, z)];
                params0_vary = params0(~fixedset);
                free_water(x, y, z) = lsqcurvefit(@(fw, data) low_wrapper(fw, data, fixedset, fixedvals),...
                    params0_vary, double(low_unique), S_low_avg_formatted, 0, 1, opts2);
            end
        end
    end
end
% --------------------------------------------------------------------------------------------------------------
function [signal] = sm_signal_low(f,perp,delta_lambda,bvalues)
D = 3.00e-3; % D is the diffusivity of free water.

signal = (1-f)*exp(-bvalues*D)+ f*(exp(-bvalues*perp).*sqrt(pi./(4*bvalues*(delta_lambda)))....
         *erf(sqrt(bvalues*(delta_lambda))));
end

function pred = low_wrapper(xvary, data, fixedset, fixedvals)
x = zeros(size(fixedset));
x(fixedset) = fixedvals;
x(~fixedset) = xvary;
pred = sm_signal_low(x,data);
end

And here is the error that I am receiving:

Array indices must be positive integers or logical values.

Error in smt_estimate>low_wrapper (line 343)
x(fixedset) = fixedvals;

Error in smt_estimate>@(fw,data)low_wrapper(fw,data,fixedset,fixedvals)

Error in lsqcurvefit (line 213)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in smt_estimate (line 279)
                free_water(x, y, z) = lsqcurvefit(@(fw, data) low_wrapper(fw, data, fixedset, fixedvals),...

Caused by:
    Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.  

I have a similar loop in a prior portion of the code which works fine, however I am not fixing any of the parameters in that fit so I believe this issue comes from an incorrect use of the low_wrapper() function. What am I doing wrong? For reference, lambda_perp and deltaLambda are 128x128x6 array, low_unique is a 5x1 array, and S_low_avg_formatted is also a 5x1 array. Neither lambda_perp nor deltaLambda have any values below 0, and I believe that any 0's would be excluded by the if/else statement.

Thank you!

-Warren

Answer 1

[Torsten]

Maybe

fixedset = logical([0 1 1]);
fixedvals = [3 4];
xvary = 33;
x = zeros(size(fixedset));
x(fixedset) = fixedvals;
x(~fixedset) = xvary;
x

x = 1×3

33 3 4

[Warren Boschen]

Works great! I think the key was setting the array to be a logical.

Answer 2

[Matt J]

I think you can set the upper bounds equal to the lower bounds to fix the parameters as the OP in the previous post originally wanted to do. lsqcurvefit in recent Matlab versions will now do the reformulation you are attempting to do automatically. It probably won't be as computationally efficient, however, as reformulating it yourself.

[John D'Errico]

In the current release, this is now true.

lsqcurvefit

x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub) defines a set of lower and upper bounds on the design variables in x, so that the solution is always in the range lb ≤ x ≤ ub. You can fix the solution component x(i) by specifying lb(i) = ub(i).

[Warren Boschen]

Do you know when this became true? I'm using MATLAB R2019a so I'm not sure.

[Matt J]

Do you know when this became true? I'm using MATLAB R2019a so I'm not sure.

Easiest is to try it and see if it works.

[Warren Boschen]

Yeah I don't think it's working. Here's the error I'm getting:

Error using erf
Input must be real and full.

Error in smt_estimate>sm_signal_low (line 342)
signal = (1-f)*exp(-bvalues*D)+ f*(exp(-bvalues*perp).*sqrt(pi./(4*bvalues*(delta_lambda))).*erf(sqrt(bvalues*(delta_lambda))));

Error in lsqcurvefit/objective (line 278)
         F = feval(funfcn_x_xdata{3},x,XDATA,varargin{:});

Error in finDiffEvalAndChkErr

Error in finitedifferences

Error in computeFinDiffGradAndJac

Error in sfdnls (line 54)
        computeFinDiffGradAndJac(x,funfcn,confcn,valx, ...

Error in snlsFixedVar>i_sfdnls (line 526)
[A,findiffevals] = sfdnls(xcurr,fvec,Jstr,group,alpha,funfcn,l,u,...

Error in snlsFixedVar (line 72)
    [JACval, initJacobFuncCount] = i_sfdnls(xstart, Jstr, fval, ...

Error in lsqncommon (line 160)
            snlsFixedVar(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...

Error in lsqcurvefit (line 271)
    lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...

Error in smt_estimate (line 281)
                free_water(x, y, z) = lsqcurvefit(@sm_signal_low, params0, double(low_unique), ...

I believe this is because delta_lambda is not being fixed during the fit. Typing the erf() section of the function into the console yields an answer though. Thoughts?

[Matt J]

We would have to see what code you ran, but it seems more likely that the sqrt() operations in your objective are generating complex values.

erf(1i)

Error using erf
Input must be real and full.

Your objective function should probably have a check in it anyway to verify that bvalues*delta_lambda is real and non-negative, since that is clearly what your model assumes.

[Matt J]

Do you know when this became true? I'm using MATLAB R2019a so I'm not sure.

It was definitely there in R2019a:

https://www.mathworks.com/help/releases/R2019a/optim/ug/lsqcurvefit.html