fmincon freezing (no NaN or Inf is involved)

Hello,

Here is a short description about the problem I'm facing and what I know so far.

The problem is associated with a master minimizing code implementing Simulating Method of Moments. This code (and its variants) uses a solver (e.g. fmincon) to minimize the error between a set of variables and the outputs of a core function. This core is MEXed (turned to c based matlab functions for higher speeds) and is used in other applications with no problem at running. However, by using this code in this application, the code freezes with no further improvement, while its status is running.

The place at which the code freezes is in a calculation and RAM intensive loop inside the core function, in a sub function that involves another minimization (this minimization is done with fmincon to be compatible wit MEX).

It seems that on average the codes pause after an average amount of time is passed. This is because by reducing the code scale, I see that the average iteration that code is stopped, increases relative to the scale down.

Also when the code itself (without using MEX) is run the code pauses at the very first iteration.

Also there is no Inf or NaN received from the function that the inner minimization is minimizing and freezes.

Here is a sketch of my code:

main function

1st optimization(core code)

core function

2nd optimization (where code freezes)

and this is the 2nd optimization and the function I try to optimize:

function [G_Max,b_pol_offg,k_pol_offg,n_pol_offg] = ...
    VFI_uni_offg1(max_iter_vfi,pr,G0,G_Max,div,b_pol_or,k_pol_or,n_pol_or,...
    b,k,nl,y,nb,nk,nn,ny,p_exit,k_2,b_2,prob_1,tol_vfi_offg,mode,ZZ)
%%%%%%%%%%%%%%%%%%%% off grid solution for uniform shock
res         = ones(1,max_iter_vfi);

G_Max(G_Max<0) = -1;

vfi_time    = zeros(1,'int16');

A           = [1,-pr.lambda,0];
B           = 0;
lb_spline   = [0,0,0];
lb_linear   = [1*min(b),1*min(k),1*min(nl)];
lb          = lb_linear*(mode==1) + lb_spline*(mode~=1);
ub_linear   = [1*max(b),1*max(k),1*max(nl)];
ub_spline   = [1*max(b),1.4*max(k),1.2*max(nl)];
ub          = ub_linear*(mode==1) + ub_spline*(mode~=1);
b_pol_offg  = b_pol_or;
k_pol_offg  = k_pol_or;
n_pol_offg  = n_pol_or;
checkpoints = linspace(0,.6,ZZ);


% ZZ = 5;

options     = optimoptions(@fmincon,'display','off','Algorithm','sqp');
tic
for ii=1:max_iter_vfi
    G0(G0<=0) = -1;
    G     = div+pr.beta.*((1-p_exit).*kron(ones(1,nn),kron((prob_1*G0'),ones(nb*nk,1)))...
        +(p_exit)*(1-pr.zeta).*(kron(ones(nb*nk*ny,1),k_2)-kron(ones(nb*nk*ny,1),b_2)));
    %    G(:,nb*nk+1:nb*nk*ns)   = pr.beta.*kron((prob_1*G0'),ones(nb*nk,1));
    %    G                       = div + G;
    fprintf('%5i \n',ii);
    for jj=1:nb*nk*ny
        x0  = [b_pol_offg(jj)+0.001,k_pol_offg(jj)*(k_pol_offg(jj)>0)+1*(k_pol_offg(jj)==0)...
            ,n_pol_offg(jj)*(n_pol_offg(jj)>0)+2*(n_pol_offg(jj)==0)];
        g   = G(jj,:);
        %         g(g==-inf)  = -1;
        g_temp      = reshape(g,nb,nk,nn);
        
        x00 = x0;
        G_Max_temp          = zeros(1,ZZ);
        b_pol_offg_temp     = zeros(1,ZZ);
        k_pol_offg_temp     = zeros(1,ZZ);
        n_pol_offg_temp     = zeros(1,ZZ);
        
        %         fprintf('%5.2f , ',jj);
        
        for zz=1:1:ZZ
            [x,fval,exitflag]    = fmincon(@(x)Vfunction(b,k,nl,g_temp,x(1),x(2),x(3),mode),x0,A,B,...
                [],[],lb,ub,@(x)nonlinconst(b,k,y,nb,nk,pr,x(1),x(2),x(3),jj),options);
            if (exitflag>=1 && -fval>=0 && -fval>=G_Max(jj))
                x0 = x00*(.7+checkpoints(zz));
                G_Max_temp(zz)       = -fval;
                b_pol_offg_temp(zz)  = x(1);
                k_pol_offg_temp(zz)  = x(2);
                n_pol_offg_temp(zz)  = x(3);
            elseif (exitflag<1 || -fval<0 || -fval<G_Max(jj))
                x0 = x00*(.7+checkpoints(zz));
                G_Max_temp(zz)       = G_Max(jj);
                b_pol_offg_temp(zz)  = b_pol_offg(jj);
                k_pol_offg_temp(zz)  = k_pol_offg(jj);
                n_pol_offg_temp(zz)  = n_pol_offg(jj);
            end
            %         if mod(jj,100)==0
            %             fprintf('%5i',jj);
            %         end
            %          fprintf('%5.4f , ',G_Max_temp(zz));
        end
        [G_Max(jj),max_ind]     = max(G_Max_temp);
        b_pol_offg(jj)          = b_pol_offg_temp(max_ind);
        k_pol_offg(jj)          = k_pol_offg_temp(max_ind);
        n_pol_offg(jj)          = n_pol_offg_temp(max_ind);
    end
    
    vfi_time = ii;
    if max(max(abs((reshape(G_Max,nb*nk,ny))-G0)))<tol_vfi_offg
        break
    end
    
    res(ii)=max(max(abs((reshape(G_Max,nb*nk,ny))-G0)));
    
    G0 = reshape(G_Max,nb*nk,ny);
    
end
time2 = toc;
fprintf('=========================================================\n');
fprintf('Offgrid Value Function converged with accuracy %6.8f in %5i iterations in %6.8f second \n',...
    res(vfi_time-1),vfi_time-1,time2);
fprintf('=========================================================\n');
end

where the code freezes at

```

[x,fval,exitflag]    = fmincon(@(x)Vfunction(b,k,nl,g_temp,x(1),x(2),x(3),mode),x0,A,B,...
    [],[],lb,ub,@(x)nonlinconst(b,k,y,nb,nk,pr,x(1),x(2),x(3),jj),options);

```

and vfunction is :

```

function [value] = Vfunction(b,k,nl,g_temp,b_var,k_var,n_var,mode)
if mode==1 
    vq      = interp3(k,b,nl,g_temp,k_var,b_var,n_var,'linear',0);
else
    vq      = interp3(k,b,nl,g_temp,k_var,b_var,n_var,'spline');
end

value   = -vq;
if isnan(value)
    value = 0;
end
if isinf(value)
    value = -100;
end
end

```

[Walter Roberson]

Could you explain more about how you know it is freezing, as compared to just taking a long time? There are some parts of fmincon that take time proportional to the square of the number of variables.

[Matt J]

If you think that fmincon is the part that is freezing, it would be advisable to dump all input variables to fmincon to a file in the iteration where the freeze happens. Then, you can re-run just the fmincon optimization in Matlab MCode in isolation from everything else. Moreoever, if you can reproduce the freeze in MCode, you can share that isolated portion with us, so that we can investigate it with you.

[Alireza Azampour]

Thank you for your response. The problem is that it doesn't freeze at a certain point and I can't record the exact point fmincon freezes, because I think I have to print down every iteration's input which makes the log file unmanageable. It would be great if you could suggest a way to pin down the inputs that freeze fmincon.

[Matt J]

Can't you tell in which iteration it freezes by printing the current iteration number to the screen? When the numbers stop increasing, you'll know that's the iteration where the freeze occurs.

because I think I have to print down every iteration's input which makes the log file unmanageable

No, don't do that. Once you know which iteration the freeze occurs in,r ecompile and rerun with lines inserted to dump only in that specific iteration.

[Alireza Azampour]

the problem is that the iteration at which the code freezes is totally random.

[Matt J]

You can give the random number generator a fixed seed so that two successive executions proceed the same.

[Alireza Azampour]

No I mean the process of the code is deterministic but each time I run the code it freezes at a different point

[Walter Roberson]

Could you explain more about how you know it is freezing, as compared to just taking a long time?

Note: when you have non-empty upper bounds or lower bounds (which is true in your case) then fmincon() uses random numbers as part of deciding which direction to move within the bounds constraints.

[Alireza Azampour]

The reason I think it freezes is that normally each iteration of the main code (which consists of VFI_uni_offg1 functions I have posted) takes about 10 minutes, but when the "freezing" happens the code is stucked at the same iteration in the line I mentioned for 80+ hours (the longest time I have tested on a server). And unfortunately when I try to pause the code on matlab on my own pc, when it's stucked, the process of pausing is also not responsive.

[Alireza Azampour]

I have understood that as you have said fmincon is not freezed but it evaluates the function many time. The number of evaluations are much more than the default number of MaxFunctionEvaluations. However, if I understand correctly if the nonlinear (or linear) constraints are not satisfied it is not counted in F-count. Is it true and if so what do you suggest cause it seems fmincon gets stuck in an area that violates the constrains and doesn't get out.

[Matt J]

@Walter Roberson

Note: when you have non-empty upper bounds or lower bounds (which is true in your case) then fmincon() uses random numbers as part of deciding which direction to move within the bounds constraints.

Does it? I don't see that anywhere in the documentation.

@Alireza Azampour

It sounds less and less tom me as though fmincon is the culprit. Is it possible you are running it over a network, or some other medium that could intermittently lose connectivity with its resources?

[Alireza Azampour]

I have run the code both on a hpc server and my pc. The "freezing" happens on both of them.

[Walter Roberson]

@Matt J

You are right, I got confused over which function was which. I had discovered the issue while working on some user code that set everything up for fmincon, but then called ga instead, and the rest of their code had been using fmincon() .

The actual random generation was in boxdirections(), called by mutationadaptfeasible() which looks to be the default mutation function for ga() and related functions. The code does some odd stuff with matrices that appear to be up to (2 * number of variables) on a side [not sure of the actual upper bound, typically number-of-variables on a side but I saw half again that large sometimes.] On large numbers of variables (over 27000 for the user!) that phase turned out to be the slowest part, and it gave me the impression that the code was "freezing" until I dug down into it.

[Alireza Azampour]

Hi

I wanted to elaborate more on the question. I have managed to pin point the exact poit that the code freezes. It indeed fezees. Here are the last instances that fmincon gets from the goal and constrain functions.

Vfunction O is -0.38535,

nonlinc O1,O2 is -1.11571, -0.21061,

Vfunction O is -0.38534,

nonlinc O1,O2 is -1.11568, -0.21038,

Vfunction O is -0.38533,

nonlinc O1,O2 is -1.11566, -0.21022,

Vfunction O is -0.38532,

nonlinc O1,O2 is -1.11565, -0.21011,

Vfunction O is -0.38532,

nonlinc O1,O2 is -1.11564, -0.21003,

Vfunction O is -0.38531,

nonlinc O1,O2 is -1.11563, -0.20998,

Vfunction O is -0.38531,

nonlinc O1,O2 is -1.11563, -0.20994,

Vfunction O is -0.38531,

nonlinc O1,O2 is -1.11562, -0.20991,

In which O1 and O2 are the nonlinear constraint function's output and O is the goal function's output.

I think the problem is with the little derivation in the goal function's output which might make the Jacobian and/or Hessian unreachable for the fmincon algorithm.

If it is true what can I do?

If not what would you suggest the problem is?

Answer 1

[Matt J]

You might try getting rid of these lines in Vfunction()

value   = -vq;
if isnan(value)
    value = 0;
end
if isinf(value)
    value = -100;
end

Because you are using the sqp algorithm, you are losing one of the advantages of the algorithm by hiding the NanNs and Infs. The sqp algorithm is supposed to be able to backtrack from NaNs and Infs when it encounters them and search elsewhere.

Because your algorithm has stochastic behavior when it shouldn't ( i.e., freezing at random, non-reproducible times) , I also think you need to trace the origin of that stochasticity. My recommendation is that you save the sequence of iteration vectors [x1,x2,x3] and make sure that they are following the same sequence every time you run. You can save them using the OutputFcn option, as illustrated in this example,

https://www.mathworks.com/help/matlab/math/output-functions.html#bsgpq6q-31

[Alireza Azampour]

The stochastic nature of the error was because of another optimizier at work in a higher level of the code which I have now put off. I have managed to get the error at a deterministic point and iteration.

I did the omit the the parts you told but it still freezes at the exact spot I told you about (where the goal function's output has very little variation).

[Matt J]

If you can now run the code deterministically, you should be able to isolate the input parameters that cause the freeze, as in my earlier advice. You should provide those to us so that we can run that particular iteration in isolation and examine it with you.

[Alireza Azampour]

Hello again

I've done as you told

here is a function that calls fmincon with the same inputs that makes MATLAB freeze in the main code

function [] = debug(b,k,nl,g_temp,mode,check,x0,A,B,lb,ub,nb,nk,pr,jj,y)
options     = optimoptions(@fmincon,'display','off','Algorithm','sqp');
[x,fval,exitflag]    = fmincon(@(x)Vfunction(b,k,nl,g_temp,x(1),x(2),x(3),mode,check),x0,A,B,...
    [],[],lb,ub,@(x)nonlinconst(b,k,y,nb,nk,pr,x(1),x(2),x(3),jj,check),options);
end

I have also made the MEX function for better compatibility between results.

Here is also the functions used in the optimization process.

function [value] = Vfunction(b,k,nl,g_temp,b_var,k_var,n_var,mode,check)
if mode==1 
    vq      = interp3(k,b,nl,g_temp,k_var,b_var,n_var,'linear');
else
    vq      = interp3(k,b,nl,g_temp,k_var,b_var,n_var,'spline');
end

value   = -vq;
% if isnan(value)
%     fprintf('NaaaN \n');
%     value = 0;
% end
% if isinf(value)
%     fprintf('Iiinf \n');
%     value = -100;
% end
if check == 1
    fprintf('Vfunction I/O is %3.5f,%3.5f,%3.5f,%3.5f, \n',k_var,b_var,n_var,value);
end
end

function [c,trash] =nonlinconst(b,k,y,nb,nk,pr,b_var,k_var,n_var,jj,check)
c       = zeros(1,2);
trash   = [];
k_st    = k(mod(floor((jj-1)/nb),nk)+1);
b_st    = b(mod((jj-1),nb)+1);
y_st    = y(floor((jj-1)/(nb*nk))+1);

del_k   = (1-pr.sig)*k_st-k_var;
del_k   = del_k*(1-pr.nu*(del_k>0));

psi     = pr.tau*abs(b_st -b_var);

shd     = -(del_k -(1+pr.r)*b_st +b_var-psi -n_var*pr.w);

if (n_var<0)||(k_st<0)
    c(1) = 1;
end

zf      = y_st*power(complex(n_var),pr.alpha)*power(complex(k_st),pr.ga);

% if ~isreal(zf)
%     fprintf('nvar is %7.4f k_st is %7.4f real part is %7.4f img is %7.4f',n_var,k_st,real(zf),imag(zf));
%     c(1) = 1;
%     return
% else
    zf = real(zf);
% end

b24ac   = (zf+pr.zeta*(k_var-b_var)).^2 - 4*zf*shd;
c(2)      = -b24ac;
if c(2)>0
    c(1)  = 1;
    return;
end
ps_shd  = 0.5*(zf +pr.zeta*(k_var-b_var) -b24ac^(0.5))*(shd>0) + shd*(shd<0);
c(1)      = -(0.5*zf-ps_shd);
if isnan(c(1))||isnan(c(2))
    fprintf('NaaaaaN');
end
if isinf(c(1))||isinf(c(2))
    fprintf('Iiiiiinf \n');
end
if check==1
    fprintf('nonlinc O1,O2 is %3.5f, %3.5f, \n',c(1),c(2));
end
end

I'm also upploading the set of inputs that makes the code freeze in the main code.

problem here is that although the first few steps of fmincon are the same in the main code and debug code, the they diverge and I don't know why this happens for the same input.

p.s. the MEX code is produced with MATLAB 2020b

[Matt J]

problem here is that although the first few steps of fmincon are the same in the main code and debug code, the they diverge and I don't know why this happens for the same input.

That's not supposed to be a problem because the x0 that I requested you isolate and include in debug_input.mat is the vector that fmincon reaches right before the freeze. In other words, I should be able to run only a single iteration of fmincon starting at x0 to make the freeze happen.

Answer 2

[Alireza Azampour]

Here is the last input that fmincon freezes on. I think there might be a problem ofcourse cause this input is stored in doubles but I don't know with what precision fmincon works.

[Alireza Azampour]

Here is the true value of the inputs to the goal function in the last iteration

[Alireza Azampour]

Sorry I am somehow confused by by how matlab passes the variables between functions. But the vfunction's inputs are double (however when I print the them it has 54 decimal precision)

[Walter Roberson]

Notice:

PI = pi;
PI_STRING = sprintf('%.999g', PI)

PI_STRING = '3.141592653589793115997963468544185161590576171875'

length(PI_STRING)

ans = 50

Is it really 48 digits (after the decimal place) precision? No! -- that is just the decimal expansion of something that is binary.

num2hex(PI)

ans = '400921fb54442d18'

next2PI = typecast(typecast(PI, 'uint64')+1, 'double')

next2PI = 3.1416

next2PI_string = sprintf('%.999g', next2PI)

next2PI_string = '3.141592653589793560087173318606801331043243408203125'

length(next2PI_string)

ans = 53

num2hex(next2PI)

ans = '400921fb54442d19'

fprintf('%.999g\n', next2PI - PI)

4.44089209850062616169452667236328125e-16

2^(-51)

ans = 4.4409e-16

next2PI is the next representable number adjacent to pi -- you can see in the num2hex displays that the 8 increments to 9 and there is no other change. Does that somehow increase the number of digits of precision by 3 (according to length())? Not really -- that's just how it comes out if you expand the binary to decimal. The actual difference between the numbers is about 4.4E-16, 2^-51 -- you only have 51 bits of fraction resolution at that point in floating point numbers.

[Alireza Azampour]

Ohh many thanks for clarifying this