dualize.m

optimization toolbox

function [Fdual,objdual,X,t,err] = dualize(F,obj,auto,extlp,extend)
% DUALIZE Create the dual of an SDP given in primal form
%
% [Fd,objd,X,t,err] = dualize(F,obj,auto)
%
% Input
%  F    : Primal constraint in form AX=b+dt, X>0, t free.
%  obj  : Primal cost CX+ct
%  auto : If set to 0, YALMIP will not automatically handle variables
%        and update variable values when the dual problem is solved.
%  extlp: If set to 0, YALMIP will not try to perform variables changes in
%         order to convert simple translated LP cones (as in x>1) to
%         standard unit cone constraints (x>0)
%
% Output
%  Fd  : Dual constraints in form C-A'y>0, c-dy==0
%  obj : Dual cost b'y (to be MAXIMIZED!)
%  X   : The detected primal cone variables
%  t   : The detected primal free variables
%  err : Error status (returns 0 if no problems)
%
% Example
%  See the HTML help.
%
% See also DUAL, SOLVESDP, PRIMALIZE

% Author Johan L鰂berg
% $Id: dualize.m,v 1.53 2006/11/21 16:10:23 joloef Exp $

% Check for unsupported problems

if isempty(F)
    F = set([]);
end

if nargin < 2
    obj = [];
end

err = 0;
p1 = ~(isreal(F) & isreal(obj));
p2 = ~(islinear(F) & islinear(obj));
p3 = any(is(F,'integer')) | any(is(F,'binary'));
if p1 | p2 | p3
    if nargout == 5
        Fdual = set([]);objdual = [];y = []; X = []; t = []; err = 1;
    else
        problems = {'Cannot dualize complex-valued problems','Cannot dualize nonlinear problems','Cannot dualize discrete problems'};
        error(problems{min(find([p1 p2 p3]))});
    end
end

if nargin<5
    extend = 1;
end

if extend
    options.dualize = 1;
    options.allowmilp = 0;
    options.solver = '';
    [F,failure,cause] = expandmodel(F,obj,options);
    if failure
        error('Failed during convexity propagation. Avoid nonlinear operators when applying dualization.');
    end
end

if nargin<3
    auto = 1;
end

if nargin<4
    extlp = 1;
end

% Cones and equalities
F_AXb = set([]);


% Shiftmatrix is a bit messy at the moment.
% We want to be able to allow cones X>shift
% by using a new variable X-shift = Z
shiftMatrix = {};
X={};

% First, get variables in initial SDP cones
% We need this to avoid adding the same variable twice
% when we add simple LP constraints (as in P>0, P(1,3)>0)
varSDP = [];
SDPset = zeros(length(F),1);
isSDP = is(F,'sdp');
for i = 1:length(F)
    if isSDP(i);
        Fi = sdpvar(F(i));
        if is(Fi,'shiftsdpcone')
            vars = getvariables(Fi);
            if isempty(findrows(varSDP,[vars(1) vars(end)]))
                SDPset(i) = 1;
                varSDP = [varSDP;vars(1) vars(end)];
                shiftMatrix{end+1} = getbasematrix(Fi,0);
                X{end+1}=Fi;
            end
        end
    end
end
F_SDP = F(find(SDPset));
F = F(find(~SDPset));

% Same thing for second order cones
% However, we must not add any SOC cones
% that we already defined as SDP cones
varSOC = [];
SOCset = zeros(length(F),1);
isSOCP = is(F,'socp');
for i = 1:length(F)
    if isSOCP(i);%is(F(i),'socp')
        Fi = sdpvar(F(i));
        if is(Fi,'socone')
            vars = getvariables(Fi);
            % Make sure these variables are not SDP cone variables
            % This can actually only happen for set(X>0) + set(Xcone((2:end,1),X(1)))
            if ~isempty(varSDP)
                inSDP = any(varSDP(:,1)<=vars(1)& vars(1) <=varSDP(:,2)) | any(varSDP(:,1)<=vars(end)& vars(end) <=varSDP(:,2));
            else
                inSDP = 0;
            end
            if ~inSDP
                SOCset(i) = 1;
                vars = getvariables(Fi);
                varSOC = [varSOC;vars(1) vars(end)];
                shiftMatrix{end+1} = getbasematrix(Fi,0);
                X{end+1}=Fi;
            end
        end
    end
end
F_SOC = F(find(SOCset));
F = F(find(~SOCset));

% Merge SDP and SOC data
varCONE = [varSDP;varSOC];
F_CONE = F_SDP + F_SOC;

% Find all LP constraints, add slacks and extract simple cones
% to speed up things, we treat LP cone somewhat different
% compared to the conceptually similiar SOCP/SDP cones
% This code is pretty messy, since there are a lot off odd
% cases to take care of (x>0 and x>1 etc etc)
elementwise = is(F,'element-wise');
elementwise_index = find(elementwise);
if ~isempty(elementwise_index)

    % Find element-wise inequalities
    Flp = F(elementwise_index);
    F = F(find(~elementwise)); % remove these LPs
    % Find LP cones
    lpconstraint = [];
    for i = 1:length(Flp)
        temp = sdpvar(Flp(i));
        if min(size(temp))>1
            temp = temp(:);
        end
        lpconstraint = [lpconstraint reshape(temp,1,length(temp))];
    end

    % Find all constraints of type a_i+x_i >0 and extract the unique and
    % most constraining inequalities (i.e. remove redundant lower bounds)
    base = getbase(lpconstraint);
    candidates = find(sum(base(:,2:end)~=0,2)==1);
    if length(candidates)>0
        % The other ones...
        alwayskeep = find(sum(base(:,2:end)~=0,2)~=1);
        w1 =  lpconstraint(alwayskeep);
        w2 =  lpconstraint(candidates);
                  
        % Find unique rows
        base = getbase(w2);
        [i,uniquerows,k] = unique(base(:,2:end)*randn(size(base,2)-1,1));   
        aUniqueRow=k(:)';
        keep = [];
        rhsLP = base(:,1);
        rr = histc(k,[1:max(k)]);
        if all(rr==1)
            lpconstraint = [w1 w2];
        else
            for i=1:length(k)
                sameRow=find(k==k(i));
                if length(sameRow)==1
                    keep = [keep sameRow];
                else
                    rhs=base(sameRow,1);
                    [val,pos]=min(rhsLP(sameRow));
                    keep = [keep sameRow(pos)];
                end
            end
            lpconstraint = [w1 w2(unique(keep))];
        end
       
    end

    % LP cone will be saved in a vector for speed
    x = [];

    % Pure cones of the type x>0
    base = getbase(lpconstraint);
    purelpcones = (base(:,1)==0) & (sum(abs(base(:,2:end)),2)==1) & (sum(base(:,2:end)==1,2)==1);
    if ~isempty(purelpcones)
        x  = [x lpconstraint(purelpcones)];      
        lpconstraint = lpconstraint(find(~purelpcones));
    end

    
    % Translated cones x>k, k positive
    % User does not want to make variable changes based on k
    % But if k>=0, we can at-least mark x as a simple LP cone variable and
    % thus avoid a free variable.
    if ~extlp & ~isempty(lpconstraint)
        base = getbase(lpconstraint);
        lpcones = (base(:,1)<0) & (sum(abs(base(:,2:end)),2)==1) & (sum(base(:,2:end)==1,2)==1);
        if ~isempty(find(lpcones))
            s = recover(getvariables(lpconstraint(find(lpcones))));
            x  = [x reshape(s,1,length(s))];
        end
    end

    % Translated cones x>k
    % Extract these and perform the associated variable change y=x-k
    if ~isempty(lpconstraint)%Avoid warning in 5.3.1
       base = getbase(lpconstraint);
       lpcones = (sum(abs(base(:,2:end)),2)==1) & (sum(base(:,2:end)==1,2)==1);
       if ~isempty(lpcones) & extlp
          x  = [x lpconstraint(find(lpcones))];
          nlp = lpconstraint(find(~lpcones));
          if ~isempty(nlp)
             s = sdpvar(1,length(nlp));
             F_AXb = F_AXb + set(nlp-s==0);
             x = [x s];
          end
       elseif length(lpconstraint) > 0
          s = sdpvar(1,length(lpconstraint));
          x = [x s]; % New LP cones
          F_AXb = F_AXb + set(lpconstraint-s==0);
       end
    end
    
    
    % Sort asccording to variable index
    % (Code below assumes x is sorted in increasing variables indicies)
    base = getbase(x);base = base(:,2:end);[i,j,k] = find(base);
    x = x(i);
    xv = getvariables(x);
        
    % For mixed LP/SDP problems, we must ensure that LP cone variables are
    % not actually an element in a SDP cone variable  
    if ~isempty(varCONE)        
        keep = zeros(length(x),1);
        for i = 1:length(xv)
            if any(varCONE(:,1)<=xv(i) & xv(i) <=varCONE(:,2))
            else
                keep(i) = 1;
            end
        end
        if ~all(keep)
            % We need to add some explicit constraints on some elements and
            % remove the x variables since they are already in a cone
            % variable
            xcone = x(find(~keep));
            s = sdpvar(1,length(xcone));
            F_AXb = F_AXb + set(xcone-s==0);
            x = x(find(keep));
            x = [x s];
        end
    end
else
    x = [];
end

% Check for mixed cones, ie terms C-A'y > 0.
keep = ones(length(F),1);
isSDP  = is(F,'sdp');
isSOCP = is(F,'socp');

% Pre-allocate all SDP slacks in one call
% This is a lot faster
if nnz(isSDP) > 0
    SDPindicies = find(isSDP)';
    for i = 1:length(SDPindicies)%find(isSDP)'
        Fi = sdpvar(F(SDPindicies(i)));
        ns(i) = size(Fi,1);
        ms(i) = ns(i);
    end
    Slacks = sdpvar(ns,ms);
    if ~isa(Slacks,'cell')
        Slacks = {Slacks};
    end
end
prei = 1;
for i = 1:length(F)
    if isSDP(i)
        % Semidefinite dual-form cone
        Fi = sdpvar(F(i));
        n  = size(Fi,1);
 %      S  = sdpvar(n,n);
        S  = Slacks{prei};prei = prei + 1;
        slack = Fi-S;
        ind = find(triu(reshape(1:n^2,n,n)));
        F_AXb =  F_AXb + set(slack(ind)==0);
        F_CONE = F_CONE + lmi(S,[],[],[],1);
        shiftMatrix{end+1} = spalloc(n,n,0);
        X{end+1}=S;
        keep(i)=0;
    elseif isSOCP(i)
        % SOC dual-form cone
        Fi = sdpvar(F(i));
        n  = size(Fi,1);
       S  = sdpvar(n,1);
%        S = Slacks{i};
        slack = Fi-S;
        F_AXb =  F_AXb + set(slack==0);
        F_CONE = F_CONE + set(cone(S(2:end),S(1)));
        shiftMatrix{end+1} = spalloc(n,1,0);
        X{end+1}=S;
        keep(i)=0;
    end
end

% Now, remove all mixed cones...
F = F(find(keep));

% Get the equalities
AXbset = is(F,'equality');
if any(AXbset)
    % Get the constraints
    F_AXb = F_AXb + F(find(AXbset));
    F = F(find(~AXbset));
end

% Is thee something we missed in our tests?
if length(F)>0
    error('DUALIZE can only treat standard SDPs (and LPs) at the moment.')
end