Too many input arguments error in objective cost function using fminsearchbnd

Hello guys,

I have been stuck with fminsearchbnd for over a day, as I keep getting too many input arguments error in the objective cost function that I used for model parameter estimation.

I am really stuck and any help would be highly appreciated. Below are some portions of the Matlab code which include the optimization, objective cost function and the model.

Best regards,

Gerald

clc;clear; close all
load PrevalenceDataHepaB.mat % Hepa B prevalence data from 2010-2023

t_data = PrevalenceDataHepaB(:,1);
I_data = PrevalenceDataHepaB(:,2);
C_data = PrevalenceDataHepaB(:,3);

%% Optimization
init_par = [2.6797 0.1]; 
LB = [0 0];
UB = [5 1];
options = optimset('Display','iter','TolX',1e-15,'TolFun',1e-15,'MaxFunEvals',1.5e+7,'MaxIter',1.5e+8);
[xopt,fval,exitflag,output] = fminsearchbnd(@model_fit,init_par,LB,UB,options);

beta = xopt(1);
q = xopt(2);

display(xopt);

%% Objective Cost Function
function Y = model_fit(x)
beta = x(1);
q = x(2);

b = [2.6501e-02,2.1984e-02]; % estimated
mu=1/44.8; % estimated
d=4.5642e-04; % estimated
epsilon=0.16;
k=4;
p=0.05;
r1=2;
r2=0.025;
r3 = 0.01;
v=0.03;
w=0.05;

I_data1 = [25784 24243 19128 13528 11028 9518 8559 6188 2874 2628 1982 2242 1745 1548]';
C_data1 = [353403 376054 377160 376131 375944 377901 418554 424359 439848 453686 436344 444763 440314 445513]';

L0 = 6446;
I0 = 25784;
C0 = 353403;
R0 = 16626604;
N0 = 49879811.5;
S0 = N0-L0-I0-C0-R0;

init = [S0; L0; I0; C0; R0];
tspan = (2010:1:2023);

options = [];
[~,HB_model] = ode45(@model,tspan,init,options,b,mu,d,beta,epsilon,k,r1,r2,r3,w,p,q,v);
Y1 = sum((HB_model(:,3) - I_data1).^2);
Y2 = sum((HB_model(:,4) - C_data1).^2);
Y = Y1+Y2;
end

%% Model
function X = model(t,X0,b,mu,d,beta,epsilon,k,r1,r2,r3,w,p,q,v)

S = X0(1);
L = X0(2);
I = X0(3);
C = X0(4);
R = X0(5);

N = S + L + I + C + R;

by = b(1);
if t > 2017
    by = b(2);
end

dS = by*w*(N - v*C) - beta*S./N*(I+epsilon*C) - (mu + r3)*S;
dL = beta*S./N*(I+epsilon*C) - (k+mu)*L;
dI = q*k*L - (mu+r1)*I; 
dC = by*w*v*C + p*r1*I + (1-q)*k*L - (mu+d+r2)*C;
dR = by*(1-w)*N + (1-p)*r1*I + r2*C + r3*S - mu*R;

X = [dS; dL; dI; dC; dR];
end

[Stephen23]

Parameterize those functions correctly:

https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Gerald Navida]

I am not sure if there's really something wrong with the parameterization of the function, as this code worked using fmincon. I also want to try using fminsearchbnd, but I got too many input arguments error in the objective function. Also, it's my first time to use fminsearchbnd and I am not so familiar with it.

[Torsten]

Your cost function must return a scalar value. At the moment, you return a vector of length 14:

Y1 = (HB_model(:,3) - I_data1).^2;
Y2 = (HB_model(:,4) - C_data1).^2;
Y = Y1+Y2;

Thus us

Y1 = (HB_model(:,3) - I_data1).^2;
Y2 = (HB_model(:,4) - C_data1).^2;
Y = [Y1;Y2]
Y = sum(Y)

instead.

And you should call ode45 in two steps: Once from 2010 to 2017 and - continuing the integration - from 2017 to 2023 in order to avoid the jump discontinuity in the parameter b.

[Gerald Navida]

Yes, I just realised that I have deleted the sum of Y earlier and so I put it back again. I already tried to use one value for b from 2010 to 2023 to prevent the jump discontinuity, but I still get too many iput arguments error in the cost function.

[Torsten]

For me, the code works (at least without the unrealistic tolerances for the optimizer).

clc;clear; close all
%load PrevalenceDataHepaB.mat % Hepa B prevalence data from 2010-2023

%t_data = PrevalenceDataHepaB(:,1);
%I_data = PrevalenceDataHepaB(:,2);
%C_data = PrevalenceDataHepaB(:,3);

%% Optimization
init_par = [2.6797 0.1]; 
LB = [0 0];
UB = [5 1];
options = optimset('Display','iter','MaxFunEvals',100000,'MaxIter',100000);%,'TolX',1e-15,'TolFun',1e-15,);
[xopt,fval,exitflag,output] = fminsearchbnd(@model_fit,init_par,LB,UB);%,options);

beta = xopt(1);
q = xopt(2);

display(xopt);

xopt = 1×2

1.1032 0.4024

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% Objective Cost Function
function Y = model_fit(x)
beta = x(1);
q = x(2);

b = [2.6501e-02,2.1984e-02]; % estimated
mu=1/44.8; % estimated
d=4.5642e-04; % estimated
epsilon=0.16;
k=4;
p=0.05;
r1=2;
r2=0.025;
r3 = 0.01;
v=0.03;
w=0.05;

I_data1 = [25784 24243 19128 13528 11028 9518 8559 6188 2874 2628 1982 2242 1745 1548]';
C_data1 = [353403 376054 377160 376131 375944 377901 418554 424359 439848 453686 436344 444763 440314 445513]';

L0 = 6446;
I0 = 25784;
C0 = 353403;
R0 = 16626604;
N0 = 49879811.5;
S0 = N0-L0-I0-C0-R0;

init = [S0; L0; I0; C0; R0];
tspan = (2010:1:2023);

options = [];
[~,HB_model] = ode45(@model,tspan,init,options,b,mu,d,beta,epsilon,k,r1,r2,r3,w,p,q,v);
Y1 = sum((HB_model(:,3) - I_data1).^2);
Y2 = sum((HB_model(:,4) - C_data1).^2);
Y = Y1+Y2;
end

%% Model
function X = model(t,X0,b,mu,d,beta,epsilon,k,r1,r2,r3,w,p,q,v)

S = X0(1);
L = X0(2);
I = X0(3);
C = X0(4);
R = X0(5);

N = S + L + I + C + R;

by = b(1);
if t > 2017
    by = b(2);
end

dS = by*w*(N - v*C) - beta*S./N*(I+epsilon*C) - (mu + r3)*S;
dL = beta*S./N*(I+epsilon*C) - (k+mu)*L;
dI = q*k*L - (mu+r1)*I; 
dC = by*w*v*C + p*r1*I + (1-q)*k*L - (mu+d+r2)*C;
dR = by*(1-w)*N + (1-p)*r1*I + r2*C + r3*S - mu*R;

X = [dS; dL; dI; dC; dR];
end
function [x,fval,exitflag,output] = fminsearchbnd(fun,x0,LB,UB,options,varargin)
% FMINSEARCHBND: FMINSEARCH, but with bound constraints by transformation
% usage: x=FMINSEARCHBND(fun,x0)
% usage: x=FMINSEARCHBND(fun,x0,LB)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB,options)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB,options,p1,p2,...)
% usage: [x,fval,exitflag,output]=FMINSEARCHBND(fun,x0,...)
% 
% arguments:
%  fun, x0, options - see the help for FMINSEARCH
%
%  LB - lower bound vector or array, must be the same size as x0
%
%       If no lower bounds exist for one of the variables, then
%       supply -inf for that variable.
%
%       If no lower bounds at all, then LB may be left empty.
%
%       Variables may be fixed in value by setting the corresponding
%       lower and upper bounds to exactly the same value.
%
%  UB - upper bound vector or array, must be the same size as x0
%
%       If no upper bounds exist for one of the variables, then
%       supply +inf for that variable.
%
%       If no upper bounds at all, then UB may be left empty.
%
%       Variables may be fixed in value by setting the corresponding
%       lower and upper bounds to exactly the same value.
%
% Notes:
%
%  If options is supplied, then TolX will apply to the transformed
%  variables. All other FMINSEARCH parameters should be unaffected.
%
%  Variables which are constrained by both a lower and an upper
%  bound will use a sin transformation. Those constrained by
%  only a lower or an upper bound will use a quadratic
%  transformation, and unconstrained variables will be left alone.
%
%  Variables may be fixed by setting their respective bounds equal.
%  In this case, the problem will be reduced in size for FMINSEARCH.
%
%  The bounds are inclusive inequalities, which admit the
%  boundary values themselves, but will not permit ANY function
%  evaluations outside the bounds. These constraints are strictly
%  followed.
%
%  If your problem has an EXCLUSIVE (strict) constraint which will
%  not admit evaluation at the bound itself, then you must provide
%  a slightly offset bound. An example of this is a function which
%  contains the log of one of its parameters. If you constrain the
%  variable to have a lower bound of zero, then FMINSEARCHBND may
%  try to evaluate the function exactly at zero.
%
%
% Example usage:
% rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2;
%
% fminsearch(rosen,[3 3])     % unconstrained
% ans =
%    1.0000    1.0000
%
% fminsearchbnd(rosen,[3 3],[2 2],[])     % constrained
% ans =
%    2.0000    4.0000
%
% See test_main.m for other examples of use.
%
%
% See also: fminsearch, fminspleas
%
%
% Author: John D'Errico
% E-mail: [e-mail]
% Release: 4
% Release date: 7/23/06

% size checks
xsize = size(x0);
x0 = x0(:);
n=length(x0);

if (nargin<3) || isempty(LB)
  LB = repmat(-inf,n,1);
else
  LB = LB(:);
end
if (nargin<4) || isempty(UB)
  UB = repmat(inf,n,1);
else
  UB = UB(:);
end

if (n~=length(LB)) || (n~=length(UB))
  error 'x0 is incompatible in size with either LB or UB.'
end

% set default options if necessary
if (nargin<5) || isempty(options)
  options = optimset('fminsearch');
end

% stuff into a struct to pass around
params.args = varargin;
params.LB = LB;
params.UB = UB;
params.fun = fun;
params.n = n;
% note that the number of parameters may actually vary if 
% a user has chosen to fix one or more parameters
params.xsize = xsize;
params.OutputFcn = [];

% 0 --> unconstrained variable
% 1 --> lower bound only
% 2 --> upper bound only
% 3 --> dual finite bounds
% 4 --> fixed variable
params.BoundClass = zeros(n,1);
for i=1:n
  k = isfinite(LB(i)) + 2*isfinite(UB(i));
  params.BoundClass(i) = k;
  if (k==3) && (LB(i)==UB(i))
    params.BoundClass(i) = 4;
  end
end

% transform starting values into their unconstrained
% surrogates. Check for infeasible starting guesses.
x0u = x0;
k=1;
for i = 1:n
  switch params.BoundClass(i)
    case 1
      % lower bound only
      if x0(i)<=LB(i)
        % infeasible starting value. Use bound.
        x0u(k) = 0;
      else
        x0u(k) = sqrt(x0(i) - LB(i));
      end
      
      % increment k
      k=k+1;
    case 2
      % upper bound only
      if x0(i)>=UB(i)
        % infeasible starting value. use bound.
        x0u(k) = 0;
      else
        x0u(k) = sqrt(UB(i) - x0(i));
      end
      
      % increment k
      k=k+1;
    case 3
      % lower and upper bounds
      if x0(i)<=LB(i)
        % infeasible starting value
        x0u(k) = -pi/2;
      elseif x0(i)>=UB(i)
        % infeasible starting value
        x0u(k) = pi/2;
      else
        x0u(k) = 2*(x0(i) - LB(i))/(UB(i)-LB(i)) - 1;
        % shift by 2*pi to avoid problems at zero in fminsearch
        % otherwise, the initial simplex is vanishingly small
        x0u(k) = 2*pi+asin(max(-1,min(1,x0u(k))));
      end
      
      % increment k
      k=k+1;
    case 0
      % unconstrained variable. x0u(i) is set.
      x0u(k) = x0(i);
      
      % increment k
      k=k+1;
    case 4
      % fixed variable. drop it before fminsearch sees it.
      % k is not incremented for this variable.
  end
  
end
% if any of the unknowns were fixed, then we need to shorten
% x0u now.
if k<=n
  x0u(k:n) = [];
end

% were all the variables fixed?
if isempty(x0u)
  % All variables were fixed. quit immediately, setting the
  % appropriate parameters, then return.
  
  % undo the variable transformations into the original space
  x = xtransform(x0u,params);
  
  % final reshape
  x = reshape(x,xsize);
  
  % stuff fval with the final value
  fval = feval(params.fun,x,params.args{:});
  
  % fminsearchbnd was not called
  exitflag = 0;
  
  output.iterations = 0;
  output.funcCount = 1;
  output.algorithm = 'fminsearch';
  output.message = 'All variables were held fixed by the applied bounds';
  
  % return with no call at all to fminsearch
  return
end

% Check for an outputfcn. If there is any, then substitute my
% own wrapper function.
if ~isempty(options.OutputFcn)
  params.OutputFcn = options.OutputFcn;
  options.OutputFcn = @outfun_wrapper;
end

% now we can call fminsearch, but with our own
% intra-objective function.
[xu,fval,exitflag,output] = fminsearch(@intrafun,x0u,options,params);

% undo the variable transformations into the original space
x = xtransform(xu,params);

% final reshape to make sure the result has the proper shape
x = reshape(x,xsize);

% Use a nested function as the OutputFcn wrapper
  function stop = outfun_wrapper(x,varargin);
    % we need to transform x first
    xtrans = xtransform(x,params);
    
    % then call the user supplied OutputFcn
    stop = params.OutputFcn(xtrans,varargin{1:(end-1)});
    
  end

end % mainline end

% ======================================
% ========= begin subfunctions =========
% ======================================
function fval = intrafun(x,params)
% transform variables, then call original function

% transform
xtrans = xtransform(x,params);

% and call fun
fval = feval(params.fun,reshape(xtrans,params.xsize),params.args{:});

end % sub function intrafun end

% ======================================
function xtrans = xtransform(x,params)
% converts unconstrained variables into their original domains

xtrans = zeros(params.xsize);
% k allows some variables to be fixed, thus dropped from the
% optimization.
k=1;
for i = 1:params.n
  switch params.BoundClass(i)
    case 1
      % lower bound only
      xtrans(i) = params.LB(i) + x(k).^2;
      
      k=k+1;
    case 2
      % upper bound only
      xtrans(i) = params.UB(i) - x(k).^2;
      
      k=k+1;
    case 3
      % lower and upper bounds
      xtrans(i) = (sin(x(k))+1)/2;
      xtrans(i) = xtrans(i)*(params.UB(i) - params.LB(i)) + params.LB(i);
      % just in case of any floating point problems
      xtrans(i) = max(params.LB(i),min(params.UB(i),xtrans(i)));
      
      k=k+1;
    case 4
      % fixed variable, bounds are equal, set it at either bound
      xtrans(i) = params.LB(i);
    case 0
      % unconstrained variable.
      xtrans(i) = x(k);
      
      k=k+1;
  end
end

end % sub function xtransform end

[Gerald Navida]

Yeah you're right. It worked! The tolerances I used might be unrealistic. Thanks

Answer 1

[John D'Errico]

Did it even work using fmincon? I can't test that out, since we lack your data. And lacking your data, I can't test fminsearchbnd on it either.

However, learn to use functions! Learn to pass in the data via a function handle, not using additional arguments to ode45. I'm pretty sure, if that still even works at all, that they will one day phase out that programming style where you pass in extra arguemnts that way. It is not even documented in MATLAB anymore as a valid way to pass arguments around.

Change your call for ode45 to this:

[~,HB_model] = ode45(@(t,x) ...
   model(t,x,b,mu,d,beta,epsilon,k,r1,r2,r3,w,p,q,v),...
   tspan,init,options);

The difference us, this call uses a function handle to pass in the parameters into your model function, and ODE45 never cares about them.

[Gerald Navida]

Yes, it worked for fmincon and so I wonder why it won't work for fminsearchbnd. Here are the data:

t_data = (2010:1:2023);

I_data = [25784 24243 19128 13528 11028 9518 8559 6188 2874 2628 1982 2242 1745 1548];

C_data = [353403 376054 377160 376131 375944 377901 418554 424359 439848 453686 436344 444763 440314 445513];

I've already changed my call for ode45 using what you've suggested, but I still get the same error.

[John D'Errico]

It is not clear what you are doing. Y is a VECTOR of size 14x1. Neither fmincon nor fminsearchbnd can minimize a VECTOR.

Y1 = (HB_model(:,3) - I_data1).^2;
Y2 = (HB_model(:,4) - C_data1).^2;
Y = Y1+Y2;

K>> size(Y)
ans =
    14     1

We see it created in those lines, and that is what you return to the optimizer. As such, fmincon would also fail to work on this code, as you have posted it. Possibly you wanted to minimize the sum of the elements of Y. But I cannot read your mind. I do know however that fmincon will NOT do what you claim it to do. I might guess that you changed things from when you had fmincon working. But we cannot know.

[Gerald Navida]

I actually took the sum of Y1 and Y2, i.e., sum(Y1) and sum(Y2), to return a scalar. I tried to delete it for fminsearchbnd but I realized that it must still be a scalar.

[John D'Errico]

And so whjat is the problem? Partly, it is this insanity:

options = optimset('Display','iter','TolX',1e-15,'TolFun',1e-15,'MaxFunEvals',1.5e+7,'MaxIter',1.5e+8);

150 million iterations? 15 milllion function evals? Of course, it will hit that limit first, since fminsearch does at least as many function evals as it does iterations.

Tolerances at 1e-15?

Let us be serious. Fminsearch related tools will NOT converge as well as gradient based tools. Not as fast. And cranking down on those tolerances will just never converge in a finite amount of time. Change it to something like this:

options = optimset('Display','final','TolX',1e-8,'TolFun',1e-8,'MaxIter',1e5);

And then return Y in the model_fit code as:

Y = sum(Y1) + sum(Y2);

And now it quickly returns a result as:

Exiting: Maximum number of function evaluations has been exceeded
- increase MaxFunEvals option.
Current function value: 2698226867.338823 

xopt =
1.1032      0.40235

Which leaves it moderately happy. It ran out of function evals before it was able to get as close as you were requesting. And that is because the objective function is VERY large. So my requested tolerance of 1e-8 was STILL way too small. It is still smaller than

eps(2698226867.338823)
ans =
4.7684e-07

Which means that I probably should have used TolFun on the order of 1 or so.

You need to understand the diference between tools like fmincon and those based on fminsearch. But just as importantly, you need to think about the tolerances you are setting. Are they meaningful, in context of double precision arithmetic?

When you required TolFun to be 1e-15, that is 8 or 9 powers of 10 smaller then eps of your objective function value! Double precision arithmetic itself fails here.

[Gerald Navida]

Yeah you're right. Those tolerances for the optimizer are unrealistic. Thanks