Re: Representing a tensor by a homogeneous (ideally, Sum Of Squares) polynomial

[Johan Löfberg]

Nothing I've seen before, and I don't even now what it means for a tensor to be psd

Perhaps you can do similarly as when you create standard scalar sos from matrix sos. To ensure M(x) is matrix sos, one way is to simply introduce a new variable v, and say that v'*M(x)*v is sos. Maybe you can do similarly on your tensor to drop it down to a 2D sos, and then apply said scalarization to drop it once again to 1D

[Johan Löfberg]

you can do manual definition of sos easily with built-in functionalities, so you should be able to do something along the lines of the stuff in the last paragraph in the first section here

https://yalmip.github.io/tutorial/sumofsquaresprogramming/

[Johan Löfberg]

I don't see that as I don't understand what you try to do when you convert to a 2D matrix. I guess the problem is I don't have any insight/knowledge in the basic problem here. Convert a given constant tensor to a polynomial? What would the 2D analog be, i.e. given a 2D matrix convert it to a polynomial?

[Johan Löfberg]

Still don't see it. What is Y , and what are the dimensions of these things, and thus p(s). In my world, if s is nx1, kron(s,s) would be n^2x1, so the only dimension that would make sense is Y nxn^2, so why talk about tensors. If p(s) is scalar, in what sense is this not a standard prolem, i.e. look for max t such that p(s)-t is sos, possibly with some positivstellensatz extras to include the domain over which you want to find a lower bound. And if you simply want to solve an optimization problem min p(s) over s in S, why start messing with SOS to begin with? 

[Johan Löfberg]

sos does not convert your cubic to a quadratic form , so I think you misunderstood the concept of sos

There is really nothing special with your objective, it is a cubic. Hence, you either use a local or global solver to solve. In yalmip,

s = sdpvar(n,1)
objective=c'*s + s'*Y*kron(s,s);
constraints = ...
optimize(constraints,-objective,sdpsettings('solver','up to you'))

[Jack1]

But, isn't that cubic term non-convex? I wanted to get guarantees on its solution; thus, had been looking to transform it. 

[Johan Löfberg]

Well of course p is defined, you created it.

The sos module fails directly as the problem makes no sense. You have a cubic, hence it trivially cannot be sos. Should be better error detection inside the sos compilator though to inform about this instead of just failing

[Johan Löfberg]

of course it is nonconvex, that's why you use a global solver if you want to solve it globally

[Johan Löfberg]

Cannot reproduce that here (for ny random sample). Works and you get for instance

sdisplay(w{1}'*Q{1}*w{1})

0.322343930628*s(1)^4+1.71525107762*s(1)^3*s(2)+3.41581827041*s(1)^2*s(2)^2+2.07266418966*s(1)*s(2)^3+0.677557960459*s(2)^4

You would have to debug to see where the crash occurs

Also, I get a different display

Initially 6 monomials in R^2

Newton polytope (0 LPs).........Keeping 3 monomials (0.015625sec)

Finding symmetries..............Found 1 symmetry  (0sec)

Partitioning using symmetry.....3x3(1) 

Old YALMIP version/Octave/OS/MATLAB version?

[Johan Löfberg]

your display is very strange, it claims you have 3 variables. I think you are messing up some versions of code or something

[Jack1]

I get the following output and error:

-------------------------------------------------------------------------
YALMIP SOS module started...
-------------------------------------------------------------------------
Detected 0 parametric variables and 2 independent variables.
Detected 0 linear inequalities, 0 equality constraints and 0 LMIs.
Using kernel representation (options.sos.model=1).

Initially 6 monomials in R^2
Newton polytope (0 LPs).........Keeping 3 monomials (0.015625sec)

Finding symmetries..............Found 1 symmetry  (0.03125sec)
Partitioning using symmetry.....3x3(1)
 
>> sol


sol =


  struct with fields:


    solvertime: 0
          info: 'No suitable solver'
       problem: -2
    yalmiptime: 2.8620


>> sdisplay(w{1}'*Q{1}*w{1})
Brace indexing is not supported for variables of this type.

The sol, w, Q seems to be empty. My Matlab version is

MATLAB Version: 9.6.0.1114505 (R2019a) Update 2

And Yalmip version is

Version 25-April-2019

I would appreciate it if you can tell me how I can get to the root of this problem. Thank you.

[Johan Löfberg]

For which model? If you mean for the models you presented above, it means you don't have any SDP solver. If you have an SDP solver, it means you have defined a problem which leads to a nonlinear sdp problem (i.e. effectively unsolvable) which you don't have any solver for.

[Johan Löfberg]

second argument to solvesos should be objective ([] in your case)

don't turn off display in solver until you know everything runs as expected

[Johan Löfberg]

If you want to minimize an objective using sos, i.e. find a lower bound, you introduce a new variable t, and then maximize t under p-t being sos. However, how you extract the minimizing argument s, if the sos relaxation is tight which it typically won't be , is not straighforward. Much easier then is to solve the dual problem formulation, i.e. the moment relaxation with solvemoment (still the sme problem that you don't necessarily have a tight relaxation, so you typically cannot extract any optimizer)

 

As I aid before, the first thing you try is to simply solve it using a global solver (such as bmibnb). sos and moment relaxations are primarily for computing lower bounds or reason and optimize over parameterized infinite-dimensional set of constraints, and not necessarily finding global minimizers to standard optimization problems (although they can be efficient for some very special problems in lucky cases)

[Johan Löfberg]

As I said, a sos decomposition does not give you any kind of solution, it gives you  certificate that a function is sos. In this case, you have a parameterized sos, and you've computed a value of t, such that p-t is sos, i.e. a certificate that p>=t for all s. Finding an s such that p(s)=t (if such a value actually exist, which you don't know) is another affair. When the sos decomposition is tight, i.e. t coincides with the global minimima of p, it can, in some cases, be recovered from the dual of the solved SDP. That's a nasty affair, and if minimization and extraction is what you want, you would start from the dual side to begin with (moment relaxation via solvemoment)

bmibnb appears to have no problems here, solving it to desired tolerance in only 37 iterations (although I guess you've solved something else than the specified model, as you have no constraints there, and that optimal value should be 0 to be consistent with the sos lower bound)

[Johan Löfberg]

and -0.5 <= s<=0.5 can not be added to the model, as that tells YALMIP that s is a parameter in the model, and thus there is nothing left to perform a sos on

If you want to minimize p(s) on g(s)>=0, you apply the positivstellensatz and find a parameterized sos polynomial q(u,s) such that p(s)-t >= q(u,s)g(s) which is relaxed to finding t and u such that p(s)-t-q(u,s)g(s) being sos in s

[Jack1]

In case the above storage problem might be difficult to address, if there might be some ways to speed up the global optimization with bmibnb say for N = 100, then, it would be greatly appreciated. For the following, I waited ~20 mins and the iterations didn't seem to have started:

M = 100;
Y = rand(M,M^2);
s = sdpvar(M,1);
objective = s'*Y*kron(s,s);
Constraints = [-0.5 <= s <= 0.5];
tic; sol = optimize(Constraints,-objective,sdpsettings('solver','bmibnb','verbose',1)); toc
solution = value(s)

It appears to be stuck in the following if I forcefully terminate:

Operation terminated by user during update_monomial_bounds (line 28)

In bmibnb (line 129)
p = update_monomial_bounds(p);

In solvesdp (line 364)
        eval(['output = ' solver.call '(interfacedata);']);

In optimize (line 31)
[varargout{1:nargout}] = solvesdp(varargin{:});

Thank you again for your assistance. 

[Johan Löfberg]

Moment relaxations exploit (if one has implemented the same tricks) exactly the same structure as sos, as they are the same thing, but interpreted from two sides of the same idea. For instance, the fourth output from solvemoment is a sos decomposition, and a solution (if relaxation is tight) if some properties hold be extracted from the dual computed when solving the sos decomposition. The solvemoment code does not implement various symmetry exploiting tricks though as they are easier to do on the the sos side. Also note that these methods are extremely expensive. It is not something you can do on anything but small polynomials, unless there is a high degree of symmetry etc which can be exploited a-priori (which typically isn't trivial to do, i.e. simply finding out how to exploit symmetry or some structure on some polynomial or class of polynomials is a full research program on its own. Just because it's been hyped as a convex optimization method, does not mean it scales well, it's pretty horrible in general

bmibnb is a spatial branch&bound, meaning it has control over both upper via nonlinear solves and lower bounds with LP relaxations, and thus has a certificate of global (sub-)optimality when it terminates.

Your problem sizes are out of scope for relaxation based methods. For N=100, your function has 171700 monomial terms so it takes time to even construct the symbolic object, and a non-symmetry exploiting sos problem will require matrices of size nchoosek(100+2,2)=5125 which you won't solve with off-the-shelf solvers

A spatial b&b with a dense quartic with 171700 monomial terms won't work either. 100 variables is doable in lucky cases, but then there has to be structure so that the number of terms that you employ bounding strategies is limited, and not on the order of 100k

As you cannot even handle the data defining the polynomial, you are massively beyond any chance of attacking this using the methods discussed here.

[Jack1]

Thank you so much for your explanation. I do realize that I cannot work with a large number of variables since I won't even be able to characterize the polynomial. However, can I still do something say for N = 100? For the problem in my last post, bmibnb seems to give  

sol = 

  struct with fields:

    yalmipversion: '20190425'

       yalmiptime: NaN

       solvertime: NaN

             info: 'Unknown problem in solver (Turn on 'debug' in sdpsettings) (Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform elementwise multiplication, use '.*'.)'

          problem: 9

with NaNs for the solution. If I try ipopt then it says

sol =

  struct with fields:

    yalmipversion: '20190425'
       yalmiptime: 1.2923
       solvertime: 4.5957
             info: 'Successfully solved (IPOPT)'
          problem: 0

however, the solution is trivial (i.e., all zeros). 

Is there another nonlinear solver that I can try? Preferably one that has a free trial version. Thank you. 

[Johan Löfberg]

Exactly what code caused that crash (very strange crash)? And turn on 'debug' in the options

100 is still of a size where you have to be lucky. For global optimization, you essentially have exponential complexity in the worst case, i.e. completely intractable for n=100

I wouldn't be surprised if the problem is easy to solve, but extremely hard to certify, i.e. a local solver will give good solutions

Supported nonlinear solvers

https://yalmip.github.io/tags/#nonlinear-programming-solver

[Jack1]

I was using the following 

clear all; clc; close all;
M = 20;

Y = rand(M,M^2);
s = sdpvar(M,1);
objective = s

'*Y*kron(s,s);
Constraints = [-0.05 <= sum(s) <= 0.05, ...
            -0.01 <= s <= 0.01];
tic; sol = optimize(Constraints,-objective,sdpsettings('solver','bmibnb','verbose',1,'debug',1)); toc
solution = value(s)
sol

and get the following report:

* Starting YALMIP global branch & bound.
* Branch-variables : 20
* Upper solver     : fmincon
* Lower solver     : SCIP
* LP solver        : SCIP
 Node       Upper      Gap(%)       Lower    Open

Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of
rows in the second matrix. To perform elementwise multiplication, use '.*'.

Error in solvelower (line 158)
            cost = output.Primal'*p_cut.Q*output.Primal + p_cut.c'*output.Primal + p.f;


Error in branch_and_bound (line 133)
        [output,cost,p,timing] = solvelower(p,options,lowersolver,x_min,upper,timing);


Error in bmibnb (line 364)
    [x_min,solved_nodes,lower,upper,lower_hist,upper_hist,timing,counter,problem] = branch_and_bound(p,x_min,upper,timing);


Error in solvesdp (line 361)

    eval(['output = ' solver.call '(interfacedata);']);

Error in optimize (line 31)

[varargout{1:nargout}] = solvesdp(varargin{:});

But, if I change the solver to ipopt, then, I get:

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************


Total number of variables............................:       20
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       20
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        2
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        2




Number of Iterations....: 4


                                   (scaled)                 (unscaled)
Objective...............:   6.9932414955296644e-54    6.9932414955296644e-54
Dual infeasibility......:   2.3822801641527197e-22    2.3822801641527197e-22
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   9.0909090909091520e-09    9.0909090909091520e-09
Overall NLP error.......:   9.0909090909091520e-09    9.0909090909091520e-09




Number of objective function evaluations             = 5
Number of objective gradient evaluations             = 5
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 5
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 5
Number of Lagrangian Hessian evaluations             = 0
Total CPU secs in IPOPT (w/o function evaluations)   =      0.161
Total CPU secs in NLP function evaluations           =      0.177


EXIT: Optimal Solution Found.
Elapsed time is 4.440863 seconds.


solution =


   1.0e-18 *


   -0.6847
   -0.6847
    0.0174
    0.0174
    0.0174
   -0.0730
   -0.0730
   -0.0730
   -0.0730
   -0.0730
   -0.0730
   -0.0730
   -0.0730
   -0.0730
   -0.0730
   -0.0730
   -0.0730
   -0.0730
   -0.0730
   -0.0730

sol =


  struct with fields:


    yalmipversion: '20190425'

       yalmiptime: 1.0451
       solvertime: 1.6549

             info: 'Successfully solved (IPOPT)'
          problem: 0

though the solution is trivial. Also, nomad seems to keep on going. I really wanted to get close to being able to tackle ~100 decision variables....

[Jack1]

though the solution is trivial. If I tried nomad then it keeps on going. I really wanted to be able to tackle something close to 100 decision variables. Also, am hoping to try out knitro and possibly baron, but, not sure if these are good approaches to proceed on? 

[Jack1]

Also, how should I interpret the following output from nomad:

------------------------------------------------------------------
 This is NOMAD v3.8.1
 Authors: M. Abramson, C. Audet, G. Couture, J. Dennis,  S. Le Digabel, C. Tribes


 Problem Properties:
 # Decision Variables:                 50
 # Number of Objectives:                1
 # Number of Nonlinear Constraints:     2
------------------------------------------------------------------
1 0.0000000000
2 -0.0000008161
3 -0.0000043706
5 -0.0000119161
8 -0.0000322226
12 -0.0000609718
110 -0.0000614974
164 -0.0000704629
165 -0.0000932122
175 -0.0000975907
185 -0.0001033187
189 -0.0001098696
193 -0.0001153159
194 -0.0001248690
196 -0.0001563574
197 -0.0001783060
200 -0.0001871463
212 -0.0002110928
213 -0.0002139896
214 -0.0002307451
219 -0.0002314382
243 -0.0002774137
244 -0.0003114578
251 -0.0003169583
299 -0.0003170891
300 -0.0003174814
301 -0.0003190501
302 -0.0003251822
303 -0.0003377322
310 -0.0003486065
344 -0.0003721223
361 -0.0004023095
403 -0.0004442718
407 -0.0004577314
408 -0.0004620652
409 -0.0004768089
442 -0.0005069975
454 -0.0005451335
492 -0.0005624675
2000 -0.0005624675
2175 -0.0005624675
2210 -0.0005624675
2230 -0.0005624675
2288 -0.0005624675
------------------------------------------------------------------
Elapsed time is 14.069672 seconds.


solution =


   -0.0100
    0.0100
   -0.0100
    0.0100
   -0.0100
   -0.0100
    0.0100
    0.0100
   -0.0100
    0.0100
    0.0100
    0.0100
   -0.0100
    0.0100
   -0.0100
    0.0100
   -0.0100
    0.0100
   -0.0100
   -0.0100
   -0.0100
   -0.0100
   -0.0100
   -0.0100
    0.0100
    0.0100
   -0.0100
   -0.0100
   -0.0100
   -0.0100
    0.0100
   -0.0100
    0.0100
    0.0100
   -0.0100
    0.0100
    0.0100
   -0.0100
    0.0100
   -0.0100
    0.0100
   -0.0100
    0.0100
   -0.0100
    0.0100
   -0.0100
    0.0100
   -0.0100
    0.0100
    0.0100

sol =


  struct with fields:


    yalmipversion: '20190425'

       yalmiptime: 1.1070
       solvertime: 11.5050
             info: 'Successfully solved (NOMAD)'
          problem: 0

Are the values next to the iteration the duality gap? Or, is this the cost value?

I am guessing that it doesn't terminate at a global optima? Is this only a local? Not sure if the solution could have improved with further iterations, which I guess could only be quantified by the value of the cost? 

[Johan Löfberg]

nomad is a solver for derivative-free optimization, thus no reason to use here

knitro is a good local nonlinear

baron is the benchmark solver for global optimization

[Johan Löfberg]

I cannot reproduce that, bmibnb easily find the global solution (up to specified tolerances)

* Starting YALMIP global branch & bound.

* Branch-variables : 20

* Upper solver     : fmincon

* Lower solver     : GUROBI

* LP solver        : GUROBI

 Node       Upper      Gap(%)       Lower    Open

    1 :   -7.135E-05     0.37     -3.813E-03   2    

* Finished.  Cost: -7.1354e-05 Gap: 0.37343

* Termination with relative gap satisfied 

* Timing: 2% spent in upper solver (2 problems solved)

*         1% spent in lower solver (41 problems solved)

*         10% spent in LP-based domain reduction (40 problems solved)

It has failed in the lower bound computation, and you use scip in the lower bound computation, so my only guess is that scip doesn't work on your machine

[Jack1]

Thank you for getting back. Yes, you are right, scip wasn't compiled on my machine....it doesn't seem to come installed with the OPTI Toolbox. 

I have the scip-6.0.2 folder downloaded, but, need to make it on a windows machine (the contents of the folder are in the attached)....couldn't find an automatic way to build it unless if you might be able to suggest. Also, am trying to get both gurobi and baron installed (at least their trial versions) 

[Johan Löfberg]

don't bother with scip. use gurobi

[Jack1]

Thank you again, I am using mosek as my lower and linear solver and get the same output as yours above! 

Should I be looking at the gap to interpret the optimality of the found optima? I.e., gap closer to 0 is a measure of the global optima or am I thinking of this incorrectly? 

And if I had another nonlinear solver like Baron for the upper in bmibnb then hopefully my gap should get closer to 0?

[Johan Löfberg]

You optimize with such a small region around the origin on a cubic, so basically any guess close to the origin will be optimal within the default tolerance settings. Hence, to see how it really performs when it actually performs a search, you would have to crank up the tolerances

>> sol = optimize(Constraints,-objective,sdpsettings('solver','bmibnb','verbose',1,'debug',1,'bmibnb.absgaptol',1e-5,'bmibnb.relgaptol',1e-5));

* Starting YALMIP global branch & bound.

* Branch-variables : 20

* Upper solver     : fmincon

* Lower solver     : CPLEX

* LP solver        : CPLEX

 Node       Upper      Gap(%)       Lower    Open

    1 :   -6.211E-05     0.38     -3.847E-03   2  Improved solution  

    2 :   -6.211E-05     0.38     -3.847E-03   3    

    3 :   -1.626E-04     0.35     -3.628E-03   4  Improved solution  

    4 :   -1.626E-04     0.37     -3.829E-03   5    

    5 :   -1.626E-04     0.37     -3.829E-03   6    

    6 :   -1.626E-04     0.37     -3.822E-03   7    

baron is a global solver just like bmibnb, so it would not make any sense to use it as the upper bound solver inside bmibnb

[Jack1]

Thank you again. Would you say that Matlab's fmincon function is pretty decent for the upper solver to use in bmibnb or can I download another solver too? 

[Johan Löfberg]

it's better than what most people think