8. The optimizers for continuous problems


The most essential part of MOSEK is the optimizers. Each optimizer is designed to solve a particular class of problems i.e. linear, conic, or general nonlinear problems. The purpose of the present chapter is to discuss which optimizers are available for the continuous problem classes and how the performance of an optimizer can be tuned, if needed.

This chapter deals with the optimizers for continuous problems with no integer variables.

8.1. How an optimizer works

When the optimizer is called, it roughly performs the following steps:

Presolve:

Preprocessing to reduce the size of the problem.

Dualizer:

Choosing whether to solve the primal or the dual form of the problem.

Scaling:

Scaling the problem for better numerical stability.

Optimize:

Solve the problem using selected method.

The first three preprocessing steps are transparent to the user, but useful to know about for tuning purposes. In general, the purpose of the preprocessing steps is to make the actual optimization more efficient and robust.

8.1.1. Presolve

Before an optimizer actually performs the optimization the problem is preprocessed using the so-called presolve. The purpose of the presolve is to

  • remove redundant constraints,
  • eliminate fixed variables,
  • remove linear dependencies,
  • substitute out free variables, and
  • reduce the size of the optimization problem in general.

After the presolved problem has been optimized the solution is automatically postsolved so that the returned solution is valid for the original problem. Hence, the presolve is completely transparent. For further details about the presolve phase, please see [10, 5].

It is possible to fine-tune the behavior of the presolve or to turn it off entirely. If presolve consumes too much time or memory compared to the reduction in problem size gained it may be disabled. This is done by setting the parameter MSK_IPAR_PRESOLVE_USE to MSK_PRESOLVE_MODE_OFF.

The two most time-consuming steps of the presolve are

  • the eliminator, and
  • the linear dependency check.

Therefore, in some cases it is worthwhile to disable one or both of these.

8.1.1.1. Eliminator

The purpose of the eliminator is to eliminate free and implied free variables from the problem using substitution. For instance, given the constraints

\begin{displaymath}\nonumber{}\begin{array}{rcl}\nonumber{}y & = & \sum _{j}x_{j},\\\nonumber{}y,x & \geq{} & 0,\end{array}\end{displaymath}

y is an implied free variable that can be substituted out of the problem, if deemed worthwhile.

If the eliminator consumes too much time or memory compared to the reduction in problem size gained it may be disabled. This can be done with the parameter MSK_IPAR_PRESOLVE_ELIMINATOR_USE to MSK_OFF.

8.1.1.2. Linear dependency checker

The purpose of the linear dependency check is to remove linear dependencies among the linear equalities. For instance, the three linear equalities

\begin{displaymath}\nonumber{}\begin{array}{rcl}\nonumber{}x_{1}+x_{2}+x_{3} & = & 1,\\\nonumber{}x_{1}+0.5x_{2} & = & 0.5,\\\nonumber{}0.5x_{2}+x_{3} & = & 0.5\end{array}\end{displaymath}

contain exactly one linear dependency. This implies that one of the constraints can be dropped without changing the set of feasible solutions. Removing linear dependencies is in general a good idea since it reduces the size of the problem. Moreover, the linear dependencies are likely to introduce numerical problems in the optimization phase.

It is best practise to build models without linear dependencies. If the linear dependencies are removed at the modeling stage, the linear dependency check can safely be disabled by setting the parameter MSK_IPAR_PRESOLVE_LINDEP_USE to MSK_OFF.

8.1.2. Dualizer

All linear, conic, and convex optimization problems have an equivalent dual problem associated with them. MOSEK has built-in heuristics to determine if it is most efficient to solve the primal or dual problem. The form (primal or dual) solved is displayed in the MOSEK log. Should the internal heuristics not choose the most efficient form of the problem it may be worthwhile to set the dualizer manually by setting the parameters:

Note that currently only linear problems may be dualized.

8.1.3. Scaling

Problems containing data with large and/or small coefficients, say 1.0e+9 or 1.0e-7, are often hard to solve. Significant digits may be truncated in calculations with finite precision, which can result in the optimizer relying on inaccurate calculations. Since computers work in finite precision, extreme coefficients should be avoided. In general, data around the same “order of magnitude” is preferred, and we will refer to a problem, satisfying this loose property, as being well-scaled. If the problem is not well scaled, MOSEK will try to scale (multiply) constraints and variables by suitable constants. MOSEK solves the scaled problem to improve the numerical properties.

The scaling process is transparent, i.e. the solution to the original problem is reported. It is important to be aware that the optimizer terminates when the termination criterion is met on the scaled problem, therefore significant primal or dual infeasibilities may occur after unscaling for badly scaled problems. The best solution to this problem is to reformulate it, making it better scaled.

By default MOSEK heuristically chooses a suitable scaling. The scaling for interior-point and simplex optimizers can be controlled with the parameters MSK_IPAR_INTPNT_SCALING and MSK_IPAR_SIM_SCALING respectively.

8.1.4. Using multiple CPU's

The interior-point optimizers in MOSEK have been parallelized. This means that if you solve linear, quadratic, conic, or general convex optimization problem using the interior-point optimizer, you can take advantage of multiple CPU's.

By default MOSEK uses one thread to solve the problem, but the number of threads (and thereby CPUs) employed can be changed by setting the parameter MSK_IPAR_INTPNT_NUM_THREADS This should never exceed the number of CPU's on the machine.

The speed-up obtained when using multiple CPUs is highly problem and hardware dependent, and consequently, it is advisable to compare single threaded and multi threaded performance for the given problem type to determine the optimal settings.

For small problems, using multiple threads will probably not be worthwhile.

8.2. Linear optimization

8.2.1. Optimizer selection

Two different types of optimizers are available for linear problems: The default is an interior-point method, and the alternatives are simplex methods. The optimizer can be selected using the parameter MSK_IPAR_OPTIMIZER.

8.2.2. The interior-point optimizer

The purpose of this section is to provide information about the algorithm employed in MOSEK interior-point optimizer.

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

\begin{math}\nonumber{}\begin{array}{lcll}\nonumber{}\mbox{minimize} & c^{T}x &  & \\\nonumber{}\mbox{subject to} & Ax & = & b,\\\nonumber{} & x\geq{}0. &  &\end{array}\end{math} (8.2.1)

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

Since it is not known beforehand whether problem (8.2.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 that MOSEK solves the so-called homogeneous model

\begin{math}\nonumber{}\begin{array}{rcl}\nonumber{}Ax-b\tau  & = & 0,\\\nonumber{}A^{T}y+s-c\tau  & = & 0,\\\nonumber{}-c^{T}x+b^{T}y-\kappa  & = & 0,\\\nonumber{}x,s,\tau ,\kappa  & \geq{} & 0,\end{array}\end{math} (8.2.2)

where y and s correspond to the dual variables in (8.2.1), and [[MathCmd 498]] and [[MathCmd 499]] are two additional scalar variables. Note that the homogeneous model (8.2.2) always has solution since

\begin{displaymath}\nonumber{}(x,y,s,\tau ,\kappa )=(0,0,0,0,0)\end{displaymath}

is a solution, although not a very interesting one.

Any solution

\begin{displaymath}\nonumber{}(x^{*},y^{*},s^{*},\tau ^{*},\kappa ^{*})\end{displaymath}

to the homogeneous model (8.2.2) satisfies

\begin{displaymath}\nonumber{}x_{j}^{*}s_{j}^{*}=0\mbox{ and }\tau ^{*}\kappa ^{*}=0.\end{displaymath}

Moreover, there is always a solution that has the property

\begin{displaymath}\nonumber{}\tau ^{*}+\kappa ^{*}>0.\end{displaymath}

First, assume that [[MathCmd 504]]. It follows that

\begin{math}\nonumber{}\begin{array}{rcl}\nonumber{}A\frac{x^{*}}{\tau ^{*}} & = & b,\\\nonumber{}A^{T}\frac{y^{*}}{\tau ^{*}}+\frac{s^{*}}{\tau *} & = & c,\\\nonumber{}-c^{T}\frac{x^{*}}{\tau ^{*}}+b^{T}\frac{y^{*}}{\tau ^{*}} & = & 0,\\\nonumber{}x^{*},s^{*},\tau ^{*},\kappa ^{*} & \geq{} & 0.\end{array}\end{math} (8.2.3)

This shows that [[MathCmd 506]] is a primal optimal solution and [[MathCmd 507]] is a dual optimal solution; this is reported as the optimal interior-point solution since

\begin{displaymath}\nonumber{}(x,y,s)=\left(\frac{x^{*}}{\tau ^{*}},\frac{y^{*}}{\tau ^{*}},\frac{s^{*}}{\tau *}\right)\end{displaymath}

is a primal-dual optimal solution.

On other hand, if [[MathCmd 509]] then

\begin{math}\nonumber{}\begin{array}{rclcl}\nonumber{}Ax^{*} & = & 0, &  & \\\nonumber{}A^{T}y^{*}+s^{*} & = & 0, &  & \\\nonumber{}-c^{T}x^{*}+b^{T}y^{*} & = & \kappa ^{*}, &  & \\\nonumber{}x^{*},s^{*},\tau ^{*},\kappa ^{*} & \geq{} & 0. &  &\end{array}\end{math} (8.2.4)

This implies that at least one of

\begin{math}\nonumber{}-c^{T}x^{*}>0\end{math} (8.2.5)

or

\begin{math}\nonumber{}b^{T}y^{*}>0\end{math} (8.2.6)

is satisfied. If (8.2.5) is satisfied then [[MathCmd 513]] is a certificate of dual infeasibility, whereas if (8.2.6) is satisfied then [[MathCmd 514]] is a certificate of dual 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 [8].

8.2.2.1. 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 every iteration, k, of the interior-point algorithm a trial solution

\begin{displaymath}\nonumber{}(x^{k},y^{k},s^{k},\tau ^{k},\kappa ^{k})\end{displaymath}

to homogeneous model is generated where

\begin{displaymath}\nonumber{}x^{k},s^{k},\tau ^{k},\kappa ^{k}>0.\end{displaymath}

Whenever the trial solution satisfies the criterion

\begin{math}\nonumber{}\begin{array}{rcl}\nonumber{}\left\|A\frac{x^{k}}{\tau ^{k}}-b\right\| & \leq{} & \varepsilon _{p}(1+\left\|b\right\|),\\\nonumber{}\left\|A^{T}\frac{y^{k}}{\tau ^{k}}+\frac{s^{k}}{\tau ^{k}}-c\right\| & \leq{} & \varepsilon _{d}(1+\left\|c\right\|),\mbox{ and}\\\nonumber{}\min \left(\frac{(x^{k})^{T} s^{k} +\tau ^{k}\kappa ^{k}}{(\tau ^{k})^{2}},\left|\frac{c^{T} x^{k}}{\tau ^{k}}-\frac{b^{T} y^{k}}{\tau ^{k}}\right|\right) & \leq{} & \varepsilon _{g}\max \left(1,\left|\frac{c^{T} x^{k}}{\tau ^{k}}\right|\right),\end{array}\end{math} (8.2.7)

the interior-point optimizer is terminated and

\begin{displaymath}\nonumber{}\frac{(x^{k},y^{k},s^{k})}{\tau ^{k}}\end{displaymath}

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

  • [[MathCmd 519]] is approximately primal feasible,
  • [[MathCmd 520]] is approximately dual feasible, and
  • the duality gap is almost zero.

On the other hand, if the trial solution satisfies

\begin{math}\nonumber{}-\varepsilon _{i}c^{T}x^{k}>\frac{\left\|c\right\|}{\max (\left\|b\right\|,1)}\left\|Ax^{k}\right\|\end{math} (8.2.8)

then the problem is declared dual infeasible and [[MathCmd 522]] is reported as a certificate of dual infeasibility. The motivation for this stopping criterion is as follows: First assume that [[MathCmd 523]]; then [[MathCmd 522]] is an exact certificate of dual infeasibility. Next assume that this is not the case, i.e.

\begin{displaymath}\nonumber{}\left\|Ax^{k}\right\|>0,\end{displaymath}

and define

\begin{displaymath}\nonumber{}\bar{x}:=\varepsilon _{i}\frac{\max (1,\left\|b\right\|) x^{k}}{\left\|A x^{k}\right\|\left\|c\right\|}.\end{displaymath}

It is easy to verify that

\begin{displaymath}\nonumber{}\left\|A\bar{x}\right\|=\varepsilon _{i}\mbox{ and }-c^{T}\bar{x}>1,\end{displaymath}

which shows [[MathCmd 528]] is an approximate certificate dual infeasibility where [[MathCmd 529]] controls the quality of the approximation. A smaller value means a better approximation.

Finally, if

\begin{math}\nonumber{}\varepsilon _{i}b^{T}y^{k}\geq{}\frac{\left\|b\right\|}{\max (1,\left\|c\right\|)}\left\|A^{T}y^{k}+s^{k}\right\|\end{math} (8.2.9)

then [[MathCmd 531]] is reported as a certificate of primal infeasibility.

It is possible to adjust the tolerances [[MathCmd 532]], [[MathCmd 533]], [[MathCmd 534]] and [[MathCmd 529]] using parameters; see table 8.1 for details.

Table 8.1: Parameters employed in termination criterion.

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 (8.2.7) reveals that quality of the solution is dependent on [[MathCmd 540]] and [[MathCmd 541]]; 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 [8]. 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, [[MathCmd 532]], [[MathCmd 533]] and [[MathCmd 534]], has to be relaxed together to achieve an effect.

The basis identification discussed in section 8.2.2.2 requires an optimal solution to work well; hence basis identification should turned off if the termination criterion is relaxed.

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

8.2.2.2. Basis identification

An interior-point optimizer does not return an optimal basic solution unless the problem has a unique primal and dual optimal solution. Therefore, the interior-point optimizer has an optional post-processing step that computes an optimal basic solution starting from the optimal interior-point solution. More information about the basis identification procedure may be found in [6].

Please note that a basic solution is often more accurate than an interior-point solution.

By default MOSEK performs a basis identification. However, if a basic solution is not needed, the basis identification procedure can be turned off. The parameters

controls when basis identification is performed.

8.2.2.3. The interior-point log

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

Optimizer  - threads                : 1
Optimizer  - solved problem         : the dual
Optimizer  - constraints            : 2                 variables              : 6
Factor     - setup time             : 0.04              order time             : 0.00
Factor     - GP order used          : no                GP order time          : 0.00
Factor     - nonzeros before factor : 3                 after factor           : 3
Factor     - offending columns      : 0                 flops                  : 1.70e+001
ITE PFEAS    DFEAS    KAP/TAU  POBJ              DOBJ              MU       TIME
0   2.0e+002 2.9e+001 2.0e+002 -0.000000000e+000 -1.204741644e+003 2.0e+002 0.44
1   2.2e+001 3.1e+000 7.3e+002 -5.885951891e+003 -5.856764353e+003 2.2e+001 0.57
2   3.8e+000 5.4e-001 9.7e+001 -7.405187479e+003 -7.413054916e+003 3.8e+000 0.58
3   4.0e-002 5.7e-003 2.6e-001 -7.664507945e+003 -7.665313396e+003 4.0e-002 0.58
4   4.2e-006 6.0e-007 2.7e-005 -7.667999629e+003 -7.667999714e+003 4.2e-006 0.59
5   4.2e-010 6.0e-011 2.7e-009 -7.667999994e+003 -7.667999994e+003 4.2e-010 0.59
The first line displays the number of threads used by the optimizer and second line tells that the optimizer choose to solve the dual problem rather the primal problem. The next line displays the problem dimensions as seen by the optimizer, and the “Factor...” lines show various statistics. This is followed by the iteration log.

Using the same notation as in section 8.2.2 the columns of the iteration log has the following meaning:

  • ITE: Iteration index.
  • PFEAS: [[MathCmd 545]]. The numbers in this column should converge monotonically towards to zero.
  • DFEAS: [[MathCmd 546]]. The numbers in this column should converge monotonically toward to zero.
  • KAP/TAU: [[MathCmd 547]]. If the numbers in this column converge toward zero then the problem has an optimal solution. Otherwise if the numbers converge towards infinity, the problem is primal or/and dual infeasible.
  • POBJ: [[MathCmd 548]]. An estimate for the primal objective value.
  • DOBJ: [[MathCmd 549]]. An estimate for the dual objective value.
  • MU: [[MathCmd 550]]. The numbers in this column should always converge monotonically to zero.
  • TIME: Time spend since the optimization started.

8.2.3. The simplex based optimizer

An alternative to the interior-point optimizer is the simplex optimizer.

The simplex optimizer uses a different method that allows exploiting an initial guess for the optimal solution to reduce the solution time. Depending on the problem it may be faster or slower to use an initial guess; see section 8.2.4 for a discussion.

MOSEK provides both a primal and a dual variant of the simplex optimizer — we will return to this later.

8.2.3.1. Simplex termination criterion

The simplex optimizer terminates when it finds an optimal basic solution or an infeasibility certificate. A basic solution is optimal when it is primal and dual feasible; see (7.1.1) and (7.1.2) for a definition of the primal and dual problem. Due the fact that to computations are performed in finite precision MOSEK allows violation of primal and dual feasibility within certain tolerances. The user can control the allowed primal and dual infeasibility with the parameters MSK_DPAR_BASIS_TOL_X and MSK_DPAR_BASIS_TOL_S.

8.2.3.2. Starting from an existing solution

When using the simplex optimizer it may be possible to reuse an existing solution and thereby reduce the solution time significantly. When a simplex optimizer starts from an existing solution it is said to perform a hot-start. If the user is solving a sequence of optimization problems by solving the problem, making modifications, and solving again, MOSEK will hot-start automatically.

Setting the parameter MSK_IPAR_OPTIMIZER to MSK_OPTIMIZER_FREE_SIMPLEX instructs MOSEK to select automatically between the primal and the dual simplex optimizers. Hence, MOSEK tries to choose the best optimizer for the given problem and the available solution.

By default MOSEK uses presolve when performing a hot-start. If the optimizer only needs very few iterations to find the optimal solution it may be better to turn off the presolve.

8.2.3.3. Numerical difficulties in the simplex optimizers

Though MOSEK is designed to minimize numerical instability, completely avoiding it is impossible when working in finite precision. MOSEK counts a “numerical unexpected behavior” event inside the optimizer as a set-back. The user can define how many set-backs the optimizer accepts; if that number is exceeded, the optimization will be aborted. Set-backs are implemented to avoid long sequences where the optimizer tries to recover from an unstable situation.

Set-backs are, for example, repeated singularities when factorizing the basis matrix, repeated loss of feasibility, degeneracy problems (no progress in objective) and other events indicating numerical difficulties. If the simplex optimizer encounters a lot of set-backs the problem is usually badly scaled; in such a situation try to reformulate into a better scaled problem. Then, if a lot of set-backs still occur, trying one or more of the following suggestions may be worthwhile:

8.2.4. The interior-point or the simplex optimizer?

Given a linear optimization problem, which optimizer is the best: The primal simplex, the dual simplex or the interior-point optimizer?

It is impossible to provide a general answer to this question, however, the interior-point optimizer behaves more predictably — it tends to use between 20 and 100 iterations, almost independently of problem size — but cannot perform hot-start, while simplex can take advantage of an initial solution, but is less predictable for cold-start. The interior-point optimizer is used by default.

8.2.5. The primal or the dual simplex variant?

MOSEK provides both a primal and a dual simplex optimizer. Predicting which simplex optimizer is faster is impossible, however, in recent years the dual optimizer has seen several algorithmic and computational improvements, which, in our experience, makes it faster on average than the primal simplex optimizer. Still, it depends much on the problem structure and size.

Setting the MSK_IPAR_OPTIMIZER parameter to MSK_OPTIMIZER_FREE_SIMPLEX instructs MOSEK to choose which simplex optimizer to use automatically.

To summarize, if you want to know which optimizer is faster for a given problem type, you should try all the optimizers.

Alternatively, use the concurrent optimizer presented in Section 8.6.3.

8.3. Linear network optimization

8.3.1. Network flow problems

Linear optimization problems with the network flow structure as specified in section 6.5 can often be solved significantly faster with a specialized version of the simplex method [7] than with the general solvers.

MOSEK includes a network simplex solver which, on avarage, solves network problems 10 to 100 times faster than the standard simplex optimizers.

To use the network simplex optimizer, do the following:

MOSEK will automatically detect the network structure and apply the specialized simplex optimizer.

8.3.2. Embedded network problems

Often problems contains both large parts with network structure and some non-network constraints or variables — such problems are said to have embedded network structure.

If the procedure described in section 8.3.1 is applied, MOSEK will attemt to exploit this structure to speed up the optimization.

This is done heuristically by detecting the largest network embedded in the problem, solving this subproblem using the network simplex optimizer, and using the solution to hot-start a normal simplex optimizer.

The MSK_IPAR_SIM_NETWORK_DETECT parameter defines how large a percentage of the problem should be a network before the specialized solver is applied. In general, it is recommended to use the network optimizer only on problems containing a substantial embedded network.

If MOSEK only finds limited network structure in a problem, consider trying to switch off presolve MSK_IPAR_PRESOLVE_USE and scaling MSK_IPAR_SIM_SCALING, since in rare cases it might disturb the network heuristic.

The network detection heuristic can also be called directly through MSK_netextraction.

8.4. Conic optimization

8.4.1. The interior-point optimizer

For conic optimization problems only an interior-point type optimizer is available. The interior-point optimizer is an implementation of the so-called homogeneous and self-dual algorithm. For a detailed description of the algorithm, please see [26].

8.4.1.1. Interior-point termination criteria

The parameters controlling when the conic interior-point optimizer terminates are shown in Table 8.2.

Parameter name Purpose
MSK_DPAR_INTPNT_CO_TOL_PFEAS Controls primal feasibility
MSK_DPAR_INTPNT_CO_TOL_DFEAS Controls dual feasibility
MSK_DPAR_INTPNT_CO_TOL_REL_GAP Controls relative gap
MSK_DPAR_INTPNT_TOL_INFEAS Controls when the problem is declared infeasible
MSK_DPAR_INTPNT_CO_TOL_MU_RED Controls when the complementarity is reduced enough
Table 8.2: Parameters employed in termination criterion.

8.5. Nonlinear convex optimization

8.5.1. The interior-point optimizer

For quadratic, quadratically constrained, and general convex optimization problems an interior-point type optimizer is available. The interior-point optimizer is an implementation of the homogeneous and self-dual algorithm. For a detailed description of the algorithm, please see [23, 22].

8.5.1.1. The convexity requirement

Continuous nonlinear problems are required to be convex. For quadratic problems MOSEK test this requirement before optimizing. Specifying a non-convex problem results in an error message.

The following parameters are available to control the convexity check:

8.5.1.2. The differentiabilty requirement

The nonlinear optimizer in MOSEK requires both first order and second order derivatives. This of course implies care should be taken when solving problems involving non-differentiable functions.

For instance, the function

\begin{displaymath}\nonumber{}f(x)=x^{2}\end{displaymath}

is differentiable everywhere whereas the function

\begin{displaymath}\nonumber{}f(x)=\sqrt{x}\end{displaymath}

is only diffrentiable for x>0. In order to make sure that MOSEK evaulates the functions at points where they are differentiable, the function domains must be defined by setting appropriate variable bounds.

In general, if a variable is not ranged MOSEK will only evaluate that variable at points strictly within the bounds. Hence, imposing the bound

\begin{displaymath}\nonumber{}x\geq{}0\end{displaymath}

in the case of [[MathCmd 554]] is sufficient to guarantee that the function will only be evaluated in points where it is differentiable.

However, if a function is differentiable on closed a range, specifying the variable bounds is not sufficient. Consider the function

\begin{math}\nonumber{}f(x)=\frac{1}{x}+\frac{1}{1-x}.\end{math} (8.5.1)

In this case the bounds

\begin{displaymath}\nonumber{}0\leq{}x\leq{}1\end{displaymath}

will not guarantee that MOSEK only evalues the function for x between 0 and 1. To force MOSEK to strictly satisfy both bounds on ranged variables set the parameter MSK_IPAR_INTPNT_STARTING_POINT to MSK_STARTING_POINT_SATISFY_BOUNDS.

For efficiency reasons it may be better to reformulate the problem than to force MOSEK to observe ranged bounds strictly. For instance, (8.5.1) can be reformulated as follows

\begin{displaymath}\nonumber{}\begin{array}{rcl}\nonumber{}f(x) & = & \frac{1}{x}+\frac{1}{y}\\\nonumber{}0 & = & 1-x-y\\\nonumber{}0 & \leq{} & x\\\nonumber{}0 & \leq{} & y.\end{array}\end{displaymath}

8.5.1.3. Interior-point termination criteria

The parameters controlling when the general convex interior-point optimizer terminates are shown in Table 8.3.

Parameter name Purpose
MSK_DPAR_INTPNT_NL_TOL_PFEAS Controls primal feasibility
MSK_DPAR_INTPNT_NL_TOL_DFEAS Controls dual feasibility
MSK_DPAR_INTPNT_NL_TOL_REL_GAP Controls relative gap
MSK_DPAR_INTPNT_TOL_INFEAS Controls when the problem is declared infeasible
MSK_DPAR_INTPNT_NL_TOL_MU_RED Controls when the complementarity is reduced enough
Table 8.3: Parameters employed in termination criteria.

8.6. Solving problems in parallel

If a computer has multiple CPUs, or has a CPU with multiple cores, it is possible for MOSEK to take advantage of this to speed up solution times.

8.6.1. Thread safety

The MOSEK API is thread-safe provided that a task is only modified or accessed from one thread at any given time — accessing two separate tasks from two separate threads at the same time is safe. Sharing an environment between threads is safe.

8.6.2. The parallelized interior-point optimizer

The interior-point optimizer is capable of using multiple CPUs or cores. This implies that whenever the MOSEK interior-point optimizer solves an optimization problem, it will try to divide the work so that each CPU gets a share of the work. The user decides how many CPUs MOSEK should exploit.

It is not always possible to divide the work equally, and often parts of the computations and the coordination of the work is processed sequentially, even if several CPUs are present. Therefore, the speed-up obtained when using multiple CPUs is highly problem dependent. However, as a rule of thumb, if the problem solves very quickly, i.e. in less than 60 seconds, it is not advantageous to use the parallel option.

The MSK_IPAR_INTPNT_NUM_THREADS parameter sets the number of threads (and therefore the number of CPUs) that the interior point optimizer will use.

8.6.3. The concurrent optimizer

An alternative to the parallel interior-point optimizer is the concurrent optimizer. The idea of the concurrent optimizer is to run multiple optimizers on the same problem concurrently, for instance, it allows you to apply the interior-point and the dual simplex optimizers to a linear optimization problem concurrently. The concurrent optimizer terminates when the first of the applied optimizers has terminated successfully, and it reports the solution of the fastest optimizer. In that way a new optimizer has been created which essentially performs as the fastest of the interior-point and the dual simplex optimizers.Hence, the concurrent optimizer is the best one to use if there are multiple optimizers available in MOSEK for the problem and you cannot say beforehand which one will be faster.

Note in particular that any solution present in the task will also be used for hot-starting the simplex algorithms. One possible scenario would therefore be running a hot-start dual simplex in parallel with interior point, taking advantage of both the stability of the interior-point method and the ability of the simplex method to use an initial solution.

By setting the

MSK_IPAR_OPTIMIZER

parameter to

MSK_OPTIMIZER_CONCURRENT

the concurrent optimizer chosen.

The number of optimizers used in parallel is determined by the

MSK_IPAR_CONCURRENT_NUM_OPTIMIZERS.

parameter. Moreover, the optimizers are selected according to a preassigned priority with optimizers having the highest priority being selected first. The default priority for each optimizer is shown in Table 8.6.3.

Table 8.4: Default priorities for optimizer selection in concurrent optimization.

For example, setting the MSK_IPAR_CONCURRENT_NUM_OPTIMIZERS parameter to 2 tells the concurrent optimizer to the apply the two optimizers with highest priorities: In the default case that means the interior-point optimizer and one of the simplex optimizers.

8.6.3.1. Concurrent optimization through the API

The following example shows how to call the concurrent optimizer through the API.

/* Copyright: Copyright (c) 1998-2011 MOSEK ApS, Denmark. All rights reserved. File: concurrent1.c Purpose: To demonstrate how to solve a problem with the concurrent optimizer. */ #include <stdio.h> #include "mosek.h" static void MSKAPI printstr(void *handle, char str[]) { printf("%s",str); } /* printstr */ int main(int argc,char *argv[]) { MSKenv_t env; MSKtask_t task; MSKintt r = MSK_RES_OK; /* Create mosek environment. */ r = MSK_makeenv(&env,NULL,NULL,NULL,NULL); if ( r==MSK_RES_OK ) MSK_linkfunctoenvstream(env,MSK_STREAM_LOG,NULL,printstr); /* Initialize the environment. */ r = MSK_initenv(env); if ( r==MSK_RES_OK ) r = MSK_maketask(env,0,0,&task); if ( r==MSK_RES_OK ) MSK_linkfunctotaskstream(task,MSK_STREAM_LOG,NULL,printstr); if (r == MSK_RES_OK) r = MSK_readdata(task,argv[1]); MSK_putintparam(task,MSK_IPAR_OPTIMIZER,MSK_OPTIMIZER_CONCURRENT); MSK_putintparam(task,MSK_IPAR_CONCURRENT_NUM_OPTIMIZERS,2); if (r == MSK_RES_OK) r = MSK_optimize(task); MSK_solutionsummary(task,MSK_STREAM_LOG); MSK_deletetask(&task); MSK_deleteenv(&env); printf("Return code: %d (0 means no error occured.)\n",r); return ( r ); } /* main */

8.6.4. A more flexible concurrent optimizer

MOSEK also provides a more flexible method of concurrent optimization by using the function MSK_optimizeconcurrent. The main advantages of this function are that it allows the calling application to assign arbitrary values to the parameters of each tasks, and that call-back functions can be attached to each task. This may be useful in the following situation: Assume that you know the primal simplex optimizer to be the best optimizer for your problem, but that you do not know which of the available selection strategies (as defined by the MSK_IPAR_SIM_PRIMAL_SELECTION parameter) is the best. In this case you can solve the problem with the primal simplex optimizer using several different selection strategies concurrently.

An example demonstrating the usage of the MSK_optimizeconcurrent function is included below. The example solves a single problem using the interior-point and primal simplex optimizers in parallel.

/* Copyright: Copyright (c) 1998-2011 MOSEK ApS, Denmark. All rights reserved. File: concurrent2.c Purpose: To demonstrate a more flexible interface for concurrent optimization. */ #include "mosek.h" static void MSKAPI printstr(void *handle, char str[]) { printf("simplex: %s",str); } /* printstr */ static void MSKAPI printstr2(void *handle, char str[]) { printf("intrpnt: %s",str); } /* printstr */ #define NUMTASKS 1 int main(int argc,char **argv) { MSKintt r=MSK_RES_OK,i; MSKenv_t env; MSKtask_t task; MSKtask_t task_list[NUMTASKS]; /* Create mosek environment. */ r = MSK_makeenv(&env,NULL,NULL,NULL,NULL); if ( r==MSK_RES_OK ) MSK_linkfunctoenvstream(env,MSK_STREAM_LOG,NULL,printstr); /* Initialize the environment. */ if ( r==MSK_RES_OK ) r = MSK_initenv(env); /* Create a task for each concurrent optimization. The 'task' is the master task that will hold the problem data. */ if ( r==MSK_RES_OK ) r = MSK_maketask(env,0,0,&task); if (r == MSK_RES_OK) r = MSK_maketask(env,0,0,&task_list[0]); if (r == MSK_RES_OK) r = MSK_readdata(task,argv[1]); /* Assign different parameter values to each task. In this case different optimizers. */ if (r == MSK_RES_OK) r = MSK_putintparam(task, MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_PRIMAL_SIMPLEX); if (r == MSK_RES_OK) r = MSK_putintparam(task_list[0], MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_INTPNT); /* Assign call-back functions to each task */ if (r == MSK_RES_OK) MSK_linkfunctotaskstream(task,MSK_STREAM_LOG,NULL,printstr); if (r == MSK_RES_OK) MSK_linkfunctotaskstream(task_list[0], MSK_STREAM_LOG, NULL, printstr2); if (r == MSK_RES_OK) r = MSK_linkfiletotaskstream(task, MSK_STREAM_LOG, "simplex.log", 0); if (r == MSK_RES_OK) r = MSK_linkfiletotaskstream(task_list[0], MSK_STREAM_LOG, "intpnt.log", 0); /* Optimize task and task_list[0] in parallel. The problem data i.e. C, A, etc. is copied from task to task_list[0]. */ if (r == MSK_RES_OK) r = MSK_optimizeconcurrent ( task, task_list, NUMTASKS); printf ("Return Code = %d\n",r); MSK_solutionsummary(task, MSK_STREAM_LOG); return r; }

8.7. Understanding solution quality

MOSEK will, in general, not produce an exact optimal solution; for efficiency reasons computations are performed in finite precision. This means that it is important to evaluate the quality of the reported solution. To evaluate the solution quality inspect the following properties:

Ideally, the primal and dual solutions should only violate the constraints of their respective problem slightly and the primal and dual objective values should be close. This should be evaluated in the context of the problem: How good is the data precision in the problem, and how exact a solution is required.

8.7.1. The solution summary

The solution summary is a small display generated by MOSEK that makes it easy to check the quality of the solution.

The function MSK_solutionsummary should be used to generate solution summary.

8.7.1.1. The optimal case

The solution summary has the format

Problem status  : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal - objective: 5.5018458883e+03   eq. infeas.: 1.20e-12 max bound infeas.: 2.31e-14
Dual   - objective: 5.5018458883e+03   eq. infeas.: 1.15e-14 max bound infeas.: 7.11e-15
i.e. it shows status information, objective values and quality measures for the primal and dual solutions.

Assumeing that we are solving a linear optimization problem and referring to the problems (7.1.1) and (7.1.2), the interpretation of the solution summary is as follows:

  • Problem status: The status of the problem.
  • Solution status: The status of the solution.
  • Primal objective: The primal objective value.
  • Primal eq. infeas: [[MathCmd 558]].
  • Primal max bound infeas.: [[MathCmd 559]].
  • Dual objective: The dual objective value.
  • Dual eq. infeas: [[MathCmd 560]].
  • Dual max bound infeas.: [[MathCmd 561]].

In the solution summary above the solution is classified as [[MathCmd 562]], meaning that the solution should be a good approximation to the true optimal solution. This seems very reasonable since the primal and dual solutions only violate their respective constraints slightly. Moreover, the duality gap is small, i.e. the primal and dual objective values are almost identical.

8.7.1.2. The primal infeasible case

For an infeasible problem the solution summary might look like this:

Problem status  : PRIMAL_INFEASIBLE
Solution status : PRIMAL_INFEASIBLE_CER
Primal - objective: 0.0000000000e+00   eq. infeas.: 0.00e+00 max bound infeas.: 0.00e+00
Dual   - objective: 1.0000000000e+02   eq. infeas.: 0.00e+00 max bound infeas.: 0.00e+00

It is known that if the problem is primal infeasible then an infeasibility certificate exists, which is a solution to the problem (7.1.3) having a positive objective value. Note that the primal solution plays no role and only the dual solution is used to specify the certificate.

Therefore, in the primal infeasible case the solution summery should report how good the dual solution is to the problem (7.1.3). The interpretation of the solution summary is as follows:

  • Problem status: The status of the problem.
  • Solution status: The status of the solution.
  • Primal objective: Should be ignored.
  • Primal eq. infeas: Should be ignored.
  • Primal max bound infeas.: Should be ignored.
  • Dual objective: [[MathCmd 563]].
  • Dual eq. infeas: [[MathCmd 564]].
  • Dual max bound infeas.: [[MathCmd 565]].

Please note that

  • any information about the primal solution should be ignored.
  • the dual objective value should be strictly positive if primal problem is minimization problem. Otherwise it should be strictly negative.
  • the bigger the ratio

    \begin{displaymath}\nonumber{}\frac{(l^{c})^{T} s_{l}^{c} - (u^{c})^{T} s_{u}^{c} + (l^{x})^{T} s_{l}^{x} - (u^{x})^{T} s_{u}^{x}}{\max (\left\|-y+s_{l}^{c}-s_{u}^{c};A^{T} y + s_{l}^{x} - s_{u}^{x} -0\right\|_{\infty },\max (-s_{l}^{c};-s_{u}^{c};-s_{l}^{x};-s_{u}^{x}))}\end{displaymath}

    is, the better the certificate is. The reason is that a certificate is a ray, and hence only the direction is important. Therefore, in principle, the certificate should be normalized before using it.

Please see Section 10.2 for more information about certificates of infeasibility.

8.7.2. Retrieving solution quality information with the API

Information about the solution quality may be retrieved in the API with the help of the following functions:

The example below shows how to use these function to determine the quality of the solution.

/* Copyright: Copyright (c) 1998-2011 MOSEK ApS, Denmark. All rights reserved. File: solutionquality.c Purpose: To demonstrate how to examine the quality of a solution. */ #include "mosek.h" #include "math.h" static void MSKAPI printstr(void *handle, char str[]) { printf("%s",str); } /* printstr */ double double_max(double arg1,double arg2) { return arg1<arg2 ? arg2 : arg1; } int main (int argc, char * argv[]) { MSKtask_t task = NULL; MSKenv_t env = NULL; MSKrescodee r = MSK_RES_OK; if (argc <= 1) { printf ("Missing argument. The syntax is:\n"); printf (" simple inputfile [ solutionfile ]\n"); } else { /* Create the mosek environment. The `NULL' arguments here, are used to specify customized memory allocators and a memory debug file. These can safely be ignored for now. */ r = MSK_makeenv(&env, NULL, NULL, NULL, NULL); /* Initialize the environment */ if ( r==MSK_RES_OK ) MSK_initenv (env); /* Create a task object linked to the environment env. Initially we create it with 0 variables and 0 columns, since we do not know the size of the problem. */ if ( r==MSK_RES_OK ) r = MSK_maketask (env, 0,0, &task); if (r == MSK_RES_OK) MSK_linkfunctotaskstream(task,MSK_STREAM_LOG,NULL,printstr); /* We assume that a problem file was given as the first command line argument (received in `argv'). */ if ( r==MSK_RES_OK ) r = MSK_readdata (task, argv[1]); /* Solve the problem */ if ( r==MSK_RES_OK ) { MSKrescodee trmcode; MSK_optimizetrm(task,&trmcode); } /* Print a summary of the solution. */ MSK_solutionsummary(task, MSK_STREAM_MSG); if (r == MSK_RES_OK) { MSKprostae prosta; MSKsolstae solsta; MSKrealt primalobj,maxpbi,maxpcni,maxpeqi,maxinti, dualobj, maxdbi, maxdcni, maxdeqi; MSKintt isdef; MSKsoltypee whichsol = MSK_SOL_BAS; int accepted = 1; MSK_getsolutioninf ( task, whichsol, &prosta, &solsta, &primalobj, &maxpbi, &maxpcni, &maxpeqi, &maxinti, &dualobj, &maxdbi, &maxdcni, &maxdeqi); switch(solsta) { case MSK_SOL_STA_OPTIMAL: case MSK_SOL_STA_NEAR_OPTIMAL: { double max_primal_infeas = 0.0; /* maximal primal infeasibility */ double max_dual_infeas = 0.0; /* maximal dual infeasibility */ double obj_gap = fabs(dualobj-primalobj); max_primal_infeas = double_max(max_primal_infeas,maxpbi); max_primal_infeas = double_max(max_primal_infeas,maxpcni); max_primal_infeas = double_max(max_primal_infeas,maxpeqi); max_dual_infeas = double_max(max_dual_infeas,maxdbi); max_dual_infeas = double_max(max_dual_infeas,maxdcni); max_dual_infeas = double_max(max_dual_infeas,maxdeqi); /* Assume the application needs the solution to be within 1e-6 ofoptimality in an absolute sense. Another approach would be looking at the relative objective gap */ printf("Objective gap: %e\n",obj_gap); if (obj_gap > 1e-6) { printf("Warning: The objective gap is too large."); accepted = 0; } printf("Max primal infeasibility: %e\n", max_primal_infeas); printf("Max dual infeasibility: %e\n" , max_dual_infeas); /* We will accept a primal infeasibility of 1e-8 and dual infeasibility of 1e-6 */ if (max_primal_infeas > 1e-8) { printf("Warning: Primal infeasibility is too large"); accepted = 0; } if (max_dual_infeas > 1e-6) { printf("Warning: Dual infeasibility is too large"); accepted = 0; } } if (accepted && r == MSK_RES_OK) { MSKintt numvar,j; MSKrealt *xx = NULL; MSK_getnumvar(task,&numvar); xx = (double *) malloc(numvar*sizeof(MSKrealt)); MSK_getsolutionslice(task, MSK_SOL_BAS, /* Request the basic solution. */ MSK_SOL_ITEM_XX,/* Which part of solution. */ 0, /* Index of first variable. */ numvar, /* Index of last variable+1. */ xx); printf("Optimal primal solution\n"); for(j=0; j<numvar; ++j) printf("x[%d]: %e\n",j,xx[j]); free(xx); } else { /* Print detailed information about the solution */ if (r == MSK_RES_OK) r = MSK_analyzesolution(task,MSK_STREAM_LOG,whichsol); } break; case MSK_SOL_STA_DUAL_INFEAS_CER: case MSK_SOL_STA_PRIM_INFEAS_CER: case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER: case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER: printf("Primal or dual infeasibility certificate found.\n"); break; case MSK_SOL_STA_UNKNOWN: printf("The status of the solution could not be determined.\n"); break; default: printf("Other solution status"); break; } } else { printf("Error while optimizing.\n"); } MSK_deletetask(&task); MSK_deleteenv(&env); } return r; }
Wed Feb 29 16:08:51 2012