msdp.m

GloptiPoly 3: moments, optimization and semidefinite programming. Gloptipoly 3 is intended to so

function [P,M] = msdp(varargin)
% @MSDP/MSDP - Constructor of class MSDP, moment semidefinite problem
%
% The instruction
%  
%   [P,M] = MSDP(VARARGIN)
%
% generates a moment SDP problem P (class MSDP) with its associated
% measures arranged in a vector M (class MEAS). The moment SDP problem
% is a finite-dimensional truncation (also called relaxation) of
% an infinite-dimensional moment optimization problem consisting of
% support constraints, moment constraints and a moment objective function.
%
% Valid input arguments that can be specified in the comma separated list
% VARARGIN are as follows:
% - SUC = support constraints (class SUPCON);
% - MOC = moment constraints (class MOMCON);
% - OBJ = objective function (class MOMCON);
% - ORD = order of the SDP relaxation (integer of class DOUBLE).
%
% All of these input arguments are optional and can be provided in
% arbitrary order. Default values are as follows:
% - SUC = no support constraints, meaning that all the measures M
%   are supported on the whole Euclidean space;
% - MOC = no moment constraints. In this case, when there is only
%   one measure (LENGTH(M) = 1), it is automatically constrained
%   to be a probability measure, i.e. such that MASS(M) == 1;
% - OBJ = no objective function. Then the sum of the traces of the
%   moment matrices of the measures M is minimized. To avoid this,
%   specify a constant (e.g. zero) objective function;
% - ORD = half the minimum degree over all polynomials defining
%   support constraints, moment constraints and objective function.
%
% When a support or moment constraint is specified with a unique
% left handside monomial, then MSDP performs explicitly moment 
% substitutions to reduce the number of variables in the moment SDP
% problem. The degree of the right handside polynomial must then be
% less than or equal to the degree of the left handside monomial.
%
% See also MSOL

% D. Henrion, 9 January 2004
% Last modified on 16 May 2007

global MMM;

if ~isfield(MMM,'yalmip')
 mset default
end

if MMM.verbose
 disp(['GloptiPoly ' gloptipolyversion]);
 disp('Define moment SDP problem')
end

% ----------------------------------------------------------------------  
% Parse input arguments
% ----------------------------------------------------------------------  
  
% Extract SDP relaxation order
% Check validity of input arguments

ord = []; % SDP relaxation order
pindmeas = []; % measure indices
pindvar = []; % active variable indices
tnmeas = max(MMM.indmeas); % total number of measures
dmax = -ones(1,tnmeas); % maximum degree wrt each measure
rankshift = ones(1,tnmeas); % rank shift for detecting global optimality

ind = 1:nargin;
for k = 1:nargin
 arg = varargin{k};
 if isa(arg,'double')
  if max(size(arg)) == 1
   if (floor(arg) ~= arg) | (arg <= 0)
    error('Relaxation order must be a positive integer')
   end
   ord = arg; % SDP relaxation order
   ind(k) = [];
  else
   error('Invalid input argument')
  end
 elseif isa(arg,'mpol')
  error('Invalid input polynomial')
 elseif ~isa(arg,'supcon') & ~isa(arg,'momcon')
  error('Invalid input argument')
 end
end
arg = {varargin{ind}}; % keep only SUPCON and MOMCON objects

% Detect objective functions (MOMCON objects of type 'min', 'max')

mobj = []; % row vector splitted wrt measures
objsign = 1; % max = +1, min = -1
for k = 1:length(arg)
 m = arg{k}; t = type(m);
 if isa(m,'momcon') & (strcmp(t,'min') | strcmp(t,'max'))
  if ~isempty(mobj)
   error('Moment objective function is not unique')
  end
  if strcmp(t,'min')
   mobj = -left(m); % moment to be minimized
   objsign = -1;
  else
   mobj = left(m); % to be maximized
  end
  if ~consistent(mobj) % check consistency
   error('Invalid objective function with inconsistent variables and measures')
  end   
  mp = indmeas(mobj); % measures
  if mp % non-constant objective function
   % store measure and variable indices
   pindmeas = [pindmeas mp];
   pindvar = [pindvar indvar(mobj)];
   % degree wrt each measure
   p = split(mobj); % polys wrt each measures
   for r = 1:length(mp)
    dmax(mp(r)) = max(dmax(mp(r)),deg(p(r)));
   end
  end
 end 
end

% Detect entrywise support constraints (SUPCON objects) and extract
% support substitutions which are support equality constraints
% with isolated left handside monomial

msupcge = []; % support inequalities (>=0)
msupceq = []; % support equalities (==0)
msups = []; % support substitutions
for k = 1:length(arg)
 m = arg{k}; 
 if isa(m,'supcon')
  for r = 1:size(m,1)
   for c = 1:size(m,2)
    p = m(r,c); t = type(p);
    lp = left(p); rp = right(p); pp = lp-rp; % polynomials
    if ~consistent(pp) % check consistency
     error('Invalid support constraint with inconsistent variables and measures')
    end
    cp = coef(lp);
    if strcmp(t, 'eq') & (length(cp)==1) & (cp(1)==1)
     % only one monomial in LHS =
     % support to be explicitly substituted
     msups = [msups struct('left',lp,'right',rp)];
    elseif strcmp(t, 'eq') % equality
     msupceq = [msupceq pp];
    elseif strcmp(t, 'ge') % inequality
     msupcge = [msupcge pp];
    else % reverse inequality
     msupcge = [msupcge -pp];
    end
    % store measure and variable indices
    newpindmeas = [indmeas(lp) indmeas(rp)];
    pindmeas = [pindmeas newpindmeas];
    pindvar = [pindvar indvar(lp) indvar(rp)];
    % degree wrt each measure
    mp = indmeas(lp);
    if mp > 0
     dp = deg(lp);
     dmax(mp) = max(dmax(mp),dp);
     rankshift(mp) = max(rankshift(mp),ceil(dp/2));
    end
    mp = indmeas(rp);
    if mp > 0
     dp = deg(rp);
     dmax(mp) = max(dmax(mp),dp);
     rankshift(mp) = max(rankshift(mp),ceil(dp/2));
    end
   end
  end
 end
end


% Detect moment constraints (MOMCON objects) and extract moment
% substitutions which are moment equality constraints
% with isolated left handside monomial

mmomcge = []; % moment inequalities (>=0)
mmomceq = []; % moment equalities (==0)
mmoms = []; % moment substitutions
for k = 1:length(arg)
 m = arg{k}; t = type(m);
 if isa(m,'momcon') & ~(strcmp(t,'min') | strcmp(t,'max'))
  for r = 1:size(m,1)
   for c = 1:size(m,2)
    p = m(r,c); t = type(p);
    lp = left(p); rp = right(p); % scalar moments
    if ~consistent(lp) | ~consistent(rp)
     error('Invalid moment constraint with inconsistent variables and measures')
    end
    lpp = split(lp); cp = coef(lpp(1));
    if strcmp(t,'eq') & (length(lpp)==1) & (length(cp)==1) & (cp(1) == 1)
     % only one monic monomial in LHS =
     % moment to be explicitly substituted
     mmoms = [mmoms struct('left',lp,'right',rp)];
    elseif strcmp(t, 'eq') % equality
     mmomceq = [mmomceq lp-rp];
    elseif strcmp(t, 'ge') % inequality
     mmomcge = [mmomcge lp-rp];
    else % reverse inequality
     mmomcge = [mmomcge -lp+rp];
    end      
    % store measure and variable indices
    newpindmeas = [indmeas(lp) indmeas(rp)];
    pindmeas = [pindmeas newpindmeas];
    pindvar = [pindvar indvar(lp) indvar(rp)];
    % degree wrt each measure
    mp = indmeas(lp);
    lp = split(lp); % vector polynomial
    for q = 1:length(lp)
     if mp(q) > 0
      dp = deg(lp(q));
      dmax(mp(q)) = max(dmax(mp(q)),dp);
     end
    end
    mp = indmeas(rp);
    rp = split(rp); % vector polynomial
    for q = 1:length(rp)
     if mp(q) > 0
      dp = deg(rp(q));
      dmax(mp(q)) = max(dmax(mp(q)),dp);
     end
    end
   end
  end
 end
end

% Some statistics

if MMM.verbose
 if isempty(mobj)
  disp('  No objective function')
 else
  disp('  Valid objective function')
 end
 fprintf('  Number of support constraints = %d',... including' ...
      length(msupcge)+length(msupceq)+length(msups));
 if length(msups) == 1
  fprintf(' including 1 substitution');
 elseif length(msups) > 1
  fprintf(' including %d substitutions',length(msups));
 end;
 fprintf('\n');
 fprintf('  Number of moment constraints = %d',... including' ...
       length(mmomcge)+length(mmomceq)+length(mmoms));
 if length(mmoms) == 1
  fprintf(' including 1 substitution');
 elseif length(mmoms) > 1
  fprintf(' including %d substitutions',length(mmoms));
 end;
 fprintf('\n');
end

% ---------------------------------
% Generate moments for each measure
% ---------------------------------

if isempty(pindmeas)
 error('No measure')
end

% Sort and remove duplicate measure indices
m = sort(pindmeas);
d = [m 0]-[0 m];
d = d(2:end-1);
i = 2:length(m);
pindmeas = m([1 i(d>0)]);

% Remove zero measure index
if (length(pindmeas) > 1) & (pindmeas(1) == 0)
 pindmeas = pindmeas(2:end);
end

if pindmeas(end) == 0
 error('No measure');
end

nmeas = length(pindmeas);

% Sort and remove duplicate variable indices
v = sort(pindvar);
d = [v 0]-[0 v];
d = d(2:end-1);
i = 2:length(v);
pindvar = v([1 i(d>0)]);

% Remove zero variable index
if (length(pindvar) > 1 ) & (pindvar(1) == 0)
 pindvar = pindvar(2:end);
end

if pindvar(end) == 0
 error('No variable');
end

% Maximum degree for each measure
dmax = dmax(pindmeas);
rankshift = rankshift(pindmeas);

% Order of SDP relaxation
if isempty(ord)
 ord = ceil(dmax/2);
else
 ord = ord*ones(size(dmax));
end

for m = 1:nmeas % relative measure index
  
  mm = pindmeas(m); % absolute measure index
  
  if MMM.verbose
   disp(['Measure ' int2str(mm)]);
  end

  if dmax(m) > 2*ord(m)
    disp('Some variables are not represented in the relaxation');
    error('Increase relaxation order'); 
  end

  % Generate moments for this measure
  nvm = genmom(mm,pindvar,ord(m));
  if nvm == 0
   error('No active variables for this measure')
  end
  
  if MMM.verbose
   disp(['  Degree = ' int2str(dmax(m))]);
   disp(['  Variables = ' int2str(MMM.M{mm}.nvar)]);
   disp(['  Moments = ' int2str(MMM.M{mm}.nm)]);
  end
 
end

disp(['Relaxation order = ' int2str(max(ord))]);

% ----------------------
% Build the SDP problem
% ----------------------

P = [];

% Algebraic constraints on moments

if isempty(mmomcge) & isempty(mmomceq) & isempty(mmoms)

 % No moment constraints so all measure masses are set to one
 for m = 1:nmeas
  mm = pindmeas(m);
  mmoms = [mmoms struct('left',mass(mm),'right',mom(1,0))];
  if MMM.verbose
   disp(['Mass of measure ' int2str(mm) ' set to one']);
  end
 end

end

% ----------------------------------------------
% Perform substitutions on the vector of moments
% in order to reduce the SDP problem

% Count monomials
for m = 1:nmeas % for each measure..
 mm = pindmeas(m);
 if m == 1
  % Starting indices in vector of all moments
  MMM.M{mm}.begind = 1;
  tnm = MMM.M{mm}.nm; % total number of moments
 else
  MMM.M{mm}.begind = tnm+1;
  tnm = tnm+MMM.M{mm}.nm;
 end
 MMM.M{mm}.indep = true(MMM.M{mm}.nm,1); % independent monomials
end

if MMM.verbose
 disp(['Total number of moments = ' int2str(tnm)]);
end

% Basis matrix of linear relations between moments
Ar = [sparse(tnm,1) speye(tnm)]; % first column = constant term
subs = zeros(tnm,1); % number of substitutions for each variable
conflict = 0; % conflicting substitutions

if ~isempty(mmoms) | ~isempty(msups)

 % --------------------
 % Support substitutions

 if ~isempty(msups)

  if MMM.verbose
   disp('Perform support substitutions');
  end

  % Scalar substitutions

  for i = 1:length(msups)
  
   % LHS
   lp = msups(i).left;
   mp = indmeas(lp); % scalar measure index
   lp = mpol(lp); % scalar polynomial
   
   % index of variable to be substituted
   lpow = locpow(lp,mp);

   % RHS
   rp = msups(i).right;

   % localization
   degp = max(deg(lp),deg(rp));
   nvar = MMM.M{mp}.nvar;
   nord = MMM.M{mp}.ord;
   % number of localization monomials
   nlm = MMM.T(nvar,nord).bin(nvar,2+2*nord-degp);
   % powers for localization
   mpow = MMM.T(nvar,nord).pow(1:nlm,:);
   % index of variables to be substituted
   lpow = reshape(lpow,1,nvar);
   lpow = reshape(repmat(lpow,nlm,1)+mpow,nlm,1,nvar);
   rind = pow2ind(lpow,mp);

   % powers for linear substitutions
   rpow = locpow(rp,mp);
   rpow = repmat(rpow,1,nlm)+...
	  repmat(reshape(mpow,1,nlm,nvar),size(rpow,1),1);
   % indices of columns for substitutions
   cind = pow2ind(rpow,mp);

   for k = 1:length(rind)
    subs(rind(k)) = subs(rind(k))+1;
    if subs(rind(k)) == 1 % new substitution
     Ar(rind(k),:) = sparse(1,tnm+1);
     Ar(rind(k),cind(:,k)+1) = coef(rp);
    end
    % otherwise discard, hoping that subs. are consistent
   end
  
   % remove linearly dependent monomials
   MMM.M{mp}.indep(rind-MMM.M{mp}.begind+1) = false;