Re: Finding an n-variable polynomial's root.

Brief recap:

Robert Israel raised the issue of difficulty (NP-completeness) of
finding e.g. a maximal independent set of vertices in graphs. Maxim
Rytin responded with some code to do this sort of thing using NMinimize
in Mathematica. I wanted to say a few words about this and related
problems.

First, Rytin used Method->{DifferentialEvolution,SearchPoints->1} and
got some fairly good results. I wanted to point out that the heavy
lifting in this case is actually done by post-processing local
optimization code. There are three methods NMinimize might try: KKT
multipliers, an interior point method, and a local optmizer with
penalties/barriers used to enforce constraints. I believe it is the
last of these that is playing the main role. DitfferentialEvolution
cannot do much of anything without at least 4 points, and in general
will not do too well with fewer than 10 or so, so that method is
basically just punting to either a different one or going right to the
local stage, I'm not sure which.

I tried several different approaches. One is a linear reformulation. In
the simplest variant we replace x[i]*x[j]==0 by a linear constraint
x[i]+x[j]<=1, then branch on any noninteger solution (choosing the one
furthest from an integer seems best), spawning a pair of subproblems by
constraining that variable to be <=floor(value) or >=ceiling(value)
respectively, until we get one with all integers. A refinement is to
constrain the system more tightly by making use of all 3-complete
subgraphs. These allow us to add constraints of the form
x[i]+x[j]+x[k]<=1 (one shows this by adding the three constraints
between pairs to get 2*(x[i]+x[j]+x[k])<=3, dividing by 2, and noting
the lhs is an integer hence <=3/2 implies <=1: a basic cutting plane
tactic).

These reformulations are effectively bullet proof provided the problem
does not crash or hang. I found that using more constraints was slower
and moreover around n=140 it could crash my memory-limited machine. The
simpler approach did not have this deficiency, and is the fastest
method I could find at n=100, but around 170 (more or less varying
jointly with the graph "density" parameter), it seemed to get lost in
the looping in trying to enforce integrality. While I will not swear my
code was bug free, all smaller problems worked fine so I think it is
somehow intrinsic to the method and not the coding thereof.

All other methods I tried have an element of imperfection in that they
are not guaranteed to get optimal results. That said, as Rytin noted
some of these do quite well.

Below is code to generate examples.

shuffle[n_Integer, m_Integer] :=
Module[{res = Range[n]}, Do[rand = Random[Integer, {j, n}];
res[[{j, rand}]] = res[[{rand, j}]];, {j, m}];
Take[res, m]]

{n, m} = {120, 25};
density = 1/4;

SeedRandom[1111];
indepset = Sort[shuffle[n, m]];

mat = Array[
Boole[#1 < #2 && Random[] < density &&
Complement[{##}, indepset] =!= {}] &, {n, n}];

Here is one version of the quadratic programming code that seems
reasonably robust. It does not always give an optimal result, but it
tends always to give an independent set near the largest size. For
reference, on my 3. GHz Pentium 4 Linux box it takes a few seconds at
n=100, and about 20 minutes at n=500. What I found to work well for
enforcing integrality is to subtract multiples of the sum of the
variables and the sum of their squares. One can regard these as
"penalty" terms.

Note that we really do not need the "m" parameter, I was just using it
as a convenience for scaling the penalty terms.

independentVerticesQP[mat_, n_, m_] := Module[
{x, vars, poly1, poly2, poly3, poly, sol},
vars = Array[x, n];
poly1 = Sum[mat[[i,j]]*x[i]*x[j], {i, n - 1}, {j, i + 1, n}];
poly2 = Total[vars];
poly3 = vars.vars;
poly = poly1 - poly2/m - poly3/m;

sol = NMinimize[{poly, Flatten[{(0 <= # <= 1 &) /@ vars}]}, vars,
Method -> {RandomSearch, SearchPoints -> 1, Method ->
InteriorPoint}];
{sol[[1]], Chop[vars /. sol[[2]]],
Flatten[Position[Round[Chop[vars /. sol[[2]]]], _?(# =!= 0 &),
Heads -> False]]}
]

Timing[independentVerticesQP[mat, n, m]]

I think that the post-processing in this class of problems was needed
in order to get good results. But having looked at quite a number of
variations, I may be getting confused.

By way of contrast, the InteriorPoint setting alone seems to do quite
well on some quadratic assignment problems, provided one tinkers with
weighting "penalties" to encourage the result to be a valid
permutation. (If not set to False, the post-processing seems to eat too
much time, at least that was my impression.)

As tests I used some of the Nugent et al problems. We are given a pair
of matrices (typically, of nonnegative values) and asked to minimize
the sum of products of corresponding elements of mat1 *
perm.mat2.Transpose[perm] (that is, we permute both rows and columns of
mat2). The minimum is taken over all possible permutations.

http://www.seas.upenn.edu/qaplib/inst.html#NVR

The code I wrote for these is as below.

QAP[mat1_, mat2_, opts___] := Module[
{n=Length[mat1], vars, x, allvars, p1, p2, p3, poly, constraints,
rng, permutedmat2, soln, res},

vars = Array[x, {n, n}];
allvars = Flatten[vars];
constraints = Map[0 <= # <= 1 &, allvars];

permutedmat2 = vars.mat2.Transpose[vars];
p1 = Apply[Plus, Flatten[mat1*permutedmat2]];
p2 = Apply[Plus,
Join[Map[(Apply[Plus, #] - 1)^2 &, vars],
Map[(Apply[Plus, #] - 1)^2 &, Transpose[vars]]]];
p3 = allvars.allvars;

poly = p1 + 5000*p2 - 5*p3;

soln = NMinimize[{poly, constraints}, allvars, opts];
soln = soln /. a_Real :> Round[a];
res = vars /. soln[[2]];
{Map[Position[#,1][[1,1]]&,res], res, res.Transpose[res] ==
IdentityMatrix[n],
p1 /. soln[[2]], Apply[Plus, Flatten[mat1*mat2]]}
]

Here is one of the smaller examples, NUG15.

mat1nug15 = {
{0,1,2,3,4,1,2,3,4,5,2,3,4,5,6},
{1,0,1,2,3,2,1,2,3,4,3,2,3,4,5},
{2,1,0,1,2,3,2,1,2,3,4,3,2,3,4},
{3,2,1,0,1,4,3,2,1,2,5,4,3,2,3},
{4,3,2,1,0,5,4,3,2,1,6,5,4,3,2},
{1,2,3,4,5,0,1,2,3,4,1,2,3,4,5},
{2,1,2,3,4,1,0,1,2,3,2,1,2,3,4},
{3,2,1,2,3,2,1,0,1,2,3,2,1,2,3},
{4,3,2,1,2,3,2,1,0,1,4,3,2,1,2},
{5,4,3,2,1,4,3,2,1,0,5,4,3,2,1},
{2,3,4,5,6,1,2,3,4,5,0,1,2,3,4},
{3,2,3,4,5,2,1,2,3,4,1,0,1,2,3},
{4,3,2,3,4,3,2,1,2,3,2,1,0,1,2},
{5,4,3,2,3,4,3,2,1,2,3,2,1,0,1},
{6,5,4,3,2,5,4,3,2,1,4,3,2,1,0}
};

mat2nug15 = {
{0,10,0,5,1,0,1,2,2,2,2,0,4,0,0},
{10,0,1,3,2,2,2,3,2,0,2,0,10,5,0},
{0,1,0,10,2,0,2,5,4,5,2,2,5,5,5},
{5,3,10,0,1,1,5,0,0,2,1,0,2,5,0},
{1,2,2,1,0,3,5,5,5,1,0,3,0,5,5},
{0,2,0,1,3,0,2,2,1,5,0,0,2,5,10},
{1,2,2,5,5,2,0,6,0,1,5,5,5,1,0},
{2,3,5,0,5,2,6,0,5,2,10,0,5,0,0},
{2,2,4,0,5,1,0,5,0,0,10,5,10,0,2},
{2,0,5,2,1,5,1,2,0,0,0,4,0,0,5},
{2,2,2,1,0,0,5,10,10,0,0,5,0,5,0},
{0,0,2,0,3,0,5,0,5,4,5,0,3,3,0},
{4,10,5,2,0,2,5,5,10,0,0,3,0,10,2},
{0,5,5,5,5,5,1,0,0,0,5,3,10,0,4},
{0,0,5,0,5,10,0,0,2,5,0,0,2,4,0}
};

The optimal objective function value is 1150. The "baseline" from using
the identity as permutation is 1492. The code shown gets 1160 in about
50 seconds.

Timing[res =
QAP[mat1nug15, mat2nug15,
Method -> {RandomSearch, SearchPoints -> 1, Method -> InteriorPoint,
PostProcess -> False}]]

For NUG25 the min is 3744 and this code gets 3752 (baseline is 4838) in
a round 7.5 minutes. There are better discrete heuristic methods e.g
evolutionary and others, usually joined with a local phase possibly
including a tabu search. All the same I think that for a simple-to-code
heuristic this is quite good, in the size range indicated (I have my
doubts about enforcing the permutation constraints via penalties should
this get much larger).

Anyway, I found these investigations to be quite interesting and wanted
to pass along some of what I've seen from them


Daniel Lichtblau
Wolfram Research

.