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)
Dual variables
The dual (dual in the sense that it is the dual related to a user
defined constraint) 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 often 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 form representation, corresponding to 471 equality
constraint, 2 semidefinite cones and 3 free variables in the
primal form. If we instead would have solved this
directly in the stated primal form, we have 3 equality
constraints, 2 semidefinite cones and no free variables,
corresponding to a dual form with 3 variables and two
semidefinite constraints. The computational effort is mainly
affected by the number of variables in the dual form and the size of the
semidefinite cones. Moreover, many solvers have robustness
problems with free variables in the primal form (equalities in the
dual form). Hence, in this case, this problem can probably be solved
much more efficiently if we could use an alternative model. 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);assign(primals{i},dual(Fd(i)));end
|
Variables are actually automatically updated, so the second
line in the code above is not needed (but can be useful to
understand what is happening). Hence, the following code is
equivalent.
The procedure 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);
solvesdp(F,obj);
double(t)
ans =
-1.9099
0.7358
[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 4th output, and
can be recovered from the dual to the equality constraints (this
is also done automatically by YALMIP in practice, see above).
solvesdp(Fd,-objd);
assign(free,dual(Fd(end)))
double(t)
ans =
-1.9099
0.7358
|
To simplify things even further, you can tell YALMIP to
dualize, solve the dual, and recover the primal variables, by
using the associated option.
solvesdp(F,obj,sdpsettings('dualize',1));
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 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.
A problem
involving translated cones X≥C where C is a
constant is automatically converted to a problem in standard primal
form, with no additional slacks variables. Hence, a lower bound on a
variable will typically reduce the size of a dualized problem, since
no free variables or slacks will be needed to model this cone.
Practice has shown that simple bound constraints of the type x≥L
where L is a large negative number can lead to problems if one tries
to perform the associated variable change in order to write it as a
simple LP cone. Essentially, the dual cost will contain large
numbers. With a primal problem {min
cTx, Ax=b, x≥-L}) will be converted to {min cTz,
Az=b+AL, z≥0}) with the dual {min (b+AL)Ty, ATy≤c}.
If you want to avoid detection of translated LP cones (and thus
treat the involved variables as free variables), set the 4th
argument in dualize.
Problems
involving second order cone constraints can also be dualized. A constraint
of the type x = sdpvar(n,1);F = set(cone(x(2:end),x(1))); is a second
order constraint in standard primal form. If your cone constraint violates this
form, slacks will be introduced, except for translated second order
cones, just as in the semidefinite case. Note
that you need a primal-dual solver that can solve second order cone
constraints natively in order to recover the original variables
(currently
SeDuMi and
SDPT3 for mixed semidefinite
second order cone problems, and
MOSEK for pure second order cone
problems).
Your solver
has to be able to return both primal and dual variables for the reconstruction of
variables to work. All SDP solvers support this, except
LMILAB.
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);
assign(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 no reason
to primalize a model.
There are however some cases where it may make sense. Consider the
problem in the KYP example
n = 50;
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 1276
equality constraints and 1326 variables is obtained. This means that
the effective number of variables is low (the degree of freedom).
[Fp,objp] = primalize(F,-obj);Fp
+++++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint| Type|
+++++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Numeric value| Matrix inequality 51x51|
| #2| Numeric value| Equality constraint 1276x1|
+++++++++++++++++++++++++++++++++++++++++++++++++++++
length(getvariables(Fp))
ans =
1326
|
For comparison, let us first solve the original problem.
solvesdp(F,obj)
ans =
yalmiptime: 0.2410
solvertime: 32.4150
info: 'No problems detected (SeDuMi)'
problem: 0
|
The primalized takes approximately the same time to solve (this
can differ between problem instances though).
solvesdp(Fp,objp)
ans =
yalmiptime: 0.3260
solvertime: 32.2530
info: 'No problems detected (SeDuMi)'
problem: 0
|
So why would we want to perform the primalization? We let YALMIP
remove the equalities constraints first!
solvesdp(Fp,objp,sdpsettings('removeequalities',1))
ans =
yalmiptime: 2.6240
solvertime: 1.1860
info: 'No problems detected (SeDuMi)'
problem: 0
|
The drastic reduction in actual solution-time of the semidefinite
program comes at a price. Removing the equality constraints and
deriving a reduced basis with a smaller number of variables requires
computation of a QR factorization of a matrix of
dimension 1326 by 1276. The cost of doing this is roughly 2 seconds.
The total saving in computation time is however still high enough to
motivate the primalization. |