How do I address "Error using barrier (line 31) Gradient at initial point contains Inf, NaN, or complex values. Fmincon cannot continue."

I have created an output function that I am using with fmincon that provides a gradient as an output. When optimizing, I keep getting the following error:

----------------------------------------------

Error using barrier (line 31)
Gradient at initial point contains Inf, NaN, or complex values. Fmincon cannot continue.

Error in fmincon (line 799)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

Error in loop_gradtest>(parfor body) (line 13)
wts_out(:,i) = fmincon(@(x)utility_gradtest(x,score_param1,
score_param2,rets),x,[],[],[mean(returns_grid_parallel)*12;ones(1,size(returns_grid_parallel,2))],[.05;1],double(limits_min(:,2)),double(limits_max(:,2)),[],fmin_options);

Error in loop_gradtest (line 9)
parfor i=1:10000;

------------------------------------------------------------

The function gradtest is used to all another function, utility_gradtest that is used to calculate maximum (min negative) utility given a vector of weights. In this function, using finite differencing, I calculate a gradient that is then supplied as part of the output for utility_gradtest.

I've used the following input in setting my options:

fmin_options = optimset('Algorithm','interior-point','GradObj','on','MaxFunEvals',10000,'MaxIter',500,'TolCon',10e-16,'TolX',10e-16,'Display','off','UseParallel',true);

Do I need to incorporate error handling capability to account for when the gradient contains an undefined point?

[Nick Hobbs]

What is your gradient function and initial point? From the error message, it is likely the case that the gradient is returning a value that fmincon cannot work with. If you can determine the value the gradient is producing for fmincon, it may be possible to handle the issue by fixing the gradient before it is returned in the function.

[Walter Roberson]

We will probably need to see utility_gradtest() and anything it calls, and knowing something about the initial x and the other variables being passed into utility_gradtest would help.

[J. Womack]

Thank you both for your responses.

My gradient function is the following:

function [gradient] = utility_gradient(x,score_param1, score_param2,returns_grid_parallel)

    %Set values for parameters
    a = score_param2; %loss aversion parameter
    p = score_param1; %risk aversion parameter

    port_annual_return = (mean(returns_grid_parallel).*12)'.*x;
    port_annual_std = (cov(returns_grid_parallel).*12)*x./sqrt(x'*(cov(returns_grid_parallel).*12)*x);

es_low_z = -(1./sqrt(2*pi).*exp(-((0-(mean(returns_grid_parallel).*12)'.*x)./((cov(returns_grid_parallel).*12)*x./sqrt(x'*(cov(returns_grid_parallel).*12)*x))).^2./2))./ ...
              (.5*(erf(((0-(mean(returns_grid_parallel).*12)'.*x)./((cov(returns_grid_parallel).*12)*x./sqrt(x'*(cov(returns_grid_parallel).*12)*x)))./sqrt(2))+1));
  es_high_z = (1./sqrt(2*pi).*exp(-((0-(mean(returns_grid_parallel).*12)'.*x)./((cov(returns_grid_parallel).*12)*x./sqrt(x'*(cov(returns_grid_parallel).*12)*x))).^2./2))./ ...
               (1-.5*(erf(((0-(mean(returns_grid_parallel).*12)'.*x)./((cov(returns_grid_parallel).*12)*x./sqrt(x'*(cov(returns_grid_parallel).*12)*x)))./sqrt(2))+1));                                                                          

    gradient = -a.*(.5*(erf(((0-(mean(returns_grid_parallel).*12)'.*x)./((cov(returns_grid_parallel).*12)*x./sqrt(x'*(cov(returns_grid_parallel).*12)*x)))./sqrt(2))+1)).*exp(-p.*(1+port_annual_return+es_low_z.*port_annual_std))+ ...
                -(1-(.5*(erf(((0-(mean(returns_grid_parallel).*12)'.*x)./((cov(returns_grid_parallel).*12)*x./sqrt(x'*(cov(returns_grid_parallel).*12)*x)))./sqrt(2))+1))).*exp(-p.*(1+port_annual_return+es_high_z.*port_annual_std));

    end

The function itself just calculates a scalar value for a utility function; I use the current weight vector, then add a tiny increment, and calculate the gradient numerically. (When I clean up the flow for production, I will eliminate running the same function twice as can be seen below). The code to calculate the scalar function is below:

    function [utility_values,grad] = utility_gradtest(x,score_param1, score_param2,returns_grid_parallel)                              

%Set values for parameters
a = score_param2; %loss aversion parameter
p = score_param1; %risk aversion parameter

port_annual_return = mean(returns_grid_parallel*x).*12;
port_annual_std = std(returns_grid_parallel*x).*sqrt(12);

p_low = normcdf((0-port_annual_return)./port_annual_std,0,1);
es_low_z = -normpdf((0-port_annual_return)./port_annual_std,0,1)./p_low;
es_high_z =pdf(9)-normpdf((0-port_annual_return)./port_annual_std,0,1))./(1-p_low);

utility_values = mean(-(-a.*p_low.*exp(-p.*(1+port_annual_return+es_low_z.*port_annual_std))+ ...
-(1-p_low).*exp(-p.*(1+port_annual_return+es_high_z.*port_annual_std))));

grad = utility_gradient(x+repmat(.0000001,size(x)),score_param1, score_param2,returns_grid_parallel)- ...
                                         utility_gradient(x,score_param1, score_param2,returns_grid_parallel);

end

My initial value for the weight vector (x input) is a vector that represents equal weights:

repmat(1/n,n,1)

I am going to work on debugging but want to make sure my function isn't bombing out when things change.

Any thoughts you can share on either error handling or approach are greatly appreciated!

[Walter Roberson]

At the command line command

dbstop if naninf

and run your program. If it stops in the debugger then you encountered a nan or inf and you will be able to explore the values in the debugger to see what is happening. If it does not stop in the debugger but gives you the barrier problem anyhow then you encountered a complex value and you can start debugging around that.

I would point out as a matter of efficiency that your code calculates some subexpressions multiple times. You should calculate them and store them in variables and use them as needed. In particular I notice mean(returns_grid_parallel).*12)'.*x) and cov(returns_grid_parallel) being repeated

My pure speculation is that you are ending up with a 0 / 0 situation.

[J. Womack]

Walter - thank you for the recommendation. I followed your advice and just focused on capturing the gradient from the parfor loop and found the problem. I've added in code to address NaNs and Infs and it's working.

I will run through the code and clean it up to try to capture additional speed and appreciate you pointing out the repetition. Where it's within the same function, I can make the change.

Thank you again for your help!

[Walter Roberson]

Did you have a 0/0 or was it an inf*0 ?

Answer 1

[J. Womack]

This question was answered.