How do I write a function for fmincon using a for loop?

Hello,

I am looking to minimize two parameters, and , with the following minimization:

I was recommended to do use the function fmincon to do this, but I'm having issues figuring out how to write the function "fun" as written on the documentation. Here is the function I have written:

function S0 = signal_min(atti, bvals, lambdas)
S0 = 0;
for b = 1:b_num_unique % For every unique b...
    % This loop determines the indices of a particular b-value.
    b_set = [];
    for i = 1:length(bvals)
        if bvals(i, :) == b_unique(b)
            b_set = [b_set; i];
        end
    end
    
    % This set of loops finds the spherical mean of the signal for a particular
    % b-value by first creating an array of the average values at each
    % point (x,y,z) and then finding the mean of that array.
    bset_avg = zeros(X, Y, Z, length(b_set));
    for i = 1:length(b_set)
        bset_avg = bset_avg + atti(:, :, slices(z), b_set(i));
    end
    S_avg = bset_avg./length(b_set);
    
    S0 = S0 + (S_avg(x, y, slices(z)) - sm_signal(bvals(b), lambdas(0), lambdas(1)))^2;
end

    function [signal] = sm_signal(bvalue, lambda_parallel, lambda_perpendicular)
        signal = exp(cast(-bvalue*lambda_perpendicular, 'double'))...
            *sqrt(pi/(cast(4*bvalue*(lambda_parallel - lambda_perpendicular), 'double')))...
            *erf(sqrt(cast(bvalue*(lambda_parallel - lambda_perpendicular), 'double')));
    end
end

And here is the error produced:

Not enough input arguments.

Error in smt_fmincon>signal_min (line 122)
for b = 1:bvals % For every b-value...

Error in fmincon (line 552)
      initVals.f = feval(funfcn{3},X,varargin{:});

Error in smt_fmincon (line 89)
               [lambda, sigval] = fmincon(@signal_min, lambda0, C, d, Ceq, deq, lb, ub);

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

From how I have seen other functions for fmincon written, they seem like they're just supposed to be one line. What I can't figure is how to include the b-dependence in just one singular expression, as the purpose of the minimization is to minimize the squared residuals across all b's. That is to say, I don't want to minimize and then take the sum which is what I think I would be doing if I were to remove the for loop from inside signal_min, along the lines of below:

for b = 1:bvals
    fmincon(@signal_min, x0, C, d, Ceq, deq, lb, ub);
end

What do you suggest I do to work around this issue?

Many thanks,

Warren

Answer 1

[Matt J]

Forget the loop. Also, use lsqcurvefit instead:

 params0=[lambda0(1), lambda0(1)-lambda0(2) ]; %re-parametrize

[params, sigval] = lsqcurvefit(@sm_signal, params0, double(bvalues),Savg(:),[0;0]);

lambda_parallel=params(1);  %undo re-parametrization
lambda_perpendicular=params(1)-params(2);

function [signal] = sm_signal( params,bvalues)
        
         bvalues=bvalues(:);
        [lambda_parallel, delta_lambda]=deal(params(1),params(2));
        
        signal = exp(-bvalues*lambda_perpendicular)...
            .*sqrt(pi./(4*bvalues*delta_lambda))...
            .*erf(sqrt(bvalues*delta_lambda));
end

[Warren Boschen]

What is the difference between lsqcurvefit and fmincon? When I posted my previous question asking how to approach this minimization with lsqcurvefit, someone recommended that I use fmincon instead.

[Matt J]

lsqcurvefit is for nonlinear least squares problems, possibly with bound constraints.

fmincon will handle more general objective functions and constraints if needed (but you don't).

[Warren Boschen]

Ah okay that makes sense, thank you. I'm having issues writing the [params, sigval] line with your adjustments to the code. Here is what I have:

[params, sigval] = lsqcurvefit(@sm_signal, lambda0, double(b_unique), double(S_avg(:, :, z, b)), 0, 3.0e-3);

And here is the error that I am getting:

Error using lsqcurvefit (line 263)
Function value and YDATA sizes are not equal.

lambda0 is a 2x1 double, b_unique is an 11x1 single, and S_avg(:, :, z, b) is a 128x128 double representing multiple voxels at a particular value of z. What is the output of sm_signal as it is written in your code? I'm looking to be able to apply this fit to every voxel of my data set. So far, I have:

for z = 1:length(slices)
    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.
                lambda_par(x, y, z) = 0;
                lambda_perp(x, y, z) = 0;
            else
                [params, sigval] = lsqcurvefit(@sm_signal, lambda0, double(b_unique),...
                double(S_avg(:, :, z, b)), 0, 3.0e-3);
                    
                lambda_par(x, y, z) = params(1);
                lambda_perp(x, y, z) = params(2);
            end
        end
    end
end

Thanks again,

Warren

[Matt J]

What is the output of sm_signal as it is written in your code?

If you run it with some input, you should be able to see.

[Warren Boschen]

That was great insight, thank you. How does the loop that you wrote account for the parameter ?

[Torsten]

You set this condition in C and d in your call to "fmincon" ?

[Matt J]

How does the loop that you wrote account for the parameter ?

Loop? There are no loops in the code I posted for you.

However, the constraint that you mention is enforced by imposing (using the 5th input argument to lsqcurvefit) a lower bound of 0 on both unknown variables. One of these variables is deltaLambda = and so deltaLambda>=0 is equivalent to your original constraint.