Duality


Problems in YALMIP are internally written in the following format (this will be referred to the dual form, or dual type representation)

The dual to this problem is (called the primal form)

Some solvers are capable of treating an extension with equality constraints in the dual form, corresponding to free (unconstrained) variables in the primal form.

The dual (in our notation) variable X can be obtained using YALMIP. Consider the following version of the Lyapunov stability example (of-course, dual variables in LP, QP and SOCP problems can also be extracted)

F = set(P > eye(n),'Normalize');
F = F + set(A'*P+P*A < 0,'Lyapunov');
solution = solvesdp(F,trace(P));

The dual variables related to the constraints P>I and ATP+PA< 0 can be obtained by using the command dual and indexing of lmi-objects.

Z1 = dual(F('Normalize'))
Z2 = dual(F('Lyapunov'))

Standard indexing can also be used.

Z1 = dual(F(1))
Z2 = dual(F(2))

Complementary slackness can easily be checked since double is overloaded on lmi-objects..

trace(dual(F(1))*double(F(1)))
trace(dual(F(2))*double(F(2)))

Notice, double(F(1)) returns double(0-(A'*P+P*A)).

Dualize

Important to note is that problems in YALMIP are modeled internally in the dual format (your primal problem is in dual form). In control theory and many other fields, this is the natural representation (we have a number of variables on which we have inequality constraints), but in some fields (combinatorial optimization), the primal form is more natural.

Due to the choice to work in the dual form, some problems are treated very inefficiently in YALMIP. Consider the following problem in YALMIP.
X = sdpvar(30,30);
Y = sdpvar(3,3);
obj = trace(X)+trace(Y);
F = set(X>0) + set(Y>0);
F = F + set(X(1,3)==9) + set(Y(1,1)==X(2,2)) + set(sum(sum(X))+sum(sum(Y)) == 20)
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                      Type|
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|   Matrix inequality 30x30|
|   #2|   Numeric value|     Matrix inequality 3x3|
|   #3|   Numeric value|   Equality constraint 1x1|
|   #4|   Numeric value|   Equality constraint 1x1|
|   #5|   Numeric value|   Equality constraint 1x1|
+++++++++++++++++++++++++++++++++++++++++++++++++++

YALMIP will explicitly parameterize the variable X with free 465 variables, Y with 6 free variables, create two semidefinite constraints and introduce 3 equality constraints in the dual representation, corresponding to  471 equality constraint, 2 semidefinite cones and 3 free variables in the primal of this dual.  If we instead would have solved this directly in the stated primal form, we would have had 3 equality constraints, 2 semidefinite cones and no free variables, corresponding to a dual problem with 3 variables and two semidefinite constraints. The computational effort is mainly affected by the number of dual variables and the size of the semidefinite cones. Moreover, many solvers have robustness problems with free variables in the primal (equalities in the dual). Hence, in this case, this problem can probably be solved much more efficiently in the alternative form instead.

The command dualize can be used to extract the primal form, and return the dual of this problem in YALMIPs preferred dual form.
[Fd,objd,primals] = dualize(F,obj);Fd
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                      Type|
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|   Matrix inequality 30x30|
|   #2|   Numeric value|     Matrix inequality 3x3|
+++++++++++++++++++++++++++++++++++++++++++++++++++

If we solve this problem in dual form, the duals to the constraints in Fd will correspond to the original variables X and Y. The optimal values of these variables can be reconstructed easily (note that the dual problem is a maximization problem)
solvesdp(Fd,-objd);
for i = 1:length(primals);setsdpvar(primals{i},dual(Fd(i)));end

The function can be applied also to problems with free variables in the primal form, corresponding to equality constraints in the dual form.

X = sdpvar(2,2);
t = sdpvar(2,1);
Y = sdpvar(3,3);
obj = trace(X)+trace(Y)+5*sum(t);

F = set(sum(X) == 6+pi*t(1)) + set(diag(Y) == -2+exp(1)*t(2))
F = F + set(Y>0) + set(X>0);
[Fd,objd,primals,free] = dualize(F,obj);Fd
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                      Type|
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|     Matrix inequality 3x3|
|   #2|   Numeric value|     Matrix inequality 2x2|
|   #3|   Numeric value|   Equality constraint 2x1|
+++++++++++++++++++++++++++++++++++++++++++++++++++

The detected free variables are returned as the 5th output. These variables can be recovered from the dual of the equality constraints in the solved problem.

solvesdp(Fd,-objd);
setsdpvar(free,dual(Fd(end)))
double(t)
ans =
   -1.9099
    0.7358

 

Mixed problems can be dualized also, i.e. problems involving constraints of both dual and primal form. Constraint in dual form S(y)≥0 are automatically changed to S(y)-X=0, X≥0, and the dualization algorithm is applied to this new problem. Note though that problems involving dual form semidefinite constraints typically not gain from being dualized, unless the dual form constraints are few and small compared to the primal form constraints.

Your solver has to be able to return dual variables for the reconstruction of variables to work. All SDP solvers support this, except LMILAB. The solver DSDP returns dual variables, but they are typically of poor quality.

Primal matrices (X and Y in the examples above) must be defined in one simple call in order to enable detection of the primal structure. In other words, a constraint set(X>0) where X is defined with the code x = sdpvar(10,1);X = [x(1) x(6);x(6) x(2)] will not be categorized as a primal matrix, but as matrix constraint in dual form with three free variables.

Primalize

For completeness, a functionality called primalize is available. This function takes an optimization problem in dual form and returns a YALMIP model in primal form. Consider the following SDP with 3 free variables, 1 equality constraint, and 1 SDP constraint of dimension 2.
C = eye(2);
A1 = randn(2,2);A1 = A1*A1';
A2 = randn(2,2);A2 = A2*A2';
A3 = randn(2,2);A3 = A3*A3';
y = sdpvar(3,1);

obj = -sum(y) % Maximize sum(y) i.e. minimize -sum(y)
F = set(C-A1*y(1)-A2*y(2)-A3*y(3) > 0) + set(y(1)+y(2)==1)

A model in primal form is (note the negation of the objective function, primalize assumes the objective function should be maximized)

[Fp,objp,free] = primalize(F,-obj);Fp
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                      Type|
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|     Matrix inequality 2x2|
|   #2|   Numeric value|   Equality constraint 3x1|
+++++++++++++++++++++++++++++++++++++++++++++++++++

The problem can now be solved in the primal form, and the original variables are reconstructed from the duals of the equality constraints (placed last). Note that the primalize function returns an objective function that should be minimized.
solvesdp(Fp,objp);
setsdpvar(free,dual(Fp(end));

Why not dualize the primalized model!

[Fd,objd,X,free] = dualize(Fp,objp);Fd
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                      Type|
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|     Matrix inequality 2x2|
|   #2|   Numeric value|   Equality constraint 1x1|
+++++++++++++++++++++++++++++++++++++++++++++++++++

The model obtained from the primalization is most often more complex than the original model, so there is typically not any reason to primalize a model.

There are however some cases where it may make sense. Consider the problem in the KYP example

n = 30;
A = randn(n);A = A - max(real(eig(A)))*eye(n)*1.5; % Stable dynamics
B = randn(n,1);
C = randn(1,n);

t = sdpvar(1,1);
P = sdpvar(n,n);

obj = t;
F = set(kyp(A,B,P,blkdiag(C'*C,-t)) < 0)

The original problem has 466 variables and one semidefinite constraint. If we primalize this problem, a new problem with 466 equality constraints and 496 variables is obtained.

[Fp,objp] = primalize(F,-obj);Fp
+++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                        Type|
+++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|     Matrix inequality 31x31|
|   #2|   Numeric value|   Equality constraint 466x1|
+++++++++++++++++++++++++++++++++++++++++++++++++++++
length(getvariables(Fp))
ans =
   496

For comparison, let us first solve the original problem.

solvesdp(F,obj)
ans = 
    yalmiptime: 0.2410
    solvertime: 52.4150
          info: 'No problems detected (SeDuMi)'
       problem: 0

The primalized problem is solved roughly twice as fast here (this can differ a lot between problem instances though).

solvesdp(Fp,objp)
ans = 
    yalmiptime: 0.3260
    solvertime: 32.2530
          info: 'No problems detected (SeDuMi)'
       problem: 0

Although we see some improvement here in terms of computational effort needed, the majority of solvers perform bad on this problem due to the equality constraints. However, if we let YALMIP remove the equalities constraints first, most solvers perform very well.

solvesdp(Fp,objp,sdpsettings('removeequalities',1))
ans = 
    yalmiptime: 10.9260
    solvertime: 2.2530
          info: 'No problems detected (SeDuMi)'
       problem: 0

The drastic reduction in actual solution-time of the semidefinite program of-course comes at a price. Removing the equality constraints and deriving a reduced basis with a smaller number of variables requires the computation of a sparse QR factorization of a matrix of dimension 496 by 466.