Working with the dual in practice

Recently we have been discussing about the use of dual problem: in our last post  we showed that building the dual can be a fairly automatized process; we also argued that there are at least two good reasons why we may want/need to work with it:

1. solving the dual may be more efficient (see this example);
2. the dual variables relate to the reduced coefficients and shadow costs (see this post);

So either we solve the primal and we would like to have the values of the dual variables, or the other way round! In both cases, we boil down to the question (remember the dual of the dual is the primal):

Given the solution of a linear program, how can we recover its dual variables?

We will discuss how to work on the dual solution for linear programs, but most of our arguments will apply to conic programs as well. Remember our primal in standard form:

$$\begin{equation}\tag{P} \begin{aligned} & \min \quad c^T x\\ & s.t.\\ & Ax -b =0\\ & x\geq 0, \end{aligned} \end{equation}$$

and the corresponding dual

$$\begin{equation}\tag{D} \begin{aligned} & \max \quad b^Ty\\ & s.t.\\ & c - A^Ty - s =0\\ & s\geq 0. \end{aligned} \end{equation}$$

For a given optimal primal solution $x^\star$ and dual $y^\star,s^\star$, duality theory tells us:

(1) Strong Duality: the optimal values of (P) and (D) coincide

(2) Complementary Slackness: it must hold  $x_i^\star s_i^\star = 0$ for all $i=1,\ldots,n.$ Thus $s_i=0$ for all $i=1,\ldots,n$ such that $x_i^\star>0$.

(3) Optimal Basis: using the partition of the primal variables induced by the optimal basis, we can define a full rank sub-matrix of the constraints coefficient matrix, and use it to compute $y^\star$...but this goes way to far for a simple blog post, you should refer for instance to Section 3.5 in [1].

If points (1) and (2) are straightforward, point (3) involves some work and attention....should we really spend time on that? The answer is NO! The reason is solvers usually provide some functionality to retrieve the dual variables after the primal have been solved:

1. Primal-dual linear solvers (such as MOSEK) solve the primal and the dual at the same time, and thus provide the dual variables directly.

2. Primal or Dual solvers based on the Simplex method have all the information to recover the dual variables easily basis (see again [1]).

So you just need to query somehow the solver that you use to retrieve the dual values. Let's play with MOSEK and a toy linear programming problem taken from [2] Section 6.2. In LP format the primal is

	\production mix - primal
	maximize
	obj: 550 x1 + 600 x2 + 350 x3 + 400 x4 + 200 x5
	subject to
	grinding: 12 x1 + 20 x2 + 25 x4 + 15 x5<= 288
	drilling: 10 x1 + 8 x2 + 16 x3 <= 192
	manpower: 20 x1 + 20 x2 + 20 x3 + 20 x4 + 20 x5<= 384
	bounds
	0<= x1
	0<= x2
	0<= x3
	0<= x4
	0<= x5
	end

view raw williams-production-mix-primal hosted with ❤ by GitHub

and its dual is

	\production mix - dual
	minimize
	obj: 288 y1 + 192 y2 + 384 y3
	subject to
	c0: 12 y1 + 10 y2 + 20 y3 >= 550
	c1: 20 y1 + 8 y2 + 20 y3 >= 600
	c2: 16 y2 + 20 y3 >= 350
	c3: 25 y1 + 20 y3 >= 400
	c4: 15 y1 + 20 y3 >= 200
	bounds
	0<= y1
	0<= y2
	0<= y3
	end

view raw williams-production-mix-dual hosted with ❤ by GitHub

Note that the two problems are formulated in a natural way. MOSEK  represents linear problems as

$$\begin{equation}\tag{MP} \begin{aligned} & \min c^T x + c_f \\ & s.t. & \\ && l^c \geq Ax \geq u^c\\ && l^x \geq x \geq u^x \end{aligned} \end{equation}$$

which is a more general and natural way for the user. After some boring but trivial calculations we obtain the dual problem in the form

$$\begin{equation}\tag{MD} \begin{aligned} & \max (l^c)^T s_l^c - (u^c)^T s_u^c + (l^x)^T s_l^x - (u^x)^T s_u^x\\ & s.t. & \\ && A^T(s_l^c - s_u^c)  + s_l^x - s_u^x=c\\ && s_l^x,s_u^x, s^c_l, s_u^c \geq 0 \end{aligned} \end{equation}$$

As you may see, there are exactly a variable for each primal constraints. Running the solver on the primal we get

	ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
	0 1.3e+01 2.7e+00 5.1e+00 1.00e+00 2.100000000e+03 0.000000000e+00 1.0e+00 0.00
	1 2.2e+00 1.7e+00 1.3e-01 0.00e+00 7.091163351e+03 6.790931921e+03 2.5e+00 0.00
	2 3.0e-01 2.3e-01 1.8e-02 7.18e-01 9.436588082e+03 9.372637344e+03 3.3e-01 0.00
	3 3.4e-02 2.6e-02 2.0e-03 8.14e-01 1.080109165e+04 1.079156482e+04 3.7e-02 0.00
	4 6.0e-04 4.7e-04 3.6e-05 9.82e-01 1.091666193e+04 1.091645458e+04 6.6e-04 0.00
	5 6.0e-08 4.7e-08 3.6e-09 1.00e+00 1.091999966e+04 1.091999964e+04 6.7e-08 0.00
	6 6.2e-12 4.7e-12 5.4e-13 1.00e+00 1.092000000e+04 1.092000000e+04 6.7e-12 0.00

view raw primal problem solver output hosted with ❤ by GitHub

and the solution looks like

	PROBLEM STATUS : PRIMAL_AND_DUAL_FEASIBLE
	SOLUTION STATUS : OPTIMAL
	OBJECTIVE NAME : obj
	PRIMAL OBJECTIVE : 1.09200000e+04
	DUAL OBJECTIVE : 1.09200000e+04
	CONSTRAINTS
	INDEX NAME AT ACTIVITY LOWER LIMIT UPPER LIMIT DUAL LOWER DUAL UPPER
	0 grinding UL 2.87999999998441e+02 NONE 2.88000000e+02 -0.00000000000000e+00 -6.24999999997661e+00
	1 drilling SB 1.77600000000379e+02 NONE 1.92000000e+02 -0.00000000000000e+00 1.76299093219222e-11
	2 manpower UL 3.83999999999842e+02 NONE 3.84000000e+02 -0.00000000000000e+00 -2.37499999999336e+01
	VARIABLES
	INDEX NAME AT ACTIVITY LOWER LIMIT UPPER LIMIT DUAL LOWER DUAL UPPER
	0 x1 SB 1.19999999999665e+01 0.00000000e+00 NONE -1.22142547778192e-10 -0.00000000000000e+00
	1 x2 SB 7.19999999992199e+00 0.00000000e+00 NONE -1.53700139339745e-10 -0.00000000000000e+00
	2 x3 LL 7.88886826796173e-11 0.00000000e+00 NONE -1.24999999999441e+02 -0.00000000000000e+00
	3 x4 LL 3.23546896699585e-12 0.00000000e+00 NONE -2.31249999999831e+02 -0.00000000000000e+00
	4 x5 LL 1.36145745561894e-11 0.00000000e+00 NONE -3.68749999999485e+02 -0.00000000000000e+00

view raw primal problem solution hosted with ❤ by GitHub

The columns named as dual represent the corresponding dual variables. Running the solver on the dual we obtain

	ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME
	0   2.7e+00  1.3e+01  5.1e+00  1.00e+00   0.000000000e+00   2.100000000e+03   1.0e+00  0.00
	1   1.7e+00  2.2e+00  1.3e-01  0.00e+00   6.790931921e+03   7.091163351e+03   2.5e+00  0.00
	2   2.3e-01  3.0e-01  1.8e-02  7.18e-01   9.372637344e+03   9.436588082e+03   3.3e-01  0.00
	3   2.6e-02  3.4e-02  2.0e-03  8.14e-01   1.079156482e+04   1.080109165e+04   3.7e-02  0.00
	4   4.7e-04  6.0e-04  3.6e-05  9.82e-01   1.091645458e+04   1.091666193e+04   6.6e-04  0.00
	5   4.7e-08  6.0e-08  3.6e-09  1.00e+00   1.091999964e+04   1.091999966e+04   6.7e-08  0.00
	6   4.7e-12  6.2e-12  5.4e-13  1.00e+00   1.092000000e+04   1.092000000e+04   6.7e-12  0.00

view raw dual problem solver output hosted with ❤ by GitHub

and the solution looks like

	PROBLEM STATUS : PRIMAL_AND_DUAL_FEASIBLE
	SOLUTION STATUS : OPTIMAL
	OBJECTIVE NAME : obj
	PRIMAL OBJECTIVE : 1.09200000e+04
	DUAL OBJECTIVE : 1.09200000e+04
	CONSTRAINTS
	INDEX NAME AT ACTIVITY LOWER LIMIT UPPER LIMIT DUAL LOWER DUAL UPPER
	0 c0 LL 5.50000000000122e+02 5.50000000e+02 NONE 1.19999999999665e+01 0.00000000000000e+00
	1 c1 LL 6.00000000000154e+02 6.00000000e+02 NONE 7.19999999992199e+00 0.00000000000000e+00
	2 c2 SB 4.74999999999441e+02 3.50000000e+02 NONE 7.88886826796173e-11 0.00000000000000e+00
	3 c3 SB 6.31249999999831e+02 4.00000000e+02 NONE 3.23546896699585e-12 0.00000000000000e+00
	4 c4 SB 5.68749999999485e+02 2.00000000e+02 NONE 1.36145745561894e-11 0.00000000000000e+00
	VARIABLES
	INDEX NAME AT ACTIVITY LOWER LIMIT UPPER LIMIT DUAL LOWER DUAL UPPER
	0 y1 SB 6.24999999997661e+00 0.00000000e+00 NONE 1.55872359232608e-09 0.00000000000000e+00
	1 y2 LL -1.76299093219222e-11 0.00000000e+00 NONE 1.43999999996206e+01 0.00000000000000e+00
	2 y3 SB 2.37499999999336e+01 0.00000000e+00 NONE 1.58100197085711e-10 0.00000000000000e+00

view raw dual problem solution hosted with ❤ by GitHub

But wait a minute...it seems MOSEK just swap primal and dual!?!?!?

True in some sense! MOSEK implements a primal-dual algorithm that solves both primal and dual at the same time. For small problems there might be no difference at all. Few comments on the solver output:

1. Feasibility and optimality are achieved at the same time: during the execution the intermediate solutions are neither primal nor dual feasible. This is shown in the columns PFEAS and DFEAS as the norm of the constraint violation.
2. Column POBJ and DOBJ show the progress towards the optimal value and, as Strong Duality Theorem told us, they are the same once optimality is reached.
   
3. Dual variables are reported in the LOWER/UPPER DUAL columns.
4. To each non active bound corresponds a zero dual variable, as for the Complementary Slackness.
5. Dual values are zero for each missing bounds: this is because a solver, as MOSEK does, actually sets a very large bound (in absolute value) in order to overcome numerical issues. A missing bound resolves in practice to a non active bound, as in the previous point.

If you want to dig more in the meaning of the dual information, follow Section 6.2 in [2]. For us,  it is important you get the message:

Dual information are readily available from your solver and can be useful!

Most of what we have said is also, and even more, true for second order conic programs: duality holds but to play with it you should ask your preferred solver.

[1] Papadimitriou, C. H., & Steiglitz, K. (1998). Combinatorial optimization: algorithms and complexity. Courier Dover Publications.
[2] Williams, H. P. (1999). Model building in mathematical programming.