YALMIP formulation of SDP for non-SDP local nonlinear solvers

[Mark L. Stone]

In the latest develop version, how does YALMIP formulate (nonlinear) SDP for solution by non-SDP local nonlinear solvers, such as FMINCON and KNITRO?

Can we say anything about the relation between (claimed) local minima found by the local nonlinear solver and local minima of the original SDP? Are additional stationary points introduced?  Are the local nonlinear solver's KKT optimality criteria valid for the original SDP if the SDP has multiple zero eigenvalues at the optimum?

I don't necessarily expect you to have answers to all these questions.  Thanks.

[Johan Löfberg]

It's a brute strategy working with sub-gradients of all eigenvalues, with a bunch of hacks to deal with the fact that computed eigenvalues are unordered etc. No barriers or extra variables or constraints anything, the central idea is that YALMIP basically does nothing to the model, and for the solver it should be like any other nonlinear constraint which returns infeasibility and jacobian. No theory, 60% of the time it works every time. Fails miserably on some trivial examples and works perfectly on some hard problems. I will hopefully get around to polish it a bit this semester and maybe write about it. The name of the solver is "Conan the slayer" (CONic And Nonlinear the Semidefinite Layer).

[Mark L. Stone]

Ah yes, you get the dynamic aspect via callbacks. I have also succeeded in quickly solving some problems with it calling KNITROL as solver, which over a range of problems, seems to work better than FMINCON. Perhaps with some polishing, you could reduce the number of failures.  Maybe pair up with a theoretician to help you do that, and then publish a nice paper?

[Mark L. Stone]

Perhaps use Conan the slayer (with specified nonlinear  solver), or some variant of it, as an upper solver option for BMIBNB? I have seen problems for which Conan the slayer produced quick good solutions, and BMIBNB wasn't coming up with any decent incumbent (with KNITRO or FMINCON as upper solver).

[Johan Löfberg]

Slayer is already completely integrated so it is the default solver/layer for bmibnb on nonlinear conic problems (it was actually the initial reason for developing it)

[Johan Löfberg]

...meaning that any nonlinear solver can be used as an upper solver for SDPs in bmibnb (although I think only fmincon and knitro are possible right now), and slayer is automatically used

[Mark L. Stone]

Here is why I thought Slayer was not used by bmibnb.

Configuration: WIN64: MATLAB R2019b. YALMIP develop 20211228. KNITRO 12.4.0. The input data, Hessian and gradient, are in attached mat file.

I realize this may not be reproducible if you don't have KNITRO (12.4) available to run. But perhaps you can figure something out from the error message. (BTW, this is just a totally made up test problem, having no relation to anything in the real world. This Hessian is indefinite.).

This example shows both a bmibnb error, as well as knitro as solver producing a better objective (not globally optimal, though) than first 2 incumbents from bmibnb with knitro as upper solver.

load 'BMI_Hessian_gradient_13by13_12-18-21.mat'
x=sdpvar(13,1);
Obj=0.5*x'*Hessian*x+gradient'*x;
BMI_matrix = sdpvar(8,8);
Con = [BMI_matrix(1,1:8) == (x(1:8).*x(8:-1:1))'];
Con = [Con,BMI_matrix(2,2:8) == x(4:10)'];
Con = [Con,BMI_matrix(3,3:8) == 2*(x(5:10).*x(6:-1:1))'];
Con = [Con,BMI_matrix(4,4:8) == [3*x(12)*x(13) -(x(1:4).*x(2:5))']];
Con = [Con,BMI_matrix(5,5:8) == .25*x(12:-1:9)'];
Con = [Con,BMI_matrix(6,6:8) == x(10:-1:8)'];
Con = [Con,BMI_matrix(7,7:8) == [x(7)*x(13) -x(2)*x(6)]];
Con = [Con,BMI_matrix(8,8) == x(4)*x(9)+4*x(5)^2-x(3)*x(4)];
optimize([Con,BMI_matrix>=0,sum(x)==26,-1<=x<=100],Obj,sdpsettings('solver','bmibnb','bmibnb.uppersolver','knitro','debug',1))

* Starting YALMIP global branch & bound.
* Upper solver     : knitro
* Lower solver     : MOSEK
* LP solver        : GUROBI
* -Extracting bounds from model
* -Perfoming root-node bound propagation
* -Calling upper solver (found a solution, objective 352.501)
* -Branch-variables : 13
* -More root-node bound-propagation
* -Performing LP-based bound-propagation
* -And some more root-node bound-propagation
* -Shifted nonconvex QP objective using SDP
* Starting the b&b process
 Node       Upper       Gap(%)       Lower     Open   Time
    1 :   2.57275E+02   199.70   -1.08990E+03    2    16s  Improved solution found by upper solver  
    2 :   2.57275E+02   199.70   -1.08990E+03    3    20s    
    3 :  -2.67629E+00   198.61   -1.05547E+03    4    34s  Improved solution found by upper solver  
    4 :  -2.67629E+00   198.61   -1.05547E+03    5    55s    
    5 :  -2.67629E+00   198.46   -9.48739E+02    6    59s    
    6 :  -2.67629E+00   198.46   -9.48739E+02    7    92s    
    7 :  -8.71489E+01   165.72   -9.39422E+02    8   106s  Improved solution found by upper solver  
    8 :  -8.71489E+01   165.72   -9.39422E+02    9   117s    
    9 :  -8.71489E+01   163.36   -8.73277E+02   10   120s    
   10 :  -8.71489E+01   163.36   -8.73277E+02   11   147s    
   11 :  -8.71489E+01   163.06   -8.65454E+02   12   152s    
   12 :  -8.71489E+01   163.06   -8.65454E+02   13   180s    
   13 :  -6.02521E+02    28.73   -8.05005E+02    6   189s  Improved solution found by upper solver | Pruned stack based on new upper bound  
   14 :  -6.02521E+02    28.73   -8.05005E+02    5   190s  Poor lower bound  
   15 :  -6.02521E+02    28.33   -8.01754E+02    6   198s  Improved solution found by upper solver  
   16 :  -6.02521E+02    28.33   -8.01754E+02    5   198s  Terminated in bound propagation  
   17 :  -6.02521E+02    24.59   -7.71709E+02    4   199s  Terminated in bound propagation  
   18 :  -6.02521E+02    24.59   -7.71709E+02    5   203s    
   19 :  -6.02521E+02    24.43   -7.70469E+02    6   210s    
   20 :  -6.02521E+02    24.43   -7.70469E+02    7   226s  Numerical problems in lower solver  
ans =
     1
'output' requires Sensor Fusion and Tracking Toolbox.
Error in solvelower (line 314)
            x(~removethese)=output.Primal;
Error in bmibnb_branch_and_bound>solvelower_safelayer (line 1107)
[output,cost,p,timing] = solvelower(p,options,lowersolver,x_min,upper,timing);
Error in bmibnb_branch_and_bound (line 346)
            [output,cost,p,timing] = solvelower_safelayer(p,options,lowersolver,x_min,upper,timing);
Error in bmibnb (line 515)
        [x_min,solved_nodes,lower,upper,lower_hist,upper_hist,solution_hist,timing,counter,problem] =
        bmibnb_branch_and_bound(p,x_min,upper,timing,solution_hist);
Error in solvesdp (line 368)
    eval(['output = ' solver.call '(interfacedata);']);
Error in optimize (line 31)
[varargout{1:nargout}] = solvesdp(varargin{:});

*********************************************
Note 2 things:

1)  'solver', 'knitro' produces
...  
 163       216  -7.309484e+001  6.284e-006  9.177e+000  0.000e+000        6

EXIT: Primal feasible solution estimate cannot be improved; desired accuracy
      in dual feasibility could not be achieved.

Final Statistics
----------------
Final objective value               = -7.30948381656738e+001
Final feasibility error (abs / rel) =   6.28e-006 / 2.73e-007
Final optimality error  (abs / rel) =   9.18e+000 / 1.28e-001

which is a better solution than produced by bmibnb with knitro as upper solver, until bmibnb iteration 7, after having already produced 2 solutions with upper solver which have worse objective than 'solver.','knitro'. Would bmibnb "accept" a knitro upper solver solution terminating with "Primal feasible solution estimate cannot be improved; desired accuracy in dual feasibility could not be achieved."? I think this corresponds to KNITRO return flag = -102 ("Current feasible solution estimate cannot be improved.").  "The returned exit flags will correspond with Knitro’s return code, rather than matching fmincon’s exit flags." Or is final feasibility error too high for bmibnb to accept?  Or maybe bmibnb is just feeding Slayer/knitro a "lousier" problem to solve than is provided directly to 'solver','knitro'?

2) bmibnb ends with error message which looks like it fell into a catch of try/catch, producing  ans = 1. I have observed this error on several runs of bmibnb solving BMIs, when using knitro as upper solver. I have not encountered it in a smaller number of runs using fmincon as upper solver, although in many of those, no incumbent was ever produced.  This error may manifest itself in different places depending on which YALMIP develop version (2011125, 20211215, or 20211228) was used, although that perhaps could be due to a "random" perturbation almost anywhere in YALMIP, changing the luck of the draw on whether/where the exact conditions to trigger the error manifest.

BTW, fmincon converges to infeasible on this, and when used as upper solver, no incumbent is produced  after 2000 iterations of bmibnb (on YALMIP develop 20211125).

The actual global optimal objective value is -642.8, which I have verified by other means.

[Johan Löfberg]

Unfortunately I think my knitro setup is broken.

It can easily be that a direct call to a nonlinear solver returns a better solution than what a particular call in bmibnb does on a nonconvex problem (since bmibnb massages the problem a lot to tighten bounds etc). This includes a call in the root, as the problem has been massaged already there.

The fact that a 1 is printed indicates it happens where I've debugged some issue (as adding display 1 is a tell-sign of debugging for me)...

I'll add this to my test-suite

[Johan Löfberg]

The issue was triggered with fmincon too with the setting below. Fixed now in develop.

-642.88 indeed :-)

optimize([Con,BMI_matrix>=0,sum(x)==26,-1<=x<=100],Obj,sdpsettings('solver','bmibnb','bmibnb.uppersolver','fmincon','debug',1,'fmincon.alg','sqp'))

* Starting YALMIP global branch & bound.

* Upper solver     : fmincon

* Lower solver     : MOSEK
* LP solver        : GUROBI
* -Extracting bounds from model
* -Perfoming root-node bound propagation

* -Calling upper solver (no solution found)

* -Branch-variables : 13
* -More root-node bound-propagation
* -Performing LP-based bound-propagation
* -And some more root-node bound-propagation
* -Shifted nonconvex QP objective using SDP
* Starting the b&b process
 Node       Upper       Gap(%)       Lower     Open   Time

    1 :   4.43077E+02   199.74   -1.10415E+03    2    24s  Solution found by upper solver  
    2 :  -6.42881E+02    52.75   -1.10415E+03    3    26s  Improved solution found by upper solver  
    3 :  -6.42881E+02    40.66   -9.71518E+02    4    39s    
    4 :  -6.42881E+02    40.66   -9.71518E+02    3    40s  Poor lower bound  
    5 :  -6.42881E+02    29.38   -8.64646E+02    4    54s    
    6 :  -6.42881E+02    29.38   -8.64646E+02    5    56s    
    7 :  -6.42881E+02    20.88   -7.92982E+02    6    75s    
    8 :  -6.42881E+02    20.88   -7.92982E+02    5    76s  Poor lower bound  
    9 :  -6.42881E+02    17.77   -7.68484E+02    6    94s    
   10 :  -6.42881E+02    17.77   -7.68484E+02    7   105s  Numerical problems in lower solver  
   11 :  -6.42881E+02    17.77   -7.68484E+02    6   105s  Terminated in bound propagation  
   12 :  -6.42881E+02    17.77   -7.68484E+02    5   107s  Infeasible node in lower solver  
   13 :  -6.42881E+02    17.13   -7.63496E+02    6   108s    
   14 :  -6.42881E+02    17.13   -7.63496E+02    7   109s    
   15 :  -6.42881E+02    15.12   -7.48203E+02    8   127s    
   16 :  -6.42881E+02    15.12   -7.48203E+02    7   127s  Terminated in bound propagation  
   17 :  -6.42881E+02    13.27   -7.34362E+02    6   128s  Terminated in bound propagation  
   18 :  -6.42881E+02    13.27   -7.34362E+02    5   128s  Infeasible node in lower solver  
   19 :  -6.42881E+02    11.73   -7.23146E+02    6   140s    
   20 :  -6.42881E+02    11.73   -7.23146E+02    7   156s    
   21 :  -6.42881E+02    11.69   -7.22805E+02    8   158s    
   22 :  -6.42881E+02    11.69   -7.22805E+02    7   158s  Poor lower bound  
   23 :  -6.42881E+02     5.95   -6.82388E+02    8   159s    
   24 :  -6.42881E+02     5.95   -6.82388E+02    7   159s  Infeasible node in lower solver  
   25 :  -6.42881E+02     4.58   -6.73039E+02    8   161s    
   26 :  -6.42881E+02     4.58   -6.73039E+02    9   162s    
   27 :  -6.42881E+02     4.43   -6.72024E+02    8   163s  Infeasible node in lower solver  
   28 :  -6.42881E+02     4.43   -6.72024E+02    7   163s  Poor lower bound  
   29 :  -6.42881E+02     2.52   -6.59320E+02    6   164s  Poor lower bound  
   30 :  -6.42881E+02     2.52   -6.59320E+02    5   164s  Infeasible node in lower solver  
   31 :  -6.42881E+02     1.84   -6.54844E+02    6   182s    
   32 :  -6.42881E+02     1.84   -6.54844E+02    5   183s  Poor lower bound  
   33 :  -6.42881E+02     1.75   -6.54245E+02    4   183s  Poor lower bound  
   34 :  -6.42881E+02     1.75   -6.54245E+02    3   184s  Infeasible node in lower solver  
   35 :  -6.42881E+02     0.95   -6.49000E+02    2   185s  Poor lower bound  
* Finished.  Cost: -642.8809 (lower bound: -648.9996, relative gap 0.94578%)
* Termination with relative gap satisfied
* Timing: 90% spent in upper solver (18 problems solved)
*         1% spent in lower solver (58 problems solved)
*         9% spent in LP-based domain reduction (2384 problems solved)
*         1% spent in upper heuristics (2372 candidates tried)

[Johan Löfberg]

this setting runs faster

optimize([Con,BMI_matrix>=0,sum(x)==26,-1<=x<=100],Obj,sdpsettings('solver','bmibnb','bmibnb.uppersolver','fmincon','debug',1,'fmincon.alg','interior-point','fmincon.InitBarrierParam',1e-5))

(that change is something I've seen in experiment that it typically is a way to get fmincons default interior-point algo to work better on the SDP reformulations)

[Mark L. Stone]

Thanks for the fix and the fmincon settings advice.

With the new develop, I did not duplicate the error, but I wouldn't take that as definitive proof of having been fixed. I'll post here if I encounter the error again.

Both fmincon settings you propose are enough to get bmibnb to produce incumbents  But using my R2019b fmincon on this bmibnb test problem,  the 1st incumbent was produced on bmibnb iteration 7 for both of your proposed fmincon, settings. And the first "decent" incumbent of objective value -603.65 was not produced until iteration 35 and 20 respectively for the 1st and 2nd of your proposed settings. I did not get the global optimum incumbent on iteration 2 as you did.

I will try your proposed fmincon settings in"combat" conditions to see how they compare to knitro. I suppose if you had knitro setup working, you could probably fine tune knitro settings. The best settings for knitro may not be the same as corresponding parameter values for fmincon due to differences between the solvers, despite having same-named algorithms and parameters.

[Johan Löfberg]

when I worked with this last spring, I used knitro and fmincon for development, and on my tests they where sort of similar, in that sometimes one worked better than the other, sometimes both failed, and sometimes both worked well. I think sqp worked best in fmincon if I recall it correctly, but then I found the trick with the barrier parameter which improved that one.

there is also a setting in slayer that can be set to convex or nonconvex which sometimes makes a huge difference (I think nonconvex worked best)

[Mark L. Stone]

Does BMIBNB now (develop version) attempt to find solutions for general Nonlinear SDPs having any type of nonlinearity within YALMIP's standard list of functions, as well as any other nonlinear constraints? This is as opposed to YALMIP's implementation for PENLAB (which when called natively addresses general Nonlinear SDPs and nonlinear constraints),  which only handles BMIs and Polynomial Matrix Inequalities which can be lifted to BMIs.

  Or is BMIBNB restricted to BMIs despite calling general nonlinear solvers such as FMINCON and KNITRO, which seem like they could support subproblems spawned by general Nonlinear SDPs and nonlinear constraints?

If this was addressed in a previous post, chalk it up to my advancing years and declining memory and mental acuity.

[Johan Löfberg]

The BMIBNB can attack nonlinear SDP solvers in two different ways, via https://yalmip.github.io/nonlinearsdpcuts or by using the Slayer framework (i.e. via  fmincon or any other nonlinear solver (although only fmincon is supported now I think))  as the local nonlinear solver