3

I'm looking at quadratic relaxation of maximum independent set problem (p.22 here), and found that FindMaximum fails for every graph I try, unless I give it optimal solution as the starting point. These quadratic programmes have 10-20 variables, so I expect them to be solvable.

  1. Is there a way to make Mathematica solve such quadratic programmes?
  2. Is there some quadratic programming package that's easy to call from within Mathematica?

Here's an example of failing FindMaximum, followed by working FindMaximum initialized at the solution

setupQuadratic[g_Graph] := (
   Ag = AdjacencyMatrix[g];
   A = IdentityMatrix[Length@VertexList@g] - Ag;
   cons = And @@ Table[0 <= x[v] <= 1, {v, VertexList@g}];
   vars = x /@ VertexList[g];
   indSet = FindIndependentVertexSet@g;
   xOpt = Array[Boole[MemberQ[indSet, #]] &, {Length@VertexList@g}];
   );

g = GraphData[{"Cubic", {10, 11}}];
setupQuadratic[g];
FindMaximum[{vars.A.vars, cons}, vars]
FindMaximum[{vars.A.vars, cons}, Thread[{vars, xOpt}]]

Here are other graphs I tried

{"DodecahedralGraph", "FruchtGraph", "TruncatedPrismGraph", \
"TruncatedTetrahedralGraph", {"Cubic", {10, 2}}, {"Cubic", {10, 
   3}}, {"Cubic", {10, 4}}, {"Cubic", {10, 6}}, {"Cubic", {10, 
   7}}, {"Cubic", {10, 11}}, {"Cubic", {10, 12}}, {"Cubic", {12, 
   5}}, {"Cubic", {12, 6}}, {"Cubic", {12, 7}}, {"Cubic", {12, 
   9}}, {"Cubic", {12, 10}}}
Yaroslav Bulatov
  • 57,332
  • 22
  • 139
  • 197

3 Answers3

2

Might try method shown in package located here. See problem 8

Daniel Lichtblau Wolfram Research

Daniel Lichtblau
  • 6,854
  • 1
  • 23
  • 30
  • interesting, thanks. BTW, I'm also getting good results by using SDP relaxation of the quadratic program -- http://mathematica-bits.blogspot.com/2011/03/semidefinite-programming-in-mathematica.html – Yaroslav Bulatov Mar 03 '11 at 10:01
1

It seems that Maximize will serve you better. Here is a modified version of your function, which returns a list of 2 results - the "manual" one and the one obtained by Maximize:

Clear[findIVSet];
findIVSet[g_Graph] :=
Module[{Ag, A, cons, vars, indSet, indSetFromMaximize, xOpt},
  Ag = AdjacencyMatrix[g];
  A = IdentityMatrix[Length@VertexList@g] - Ag;
  cons = And @@ Table[0 <= x[v] <= 1, {v, VertexList@g}];
  vars = x /@ VertexList[g];
  indSet = FindIndependentVertexSet@g;
  xOpt = Array[Boole[MemberQ[indSet, #]] &, {Length@VertexList@g}];
  {indSet, DeleteCases[vars /. (Last@
    Maximize[{vars.A.vars, cons}, vars,Integers] /. (x[i_] -> 1) :> (x[i] -> i)), 0]}];

Here are the results:

In[32]:= graphs = GraphData /@ {"DodecahedralGraph", "FruchtGraph", 
"TruncatedPrismGraph", "TruncatedTetrahedralGraph", {"Cubic", {10, 2}}, {"Cubic", {10, 
  3}}, {"Cubic", {10, 4}}, {"Cubic", {10, 6}}, {"Cubic", {10, 
  7}}, {"Cubic", {10, 11}}, {"Cubic", {10, 12}}, {"Cubic", {12, 
  5}}, {"Cubic", {12, 6}}, {"Cubic", {12, 7}}, {"Cubic", {12, 
  9}}, {"Cubic", {12, 10}}};


In[33]:= sets = findIVSet /@ graphs

Out[33]= {{{1, 2, 3, 8, 10, 11, 17, 20}, {5, 6, 7, 8, 14, 15, 17, 18}},
{{2, 4, 6, 11, 12}, {2, 4, 6, 11, 12}}, {{2, 7, 10, 12, 16, 18}, {8, 11, 13, 16, 17, 18}}, 
{{1, 4, 7, 12}, {4, 7, 9, 12}}, {{2,3, 8, 9}, {2, 3, 8, 9}}, {{1, 4, 7, 10}, {2, 5, 8, 9}}, 
{{1, 4, 7, 10}, {2, 4, 7, 9}}, {{2, 4, 5, 8}, {3, 6, 7, 9}}, {{2, 5, 8, 9}, {2, 5, 8, 9}}, 
{{1, 3, 7, 10}, {4, 5, 8, 9}}, {{1, 6, 8, 9}, {2, 3, 6, 10}}, {{1, 6, 7, 12}, {4, 5, 9, 10}}, 
{{3, 4, 7, 8, 12}, {3, 4, 7, 8, 12}}, {{1, 5, 8, 9}, {4, 5, 10, 11}}, 
{{1, 5, 6, 9, 10}, {3, 4, 7, 8, 12}}, {{3, 4, 7, 9, 10}, {3, 4, 7, 9, 10}}}

They are not always the same for "manual" ones and those from Maximize, but then there is more than one solution for an independent set. The results from Maximize are all independent sets, which is easily verified: 

In[34]:= MapThread[IndependentVertexSetQ, {graphs, sets[[All, 2]]}]

Out[34]= {True, True, True, True, True, True, True, True, True, True, True, True, True, 
True, True,True}
Leonid Shifrin
  • 22,449
  • 4
  • 68
  • 100
  • Well, I already have the independent set from FindIndependentVertexSet, which might internally use something like what you propose, I'm actually interested in finding maximum of the original quadratic, Real domain – Yaroslav Bulatov Mar 01 '11 at 17:48
  • @Yaroslav Why are you interested in the real domain? As far as I could see, it is clearly stated in the paper you linked, that all variables can only have values 0 or 1 (that is, Integer domain). And, for that matter, what makes you think that the results for the real domain will have anything to do with the maximal independent set problem? It seems quite clear that we are not searching for local extrema here, so the constraints are important. Extrema found with continuous constraints can be quite different from those with discrete constraints. Or did I totally miss the point? – Leonid Shifrin Mar 01 '11 at 18:04
  • Because solving it in Integer domain is NP-complete, whereas real-valued QP programming is tractable. Discarding integer constraints is known as a "relaxation", and it works surprisingly well for certain classes of graphs -- http://cstheory.stackexchange.com/questions/5170/lp-relaxation-of-independent-set – Yaroslav Bulatov Mar 01 '11 at 18:15
  • @Yaroslav It is tractable but maximizing is nonconvex (minimizing in this case is convex). So no guarantees of globality. – Daniel Lichtblau Mar 01 '11 at 18:29
  • Interesting, I guess this QP is not as easy as I thought. Apparently such non-convex QP problems are known as "Box-QP" – Yaroslav Bulatov Mar 01 '11 at 19:37
0

IMO, the reason FindMaximum doesn't work here is the wild nature of your function. I tried a grid with 1,048,576 samples in variable space and none achieve a higher value than zero. Your optimum starting value gets -20.

In[10]:= (x[1]^2 + x[2]^2 + x[3]^2 - 2 x[3] x[4] + x[4]^2 - 
  2 x[2] (x[3] + x[4]) + x[5]^2 - 2 x[3] x[6] - 2 x[5] x[6] + 
  x[6]^2 - 2 x[5] x[7] + x[7]^2 - 2 x[6] x[8] - 2 x[7] x[8] + 
  x[8]^2 - 2 x[7] x[9] + x[9]^2 - 2 x[1] (x[2] + x[5] + x[9]) - 
  2 x[4] x[10] - 2 x[8] x[10] - 2 x[9] x[10] + x[10]^2 /. 
 Thread[vars -> #]) & @@@ Tuples[{0.0, 0.333, 0.667, 1.0}, 10] // Max

Out[10]= 0.

In[11]:= (x[1]^2 + x[2]^2 + x[3]^2 - 2 x[3] x[4] + x[4]^2 - 
 2 x[2] (x[3] + x[4]) + x[5]^2 - 2 x[3] x[6] - 2 x[5] x[6] + 
 x[6]^2 - 2 x[5] x[7] + x[7]^2 - 2 x[6] x[8] - 2 x[7] x[8] + 
 x[8]^2 - 2 x[7] x[9] + x[9]^2 - 2 x[1] (x[2] + x[5] + x[9]) - 
 2 x[4] x[10] - 2 x[8] x[10] - 2 x[9] x[10] + x[10]^2 /. 
Thread[vars -> #]) & @@@ {xOpt}

Out[11]= {-20}
Sjoerd C. de Vries
  • 16,122
  • 3
  • 42
  • 94
  • The objective above simplifies to `Times[-20,Power[Slot[1],2]]` after `ReplaceAll`. One way to evaluate original objective is `vars.A.vars /. Thread[vars -> xOpt]` – Yaroslav Bulatov Mar 01 '11 at 18:00
  • btw, above can be fixed if you replace `@@@` with `/@` or `#` with `{##}` (but not both at the same time!) – Yaroslav Bulatov Mar 01 '11 at 18:09
  • I expanded vars.A.vars because it saves a factor of 2.3 in time on my machine. I am aware the 2nd part could have been written differently. I just did it this way to keep the structure of the code the same. – Sjoerd C. de Vries Mar 01 '11 at 23:40
  • my point is that there's a bug, objective should be `4` at `xOpt`, not `-20` – Yaroslav Bulatov Mar 02 '11 at 00:59