msdp.m

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

  end % for i ..

 end % if ~isempty(msups)
 
 % --------------------
 % Moment substitutions

 if ~isempty(mmoms)

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

  % Scalar substitutions
  
  for i = 1:length(mmoms)

   % LHS
   lp = mmoms(i).left;
   mlp = indmeas(lp); % scalar measure index
   lp = split(lp); % scalar polynomial
   
   % index of variable to be substituted
   lpow = locpow(lp,mlp); rind = pow2ind(lpow,mlp);

   % RHS
   rp = mmoms(i).right;
   mrp = indmeas(rp); % vector measure indices
   rp = split(rp); % vector polynomial
   
   subs(rind) = subs(rind)+1; 
   if subs(rind) == 1 % new substitution
    Ar(rind,:) = sparse(1,tnm+1);
    for k = 1:length(rp) % linear combination of moments
     if mrp(k) == 0 % constant term     
      Ar(rind,1) = Ar(rind,1)+coef(rp(k));
     else % monomial
      % index for substitution
      rpow = locpow(rp(k),mrp(k)); cind = pow2ind(rpow,mrp(k));
      cp = coef(rp(k))';
      Ar(rind,cind+1) = Ar(rind,cind+1)+cp;
     end
    end
   else % old subs.
    % transform potentially conflicting subs. into
    % explicit moment equality constraint
    mmomceq = [mmomceq mmoms(i).left-mmoms(i).right];
    conflict = conflict+1;   
   end % if
   
  end % for i = ..

 end % if ~isempty(mmoms)
 
% char(logical(full(Ar))+char('0'))

 % Detect triangular structure
 if ~any(any(triu(Ar(:,2:end),1)))
  % lower triangular
  % propagate linear dependence relations
  for c = 1:tnm
   if Ar(c,c+1) == 0
    ind = find(Ar(:,c+1))';
    if ~isempty(ind)
     for r = ind
      k = Ar(r,c+1); Ar(r,c+1) = 0;
      Ar(r,:) = Ar(r,:) + k*Ar(c,:);
     end
    end
   end
  end
 elseif ~any(any(tril(Ar(:,2:end),-1)))
  % upper triangular
  % propagate linear dependence relations
  for r = 1:tnm
   if Ar(r,r+1) == 0
    ind = find(Ar(r,:))-1;
    if ~isempty(ind)
     for c = ind
      k = Ar(r,c+1); Ar(r,c+1) = 0;
      Ar(:,c+1) = Ar(:,c+1) + k*Ar(:,r+1);
     end
    end
   end
  end
 else
  error(['Monomial substitutions should have triangular structure' 13 ...
  'Modify accordingly the sequence of substitutions' 13 ...
  'or replace substitutions constraints with equality constraints']);
 end

 if MMM.verbose
  if conflict > 0
   warning([int2str(conflict) ' potentially conflicting moment substitutions' 13 ...
           'replaced with explicit moment equality constraints'])
  end
 end

end % if moment substitutions

% Global mask of independent monomials
indep = 1+find(sum(abs(Ar(:,2:tnm+1))));
nvarlmi = length(indep);

% Linear dependence patterns of monomials for each measure
for m = 1:nmeas % relative measure index
 mm = pindmeas(m); % absolute measure index
 MMM.M{mm}.Ar = Ar(~MMM.M{mm}.indep,:);
end

if MMM.verbose
 if tnm ~= nvarlmi
  disp(['Number of moments after substitutions = ' ...
       int2str(nvarlmi)]);
 end
end

% -----------
% Constraints

% first column = constant term
Af = sparse(0,tnm+1); % equality constraints
Al = sparse(0,tnm+1); % scalar inequality constraints
As = sparse(0,tnm+1); % SDP constraints
K = struct('f',0,'l',0,'q',0,'s',[]); % cone structure

% indices for retrieving dual variables
indmomeq = []; % moment equality constraints
indmomge = []; % moment inequality constraints

if MMM.verbose
 disp('Generate moment and support constraints')
end

% indices in moment matrix
indmmat = zeros(nmeas,2); % indices of moment matrices in z

% Moment matrix SDP constraints

for m = 1:nmeas % relative measure index
 
 mm = pindmeas(m); % absolute measure index
 
 nord = MMM.M{mm}.ord;
 nvar = MMM.M{mm}.nvar;
 nmm = MMM.T(nvar,nord).bin(end,nord+2); 
 ind = MMM.M{mm}.indep(1:nmm); % only independent monomials
 nmm = sum(ind); % size of moment matrix
 nmm2 = nmm^2; % number of entries in moment matrix
 mpow = MMM.T(nvar,nord).bas(ind,ind,:); % powers in 3D matrix
 cind = pow2ind(mpow,mm); % indices in vector of moments
 % update SDP constraints
 As = [As;sparse(1:nmm2,cind(:)+1,ones(1,nmm2),nmm2,tnm+1)];
 K.s = [K.s nmm];
 indmmat(m,:) = size(As,1)-[nmm2-1 0];

end

% Bound on the sum of all moments

% if MMM.bound
%  disp(['Bound on the sum of all moments = ' num2str(MMM.bound)]);
%  Alrow = [MMM.bound -ones(1,tnm)];
%  % update inequality constraints
%  indmomeq = [indmomeq; size(Al,1)+[1 size(Alrow,1)]];
%  Al = [Al; Alrow];
%  K.l = K.l + 1;   
% end
 
% Equality moment constraints
 
if ~isempty(mmomceq)
   
  for i = 1:length(mmomceq)
     
   p = split(mmomceq(i)); % scalar moment to vector polynomial
   mp = indmeas(mmomceq(i)); % measure indices
   Afrow = sparse(1,tnm+1); % new equality constraint in A

   for k = 1:length(p) % linear combination of moments
   
    dp = deg(p(k));
    cp = coef(p(k))'; % row vector
    if mp(k) == 0 % scalar constant
     Afrow(1) = Afrow(1) + cp;
    else % polynomial
     nmon = length(cp);
     ppow = locpow(p(k),mp(k)); % powers
     cind = pow2ind(ppow,mp(k)); % indices
     Afrow(cind+1) = Afrow(cind+1) + cp; % constraints
    end;
    
   end

   % update equality constraints
   indmomeq = [indmomeq; size(Af,1)+[1 size(Afrow,1)]];
   Af = [Af; Afrow];
   K.f = K.f + 1;
  end
   
end
 
% Inequality moment constraints
 
if ~isempty(mmomcge)
  
  for i = 1:length(mmomcge)
     
   p = split(mmomcge(i)); % scalar moment to vector polynomial
   mp = indmeas(mmomcge(i)); % measure indices
   Alrow = sparse(1,tnm+1); % new inequality constraint in A
   
   for k = 1:length(p) % linear combination of moments
   
    dp = deg(p(k));
    cp = coef(p(k))'; % row vector
    if mp(k) == 0 % scalar constant
     Alrow(1) = Alrow(1) + cp;
    else
     nmon = length(cp);
     ppow = locpow(p(k),mp(k)); % powers
     cind = pow2ind(ppow,mp(k)); % indices
     Alrow(cind+1) = Alrow(cind+1) + cp; % constraints
    end
    
   end
   
   % update inequality constraints
   indmomge = [indmomge; size(Al,1)+[1 size(Alrow,1)]];
   Al = [Al; Alrow];
   K.l = K.l + 1;
   
  end
   
end
 
% Localize equality support constraints

if ~isempty(msupceq)
  
  for i = 1:length(msupceq)
    
   p = msupceq(i); % scalar polynomial

   dp = deg(p);
   mp = indmeas(p);
   cp = coef(p);
   if (mp == 0) & (cp ~= 0)
    error('Inconsistent constant support constraint')
   end
   nmon = size(cp,1);
   ppow = locpow(p,mp);
   nvar = MMM.M{mp}.nvar;
   nord = MMM.M{mp}.ord;
   nlm = MMM.T(nvar,nord).bin(nvar,2+2*nord-dp); % size of loc. mat.
   ind = MMM.M{mp}.indep(1:nlm); % keep independent monomials only
   nlm = sum(ind); % number of localization monomials
   % 2D powers for localization
   mpow = MMM.T(nvar,nord).pow(ind,:);
   mpow = reshape(mpow,1,nlm,nvar); % 3D columnwise
   ppow = repmat(mpow,nmon,1)+repmat(ppow,1,nlm); % localize
   cind = pow2ind(ppow,mp); % indices in vector of moments
   cind = reshape(cind,nmon*nlm,1); % columwise
   rind = repmat(1:nlm,nmon,1); % row indices
   cp = repmat(cp,1,nlm); % coefs
   % update equality constraints
   Af = [Af;sparse(rind(:),cind+1,cp(:),nlm,tnm+1)];
   K.f = K.f+nlm;

  end 

end

% Localize inequality support constraints
 
if ~isempty(msupcge)
 
  for i = 1:length(msupcge)

   p = msupcge(i); % scalar polynomial
   
   dp = deg(p);
   mp = indmeas(p);
   cp = coef(p);
   if (mp == 0) & (cp < 0)
    error('Inconsistent constant support constraint')
   end
   nmon = size(cp,1);
   ppow = locpow(p,mp);
   nvar = MMM.M{mp}.nvar;
   nord = MMM.M{mp}.ord;
   nlm = MMM.T(nvar,nord).bin(nvar,2+(nord-ceil(dp/2))); % size of loc. mat.
   ind = MMM.M{mp}.indep(1:nlm); % keep independent monomials
   nlm = sum(ind); % number of localization monomials
   mpow = MMM.T(nvar,nord).bas(ind,ind,:); % 3D powers for localization
   nlm2 = nlm^2; % total number of elements (with repeated entries)
   mpow = reshape(mpow,1,nlm2,nvar); % 3D columnwise
   ppow = repmat(mpow,nmon,1)+repmat(ppow,1,nlm2); % localize
   cind = pow2ind(ppow,mp); % indices in vector of moments
   cind = reshape(cind,nmon*nlm2,1); % columwise
   rind = repmat(1:nlm2,nmon,1); % row indices
   cp = repmat(cp,1,nlm2); % coefs
   % update SDP constraints
   if nlm2 == 1 % scalar SDP = linear constraint
    Al = [Al;sparse(rind(:),cind+1,cp(:),nlm2,tnm+1)];
    K.l = K.l+nlm2;
   else % matrix SDP
    As = [As;sparse(rind(:),cind+1,cp(:),nlm2,tnm+1)];
    K.s = [K.s nlm];
   end 
  
  end
  
end

% Objective function

objshift = 0; % constant shift
obj = sparse(1,tnm);

if ~isempty(mobj)

  if indvar(mobj) | indmeas(mobj) % not a constant
    
   mo = indmeas(mobj); % vector measure indices
   po = split(mobj); % vector polynomial
  
   for k = 1:length(po) % linear combination of moments
    if mo(k) == 0 % constant term
     objshift = objshift + double(po(k));
    else % monomial
     % powers to indices
     powo = locpow(po(k),mo(k)); cind = pow2ind(powo,mo(k));
     % store coeffs
     obj(cind) = coef(po(k));
    end
   end
  
  else % constant
    
   if MMM.verbose
    disp('Constant objective function: generate feasibility problem')
   end
    
   objshift = coef(split(mobj));

  end
 
else % no objective function
     % so minimize trace of moment matrix of first measure 

 if MMM.verbose
  disp('No objective function')
  disp(['Minimize trace of moment matrix of measure ' ...
	int2str(pindmeas(1))])
 end

 mo = pindmeas(1); % first measure
 powb = MMM.T(MMM.M{mo}.nvar,MMM.M{mo}.ord).bas; % 3D powers
 nr = size(powb,1);
 powo = zeros(nr,1,size(powb,3)); 
 for r = 1:nr % extract diagonal entries
  powo(r,1,:) = powb(r,r,:); % in 3D power vector
 end
 cind = pow2ind(powo,mo); % indices
 obj(cind) = 1; % all ones
 
end


% Reduce constraints
A = [Af;Al;As];
c = A(:,1)+A(:,2:end)*Ar(:,1);
A = A(:,2:end)*Ar(:,indep);
b = obj*Ar(:,indep);
objshift = objshift + obj*Ar(:,1);

% Update indices
indmmat = indmmat + K.f + K.l + K.q;
indmomge = indmomge + K.f;

% Display some info on SDP problem

if MMM.verbose
 disp('Generate moment SDP problem')
end

% --------------------------------------
% Parameters relative to the SDP problem
% --------------------------------------

% SDP problem data in SeDuMi dual format
% max b'*y s.t. z = c-A*y \in K
P.A = -A; P.b = b; P.c = c; P.K = K;

% Linearly independent monomials
P.indep = indep;
% Dependence pattern
P.Ar = Ar;

% Indices of measures
P.indmeas = pindmeas;
% Indices of moment matrices in z
P.indmmat = indmmat;
% Indices of moment equality constraints in z
P.indmomeq = indmomeq;
% Indices of moment inequality constraints in z
P.indmomge = indmomge;

% Support constraints used to check feasibility in MSOL
P.msupceq = msupceq;
P.msupcge = msupcge;
% Moment constraints used when retrieving multipliers
P.mmomceq = mmomceq;
P.mmomcge = mmomcge;
% Data used to check optimality in MSOL
P.mobj = mobj; % objective function
P.objshift = objshift; % constant shift
P.objsign = objsign; % max = +1, min = -1
% Rank shift to detect global optimality in MSOL
P.rankshift = rankshift;
% SDP relaxation orders
P.order = ord;

% Numerical value of the objective function
P.sdpobj = [];

% Measures in vectorized MEAS format
M = [];
for k = 1:nmeas
 M = [M meas(pindmeas(k))];
end

% Constructor
P = class(P,'msdp');

% End of main function MDEF
% -----------------------------------------------------------------

% Auxillary functions