Using 'fminunc', I receive Optimization completed because the size of the gradient is less than the value of the optimality tolerance.

I use 'fminunc' to find the value of B (2*4 matrix) that minimzes the difference between the corresponding elements in two vectors as indicated in the attached code. In other words, I want to find the B (2*4 matrix) that makes the elements of beta_d (1*4 vector) which is a function of B matrix, equal to the corresponding ones of a "given" beta_u (1*4 vector), i.e. beta_d(1,1) = beta_u(1,1) && beta_d(1,2) = beta_u(1,2) && beta_d(1,3) = beta_u(1,3) && beta_d(1,4) = beta_u(1,4).

However, I usually receive the following message without getting any result and the program seems to go on an infinite loop!

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

<stopping criteria details>

the code:

clc;
clear;

% System paramters:
N = 2;
K = 4;
C_l = 4;
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
B = zeros(2,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);

%fminunc optimization
while (sum(abs(beta_u-beta_d))>=0.1)
    initial_guess = randn(2,4);
    OLS = @(B_d,input_vars)sum((beta_u-myfun(B_d,input_vars)).^2);%ordinary least squares cost function
    opts = optimoptions(@fminunc,'MaxIterations',10000,'MaxFunctionEvaluations',50000,'CheckGradients',true);
    B = fminunc(OLS, initial_guess, opts,s);
    input_vars = s;
    [beta_d, D_d]=myfun(B,input_vars);      
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

Answer 1

[Matt J]

If the problem doesn't have a solution satisfying the tolerance condition in your while loop

(sum(abs(beta_u-beta_d))<0.1)

then the loop has no way of stopping.

[Ahmed Nasr]

But, when I change the values of the given H, A, C, P matrices and beta_u vector, the fminunc sometimes gives me a result matched with the while function! .. I think this means that the math problem is solvable such as:

clc;
clear;

% System paramters:
N = 2;
K = 4;
C_l = 4;
% >>> New values 
H = [-0.5119 -0.8390; -0.0719  0.6697; 1.4303  2.0519; 1.4391 0.6676]
A = [-1 1; 1 0];
C = [-0.20  0.4  0.2 -0.6; -0.2  0.4  0.2 -0.6; -0.0154 0.1077 -0.0615 -0.2769; -0.0154 0.1077 -0.0615 -0.2769];
P = [2500   0   0   0; 0  2500  0   0; 0   0  2500  0; 0  0  0  2500];
beta_u = [22.6467   22.6467   17.5767   17.5767];
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)
    initial_guess = randn(2,4);
    OLS = @(B_d,input_vars)sum((beta_u-myfun(B_d,input_vars)).^2);%ordinary least squares cost function
    opts = optimoptions(@fminunc,'MaxIterations',10000,'MaxFunctionEvaluations',50000,'CheckGradients',true);
    B = fminunc(OLS, initial_guess, opts,s);
    input_vars = s;
    [beta_d, D_d]=myfun(B,input_vars);      
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

The optimization is solved correctly, and the obtained beta_d is equal to beta_u as required as:

beta_d =

   22.6467   22.6467   17.5767   17.5767

I want to know why when I use the values of H, P, A, C, and beta_u indicated in the question, the optimization is not solved?!

[Matt J]

Because a solution may not exist for certain choices of the parameters. In those cases, your while loop has no way of deciding when to stop trying to find one.

[Ahmed Nasr]

But, If I remove the while loop in the code indicated in the question, the result of the fminunc is not the desired one as I want the elements of beta_d to be equal that in beta_u, and I have beta_u is equal

beta_u =

   50.2207   50.2207   20.3433   20.3433

wheras, the obtained beta_d from the optimization after removing while function is:

   beta_d =

   50.0010   50.1870   20.3158    3.8982

which is not the desired output. How can I overcome this problem without using while loop and enter on an infinite loop?

[Matt J]

the result of the fminunc is not the desired one as I want the elements of beta_d to be equal that in beta_u

But again, you cannot be sure that such a thing is possible for all choices of H, P, A, C, and beta_u. To handle such cases, you need to rewrite the while loop so that it gives up the search at some point and reports that no solution could be found.

[Ahmed Nasr]

How can I do this?

[Matt J]

How about this?

k=0;
while (sum(abs(beta_u-beta_d))>=0.1)
    k=k+1;
    
    initial_guess = randn(2,4);
    OLS = @(B_d,input_vars)sum((beta_u-myfun(B_d,input_vars)).^2);%ordinary least squares cost function
    opts = optimoptions(@fminunc,'MaxIterations',10000,'MaxFunctionEvaluations',50000,'CheckGradients',true);
    B = fminunc(OLS, initial_guess, opts,s);
    input_vars = s;
    [beta_d, D_d]=myfun(B,input_vars);      
    
    if k>100
        warning( 'Search for a solution failed too many times. Bailing out.' )
        break;
    end 
    
end

[Ahmed Nasr]

Thanks for your help, the code now is working well. The last question is during the running of fminunc on the same code, sometimes the code is running correctly and gives me a result and sometimes I receive the following error message !! .. What's the problem?

Error using lineSearch
Search direction is not a descent direction; roundoff errors may be affecting convergence.

Error in fminusub (line 186)
          lineSearch(funfcn,x,dir,f,dirDerivative, ...

Error in fminunc (line 472)
   [x,FVAL,GRAD,HESSIAN,EXITFLAG,OUTPUT] = fminusub(funfcn,x, ...

Error in random_algorithm (line 26)
    B = fminunc(OLS, initial_guess, opts,s);

[Matt J]

I would need to have the code that procudes the error.

[Ahmed Nasr]

This code sometimes gives me a result, and sometimes gives me the previous mentioned error!

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]

fminunc is not succeeding in approximating the gradient of your objective function accurately using its default finite difference method. You could try supplying your own gradient calculation using the SpecifyObjectiveGradient option,

https://www.mathworks.com/help/optim/ug/writing-scalar-objective-functions.html#bu2w6h1