Home > examples > elliptope_SDP.m

elliptope_SDP

PURPOSE ^

Solver for semidefinite programs (SDP's) with unit diagonal constraints.

SYNOPSIS ^

function [Y, problem, S] = elliptope_SDP(A, p, Y0)

DESCRIPTION ^

 Solver for semidefinite programs (SDP's) with unit diagonal constraints.
 
 function [Y, problem, S] = elliptope_SDP(A)
 function [Y, problem, S] = elliptope_SDP(A, p)
 function [Y, problem, S] = elliptope_SDP(A, p, Y0)

 A is a real, symmetric matrix of size n.

 This function uses a local optimization method in Manopt to solve the SDP

   min_X  trace(A*X)  s.t.  diag(X) = 1 and X is positive semidefinite.

 In practice, the symmetric matrix X of size n is parameterized
 as X = Y*Y', where Y has size n x p. By default, p is taken large enough
 (about sqrt(2n)) to ensure that there exists an optimal X whose rank is
 smaller than p. This ensures that the SDP is equivalent to the new
 problem in Y:

   min_Y  trace(Y'*A*Y)  s.t.  diag(Y*Y') = 1.

 The constraints on Y require each row of Y to have unit norm, which is
 why Manopt is appropriate software to solve this problem. An optional
 initial guess can be specified via the input Y0.

 See the paper below for theory, specifically, for a proof that, for
 almost all A, second-order critical points of the problem in Y are
 globally optimal. In other words: there are no local traps in Y, despite
 non-convexity.

 Outputs:

       Y: is the best point found (an nxp matrix with unit norm rows.)
          To find X, form Y*Y' (or, more efficiently, study X through Y.)
 
       problem: is the Manopt problem structure used to produce Y.
 
       S: is a dual optimality certificate (a symmetric matrix of size n,
          sparse if A is sparse). The optimality gap (in the cost
          function) is at most n*min(eig(S)), for both Y and X = Y*Y'.
          Hence, if min(eig(S)) is close to zero, Y is close to globally
          optimal. This can be computed via eigs(S, 1, 'SR').
 
 Paper: https://arxiv.org/abs/1606.04970

 @inproceedings{boumal2016bmapproach,
   author  = {Boumal, N. and Voroninski, V. and Bandeira, A.S.},
   title   = {The non-convex {B}urer-{M}onteiro approach works on smooth semidefinite programs},
   booktitle={Neural Information Processing Systems (NIPS 2016)},
   year    = {2016}
 }
 
 See also: maxcut elliptope_SDP_complex

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [Y, problem, S] = elliptope_SDP(A, p, Y0)
0002 % Solver for semidefinite programs (SDP's) with unit diagonal constraints.
0003 %
0004 % function [Y, problem, S] = elliptope_SDP(A)
0005 % function [Y, problem, S] = elliptope_SDP(A, p)
0006 % function [Y, problem, S] = elliptope_SDP(A, p, Y0)
0007 %
0008 % A is a real, symmetric matrix of size n.
0009 %
0010 % This function uses a local optimization method in Manopt to solve the SDP
0011 %
0012 %   min_X  trace(A*X)  s.t.  diag(X) = 1 and X is positive semidefinite.
0013 %
0014 % In practice, the symmetric matrix X of size n is parameterized
0015 % as X = Y*Y', where Y has size n x p. By default, p is taken large enough
0016 % (about sqrt(2n)) to ensure that there exists an optimal X whose rank is
0017 % smaller than p. This ensures that the SDP is equivalent to the new
0018 % problem in Y:
0019 %
0020 %   min_Y  trace(Y'*A*Y)  s.t.  diag(Y*Y') = 1.
0021 %
0022 % The constraints on Y require each row of Y to have unit norm, which is
0023 % why Manopt is appropriate software to solve this problem. An optional
0024 % initial guess can be specified via the input Y0.
0025 %
0026 % See the paper below for theory, specifically, for a proof that, for
0027 % almost all A, second-order critical points of the problem in Y are
0028 % globally optimal. In other words: there are no local traps in Y, despite
0029 % non-convexity.
0030 %
0031 % Outputs:
0032 %
0033 %       Y: is the best point found (an nxp matrix with unit norm rows.)
0034 %          To find X, form Y*Y' (or, more efficiently, study X through Y.)
0035 %
0036 %       problem: is the Manopt problem structure used to produce Y.
0037 %
0038 %       S: is a dual optimality certificate (a symmetric matrix of size n,
0039 %          sparse if A is sparse). The optimality gap (in the cost
0040 %          function) is at most n*min(eig(S)), for both Y and X = Y*Y'.
0041 %          Hence, if min(eig(S)) is close to zero, Y is close to globally
0042 %          optimal. This can be computed via eigs(S, 1, 'SR').
0043 %
0044 % Paper: https://arxiv.org/abs/1606.04970
0045 %
0046 % @inproceedings{boumal2016bmapproach,
0047 %   author  = {Boumal, N. and Voroninski, V. and Bandeira, A.S.},
0048 %   title   = {The non-convex {B}urer-{M}onteiro approach works on smooth semidefinite programs},
0049 %   booktitle={Neural Information Processing Systems (NIPS 2016)},
0050 %   year    = {2016}
0051 % }
0052 %
0053 % See also: maxcut elliptope_SDP_complex
0054 
0055 % This file is part of Manopt: www.manopt.org.
0056 % Original author: Nicolas Boumal, June 28, 2016
0057 % Contributors:
0058 % Change log:
0059 
0060 
0061     % If no inputs are provided, since this is an example file, generate
0062     % a random Erdos-Renyi graph. This is for illustration purposes only.
0063     if ~exist('A', 'var') || isempty(A)
0064         n = 100;
0065         A = triu(rand(n) <= .1, 1);
0066         A = (A+A.')/(2*n);
0067     end
0068 
0069     n = size(A, 1);
0070     assert(n >= 2, 'A must be at least 2x2.');
0071     assert(isreal(A), 'A must be real.');
0072     assert(size(A, 2) == n, 'A must be square.');
0073     
0074     % Force A to be symmetric
0075     A = (A+A.')/2;
0076     
0077     % By default, pick a sufficiently large p (number of columns of Y).
0078     if ~exist('p', 'var') || isempty(p)
0079         p = ceil(sqrt(8*n+1)/2);
0080     end
0081     
0082     assert(p >= 2 && p == round(p), 'p must be an integer >= 2.');
0083 
0084     % Pick the manifold of n-by-p matrices with unit norm rows.
0085     manifold = obliquefactory(p, n, true);
0086     
0087     problem.M = manifold;
0088     
0089     
0090     % These three, quick commented lines of code are sufficient to define
0091     % the cost function and its derivatives. This is good code to write
0092     % when prototyping. Below, a more advanced use of Manopt is shown,
0093     % where the redundant computation A*Y is avoided between the gradient
0094     % and the cost evaluation.
0095     % % problem.cost  = @(Y) .5*sum(sum((A*Y).*Y));
0096     % % problem.egrad = @(Y) A*Y;
0097     % % problem.ehess = @(Y, Ydot) A*Ydot;
0098     
0099     % Products with A dominate the cost, hence we store the result.
0100     % This allows to share the results among cost, grad and hess.
0101     % This is completely optional.
0102     function store = prepare(Y, store)
0103         if ~isfield(store, 'AY')
0104             AY = A*Y;
0105             store.AY = AY;
0106             store.diagAYYt = sum(AY .* Y, 2);
0107         end
0108     end
0109     
0110     % Define the cost function to be /minimized/.
0111     problem.cost = @cost;
0112     function [f, store] = cost(Y, store)
0113         store = prepare(Y, store);
0114         f = .5*sum(store.diagAYYt);
0115     end
0116 
0117     % Define the Riemannian gradient.
0118     problem.grad = @grad;
0119     function [G, store] = grad(Y, store)
0120         store = prepare(Y, store);
0121         G = store.AY - bsxfun(@times, Y, store.diagAYYt);
0122     end
0123 
0124     % If you want to, you can specify the Riemannian Hessian as well.
0125     problem.hess = @hess;
0126     function [H, store] = hess(Y, Ydot, store)
0127         store = prepare(Y, store);
0128         SYdot = A*Ydot - bsxfun(@times, Ydot, store.diagAYYt);
0129         H = manifold.proj(Y, SYdot);
0130     end
0131 
0132 
0133     % If no initial guess is available, tell Manopt to use a random one.
0134     if ~exist('Y0', 'var') || isempty(Y0)
0135         Y0 = [];
0136     end
0137 
0138     % Call your favorite solver.
0139     opts = struct();
0140     opts.verbosity = 0;      % Set to 0 for no output, 2 for normal output
0141     opts.maxinner = 500;     % maximum Hessian calls per iteration
0142     opts.tolgradnorm = 1e-6; % tolerance on gradient norm
0143     Y = trustregions(problem, Y0, opts);
0144     
0145     % If required, produce an optimality certificate.
0146     if nargout >= 3
0147         S = A - spdiags(sum((A*Y).*Y, 2), 0, n, n);
0148     end
0149 
0150 end

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005