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

Question

Warren Boschen 2022-10-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1822728-how-do-i-write-a-function-for-fmincon-using-a-for-loop

编辑： Matt J 2022-10-21

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2022-10-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1822728-how-do-i-write-a-function-for-fmincon-using-a-for-loop#answer_1071843

编辑：Matt J 2022-10-21

在 MATLAB Online 中打开

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

7 个评论
显示 5更早的评论隐藏 5更早的评论

Warren Boschen 2022-10-13

在 MATLAB Online 中打开

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

Torsten 2022-10-21

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

Matt J 2022-10-21

编辑：Matt J 2022-10-21

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.

请先登录，再进行评论。