maxcut.gms : Goemans/Williamson Randomized Approximation Algorithm for MaxCut

Description

```Let G(N, E) denote a graph. A cut is a partition of the vertices N
into two sets S and T. Any edge (u,v) in E with u in S and v in T is
said to be crossing the cut and is a cut edge. The size of the cut is
defined to be sum of weights of the edges crossing the cut.

This model presents a simple MIP formulation of the problem that is
seeded with a solution from the Goemans/Williamson randomized
approximation algorithm based on a semidefinite programming
relaxation. By default CSDP is used to solve the SDP.
Use --SDPSOLVER=MOSEK to switch to Mosek.

The MaxCut instance tg20_7777 is available from the Biq Mac Library
and comes from applications in statistical physics.
```

Large Model of Type : MIP

Category : GAMS Model library

Main file : maxcut.gms   includes :  runcsdp.inc  tg207777.inc

``````\$title Goemans/Williamson Randomized Approximation Algorithm for MaxCut (MAXCUT,SEQ=338)

\$onText
Let G(N, E) denote a graph. A cut is a partition of the vertices N
into two sets S and T. Any edge (u,v) in E with u in S and v in T is
said to be crossing the cut and is a cut edge. The size of the cut is
defined to be sum of weights of the edges crossing the cut.

This model presents a simple MIP formulation of the problem that is
seeded with a solution from the Goemans/Williamson randomized
approximation algorithm based on a semidefinite programming
relaxation. By default CSDP is used to solve the SDP.
Use --SDPSOLVER=MOSEK to switch to Mosek.

The MaxCut instance tg20_7777 is available from the Biq Mac Library
and comes from applications in statistical physics.

Wiegele A., Biq Mac Library - Binary Quadratic and Max Cut Library.
http://biqmac.uni-klu.ac.at/biqmaclib.html

Goemans M.X., and Williamson, D.P., Improved Approximation Algorithms
for Maximum Cut and Satisfiability Problems Using Semidefinite
Programming. Journal of the ACM 42 (1995), 1115-1145.
http://www-math.mit.edu/~goemans/PAPERS/maxcut-jacm.pdf

Keywords: mixed integer linear programming, approximation algorithms,
convex optimization, randomized algorithms, maximum cut problem,
mathematics
\$offText

\$if not set SDPSOLVER \$set SDPSOLVER CSDP

Set n 'nodes';

Alias (n,i,j);

Parameter w(i,j) 'edge weights';

Set e(i,j) 'edges';

\$if not set instance \$set instance tg207777.inc
\$onEmbeddedCode Python:
with open("%instance%", "r") as f:
n, _ = [int(i) for i in f.readline().split()]
gams.set('n', [str(i+1) for i in range(n)])
gams.set('w', [(i,j,int(v)) for i,j,v in [l.split() for l in f.readlines()]])
\$offEmbeddedCode n w

* We want all edges to be i-j with i<j;
e(i,j)    = ord(i) < ord(j);
w(e(i,j)) = w(i,j) + w(j,i);
w(i,j)\$(not e(i,j)) = 0;

option e < w;

* Simple MIP model
Variable
x(n)     'decides on what side of the cut'
cut(i,j) 'edge is in the cut'
z        'objective';

Binary Variable x;

Equation obj, xor1(i,j), xor2(i,j), xor3(i,j), xor4(i,j);

obj..          z      =e= sum(e, w(e)*cut(e));

xor1(e(i,j)).. cut(e) =l= x(i) + x(j);

xor2(e(i,j)).. cut(e) =l= 2 - x(i) - x(j);

xor3(e(i,j)).. cut(e) =g= x(i) - x(j);

xor4(e(i,j)).. cut(e) =g= x(j) - x(i);

Model maxcut / all /;

\$onText
Set up the SDP
max W*Y s.t. Y_ii = 1, Y positive semidefinite (psd)
We need to pass on the dual to csdp
min x1 + x2 + ... + xn s.t. X = F1*x1 + F2*x2 + ... + Fn*xn - W, X psd
with F_i = 1 for F_ii and 0 otherwise
\$offText

Scalar SDPRelaxation;
Parameter L(i,j) 'Cholesky factor of Y';

\$ifthen %SDPSOLVER% == CSDP
Parameter
c(n)     'cost coefficients'
F(n,i,j) 'constraint matrix'
F0(i,j)  'constant term'
Y(i,j)   'dual solution';

c(n)     =  1;
F(n,n,n) =  1;
F0(i,j)  = -w(i,j);

\$set ISYM n
\$set ASYM Y
\$set LSYM L

execute_unload 'csdpin.gdx' n = m, n, c, F, F0;
execute.checkErrorLevel 'gams runcsdp.inc lo=%gams.lo% --strict=1'
\$libInclude linalg cholesky n Y L

* Check if Cholesky factorization is correct
Parameter Y_, Ydiff;
Y_(i,j)    = sum(n, L(i,n)*L(j,n));
Ydiff(i,j) = round(Y(i,j) - Y_(i,j),1e-6);

option Ydiff:8:0:1;
abort\$card(Ydiff) Ydiff;

SDPRelaxation = 0.5*sum(e, w(e)*(1 - Y(e)));

\$else

Variable Y(i,j)    'PSDMATRIX';
Variable sdpobj    'objective function variable';
Equation sdpobjdef 'objective function W*Y';
sdpobjdef.. sum(e(i,j),w(i,j)*(Y(i,j)+Y(j,i))/2.0) + sum((i,j),eps*Y(i,j)) =E= sdpobj;
Y.fx(i,i) = 1.0;
Model sdp / sdpobjdef /;

option lp = %SDPSOLVER%;
option limrow = 0;
sdp.solprint = 0;
Solve sdp min sdpobj using lp;

Parameter Yl(i,j)   'level values of Y as parameter';
Yl(i,j) = Y.l(i,j);
\$libInclude linalg cholesky n Yl L

* Check if Cholesky factorization is correct
Parameter Y_, Ydiff;
Y_(i,j)    = sum(n, L(i,n)*L(j,n));
Ydiff(i,j) = round(Y.l(i,j) - Y_(i,j),1e-6);
option Ydiff:8:0:1;
abort\$card(Ydiff) Ydiff;

SDPRelaxation = 0.5*sum(e, w(e)*(1 - Y.l(e)));

\$endif

display SDPRelaxation;

* Now do the random hyperplane r
Parameter r(n);

Set S(n), T(n), bestS(n);

Scalar
wS    'weight of cut S'
maxwS 'best weight' / -inf /
cnt;

for(cnt = 1 to 10,
r(n) = uniform(-1,1);
S(n) = sum(i, L(n,i)*r(i)) < 0;
T(n) = yes;
T(S) =  no;
wS   = sum(e(S,T), w(S,T)) + sum(e(T,S), + w(T,S));
if(wS > maxwS, maxwS = wS; bestS(n) = S(n););
);
display maxwS;

* use computed feasible solution as starting point for MIP solve
x.l(bestS)    = 1;
cut.l(e(i,j)) = x.l(i) xor x.l(j);

* SCIP and COPT do this by default, for other solvers we need to enable it
\$set MIPSTART
\$if %gams.mip% == cplex  \$set MIPSTART mipStart
\$if %gams.mip% == cbc    \$set MIPSTART mipStart
\$if %gams.mip% == gurobi \$set MIPSTART mipStart
\$if %gams.mip% == highs  \$set MIPSTART mipStart
\$if %gams.mip% == xpress \$set MIPSTART loadmipsol
\$ifThen not x%MIPSTART% == x
\$  echo %MIPSTART% 1 > %gams.mip%.opt
maxcut.optFile = 1;
\$endIf
solve maxcut max z using mip;
``````
GAMS Development Corp.
GAMS Software GmbH

General Information and Sales
U.S. (+1) 202 342-0180
Europe: (+49) 221 949-9170