Home > examples > elliptope_SDP_complex.m

elliptope_SDP_complex

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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

 A is a Hermitian 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, X is complex, positive semidefinite.

 In practice, the Hermitian 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
 (that is, sqrt(n)) 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, Y complex

 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 Hermitian 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [Y, problem, S] = elliptope_SDP_complex(A, p, Y0)
0002 % Solver for complex semidefinite programs (SDP's) with unit diagonal.
0003 %
0004 % function [Y, problem, S] = elliptope_SDP_complex(A)
0005 % function [Y, problem, S] = elliptope_SDP_complex(A, p)
0006 % function [Y, problem, S] = elliptope_SDP_complex(A, p, Y0)
0007 %
0008 % A is a Hermitian 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, X is complex, positive semidefinite.
0013 %
0014 % In practice, the Hermitian matrix X of size n is parameterized as
0015 % X = Y*Y', where Y has size n x p. By default, p is taken large enough
0016 % (that is, sqrt(n)) 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, Y complex
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 Hermitian 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
0054 
0055 % This file is part of Manopt: www.manopt.org.
0056 % Original author: Nicolas Boumal, Oct. 21, 2016
0057 % Contributors:
0058 % Change log:
0059 
0060 
0061     % If no inputs are provided, since this is an example file, generate
0062     % a random complex matrix. This is for illustration purposes only.
0063     if ~exist('A', 'var') || isempty(A)
0064         n = 100;
0065         A = randn(n) + 1i*randn(n);
0066         A = (A+A')/sqrt(2*n);
0067     end
0068 
0069     n = size(A, 1);
0070     assert(n >= 2, 'A must be at least 2x2.');
0071     assert(size(A, 2) == n, 'A must be square.');
0072     
0073     % Force A to be Hermitian
0074     A = (A+A')/2;
0075     
0076     % By default, pick a sufficiently large p (number of columns of Y).
0077     if ~exist('p', 'var') || isempty(p)
0078         p = floor(sqrt(n)+1);
0079     end
0080     
0081     assert(p >= 1 && p == round(p), 'p must be an integer >= 1.');
0082 
0083     % Pick the manifold of complex n-by-p matrices with unit norm rows.
0084     manifold = obliquecomplexfactory(p, n, true);
0085     
0086     problem.M = manifold;
0087     
0088     
0089     % These three, quick commented lines of code are sufficient to define
0090     % the cost function and its derivatives. This is good code to write
0091     % when prototyping. Below, a more advanced use of Manopt is shown,
0092     % where the redundant computation A*Y is avoided between the gradient
0093     % and the cost evaluation.
0094     % % problem.cost  = @(Y) .5*sum(sum(real((A*Y).*conj(Y))));
0095     % % problem.egrad = @(Y) A*Y;
0096     % % problem.ehess = @(Y, Ydot) A*Ydot;
0097     
0098     % Products with A dominate the cost, hence we store the result.
0099     % This allows to share the results among cost, grad and hess.
0100     % This is completely optional.
0101     function store = prepare(Y, store)
0102         if ~isfield(store, 'AY')
0103             AY = A*Y;
0104             store.AY = AY;
0105             store.diagAYYt = sum(real(AY .* conj(Y)), 2);
0106         end
0107     end
0108     
0109     % Define the cost function to be /minimized/.
0110     problem.cost = @cost;
0111     function [f, store] = cost(Y, store)
0112         store = prepare(Y, store);
0113         f = .5*sum(store.diagAYYt);
0114     end
0115 
0116     % Define the Riemannian gradient.
0117     problem.grad = @grad;
0118     function [G, store] = grad(Y, store)
0119         store = prepare(Y, store);
0120         G = store.AY - bsxfun(@times, Y, store.diagAYYt);
0121     end
0122 
0123     % If you want to, you can specify the Riemannian Hessian as well.
0124     problem.hess = @hess;
0125     function [H, store] = hess(Y, Ydot, store)
0126         store = prepare(Y, store);
0127         SYdot = A*Ydot - bsxfun(@times, Ydot, store.diagAYYt);
0128         H = manifold.proj(Y, SYdot);
0129     end
0130 
0131 
0132     % If no initial guess is available, tell Manopt to use a random one.
0133     if ~exist('Y0', 'var') || isempty(Y0)
0134         Y0 = [];
0135     end
0136 
0137     % Call your favorite solver.
0138     opts = struct();
0139     opts.verbosity = 0;      % Set to 0 for no output, 2 for normal output
0140     opts.maxinner = 500;     % maximum Hessian calls per iteration
0141     opts.tolgradnorm = 1e-6; % tolerance on gradient norm
0142     Y = trustregions(problem, Y0, opts);
0143     
0144     % If required, produce an optimality certificate.
0145     if nargout >= 3
0146         S = A - spdiags(sum(real((A*Y).*conj(Y)), 2), 0, n, n);
0147     end
0148 
0149 end

Generated on Mon 10-Sep-2018 11:48:06 by m2html © 2005