Home > examples > maxcut.m

maxcut

PURPOSE ^

Algorithm to (try to) compute a maximum cut of a graph, via SDP approach.

SYNOPSIS ^

function [x, cutvalue, cutvalue_upperbound, Y] = maxcut(L, r)

DESCRIPTION ^

 Algorithm to (try to) compute a maximum cut of a graph, via SDP approach.
 
 function x = maxcut(L)
 function [x, cutvalue, cutvalue_upperbound, Y] = maxcut(L, r)

 L is the Laplacian matrix describing the graph to cut. The Laplacian of a
 graph is L = D - A, where D is the diagonal degree matrix (D(i, i) is the
 sum of the weights of the edges adjacent to node i) and A is the
 symmetric adjacency matrix of the graph (A(i, j) = A(j, i) is the weight
 of the edge joining nodes i and j). If L is sparse, this will be
 exploited.

 If the graph has n nodes, then L is nxn and the output x is a vector of
 length n such that x(i) is +1 or -1. This partitions the nodes of the
 graph in two classes, in an attempt to maximize the sum of the weights of
 the edges that go from one class to the other (MAX CUT problem).

 cutvalue is the sum of the weights of the edges 'cut' by the partition x.

 If the algorithm reached the global optimum of the underlying SDP
 problem, then it produces an upperbound on the maximum cut value. This
 value is returned in cutvalue_upperbound if it is found. Otherwise, that
 output is set to NaN.

 If r is specified (by default, r = n), the algorithm will stop at rank r.
 This may prevent the algorithm from reaching a globally optimal solution
 for the underlying SDP problem (but can greatly help in keeping the
 execution time under control). If a global optimum of the SDP is reached
 before rank r, the algorithm will stop of course.

 Y is a matrix of size nxp, with p <= r, such that X = Y*Y' is the best
 solution found for the underlying SDP problem. If cutvalue_upperbound is
 not NaN, then Y*Y' is optimal for the SDP and cutvalue_upperbound is its
 cut value.
 
 By Goemans and Williamson 1995, it is known that if the optimal value of
 the SDP is reached, then the returned cut, in expectation, is at most at
 a fraction 0.878 of the optimal cut. (This is not exactly valid because
 we do not use random projection here; sign(Y*randn(size(Y, 2), 1)) will
 give a cut that respects this statement -- it's usually worse though).

 The algorithm is essentially that of:
 Journee, Bach, Absil and Sepulchre, SIAM 2010
 Low-rank optimization on the cone of positive semidefinite matrices.

 It is itself based on the famous SDP relaxation of MAX CUT:
 Goemans and Williamson, 1995
 Improved approximation algorithms for maximum cut and satisfiability
 problems using semidefinite programming.
 
 See also: elliptope_SDP

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [x, cutvalue, cutvalue_upperbound, Y] = maxcut(L, r)
0002 % Algorithm to (try to) compute a maximum cut of a graph, via SDP approach.
0003 %
0004 % function x = maxcut(L)
0005 % function [x, cutvalue, cutvalue_upperbound, Y] = maxcut(L, r)
0006 %
0007 % L is the Laplacian matrix describing the graph to cut. The Laplacian of a
0008 % graph is L = D - A, where D is the diagonal degree matrix (D(i, i) is the
0009 % sum of the weights of the edges adjacent to node i) and A is the
0010 % symmetric adjacency matrix of the graph (A(i, j) = A(j, i) is the weight
0011 % of the edge joining nodes i and j). If L is sparse, this will be
0012 % exploited.
0013 %
0014 % If the graph has n nodes, then L is nxn and the output x is a vector of
0015 % length n such that x(i) is +1 or -1. This partitions the nodes of the
0016 % graph in two classes, in an attempt to maximize the sum of the weights of
0017 % the edges that go from one class to the other (MAX CUT problem).
0018 %
0019 % cutvalue is the sum of the weights of the edges 'cut' by the partition x.
0020 %
0021 % If the algorithm reached the global optimum of the underlying SDP
0022 % problem, then it produces an upperbound on the maximum cut value. This
0023 % value is returned in cutvalue_upperbound if it is found. Otherwise, that
0024 % output is set to NaN.
0025 %
0026 % If r is specified (by default, r = n), the algorithm will stop at rank r.
0027 % This may prevent the algorithm from reaching a globally optimal solution
0028 % for the underlying SDP problem (but can greatly help in keeping the
0029 % execution time under control). If a global optimum of the SDP is reached
0030 % before rank r, the algorithm will stop of course.
0031 %
0032 % Y is a matrix of size nxp, with p <= r, such that X = Y*Y' is the best
0033 % solution found for the underlying SDP problem. If cutvalue_upperbound is
0034 % not NaN, then Y*Y' is optimal for the SDP and cutvalue_upperbound is its
0035 % cut value.
0036 %
0037 % By Goemans and Williamson 1995, it is known that if the optimal value of
0038 % the SDP is reached, then the returned cut, in expectation, is at most at
0039 % a fraction 0.878 of the optimal cut. (This is not exactly valid because
0040 % we do not use random projection here; sign(Y*randn(size(Y, 2), 1)) will
0041 % give a cut that respects this statement -- it's usually worse though).
0042 %
0043 % The algorithm is essentially that of:
0044 % Journee, Bach, Absil and Sepulchre, SIAM 2010
0045 % Low-rank optimization on the cone of positive semidefinite matrices.
0046 %
0047 % It is itself based on the famous SDP relaxation of MAX CUT:
0048 % Goemans and Williamson, 1995
0049 % Improved approximation algorithms for maximum cut and satisfiability
0050 % problems using semidefinite programming.
0051 %
0052 % See also: elliptope_SDP
0053 
0054 % This file is part of Manopt: www.manopt.org.
0055 % Original author: Nicolas Boumal, July 18, 2013
0056 % Contributors:
0057 % Change log:
0058 %
0059 %   April 3, 2015 (NB):
0060 %       L products now counted with the new shared memory system. This is
0061 %       more reliable and more flexible than using a global variable.
0062 
0063 
0064     % If no inputs are provided, generate a random graph Laplacian.
0065     % This is for illustration purposes only.
0066     if ~exist('L', 'var') || isempty(L)
0067         n = 20;
0068         A = triu(randn(n) <= .4, 1);
0069         A = A+A';
0070         D = diag(sum(A, 2));
0071         L = D-A;
0072     end
0073 
0074 
0075     n = size(L, 1);
0076     assert(size(L, 2) == n, 'L must be square.');
0077 
0078     if ~exist('r', 'var') || isempty(r) || r > n
0079         r = n;
0080     end
0081     
0082     % We will let the rank increase. Each rank value will generate a cut.
0083     % We have to go up in the rank to eventually find a certificate of SDP
0084     % optimality. This in turn will provide an upperbound on the MAX CUT
0085     % value and ensure that we're doing well, according to Goemans and
0086     % Williamson's argument. In practice though, the good cuts often come
0087     % up for low rank values, so we better keep track of the best one.
0088     best_x = ones(n, 1);
0089     best_cutvalue = 0;
0090     cutvalue_upperbound = NaN;
0091     
0092     time = [];
0093     cost = [];
0094     
0095     for rr = 2 : r
0096         
0097         manifold = elliptopefactory(n, rr);
0098         
0099         if rr == 2
0100             
0101             % At first, for rank 2, generate a random point.
0102             Y0 = manifold.rand();
0103              
0104         else
0105             
0106             % To increase the rank, we could just add a column of zeros to
0107             % the Y matrix. Unfortunately, this lands us in a saddle point.
0108             % To escape from the saddle, we may compute an eigenvector of
0109             % Sy associated to a negative eigenvalue: that will yield a
0110             % (second order) descent direction Z. See Journee et al ; Sy is
0111             % linked to dual certificates for the SDP.
0112             Y0 = [Y zeros(n, 1)];
0113             LY0 = L*Y0;
0114             Dy = spdiags(sum(LY0.*Y0, 2), 0, n, n);
0115             Sy = (Dy - L)/4;
0116             % Find the smallest (the "most negative") eigenvalue of Sy.
0117             eigsopts.issym = true;
0118             eigsopts.isreal = true;
0119             [v, s] = eigs(Sy, 1, 'SA', eigsopts);
0120             % If there is no negative eigenvalue for Sy, than we are not at
0121             % a saddle point: we're actually done!
0122             if s >= -1e-8
0123                 % We can stop here: we found the global optimum of the SDP,
0124                 % and hence the reached cost is a valid upper bound on the
0125                 % maximum cut value.
0126                 cutvalue_upperbound = max(-[info.cost]);
0127                 break;
0128             end
0129             
0130             % This is our escape direction.
0131             Z = manifold.proj(Y0, [zeros(n, rr-1) v]);
0132             
0133             % % These instructions can be uncommented to see what the cost
0134             % % function looks like at a saddle point. But will require the
0135             % % problem structure which is not defined here: see the helper
0136             % % function.
0137             % plotprofile(problem, Y0, Z, linspace(-1, 1, 101));
0138             % drawnow; pause;
0139             
0140             % Now make a step in the Z direction to escape from the saddle.
0141             % It is not obvious that it is ok to do a unit step ... perhaps
0142             % need to be cautious here with the stepsize. It's not too
0143             % critical though: the important point is to leave the saddle
0144             % point. But it's nice to guarantee monotone decrease of the
0145             % cost, and we can't do that with a constant step (at least,
0146             % not without a proper argument to back it up).
0147             stepsize = 1;
0148             Y0 = manifold.retr(Y0, Z, stepsize);
0149             
0150         end
0151         
0152         % Use the Riemannian optimization based algorithm lower in this
0153         % file to reach a critical point (typically a local optimizer) of
0154         % the max cut cost with fixed rank, starting from Y0.
0155         [Y, info] = maxcut_fixedrank(L, Y0);
0156         
0157         % Some info logging.
0158         thistime = [info.time];
0159         if ~isempty(time)
0160             thistime = time(end) + thistime;
0161         end
0162         time = [time thistime]; %#ok<AGROW>
0163         cost = [cost [info.cost]]; %#ok<AGROW>
0164 
0165         % Time to turn the matrix Y into a cut.
0166         % We can either do the random rounding as follows:
0167         % x = sign(Y*randn(rr, 1));
0168         % or extract the "PCA direction" of the points in Y and cut
0169         % orthogonally to that direction, as follows (seems faster than
0170         % calling svds):
0171         [U, ~, ~] = svd(Y, 0);
0172         u = U(:, 1);
0173         x = sign(u);
0174 
0175         cutvalue = (x'*L*x)/4;
0176         if cutvalue > best_cutvalue
0177             best_x = x;
0178             best_cutvalue = cutvalue;
0179         end
0180         
0181     end
0182     
0183     x = best_x;
0184     cutvalue = best_cutvalue;
0185     
0186     plot(time, -cost, '.-');
0187     xlabel('Time [s]');
0188     ylabel('Relaxed cut value');
0189     title('The relaxed cut value is an upper bound on the optimal cut value.');
0190 
0191 end
0192 
0193 
0194 function [Y, info] = maxcut_fixedrank(L, Y)
0195 % Try to solve the (fixed) rank r relaxed max cut program, based on the
0196 % Laplacian of the graph L and an initial guess Y. L is nxn and Y is nxr.
0197 
0198     [n, r] = size(Y);
0199     assert(all(size(L) == n));
0200     
0201     % The fixed rank elliptope geometry describes symmetric, positive
0202     % semidefinite matrices of size n with rank r and all diagonal entries
0203     % are 1.
0204     manifold = elliptopefactory(n, r);
0205     
0206     % % If you want to compare the performance of the elliptope geometry
0207     % % against the (conceptually simpler) oblique manifold geometry,
0208     % % uncomment this line.
0209     % manifold = obliquefactory(r, n, true);
0210     
0211     problem.M = manifold;
0212     
0213     % % For rapid prototyping, these lines suffice to describe the cost
0214     % % function and its gradient and Hessian (here expressed using the
0215     % % Euclidean gradient and Hessian).
0216     % problem.cost  = @(Y)  -trace(Y'*L*Y)/4;
0217     % problem.egrad = @(Y) -(L*Y)/2;
0218     % problem.ehess = @(Y, U) -(L*U)/2;
0219     
0220     % Instead of the prototyping version, the functions below describe the
0221     % cost, gradient and Hessian using the caching system (the store
0222     % structure). This alows to execute exactly the required number of
0223     % multiplications with the matrix L. These multiplications are counted
0224     % using the shared memory in the store structure: that memory is
0225     % shared , so we get access to the same data, regardless of the
0226     % point Y currently visited.
0227 
0228     % For every visited point Y, we will need L*Y. This function makes sure
0229     % the quantity L*Y is available, but only computes it if it wasn't
0230     % already computed.
0231     function store = prepare(Y, store)
0232         if ~isfield(store, 'LY')
0233             % Compute and store the product for the current point Y.
0234             store.LY = L*Y;
0235             % Create / increment the shared counter (independent of Y).
0236             if isfield(store.shared, 'counter')
0237                 store.shared.counter = store.shared.counter + 1;
0238             else
0239                 store.shared.counter = 1;
0240             end
0241         end
0242     end
0243 
0244     problem.cost = @cost;
0245     function [f, store] = cost(Y, store)
0246         store = prepare(Y, store);
0247         LY = store.LY;
0248         f = -(Y(:)'*LY(:))/4; % = -trace(Y'*LY)/4; but faster
0249     end
0250 
0251     problem.egrad = @egrad;
0252     function [g, store] = egrad(Y, store)
0253         store = prepare(Y, store);
0254         LY = store.LY;
0255         g = -LY/2;
0256     end
0257 
0258     problem.ehess = @ehess;
0259     function [h, store] = ehess(Y, U, store)
0260         store = prepare(Y, store); % this line is not strictly necessary
0261         LU = L*U;
0262         store.shared.counter = store.shared.counter + 1;
0263         h = -LU/2;
0264     end
0265 
0266     % statsfun is called exactly once after each iteration (including after
0267     % the evaluation of the cost at the initial guess). We then register
0268     % the value of the L-products counter (which counts how many products
0269     % were needed so far).
0270     % options.statsfun = @statsfun;
0271     % function stats = statsfun(problem, Y, stats, store) %#ok
0272     %     stats.Lproducts = store.shared.counter;
0273     % end
0274     % Equivalent, but simpler syntax:
0275     options.statsfun = statsfunhelper('Lproducts', ...
0276                      @(problem, Y, stats, store) store.shared.counter );
0277     
0278 
0279     % % Diagnostics tools: to make sure the gradient and Hessian are
0280     % % correct during the prototyping stage.
0281     % checkgradient(problem); pause;
0282     % checkhessian(problem); pause;
0283     
0284     % % To investigate the effect of the rotational invariance when using
0285     % % the oblique or the elliptope geometry, or to study the saddle point
0286     % % issue mentioned above, it is sometimes interesting to look at the
0287     % % spectrum of the Hessian. For large dimensions, this is slow!
0288     % stairs(sort(hessianspectrum(problem, Y)));
0289     % drawnow; pause;
0290     
0291     
0292     % % When facing a saddle point issue as described in the master
0293     % % function, and when no sure mechanism exists to find an escape
0294     % % direction, it may be helpful to set useRand to true and raise
0295     % % miniter to more than 1, when using trustregions. This will tell the
0296     % % solver to not stop before at least miniter iterations were
0297     % % accomplished (thus disregarding the zero gradient at the saddle
0298     % % point) and to use random search directions to kick start the inner
0299     % % solve (tCG) step. It is not as efficient as finding a sure escape
0300     % % direction, but sometimes it's the best we have.
0301     % options.useRand = true;
0302     % options.miniter = 5;
0303     
0304     options.verbosity = 2;
0305     [Y, Ycost, info] = trustregions(problem, Y, options); %#ok<ASGLU>
0306     
0307     fprintf('Products with L: %d\n', max([info.Lproducts]));
0308 
0309 end

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005