How can I use 'Least Squares Solver and Jacobian' instead of 'fminunc'?

I want to use Least Squares Solver and Jacobian instead of fminunc, but I don't know what is the required changes must be done in the attached code.

Can anyone help me?

clc;
clear;
% System paramters:
N = 2;
K = 4;
C_l = 4;
counter = 0;
H = [-0.3208  -0.9784; -1.5994   -1.4689; -1.5197   -0.4568; -0.0993   -0.7667]; % 4*2 matrix
A = [-1 1; 0 1]; % 2*2 matrix
C = [-0.20  0.4  0.6 -0.2; -0.2  0.4  0.6 -0.2; 0.4  0.2 -0.2  0.4; 0.4  0.2 -0.2  0.4]; % 4*4 matrix
P = [250000   0   0   0; 0  250000  0   0; 0   0  250000  0; 0  0  0  250000]; % 4*4 matrix
beta_u = [50.2207   50.2207   20.3433   20.3433]; % 1*4 vector
beta_d = zeros(1,4); % intial value
%store inputs to a struct for shorter syntax
s=struct;[s.H,s.A,s.C,s.P,s.C_l,s.N,s.K]=deal(H,A,C,P,C_l,N,K);
while (sum(abs(beta_u-beta_d))>0.1 && counter< 500)
    initial_guess = randn(2,4);
    OLS = @(B_d,input_vars)sum(abs(beta_u-myfun(B_d,input_vars)).^2);%ordinary least squares cost function
    opts = optimoptions(@fminunc,'MaxIterations',10000,'MaxFunctionEvaluations',50000);
    [B,FVAL,grad] = fminunc(OLS, initial_guess, opts,s);
    input_vars = s;
    [beta_d, D_d]=myfun(B,input_vars);
    counter = counter+1;
end
%calculate beta_d from B and the other inputs
function [beta_d, D_d]=myfun(B,input_vars)
%load parameters
s=input_vars;[H,A,C,P,C_l,N,K]=deal(s.H,s.A,s.C,s.P,s.C_l,s.N,s.K);
for j=1:1:N
    d(j) = (B(j,:)*P*B(j,:)')/((2^(2*C_l))-(norm(A(:,j))^2));
end
D_d = diag(d);
for i=1:1:K   
    V_d(i)=C(i,:)*P*B'*H(i,:)'*inv(1+H(i,:)*(A'*D_d*A+B*P*B')*H(i,:)');
    sigma_d(i)=norm((V_d(i)*H(i,:)*B-C(i,:))*(P^(1/2)))^2+(V_d(i)^2)*(1+H(i,:)*A'*D_d*A*H(i,:)');
    beta_d(i)=((P(i,i))/sigma_d(:,i));   
end
end

[Matt J]

Is D_d the Jacobian of myfun?

[Ahmed Nasr]

Yes

Answer 1

[Matt J]

Assuming D_d is the Jacobian of myfun, it would be,

while (sum(abs(beta_u-beta_d))>0.1 && counter< 500)
    
    initial_guess = randn(2,4);

    opts = optimoptions(@lsqcurvefit,'MaxIterations',10000,'MaxFunctionEvaluations',50000,...
                        'SpecifyObjectiveGradient',true); 
    
    [B,FVAL,residual,exitflag,output,lambda,jacobian]= ...
     lsqcurvefit(@myfun, initial_guess,input_vars,beta_u,[],[],opts);
    
    grad=jacobian.'*residual(:);
 
    input_vars = s;
    [beta_d, D_d]=myfun(B,input_vars);
    
    counter = counter+1;
end

[Ahmed Nasr]

I receive the following error when I am running the proposed code. How can I solve it?

Error using lsqncommon (line 45)
Supplied Jacobian is not the correct size:
the Jacobian matrix should be 4-by-8.

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

Error in LSalgorithm (line 32)
     lsqcurvefit(@myfun, initial_guess,s,beta_u,[],[],opts);

[Matt J]

The error message has told you that your Jacobian calculation D_d should be 4-by-8. Why isn't it?

[Ahmed Nasr]

But, D_d must be 2*2 matrix.

[Matt J]

You told us above that D_d is the Jacobian of myfun. That means D_d must have as many columns as there are unknown variables, in this case 8.