Error using fmincon MATLAB function

I would like to obtain the minimum of the function using MATLAB 'fmincon' function where the first eigenvalue derived from the given matrix is my objective function. I was to write/supply its objective function, but it is too long to put it inside the code. I wrote the following code, but I got stuck to continue. The error I found is:

Error using symengine
Unable to convert expression into double array.

I was wondering if I could get help? Thanks.

options = optimoptions('fmincon','GradObj','on');
x = fmincon(@myfun,[1,1],[],[],[],[],[1,1],[2,2],[],options)
function [f,g] = myfun(x)
syms x1 x2;
% 4 by 4 matrix with two variables inside
C=[1    2   x1  4;...
   1/2    1   3  x2;...
  1/x1  5   1  3/2;...
  1/4    1/x2 2  1];
f1 = eig(C); % calculates eigenvalues 
f=f1(1); % the first eigenvalue
if nargout > 1 % gradient required
% Calculate partial derivatives of f
    g(1) = diff(f,x1);
    g(2) = diff(f,x2);
    g = double([g(1);g(2)]);% data type double    
end
end

[Walter Roberson]

When you say "first" do you mean the one with largest absolute value?

[Walter Roberson]

using the eigenvalue with the largest absolute value is potentially discontinuous in the variables. fmincon requires differentiable jacobian but anything equivalent to a max() operation is not differentiateable at the boundary (and the jacobian potentially turns into zero with respect to some of the variables.)

Answer 1

[Matt J]

Because of the discontinuity at x1=0 and x2=0, you need to optimize over each quadrant separately.

 s=[+1,+1; %quadrant signs
    -1,+1;
    +1,-1;
    -1,-1];

[X1,X2] =ndgrid(linspace(1e-5,2,100)); 

clear Q xcell
for i=4:-1:1 %optimize over each quadrant
 
  F=arrayfun(@myfun,X1*s(i,1),X2*s(i,2)); %initial discrete search for initial guess
  [r,c]=find(F==min(F(:)),1);
 
   x0=[X1(r,c)*s(i,1),X2(r,c)*s(i,2)];
 
[xcell{i},Q(i)]= fminsearch(@(x)myfun(x(1),x(2)), x0 ,...
                           optimset('MaxFunEvals',1e8,'MaxIter',1e6));

end

[~,j]=min(Q); %pick the best quadrant


[fOptimal,xOptimal]=deal(Q(j),xcell{j}) %final solution

fOptimal = 4.0000

xOptimal = 1×2

1.5995 -0.1835

function fval = myfun(x1,x2)

 if abs(abs(x1)-1)>=1, fval=inf; return; end
 if abs(abs(x2)-1)>=1, fval=inf; return; end

  fval= max(abs(eig([1    2   x1  4;...
                    1/2    1   3  x2;...
                    1/x1  5   1  3/2;...
                     1/4    1/x2 2  1])));

end

[Matt J]

No, the solution I gave you is constrained. The bounds given in your post are -2<=x1<=2 and -2<=x2<=2 and that's what the code implements via the lines,

 if abs(abs(x1)-1)>=1, fval=inf; return; end
 if abs(abs(x2)-1)>=1, fval=inf; return; end

[Matt J]

 s=[+1,+1]; %quadrant signs


[X1,X2] =ndgrid(linspace(1,2,100)); 

clear Q xcell
for i=size(s,1):-1:1 %optimize over each quadrant
 
  F=arrayfun(@myfun,X1*s(i,1),X2*s(i,2)); %initial discrete search for initial guess
  [r,c]=find(F==min(F(:)),1);
 
   x0=[X1(r,c)*s(i,1),X2(r,c)*s(i,2)];
 
[xcell{i},Q(i)]= fminsearch(@(x)myfun(x(1),x(2)), x0 ,...
                           optimset('MaxFunEvals',1e8,'MaxIter',1e6));

end

[~,j]=min(Q); %pick the best quadrant


[fOptimal,xOptimal]=deal(Q(j),xcell{j}) %final solution

fOptimal = 6.3056

xOptimal = 1×2

2.0000 1.0001

function fval = myfun(x1,x2)

 if abs(abs(x1)-3/2)>=1/2, fval=inf; return; end
 if abs(abs(x2)-3/2)>=1/2, fval=inf; return; end

  fval= max(abs(eig([1    2   x1  4;...
                    1/2    1   3  x2;...
                    1/x1  5   1  3/2;...
                     1/4    1/x2 2  1])));

end

[Matt J]

As Walter told you, fmincon is not applicable to this problem because you do not have a differentiable cost function.

[Matt J]

On what intervals are x1,x2 bounded in this case?

[HAT]

@Matt J On the same intervals. i.e., both in [1,2]. Thanks.

[Matt J]

Well, I don't know why you think you need to supply your own gradient with fmincon. Using fmincon with its default finite difference approximation produces a minimum in close agreement both with fminbnd and with a surf plot of the objective,

%% Solution via fmincon

  [X1,X2]=ndgrid(linspace(1,2,101));

  F=arrayfun(@objfunc,X1,X2);
 [i,j]=find(F==min(F(:)),1);
 x0=[X1(i,j),X2(i,j)];
 
 fun=@(x)objfunc(x(1),x(2));
xfmincon=fmincon(fun,x0, [],[],[],[],[1,1],[+2,+2])

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

xfmincon = 1×2

1.4804 1.0000

%% Solution via fminbnd
x1=fminbnd(@fun1,1,2);
x2=fminbnd(  @(x2)objfunc(x1,x2),  1,2);

xfminbnd=[x1,x2]

xfminbnd = 1×2

1.4804 1.0001

%% Visual check
surf(X1,X2,F,'FaceAlpha',.4), xlabel x1, ylabel x2

function f=fun1(x1)
 [~,f]=fminbnd(  @(x2)objfunc(x1,x2), 1,2);
end

function  fval=objfunc(x1,x2)
 
    C=[1      2   exp(x1)  4;...
       1/2    1   3   exp(x2);...
      1/exp(x1)  1/3   1  1/5;...
      4    1/exp(x2)   5  1];
 
    fval=max(abs(eig(C)));

end

[HAT]

Ok, thanks!

[Matt J]

You're welcome, but please Accept-click the answer if you consider the question resolved.

Answer 2

[Walter Roberson]

You can use the symbolic toolbox to generate eig() on the matrix with Symbolic x1 and x2, and then you can matlabFunction that and pass the handle into the objective function. You would also calculate the g functions with respect to each of the four eigenvalues, and pass those handles as well.

The objective would then not have any symbolic operations. It would call the stored handle upon the input x values, getting back a vector of four eigenvalues. It would figure out which had the largest absolute value and store in f. If nargout requires, it would use the index to select the jacobian handle and evaluate it on the current inputs and return those.

[Walter Roberson]

% fmincon with gradient vector
syms x1 x2
vars = [x1, x2];
C=[1      2   x1  4;...
   1/2    1   3   x2;...
  1/x1  1/3   1  1/5;...
  4    1/x2   5  1];
f=eig(C) % generate eigenvalues on the matrix with Symbolic x1 and x2
eigFunction = matlabFunction(f, 'vars', {vars}); % matlabFunction that and pass the handle into the objective function

% % calculate the g functions with respect to each of the four eigenvalues
g = jacobian([f(1), f(2), f(3),f(4)], vars); % jacobian
gradFunction = matlabFunction(g, 'vars', {vars});

    x=ones(2,1); % initial value

% MATLAB fmincon
 options = optimoptions('fmincon','GradObj','on');
 x = fmincon(@(x)myfun(x, eigFunction, gradFunction), x, [], [], [], [], [1,1], [2,2], [], options)

 function [objFunction2 grad2] = myfun(x, eigFunction, gradFunction)
    
  t = eigFunction(x);
  [~, maxidx] = max(abs(t));
  objFunction2 = t(maxidx);
  
  % If nargout requires, it would use the index 
  %to select the jacobian handle and evaluate it on the current inputs and return those.
  if nargout > 1
       gradvalues = gradFunction(x);
       grad2 = gradvalues(maxidx);
  end
end