13.2 Linear Interior-point optimizer

The purpose of this section is to provide information about the algorithm employed in the MOSEK interior-point optimizer for linear problems and about its termination criteria.

13.2.1 The homogeneous primal-dual problem

In order to keep the discussion simple it is assumed that MOSEK solves linear optimization problems of standard form

(13.1)\[\begin{split}\begin{array}{lcll} \mbox{minimize} & c^T x & & \\ \mbox{subject to} & A x & = & b, \\ & x\geq 0. & & \end{array}\end{split}\]

This is in fact what happens inside MOSEK; for efficiency reasons MOSEK converts the problem to standard form before solving, then converts it back to the input form when reporting the solution.

Since it is not known beforehand whether problem (13.1) has an optimal solution, is primal infeasible or is dual infeasible, the optimization algorithm must deal with all three situations. This is the reason why MOSEK solves the so-called homogeneous model

(13.2)\[\begin{split}\begin{array}{rcl} A x - b \tau & = & 0, \\ A^T y + s - c \tau & = & 0, \\ -c^T x + b^T y - \kappa & = & 0, \\ x,s,\tau ,\kappa & \geq & 0, \end{array}\end{split}\]

where \(y\) and \(s\) correspond to the dual variables in (13.1), and \(\tau\) and \(\kappa\) are two additional scalar variables. Note that the homogeneous model (13.2) always has solution since

\[(x,y,s,\tau ,\kappa ) = (0,0,0,0,0)\]

is a solution, although not a very interesting one. Any solution

\[(x^*,y^*,s^*,\tau^*,\kappa^*)\]

to the homogeneous model (13.2) satisfies

\[x_{j}^* s_{j}^* = 0 \mbox{ and } \tau^* \kappa^* = 0.\]

Moreover, there is always a solution that has the property \(\tau^* + \kappa^* > 0\).

First, assume that \(\tau^*>0\) . It follows that

\[\begin{split}\begin{array}{rcl} A \frac{x^*}{\tau^*} & = & b, \\ A^T \frac{y^*}{\tau^*} + \frac{s^*}{\tau^*} & = & c, \\ - c^T \frac{x^*}{\tau^*} + b^T\frac{y^*}{\tau^*} & = & 0, \\ x^*,s^*,\tau^*,\kappa^* & \geq & 0. \end{array}\end{split}\]

This shows that \(\frac{x^*}{\tau^*}\) is a primal optimal solution and \((\frac{y^*}{\tau^*},\frac{s^*}{\tau *})\) is a dual optimal solution; this is reported as the optimal interior-point solution since

\[(x,y,s) = \left\lbrace \frac{x^*}{\tau^*},\frac{y^*}{\tau^*},\frac{s^*}{\tau^*} \right\rbrace\]

is a primal-dual optimal solution (see Sec. 12.1 (Linear Optimization) for the mathematical background on duality and optimality).

On other hand, if \(\kappa^* >0\) then

\[\begin{split}\begin{array}{rcl} A x^* & = & 0, \\ A^T y^* + s^* & = & 0, \\ -c^T x^* + b^T y^* & = & \kappa^*, \\ x^*, s^*, \tau^* ,\kappa^* & \geq & 0. \end{array}\end{split}\]

This implies that at least one of

(13.3)\[c^T x^* < 0\]

or

(13.4)\[b^T y^* > 0\]

is satisfied. If (13.3) is satisfied then \(x^*\) is a certificate of dual infeasibility, whereas if (13.4) is satisfied then \(y^*\) is a certificate of primal infeasibility.

In summary, by computing an appropriate solution to the homogeneous model, all information required for a solution to the original problem is obtained. A solution to the homogeneous model can be computed using a primal-dual interior-point algorithm [And09].

13.2.2 Interior-point Termination Criterion

For efficiency reasons it is not practical to solve the homogeneous model exactly. Hence, an exact optimal solution or an exact infeasibility certificate cannot be computed and a reasonable termination criterion has to be employed.

In the \(k\)-th iteration of the interior-point algorithm a trial solution

\[(x^k,y^k,s^k,\tau^k,\kappa^k)\]

to homogeneous model is generated, where

\[x^k,s^k,\tau^k,\kappa^k > 0.\]

Optimal case

Whenever the trial solution satisfies the criterion

(13.5)\[\begin{split}\begin{array}{rcl} \left\| A \frac{x^k}{\tau^k} - b \right\|_\infty & \leq & \epsilon_p (1+\left\| b \right\|_\infty ), \\ \left\| A^T \frac{y^k}{\tau^k} + \frac{s^k}{\tau^k}- c \right\|_\infty & \leq & \epsilon_d (1+\left\| c \right\|_\infty ), \mbox{ and} \\ \min \left( \frac{(x^k)^T s^k}{ (\tau^k)^2 }, \vert \frac{c^T x^k}{\tau^k} - \frac{ b^T y^k}{\tau^k} \vert \right) & \leq & \epsilon_g \max \left( 1, \frac{ \min\left( \left\vert c^T x^k \right\vert, \left\vert b^T y^k \right\vert \right) }{\tau^k} \right), \end{array}\end{split}\]

the interior-point optimizer is terminated and

\[\frac{(x^k,y^k,s^k)}{\tau^k}\]

is reported as the primal-dual optimal solution. The interpretation of (13.5) is that the optimizer is terminated if

  • \(\frac{x^k}{\tau^k}\) is approximately primal feasible,

  • \(\left\lbrace \frac{y^k}{\tau^k},\frac{s^k}{\tau^k} \right\rbrace\) is approximately dual feasible, and

  • the duality gap is almost zero.

Dual infeasibility certificate

On the other hand, if the trial solution satisfies

\[-\epsilon_i c^T x^k > \frac{ \left\| c \right\|_\infty }{ \max\left( 1, \left\| b \right\|_\infty \right) } \left\| A x^k \right\|_\infty\]

then the problem is declared dual infeasible and \(x^k\) is reported as a certificate of dual infeasibility. The motivation for this stopping criterion is as follows: First assume that \(\left\| A x^k \right\|_\infty = 0\) ; then \(x^k\) is an exact certificate of dual infeasibility. Next assume that this is not the case, i.e.

\[\left\| A x^{k} \right\|_\infty > 0,\]

and define

\[\bar{x} := \epsilon_i \frac{ \max \left( 1, \left\| b \right\|_\infty \right) } { \left\| A x^k \right\|_\infty \left\| c \right\|_\infty } x^k.\]

It is easy to verify that

\[\left\| A \bar x \right\|_\infty = \epsilon_i \frac{\max \left( 1,\left\| b \right\|_\infty\right)} {\left\| c \right\|_\infty} \mbox{ and } -c^T \bar{x} > 1,\]

which shows \(\bar{x}\) is an approximate certificate of dual infeasibility, where \(\varepsilon_{i}\) controls the quality of the approximation. A smaller value means a better approximation.

Primal infeasibility certificate

Finally, if

\[\epsilon_i b^T y^k > \frac{ \left\| b \right\|_\infty } { \max \left( 1,\left\| c \right\|_\infty \right) } \left\| A^T y^k + s^k \right\|_\infty\]

then \(y^k\) is reported as a certificate of primal infeasibility.

13.2.3 Adjusting optimality criteria

It is possible to adjust the tolerances \(\varepsilon_{p}\), \(\varepsilon_{d}\), \(\varepsilon_{g}\) and \(\varepsilon_{i}\) using parameters; see table for details.

Table 13.1 Parameters employed in termination criterion

ToleranceParameter

name

\(\varepsilon_{p}\)

MSK_DPAR_INTPNT_TOL_PFEAS

\(\varepsilon_{d}\)

MSK_DPAR_INTPNT_TOL_DFEAS

\(\varepsilon_{g}\)

MSK_DPAR_INTPNT_TOL_REL_GAP

\(\varepsilon_{i}\)

MSK_DPAR_INTPNT_TOL_INFEAS

The default values of the termination tolerances are chosen such that for a majority of problems appearing in practice it is not possible to achieve much better accuracy. Therefore, tightening the tolerances usually is not worthwhile. However, an inspection of (13.5) reveals that the quality of the solution depends on \(\left\| b \right\|_{\infty}\) and \(\left\| c \right\|_{\infty}\); the smaller the norms are, the better the solution accuracy.

The interior-point method as implemented by MOSEK will converge toward optimality and primal and dual feasibility at the same rate [And09]. This means that if the optimizer is stopped prematurely then it is very unlikely that either the primal or dual solution is feasible. Another consequence is that in most cases all the tolerances, \(\varepsilon_p\), \(\varepsilon_d\), \(\varepsilon_g\) and \(\varepsilon_i\), have to be relaxed together to achieve an effect.

The basis identification discussed in Sec. 13.6 (Basis Identification) requires an optimal solution to work well; hence basis identification should be turned off if the termination criterion is relaxed.

To conclude the discussion in this section, relaxing the termination criterion is usually not worthwhile.

13.2.4 The Interior-point Log

Below is a typical log output from the interior-point optimizer:

Optimizer  - threads                : 32
Optimizer  - solved problem         : the primal
Optimizer  - Constraints            : 699
Optimizer  - Cones                  : 0
Optimizer  - Scalar variables       : 1749              conic                  : 0
Optimizer  - Semi-definite variables: 0                 scalarized             : 0
Factor     - GP order time          : 0.01
Factor     - nonzeros before factor : 1.14e+04          after factor           : 3.00e+04
Factor     - flops                  : 2.13e+06
Interior-point optimizer setup terminated. Time: 0.02
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME
0   1.4e+03  9.0e+01  9.0e+03  0.00e+00   4.127541176e+03   -4.825341578e+03  3.4e+00  0.02
1   7.1e+02  4.5e+01  4.5e+03  -6.92e-01  4.408986445e+03   -2.847314813e+03  1.7e+00  0.02
2   3.2e+02  2.0e+01  2.0e+03  -2.95e-01  5.251854912e+03   4.300306938e+02   7.6e-01  0.02
...
13  3.0e-07  1.9e-08  1.9e-06  1.00e+00   5.501845893e+03   5.501845887e+03   7.1e-10  0.04
14  2.9e-11  1.9e-12  2.2e-10  1.00e+00   5.501845888e+03   5.501845888e+03   7.2e-14  0.05
Interior-point optimizer terminated. Time: 0.05

The first lines summarize the problem the optimizer is solving and various settings. This is followed by the iteration log, with the following meaning:

  • ITE: Iteration index \(k\).

  • PFEAS: \(\left\| Ax^k - b \tau^k \right\|_{\infty}\) . The numbers in this column should converge monotonically towards zero but may stall at low level due to rounding errors.

  • DFEAS: \(\left\| A^T y^{k} + s^k - c \tau^k \right\|_{\infty}\) . The numbers in this column should converge monotonically towards zero but may stall at low level due to rounding errors.

  • GFEAS: \(| - c^T x^k + b^T y^{k} - \kappa^k|\) . The numbers in this column should converge monotonically towards zero but may stall at low level due to rounding errors.

  • PRSTATUS: This number converges to \(1\) if the problem has an optimal solution whereas it converges to \(-1\) if that is not the case.

  • POBJ: \(c^T x^k / \tau^k\). An estimate for the primal objective value.

  • DOBJ: \(b^T y^k / \tau^k\). An estimate for the dual objective value.

  • MU: \(\frac{ (x^k)^T s^k + \tau^k \kappa^k }{n+1}\) . The numbers in this column should always converge to zero.

  • TIME: Time spent since the optimization started.