Search code examples
matlableast-squaresminimization

Least squares minimization for multiple variables matlab


I need to find the value of tree variables: a, b and c, by finding a global minimum for least squares method. My function is as follows:

f = (1/a)*(asinh((Z(i)/b)^(1/c))^(-1)

where i is the index of vector Z. The vector Z has 9 values given in the task. I also have a vector with 9 values for s. The least squares method need to sum differences between values computed by function f with the values from vector s. It should look something like this:

((s(i)-f(i))/s(i))^2

I also have the boundaries for value a, b and c:

10e10>a>10e19, 10e-7>b>50, 10e-15>c>10.

I tried to use the lsqnonlin but I have no idea how to do it. I will appreciate any of your help!

I tried to do something like this:

function f=Fsigma(x, Z, sigma)
f=0;
for i=1:length(sigma)
    f=f+((sigma(i)-((1/x(1,:))*(asinh((Z(i)/x(2,:))^(1/x(3,:)))^(-1))))/sigma(i))^2
end
end

and call this function in lsqnonlin like this:

Z= [1.49E+18 1.49E+19 1.49E+20 1.99E+15 1.99E+16 1.99E+17 1.49E+13 1.49E+14 1.49E+15];
sigma = [55.1705 79.1016 105.636 25.4809 40.8572 61.7238 12.8147 21.4054 34.8319];

a=linspace(10e10,10e19);
b=linspace(10e-7,50);
c=linspace(10e-15,10);
x=[a; b; c];
p=lsqnonlin(Fsigma(x,Z,sigma));

Solution

  • OK @Agata, I'll go through it in detail, so you can learn some basics.

    First, how to hand functions to functions: function handles aka. the @-operator

    fnc = @(x) Fsigma(x,Z,sigma);
    

    fnc is an object -- in fact, it is a something that points at the function Fsigma. However, it even overshadows the additional inputs of Fsigma by declaring that its only input is x (this is called an anonymous function handle, indicated by the (), in between which you define the inputs that someone calling fnc may use. The other inputs to Fsigma in this line are the values of the variables that they have in this line.

    boundaries Boundaries should be provided as vectors:

    % bounds
    lb = [  10e10;
            10e-7;
            10e-15];
    ub = [  10e19;
            50;
            10];
    

    initial guess + optimization call if you read the docs of lsqnonlin, it requires an initial guess

    % initial guess
    x0 = ones(3,1);
    % optimization call: x = lsqnonlin(fun,x0,lb,ub)
    [p,fval] = lsqnonlin(fnc,x0,lb,ub)
    

    Note to the cost function lsqnonlin is better if you do not provide a summed cost but an array of errors (see the docs, again). So I adjust you Fsigma-function

    function f = Fsigma(x, Z, sigma)
    f = ((sigma-((1/x(1))./ asinh((Z./x(2)).^(1/x(3)))) )./sigma).^2;
    end
    

    If you want to have use a single output -- as most optimization algorithms require from their cost function -- you can use fmincon

    Have a look at the full code

    Z = [1.49E+18 1.49E+19 1.49E+20 1.99E+15 1.99E+16 1.99E+17 1.49E+13 1.49E+14 1.49E+15];
    sigma = [55.1705 79.1016 105.636 25.4809 40.8572 61.7238 12.8147 21.4054 34.8319];
    
    % bounds
    lb = [  10e10;
            10e-7;
            10e-15];
    ub = [  10e19;
            50;
            10];
    % initial guess
    x0 = ones(3,1);
    
    % create an anyonymous function handle (using @(x)
    fnc = @(x) Fsigma(x,Z,sigma);
    % optimization call: x = lsqnonlin(fun,x0,lb,ub)
    [p,fval] = lsqnonlin(fnc,x0,lb,ub);
    
    % optimization call: x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub)
    fnc2 = @(x) sqrt(sum(fnc(x).^2));
    [p2,fval2] = fmincon(fnc2,x0,[],[],[],[],lb,ub);
    
    
    % cost function
    function f = Fsigma(x, Z, sigma)
    f = ((sigma-((1/x(1))./ asinh((Z./x(2)).^(1/x(3)))) )./sigma).^2;
    end
    

    PS: note that this community is not a coding service, so do your readings the next time!