e04pt: FL CL CPP AD

NAG FL Interface
e04ptf (handle_solve_socp_ipm)

Note: this routine uses optional parameters to define choices in the problem specification and in the details of the algorithm. If you wish to use default settings for all of the optional parameters, you need only read Sections 1 to 10 of this document. If, however, you wish to reset some or all of the settings please refer to Section 11 for a detailed description of the algorithm and to Section 12 for a detailed description of the specification of the optional parameters.

Keyword Search:

NAG Library Manual, Mark 27

Interfaces: FL CL CPP AD

NAG FL Interface Introduction

E04 (Opt) Chapter Contents

E04 (Opt) Chapter Introduction

e04pt: FL CL CPP AD

▸▿ Contents

1 Purpose

2 Specification

▸▿ 3 Description

3.1 Structure of the Lagrangian Multipliers

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

▸▿ 9 Further Comments

▸▿ 9.1 Formulating Problems as SOCPs

9.1.1 Quadratically Constrained Quadratic Programming

9.1.2 Norm Minimization Problems

9.1.3 Robust Linear Programming

9.2 Description of the Printed Output

▸▿ 10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

▸▿ 11 Algorithmic Details

11.1 Homogeneous Self-Dual Algorithm

11.2 The Nesterov–Todd Search Direction

11.3 Mehrotra's Predictor-Corrector Method

11.4 Solving the KKT System

▸▿ 11.5 Stopping Criteria

11.5.1 Convergence-optimal termination

11.5.2 Infeasibility/Unboundedness Detection

11.6 Further Details

▸▿ 12 Optional Parameters

12.1 Description of the Optional Parameters

1 Purpose

e04ptf is a solver from the NAG optimization modelling suite for large-scale second-order cone programming (SOCP) problems. It is based on an interior point method (IPM).

2 Specification

Fortran Interface

Subroutine e04ptf (

handle, nvar, x, nnzu, u, nnzuc, uc, rinfo, stats, monit, iuser, ruser, cpuser, ifail)

Integer, Intent (In)	::	nvar, nnzu, nnzuc
Integer, Intent (Inout)	::	iuser(*), ifail
Real (Kind=nag_wp), Intent (Inout)	::	x(nvar), u(nnzu), uc(nnzuc), ruser(*)
Real (Kind=nag_wp), Intent (Out)	::	rinfo(100), stats(100)
Type (c_ptr), Intent (In)	::	handle, cpuser
External	::	monit

C Header Interface

#include <nag.h>

void

e04ptf_ (void **handle, const Integer *nvar, double x[], const Integer *nnzu, double u[], const Integer *nnzuc, double uc[], double rinfo[], double stats[],
void (NAG_CALL *monit)(void **handle, const double rinfo[], const double stats[], Integer iuser[], double ruser[], void **cpuser, Integer *inform),
Integer iuser[], double ruser[], void **cpuser, Integer *ifail)

The routine may be called by the names e04ptf or nagf_opt_handle_solve_socp_ipm.

3 Description

e04ptf solves a large-scale SOCP optimization problem in the following form

\begin{array}{l} \underset{x \in ℝ^{n}}{minimize} & c^{T} x \\ subject to & l_{A} \leq A x \leq u_{A}, \\ l_{x} \leq x \leq u_{x}, \\ x \in K, \end{array}

(1)

where

K = K^{n_{1}} \times \dots \times K^{n_{r}} \times ℝ^{n_{l}}

is a Cartesian product of

r

quadratic (second-order type) cones and

n_{l}

-dimensional real space, and

n = \sum_{i = 1}^{r} n_{i} + n_{l}

is the number of decision variables. Here

c

x

l_{x}

and

u_{x}

are

n

-dimensional vectors,

A

is an

m

n

sparse matrix, and

l_{A}

and

u_{A}

are

m

-dimensional vectors. Note that

x \in K

partitions subsets of variables into quadratic cones and each

K^{n_{i}}

can be either a quadratic cone or a rotated quadratic cone. These are defined as follows:

•Quadratic cone:
$K_{q}^{n_{i}} ≔ \{z = (z_{1}, z_{2}, \dots, z_{n_{i}}) \in ℝ^{n_{i}} : z_{1}^{2} \geq \sum_{j = 2}^{n_{i}} z_{j}^{2}, z_{1} \geq 0\} .$ (2)

•Rotated quadratic cone:

K_{r}^{n_{i}} ≔ \{z = (z_{1}, z_{2}, \dots, z_{n_{i}}) \in ℝ^{n_{i}} : 2 z_{1} z_{2} \geq \sum_{j = 3}^{n_{i}} z_{j}^{2}, z_{1} \geq 0, z_{2} \geq 0\} .

(3)

e04ptf solves SOCP problems stored as a handle. The handle points to an internal data structure which defines the problem and serves as a means of communication for routines in the NAG optimization modelling suite. First, the problem handle is initialized by calling e04raf. Then some of the routines e04rbf, e04ref, e04rff, e04rhf or e04rjf may be called to formulate the quadratic cones, linear objective function, bounds of the variables, and the block of linear constraints, respectively. Alternatively, the whole model can be loaded from a file by e04saf. When the handle is no longer needed, e04rzf should be called to destroy it and deallocate the memory held within. See Section 3.1 in the E04 Chapter Introduction for more details about the NAG optimization modelling suite.

The solver method can be modified by various optional parameters (see Section 12) which can be set by e04zmf and e04zpf anytime between the initialization of the handle and a call to the solver. Once the solver has finished, options may be modified for the next solve. The solver may be called repeatedly with various optional parameters.

The optional parameter Task may be used to switch the problem to maximization or to ignore the objective function and find only a feasible point.

Several options may have significant impact on the performance of the solver. Even if the defaults were chosen to suit the majority of problems, it is recommended that you experiment in order to find the most suitable set of options for a particular problem, see Sections 11 and 12 for further details.

3.1 Structure of the Lagrangian Multipliers

The algorithm works internally with estimates of both the decision variables, denoted by

x

, and the Lagrangian multipliers (dual variables), denoted by

u

for bound and linear constraints, and

u c

for quadratic cone constraints.

If the simple bounds have been defined (e04rhf was successfully called), the first

2 n

elements of

u

belong to the corresponding Lagrangian multipliers, interleaving a multiplier for the lower and the upper bound for each

x_{i}

. If any of the bounds were set to infinity, the corresponding Lagrangian multipliers are set to

0

and may be ignored.

Similarly, the following

2 m

elements of

u

belong to multipliers for the linear constraints (if e04rjf has been successfully called). The organization is the same, i.e., the multipliers for each constraint for the lower and upper bounds are alternated and zeros are used for any missing (infinite bound) constraints.

Some solvers merge multipliers for both lower and upper inequality into one element whose sign determines the inequality. Negative multipliers are associated with the upper bounds and positive with the lower bounds. An equivalent result can be achieved with this storage scheme by subtracting the upper bound multiplier from the lower one. This is also consistent with equality constraints.

Finally, the elements of

u c

are the corresponding Lagrangian multipliers for the variables in the quadratic cone constraints that have been defined by e04rbf. All multipliers are stored next to each other in array uc in the same order as the cone constraints were defined by e04rbf. For example, if the first cone constraint contains variables

x_{4}

x_{2}

x_{3}

and the second cone constraint contains variables

x_{1}

x_{7}

x_{6}

x_{5}

, then the dimension of array uc must be

7

and the first

3

elements are the corresponding Lagrangian multipliers for the cone composed of

x_{4}

x_{2}

x_{3}

, followed by

4

elements that are the corresponding Lagrangian multipliers for the cone of

x_{1}

x_{7}

x_{6}

x_{5}

4 References

Alizadeh F and Goldfarb D (2003) Second-order cone programming Mathematical programming 95(1) 3–51

Andersen E D, Roos C and Terlaky T (2003) On implementing a primal-dual interior-point method for conic quadratic optimization Mathematical programming 95(2) 249–277

Goldfarb D and Scheinberg K (2005) Product-form Cholesky factorization in interior point methods for second-order cone programming Mathematical programming 103(1) 153–179

Goldman A J and Tucker A W (1956) Theory of linear programming Linear inequalities and related systems 38 53–97

Hogg J D and Scott J A (2011) HSL MA97: a bit-compatible multifrontal code for sparse symmetric systems RAL Technical Report. RAL-TR-2011-024

HSL (2011) A collection of Fortran codes for large-scale scientific computation http://www.hsl.rl.ac.uk/

Karypis G and Kumar V (1998) A fast and high quality multilevel scheme for partitioning irregular graphs SIAM J. Sci. Comput. 20(1) 359–392

Lobo M S, Vandenberghe L, Boyd S and Levret H (1998) Applications of second-order cone programming Linear Algebra and its Applications 284(1-3) 193–228

Lustig I J, Marsten R E and Shanno D F (1992) On implementing Mehrotra's predictor–corrector interior-point method for linear programming SIAM J. Optim. 2(3) 435–449

Mehrotra S (1992) On the implementation of a primal-dual interior point method SIAM J. Optim. 2 575–601

Nesterov Y E and Todd M J (1997) Self-scaled barriers and interior-point methods for convex programming Mathematics of Operations research 22(1) 1–42

Nesterov Y E and Todd M J (1998) Primal-dual interior-point methods for self-scaled cones SIAM J. Optim. 8(2) 324–364

Nocedal J and Wright S J (2006) Numerical Optimization (2nd Edition) Springer Series in Operations Research, Springer, New York

Sturm J F (2002) Implementation of Interior Point Methods for Mixed Semidefinite and Second Order Cone Optimization Problems Optimization Methods and Software 17(6) 151–171

Xu X, Hung P-F and Ye Y (1996) A simplified homogeneous and self-dual linear programming algorithm and its implementation Annals of Operations Research 62(1) 151–171

5 Arguments

1: $handle$ – Type (c_ptr) Input

On entry: the handle to the problem. It needs to be initialized by e04raf and the problem formulated by some of the routines e04rbf, e04ref, e04rff, e04rhf and e04rjf. It must not be changed between calls to the NAG optimization modelling suite.

2: $nvar$ – Integer Input

On entry:

n

, the number of variables in the problem.

3: $x (nvar)$ – Real (Kind=nag_wp) array Input/Output

On entry: the input of x is reserved for future releases of the NAG Library and it is ignored at the moment.

On exit: the final values of the variables

x

4: $nnzu$ – Integer Input

On entry: the dimension of array u.

nnzu = 0

, u will not be referenced; otherwise it needs to match the dimension of constraints defined by e04rhf and e04rjf as explained in Section 3.1.

Constraint:

nnzu \geq 0

5: $u (nnzu)$ – Real (Kind=nag_wp) array Input/Output

Note: if

nnzu > 0

, u holds Lagrange multipliers (dual variables) for the bound constraints and linear constraints. If

nnzu = 0

, u will not be referenced.

On entry: the input of u is reserved for future releases of the NAG Library and it is ignored at the moment.

On exit: the final values of the variables

u

6: $nnzuc$ – Integer Input

On entry: the dimension of array uc.

nnzuc = 0

, uc will not be referenced; otherwise it needs to match the total number of cone variables defined by e04rbf as explained in Section 3.1.

Constraint:

nnzuc \geq 0

7: $uc (nnzuc)$ – Real (Kind=nag_wp) array Input/Output

Note: if

nnzuc > 0

, uc holds Lagrange multipliers (dual variables) for second-order cones as defined by e04rbf. If

nnzuc = 0

, uc will not be referenced.

On entry: the input of uc is reserved for future releases of the NAG Library and it is ignored at the moment.

On exit: the final values of the variables

uc

8: $rinfo (100)$ – Real (Kind=nag_wp) array Output

On exit: error measures and various indicators of the algorithm (see Section 11 for details) as given in the table below:

$1$	Value of the primal objective.
$2$	Value of the dual objective.
$3$	Flag indicating the system formulation used by the solver, $0$ : augmented system, $1$ : normal equation.
$4$	Factorization type, $3$ : Cholesky, $4$ : Bunch–Parlett.
$5$ – $14$	Not referenced in this solver.
$15$	Relative primal infeasibility, see Section 11.5.1.
$16$	Relative duality gap, see Section 11.5.1.
$17$	Relative dual infeasibility, see Section 11.5.1.
$18$	Accuracy, see Section 11.5.1.
$19$	$τ$ , see (23).
$20$	$κ$ , see (23).
$21$	Step length.
$22$ – $100$	Reserved for future use.

9: $stats (100)$ – Real (Kind=nag_wp) array Output

On exit: solver statistics as given in the table below. Note that times are measured in seconds, see optional parameter Stats Time.

$1$	Number of iterations.
$2$	Not referenced.
$3$	Total number of iterative refinements performed.
$4$	Value of the perturbation added to the diagonal in the normal equation formulation or the augmented system formulation.
$5$	Total number of factorizations performed.
$6$	Total time spent in the solver.
$7$	Time spent in the presolve phase.
$8$	Time spent in the last iteration.
$9$	Total time spent factorizing the system matrix.
$10$	Total time spent backsolving the system matrix.
$11$	Not referenced.
$12$	Time spent in the initialization phase.
$13$	Number of nonzeros in the system matrix.
$14$	Number of nonzeros in the system matrix factor.
$15$	Maximum error of the backsolve.
$16$	Number of columns in $A$ considered dense by the solver.
$17$	Number of conic constraints considered dense by the solver.
$18$ – $100$	Reserved for future use.

10: $monit$ – Subroutine, supplied by the NAG Library or the user. External Procedure

monit is provided to enable you to monitor the progress of the optimization. It is invoked at the end of every

i

th iteration where

i

is given by the optional parameter SOCP Monitor Frequency (the default is

0

, monit is not called).

monit may be the dummy subroutine e04ptu (e04ptu is included in the NAG Library).

The specification of monit is:

Fortran Interface

Subroutine monit (

handle, rinfo, stats, iuser, ruser, cpuser, inform)

Integer, Intent (Inout)	::	iuser(*), inform
Real (Kind=nag_wp), Intent (In)	::	rinfo(100), stats(100)
Real (Kind=nag_wp), Intent (Inout)	::	ruser(*)
Type (c_ptr), Intent (In)	::	handle, cpuser

C Header Interface

void	monit_ (void handle, const double rinfo[], const double stats[], Integer iuser[], double ruser[], void cpuser, Integer *inform)

1: $handle$ – Type (c_ptr) Input: On entry: the handle to the problem. It needs to be initialized by e04raf and must not be changed between calls to the NAG optimization modelling suite.
2: $rinfo (100)$ – Real (Kind=nag_wp) array Input: On entry: error measures and various indicators at the end of the current iteration as described in rinfo.
3: $stats (100)$ – Real (Kind=nag_wp) array Input: On entry: solver statistics at the end of the current iteration as described in stats, however, elements $3$ , $5$ , $9$ , $10$ and $15$ refer to the quantities in the last iteration rather than accumulated over all iterations through the whole algorithm run.
4: $iuser (*)$ – Integer array User Workspace
5: $ruser (*)$ – Real (Kind=nag_wp) array User Workspace
6: $cpuser$ – Type (c_ptr) User Workspace: monit is called with the arguments iuser, ruser and cpuser as supplied to e04ptf. You should use the arrays iuser and ruser, and the data handle cpuser to supply information to monit.
7: $inform$ – Integer Input/Output: On entry: a non-negative value.

On exit: must be set to a value describing the action to be taken by the solver on return from monit. Specifically, if the value is negative the solution of the current problem will terminate immediately with $ifail = 20$ ; otherwise, computations will continue.

monit must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e04ptf is called. Arguments denoted as Input must not be changed by this procedure.

11: $iuser (*)$ – Integer array User Workspace

12: $ruser (*)$ – Real (Kind=nag_wp) array User Workspace

13: $cpuser$ – Type (c_ptr) User Workspace

iuser, ruser and cpuser are not used by e04ptf, but are passed directly to monit and may be used to pass information to this routine. If you do not need to reference cpuser, it should be initialized to c_null_ptr.

14: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

- 1

1

to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of

0

causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of

- 1

means that an error message is printed while a value of

1

means that it is not.

If halting is not appropriate, the value

- 1

1

is recommended. If message printing is undesirable, then the value

1

is recommended. Otherwise, the value

- 1

is recommended since useful values can be provided in some output arguments even when

ifail \neq 0

on exit. When the value $- 1$ or $1$ is used it is essential to test the value of ifail on exit.

On exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

ifail = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

Note: in some cases e04ptf may return useful information.

$ifail = 1$: The supplied handle does not define a valid handle to the data structure for the NAG optimization modelling suite. It has not been initialized by e04raf or it has been corrupted.

$ifail = 2$: The problem is already being solved.

This solver does not support the model defined in the handle.

$ifail = 4$: On entry, $nvar = 〈value〉$ , expected $value = 〈value〉$ .
Constraint: nvar must match the value given during initialization of handle.

$ifail = 5$: On entry, $nnzu = 〈value〉$ .
nnzu does not match the size of the Lagrangian multipliers for constraints.
The correct value is $0$ for no constraints.

On entry, $nnzu = 〈value〉$ .
nnzu does not match the size of the Lagrangian multipliers for constraints.
The correct value is either $0$ or $〈value〉$ .

On entry, $nnzuc = 〈value〉$ .
nnzuc does not match the size of the Lagrangian multipliers for second-order cone constraints.
$nnzuc = 0$ when there are no second-order cone constraints.

On entry, $nnzuc = 〈value〉$ .
nnzuc does not match the size of the Lagrangian multipliers for second-order cone constraints.
The correct value is either $0$ or $〈value〉$ .

$ifail = 20$: User requested termination during a monitoring step.

$ifail = 22$: Maximum number of iterations exceeded.

$ifail = 23$: The solver terminated after the maximum time allowed was exceeded.
Maximum number of seconds exceeded. Use optional parameter Time Limit to reset the limit.

$ifail = 24$: No progress, stopping early.
The solver predicted that it is unable to make further progress and stopped prematurely. This might be due to the scaling of the problem, its conditioning or numerical difficulties.

$ifail = 50$: Suboptimal solution.
The solver predicted that it is unable to reach a better estimate of the solution. However, the error measures indicate that the point is a reasonable approximation.

$ifail = 51$: The problem was found to be primal infeasible.

$ifail = 52$: The problem was found to be dual infeasible.
This error indicates that the primal problem is unbounded or infeasible.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

The accuracy of the solution is determined by optional parameters SOCP Stop Tolerance and SOCP Stop Tolerance 2

ifail = 0

on the final exit, the returned point satisfies Karush–Kuhn–Tucker (KKT) conditions to the requested accuracy (under the default settings close to

\sqrt{ε}

) and thus it is a good estimate of the solution. If

ifail = 50

, some of the convergence conditions were not fully satisfied but the point is a reasonable estimate and still usable. Please refer to Section 11.5 and the description of the particular options.

8 Parallelism and Performance

e04ptf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

e04ptf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.

Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

9.1 Formulating Problems as SOCPs

This SOCP solver can solve several common convex problems covering a large variety of applications. However, in this release, the cone constraints need to be specified explicitly which means that in certain cases a reformulation is needed. In this section, we cover QCQP, norm minimization problems and robust linear programming, see Alizadeh and Goldfarb (2003) and Lobo et al. (1998) for further details.

9.1.1 Quadratically Constrained Quadratic Programming

Convex quadratically constrained quadratic programming (QCQP) appears in applications such as modern portfolio theory and wireless sensor network localization. The general convex QCQP problem has the following form

\begin{array}{l} \underset{x \in ℝ^{n}}{minimize} & \frac{1}{2} x^{T} P_{0} x + q_{0}^{T} x + r_{0} \\ subject to & \frac{1}{2} x^{T} P_{i} x + q_{i}^{T} x + r_{i} \leq 0, i = 1, \dots, p, \end{array}

(4)

where

P_{i} \in ℝ^{n \times n}

, for

i = 0, \dots, p

, are symmetric and positive semidefinite matrices, hence there exist matrices

F_{i} \in ℝ^{k_{i} \times n}

such that

P_{i} = F_{i}^{T} F_{i}, i = 0, 1, \dots, p .

In many practical problems this decomposition is already available. Otherwise it needs to be computed, for example, via Cholesky or eigenvalue decomposition, such as f07fdf for positive definite matrices, and f07kdf or f08fcf for positive semidefinite matrices. Let's introduce new artificial variables

t_{i}

such that

t_{i} + q_{i}^{T} x + r_{i} = 0

, then we have an equivalent characterization of cone constraints as

\frac{1}{2} x^{T} P_{i} x + q_{i}^{T} x + r_{i} \leq 0 ⟷ t_{i} + q_{i}^{T} x + r_{i} = 0 and \frac{1}{2} {‖F_{i} x‖}_{2} \leq t_{i} .

By the definition of rotated quadratic cone (3) we have

\frac{1}{2} {‖F_{i}‖}_{2} \leq t_{i} ⟷ (t_{i}, 1, F_{i} x) \in K_{r}^{k i + 2} .

Therefore, model (4) can be transformed equivalently to the following SOCP problem

\begin{array}{l} \underset{x \in ℝ^{n}, t \in ℝ^{p + 1}}{minimize} & t_{0} + q_{0}^{T} x + r_{0} \\ subject to & t_{i} + q_{i}^{T} x + r_{i} = 0, i = 1, \dots, p, \\ (t_{i}, 1, F_{i} x) \in K_{r}^{k_{i} + 2}, i = 0, \dots, p . \end{array}

(5)

See Section 10 for an example of how to formulate and solve a QCQP problem.

9.1.2 Norm Minimization Problems

Consider the following problem that minimizes the sum of Euclidean norms

\begin{array}{l} \underset{x \in ℝ^{n}}{minimize} & \sum_{i = 1}^{r} {‖A_{i} x + b_{i}‖}_{2}, \end{array}

(6)

where

A_{i} \in ℝ^{n_{i} \times n}

and

b_{i} \in ℝ^{n_{i}}

. Problem (6) can be formulated as SOCP by introducing

r

auxiliary variables

t_{i}

, for

i = 1, \dots, r

, and adding constraints

{‖A_{i} x + b_{i}‖}_{2} \leq t_{i} ⟷ (t_{i}, A_{i} x + b_{i}) \in K_{q}^{n_{i} + 1}, i = 1, \dots, r .

Then the resulting SOCP is

\begin{array}{l} \underset{x \in ℝ^{n}, t \in ℝ^{r}}{minimize} & \sum_{i = 1}^{r} t_{i} \\ subject to & (t_{i}, A_{i} x + b_{i}) \in K_{q}^{n_{i} + 1}, i = 1, \dots, r . \end{array}

(7)

Observe that if (6) had non-negative weights in the sum, the problem would still be an SOCP.

Similarly, minimizing the maximum of Euclidean norms can be expressed as SOCPs. By introducing one auxiliary variable

t \in ℝ

, the problem

\begin{array}{l} \underset{x \in ℝ^{n}}{minimize} & \max_{i = 1, \dots, r} {‖A_{i} x + b_{i}‖}_{2}, \end{array}

(8)

is equivalent to

\begin{array}{l} \underset{x \in ℝ^{n}, t \in ℝ}{minimize} & t \\ subject to & {‖A_{i} x + b_{i}‖}_{2} \leq t, i = 1, \dots, r . \end{array}

(9)

Hence, problem (8) can be cast as the following SOCP

\begin{array}{l} \underset{x \in ℝ^{n}, t \in ℝ}{minimize} & t \\ subject to & (t, A_{i} x + b_{i}) \in K_{q}^{n_{i} + 1}, i = 1, \dots, r . \end{array}

(10)

As an interesting special case, an

l_{1}

-norm minimization problem can also be solved by SOCP. Consider the following unconstrained problem

\begin{array}{l} \underset{x \in ℝ^{n}}{minimize} & {‖A x + b‖}_{1}, \end{array}

(11)

where

A \in ℝ^{m \times n}

and

b \in ℝ^{m}

, introduce an auxiliary variable

u \in ℝ^{m}

such that

A x + b = u

, then problem (11) is transformed to

\begin{array}{l} \underset{x \in ℝ^{n}, u \in ℝ^{m}}{minimize} & {‖u‖}_{1} \\ subject to & A x + b = u . \end{array}

(12)

By adding an auxiliary variable

t \in ℝ^{m}

, the above problem is equivalent to

\begin{array}{l} \underset{x \in ℝ^{n}, u \in ℝ^{m}, t \in ℝ^{m}}{minimize} & \sum_{i = 1}^{m} t_{i} \\ subject to & A x + b = u, \\ |u_{i}| \leq t_{i}, i = 1, \dots, m . \end{array}

(13)

Note that the inequality

|u_{i}| \leq t_{i}

is equivalent to

(t_{i}, u_{i}) \in K_{q}^{2}

, therefore the final SOCP is

\begin{array}{l} \underset{x \in ℝ^{n}, u \in ℝ^{m}, t \in ℝ^{m}}{minimize} & \sum_{i = 1}^{m} t_{i} \\ subject to & A x + b = u, & (t_{i}, u_{i}) \in K_{q}^{2}, i = 1, \dots, m . \end{array}

(14)

9.1.3 Robust Linear Programming

Consider a linear programming problem

\begin{array}{l} \underset{x \in ℝ^{n}}{minimize} & c^{T} x \\ subject to & a_{i}^{T} x \leq b_{i}, i = 1, \dots, r, \end{array}

(15)

where

c

and

b_{i}

are given but there is some uncertainty in parameter

a_{i}

. In such a situation you might want to solve problem (15) in the worst-case sense, i.e., find the best solution

x

with respect to the most adverse choice of

a_{i}

. Introducing uncertainty set to some or all your data and solving the problem in the worst-case scenario helps to avoid high sensitivity of your results even for a small perturbation in the input data. Assume

a_{i}

are known to lie in given ellipsoids around its known centre

{\bar{a}}_{i}

a_{i} \in E_{i} ≔ \{{\bar{a}}_{i} + P_{i} u | {‖u‖}_{2} \leq 1\},

where

P_{i} \in ℝ^{n \times n}

are positive semidefinite matrices, this problem is also known as robust linear programming which can be modelled as

\begin{array}{l} \underset{x \in ℝ^{n}}{minimize} & c^{T} x \\ subject to & a_{i}^{T} x \leq b_{i},   for all a_{i} \in E_{i}, i = 1, \dots, r . \end{array}

(16)

Constraints

a_{i}^{T} x \leq b_{i},   for all a_{i} \in E_{i}, i = 1, \dots, r

are equivalent to

\underset{u \in ℝ^{n}}{maximize} \{{\bar{a}}_{i}^{T} x + {(P_{i} x)}^{T} u | ‖u‖ \leq 1\} \leq b_{i},

{\bar{a}}_{i}^{T} x + \underset{u \in ℝ^{n}}{maximize} \{{(P_{i} x)}^{T} u | ‖u‖ \leq 1\} \leq b_{i} .

Using the definition of the dual norm of the Euclidean norm we can write down the equivalent for problem (16) as

\begin{array}{l} \underset{x \in ℝ^{n}}{minimize} & c^{T} x \\ subject to & {\bar{a}}_{i}^{T} x + ‖P_{i} x‖ \leq b_{i}, i = 1, \dots, r . \end{array}

(17)

By adding auxiliary variables

t_{i}

i = 1, \dots, r

that

{\bar{a}}_{i}^{T} x + t_{i} = b_{i}

, we have the final equivalent SOCP as

\begin{array}{l} \underset{x \in ℝ^{n}, t \in ℝ^{r}}{minimize} & c^{T} x \\ subject to & {\bar{a}}_{i}^{T} x + t_{i} = b_{i}, i = 1, \dots, r, \\ (t_{i}, P_{i} x) \in K_{q}^{n + 1}, i = 1, \dots, r . \end{array}

(18)

Note we can also get SOCP formulation if there is some uncertainty or variation in the parameters

c

and

b_{i}

following a similar procedure.

9.2 Description of the Printed Output

The solver can print information to give an overview of the problem and of the progress of the computation. The output may be sent to two independent streams (files) which are set by optional parameters Print File and Monitoring File. Optional parameters Print Level, Monitoring Level, Print Solution and Print Options determine the exposed level of detail. This allows, for example, a detailed log file to be generated while the condensed information is displayed on the screen.

By default (

Print File = 6

Print Level = 2

), six sections are printed to the standard output:

Header
Optional parameters list (if $Print Options = YES$ )
Problem statistics
Iteration log
Summary
Solution (if $Print Solution = YES$ )

Header

The header is a message indicating the start of the solver. It should look like:

------------------------------------------------
 E04PT, Interior point method for SOCP problems
------------------------------------------------

Optional parameters list

The list shows all options of the solver, each displayed on one line. The output contains the option name, its current value and an indicator for how it was set. The options unchanged from the default setting are noted by ‘d’, options you set are noted by ‘U’, and options reset by the solver are noted by ‘S’. Note that the output format is compatible with the file format expected by e04zpf. The output might look as follows:

Stats Time                    =                 Yes     * U
Task                          =            Minimize     * d
Socp Iteration Limit          =                 100     * d
Socp Presolve                 =                 Yes     * d

Problem statistics

Print Level \geq 2

, statistics on the original and the presolved problems are printed. More detailed statistics, as well as a list of the presolve operations, are also printed for Print Level

3

or above, for example:

Original Problem Statistics

  Number of variables                        2111
  Number of linear constraints               2110
  Number of nonzeros                         3428
  Number of cones                              20


Presolved Problem Statistics

  Number of variables                        1268
  Number of linear constraints               1199
  Number of nonzeros                         1676
  Number of cones                              20

Iteration log

Print Level \geq 2

, the solver prints the status of each iteration.

Print Level = 2

, the output shows the iteration number (

0

represents the starting point), the current primal and dual objective value, convergence measures (primal infeasibility, dual infeasibility and duality gap defined in Section 11.5.1) and the value of the additional variable

τ

(see Section 11.1). The output might look as follows:

------------------------------------------------------------------------
  it|    pobj    |    dobj    |  p.inf  |  d.inf  |  d.gap  |   tau  | I
------------------------------------------------------------------------
   0  2.34871E+00  0.00000E+00  8.89E-01  1.09E-01  1.80E-01  1.0E+00
   1  5.95233E+00  7.97442E+00  1.49E-01  1.83E-02  3.00E-02  1.9E-01
   2  1.71247E+01  1.59748E+01  1.10E-01  1.35E-02  2.22E-02  3.0E-01
   3  2.55291E+01  2.53467E+01  1.61E-02  1.98E-03  3.25E-03  2.9E-01

Print Level = 3

, the solver also prints for each iteration

ρ_{A}

(defined in Section 11.5.1), the value of the variable

κ

(see Section 11.1), the step size, the maximum error of the backsolves performed as well as the total number of iterative refinements performed. The output takes the following form:

----------------------------------------------------------------------------------------------------------------------
  it|    pobj    |    dobj    |  p.inf  |  d.inf  |  d.gap  |  rhoa  |   tau  |  kappa |   step  |  errbs  | nrefi | I
----------------------------------------------------------------------------------------------------------------------
   0  2.34871E+00  0.00000E+00  8.89E-01  1.09E-01  1.80E-01  2.3E+00
   1  5.95233E+00  7.97442E+00  1.49E-01  1.83E-02  3.00E-02  2.3E-01  1.9E-01  9.5E-01  8.44E-01  1.07E-14   9
   2  1.71247E+01  1.59748E+01  1.10E-01  1.35E-02  2.22E-02  6.8E-02  3.0E-01  6.5E-02  4.30E-01  9.24E-15   9
   3  2.55291E+01  2.53467E+01  1.61E-02  1.98E-03  3.25E-03  6.9E-03  2.9E-01  7.9E-03  8.67E-01  3.04E-14   7

Occasionally, when numerical instabilities are too big, the solver will restart the iteration and switch to an augmented system formulation. In such cases the letters RS will be printed in the information column (I).

Print Level > 3

, each iteration produces more information that expands over several lines. This additional information contains:

the method used (normal equation, augmented system);
the number of factorizations performed at the current iteration;
the type of factorization performed (Cholesky, Bunch–Parlett);
the value of the perturbation added to the diagonal in the normal equation formulation or on the zero block in the augmented system formulation;
the total time spent in the iteration if Stats Time is not set to $NO$ .

The output might look as follows:

------------ Details of Iteration   1 ------------
method                           Augmented System
iterative refinements                           9
factorizations                                  1
matrix type                         Bunch-Parlett
diagonal perturbation                    7.00E-08
time iteration                           0.05 sec
--------------------------------------------------

Summary

Once the solver finishes, a detailed summary is produced:

-------------------------------------------------
Status: converged, an optimal solution found
-------------------------------------------------
Final primal objective value         2.688878E+01
Final dual objective value           2.688878E+01
Absolute primal infeasibility        2.264154E-07
Relative primal infeasibility        6.788104E-09
Absolute dual infeasibility          7.639479E-09
Relative dual infeasibility          1.371539E-09
Absolute complementarity gap         2.558237E-08
Relative complementarity gap         8.342957E-10
Iterations                                      8

It starts with the status line of the overall result which matches the ifail value and is followed by the final primal and dual objective values as well as the error measures and iteration count.

Optionally, if

Stats Time = YES

, the timings of the different parts of the algorithm are displayed. It might look as follows:

Timing
  Total time                             0.16 sec
  Presolver                              0.00 sec   (  1.3%)
  Core                                   0.15 sec   ( 98.7%)
    Initialization                       0.00 sec   (  1.4%)
    Factorization                        0.13 sec   ( 88.2%)
    Compute directions                   0.02 sec   ( 10.4%)
  Iterative refinement                   0.01 sec   (  9.7%)

Solution

Print Solution = X

, the values of the primal variables and their bounds on the primary and secondary outputs. It might look as follows:

Primal variables:
x_idx   Lower bound        Value      Upper bound
    1   0.00000E+00    1.02411E-08         inf
    2   0.00000E+00    1.43619E-08         inf
    3   4.00000E+00    1.00000E+01    1.00000E+01
    4   0.00000E+00    2.05523E+00    4.00000E+00
    5  -1.00000E+01   -6.28719E+00    1.00000E+01
    6  -8.00000E+00   -7.49982E+00    8.00000E+00
    7   1.00000E+00    2.08866E+00    3.00000E+00
    8   5.00000E-01    2.52602E+00    5.00000E+00

Print Solution = YES

ALL

, the values of the dual variables are also printed. It should look as follows:

Box bounds dual variables:
x_idx   Lower bound        Value      Upper bound        Value
    1   0.00000E+00    1.03294E+01         inf       0.00000E+00
    2   0.00000E+00    4.77419E+00         inf       0.00000E+00
    3   4.00000E+00    0.00000E+00    1.00000E+01    4.00326E+00
    4   0.00000E+00    0.00000E+00    4.00000E+00    1.88512E-08
    5  -1.00000E+01    9.77434E-09    1.00000E+01    0.00000E+00
    6  -8.00000E+00    1.18996E-07    8.00000E+00    0.00000E+00
    7   1.00000E+00    0.00000E+00    3.00000E+00    2.13077E-08
    8   5.00000E-01    2.00243E-09    5.00000E+00    0.00000E+00

Linear constraints dual variables:
  idx   Lower bound        Value      Upper bound        Value
    1   7.00000E+00    0.00000E+00    9.00000E+00    1.73118E+00
    2  -1.00000E+01    0.00000E+00   -8.00000E+00    1.20039E+00
    3  -1.50000E+01    0.00000E+00   -1.10000E+01    4.30107E-02

Cone constraints dual variables:
idgroup       x_idx       Value
    1             6    2.02570E+00
                  5   -2.02453E-01
                  4   -2.50000E-01
                  7    1.99999E+00

    2             8    7.11750E+00
                  2   -7.11749E+00

10 Example

As described in Section 9.1, second-order cone programming has many applications, however, a reformulation might be needed. This example demonstrates how to formulate and solve a quadratically constrained quadratic programming problem via SOCP. See also Section 10 in e04rbf for a simpler example where the cone constraints are explicitly defined, and Section 10 in e04saf where the input is read from a file.

This example solves a quadratically constrained quadratic programming problem

\begin{array}{l} \underset{x \in ℝ^{3}}{minimize} & \frac{1}{2} x^{T} P_{0} x + q_{0}^{T} x \\ subject to & \frac{1}{2} x^{T} P_{1} x + q_{1}^{T} x + r_{1} \leq 0, \end{array}

where

\begin{matrix} P_{0} = (\begin{matrix} 0.493 & 0.382 & 0.270 \\ 0.382 & 0.475 & 0.448 \\ 0.270 & 0.448 & 0.515 \end{matrix}) & P_{1} = (\begin{matrix} 0.737 & 0.453 & 1.002 \\ 0.453 & 0.316 & 0.635 \\ 1.002 & 0.635 & 1.590 \end{matrix}), \\ q_{0} = (\begin{matrix} 0.847 \\ 0.08 \\ 0.505 \end{matrix}) & q_{1} = (\begin{matrix} 0.065 \\ 0.428 \\ 0.097 \end{matrix}) \end{matrix}

and

r_{1} = 1.276

. We factorize

P_{0} = F_{0}^{T} F_{0}

and

P_{1} = F_{1}^{T} F_{1}

using an eigenvalue decomposition, then we have the following second-order cone programming problem

\begin{array}{l} \underset{x \in ℝ^{3}, t \in ℝ^{2}}{minimize} & t_{0} + q_{0}^{T} x + r_{0} \\ subject to & t_{1} + q_{1}^{T} x + r_{1} = 0, \\ (t_{0}, 1, F_{0} x) \in K_{r}^{5}, \\ (t_{1}, 1, F_{1} x) \in K_{r}^{5} \end{array}

as in Section 9.1.1. By adding artificial variables

u = F_{0} x

v = F_{1} x

y_{0} = 1

, and

y_{1} = 1

to the above formulation, we have the following standard SOCP problem

\begin{array}{l} \underset{x \in ℝ^{3}, t \in ℝ^{2}, u \in ℝ^{3}, v \in ℝ^{3}, y \in ℝ^{2}}{minimize} & t_{0} + q_{0}^{T} x \\ subject to & t_{1} + q_{1}^{T} x + r_{1} = 0, \\ F_{0} x = u, \\ F_{1} x = v, \\ y_{0} = 1, \\ y_{1} = 1, \\ (t_{0}, y_{0}, u) \in K_{r}^{5}, \\ (t_{1}, y_{1}, v) \in K_{r}^{5} . \end{array}

11 Algorithmic Details

This section contains the description of the underlying algorithms used in e04ptf, which implements the standard primal-dual path-following interior point method with Nesterov–Todd scaling and self-dual embedding. For further details, see Nesterov and Todd (1998), Nesterov and Todd (1997) and Andersen et al. (2003).

For simplicity, we consider the following primal second-order cone programming formulation

\begin{array}{l} \underset{x \in ℝ^{n}}{minimize} & c^{T} x \\ subject to & A x = b, \\ x \in \bar{K}, \end{array}

(19)

where

c

x \in ℝ^{n}

b \in ℝ^{m}

A \in ℝ^{m \times n}

with full row rank, and

\bar{K} = K^{n_{1}} \times \dots \times K^{n_{r}} \times ℝ_{+}^{n_{l}}

. The dual formulation for problem (19) is given by

\begin{array}{l} \underset{y \in ℝ^{m}, z \in ℝ^{n}}{maximize} & b^{T} y \\ subject to & A^{T} y + z = c, \\ z \in \bar{K}, \end{array}

(20)

where

y

and

z

denote the dual variables and

\bar{K}

is as defined above (it is a self-dual cone). Solutions of the primal (19) and dual (20) problem are connected by the strong duality theory (see, for example, Nocedal and Wright (2006)) and are characterized by the first-order optimality conditions, the so-called Karush–Kuhn–Tucker (KKT) conditions, which are stated as follows:

\begin{array}{l} A x = b, x \in \bar{K} & (primal feasibility) \\ A^{T} y + z = c, z \in \bar{K} & (dual feasibility) \\ x \circ z = 0 & (complementarity), \end{array}

(21)

where

\circ

is the multiplication operator defined in a special case of a so-called Euclidean Jordan algebra

(ℝ^{n}, \circ)

with the following definition

x \circ y ≔ (\begin{matrix} x^{T} y \\ x_{0} y_{1} + y_{0} x_{1} \\ ⋮ \\ x_{0} y_{n} + y_{0} x_{n} \end{matrix}) .

(22)

If (19) and (20) have a strictly feasible solution (i.e., there is a feasible solution

(\hat{x}, (\hat{y}, \hat{z}))

such that

\hat{x} \in int \bar{K}

and

\hat{z} \in int \bar{K}

), then they both have optimal solutions and the duality gap is zero. Moreover, a feasible solution pair

(x^{*}, y^{*}, z^{*})

is optimal if, and only if, the KKT conditions (21) hold at this point, see Alizadeh and Goldfarb (2003) for more details.

The underlying algorithm applies an iterative method to find an optimal solution

(x^{*}, y^{*}, z^{*})

of the system (21) employing variants of Newton's method and modifying the search direction and step length so that the cone constraints are preserved at every iteration.

11.1 Homogeneous Self-Dual Algorithm

The homogeneous and self-dual (HSD) model was first studied by Goldman and Tucker (1956) for linear programming and simplified by Xu et al. (1996). Then a generalization of HSD was employed to solve SOCP problems by Andersen et al. (2003) and Sturm (2002). As its name suggests, the HSD model and its dual are equivalent. Self-dual formulations embed the original problem (19) in a larger conic optimization problem such that the latter is primal and dual feasible, with known feasible points, and from which solution we can extract optimal solutions or certificates of infeasibility of the original problem.

We define the homogeneous and self-dual model for problem (19) as follows:

\begin{array}{l} A x - b τ = 0, \\ A^{T} y + z - c τ = 0, \\ - c^{T} x + b^{T} y - κ = 0, \\ (x; τ) \in \tilde{K}, (z; κ) \in \tilde{K} . \end{array}

(23)

Here

τ

and

κ

are two additional variables and we use the notation that

\tilde{K} ≔ \bar{K} \times ℝ_{+} .

The model (23) can be viewed as a self-dual optimization problem with a zero objective function. If

(\hat{x}, \hat{τ}, \hat{y}, \hat{z}, \hat{κ})

is any feasible solution to (23), then if

\hat{τ} > 0

, a primal-dual optimal solution to (19) and (20) is given by

(x^{*}, y^{*}, z^{*}) = (\hat{x}, \hat{y}, \hat{z}) / \hat{τ},

and the duality gap is given by

c^{T} x^{*} - b^{T} y^{*} = \hat{κ} / \hat{τ} = 0

. The homogeneous algorithm is an application of the primal-dual method for the computation of a feasible solution to (23). In order to achieve this, we follow the guideline of path-following interior point method and define a central path that is a smooth curve connecting an initial interior point and a complementary solution. So the set of nonlinear equations

\begin{matrix} A x - b τ & = & γ (A x^{0} - b τ^{0}), \\ A^{T} y + z - c τ & = & γ (A^{T} y^{0} + z^{0} - c τ^{0}), \\ - c^{T} x + b^{T} y - κ & = & γ (- c^{T} x^{0} + b^{T} y^{0} - κ^{0}), \\ x \circ z & = & γ μ^{0} e, \\ τ κ & = & γ μ^{0}, \end{matrix}

(24)

defines the central path of the homogeneous model parameterized by

γ \in [0, 1]

(x^{0}, z^{0}, y^{0}, τ^{0}, κ^{0})

is an initial feasible point and

μ

has the expression

μ ≔ \frac{x^{T} s + τ κ}{r + 1}

where

r

is the number of cones.

11.2 The Nesterov–Todd Search Direction

The Newton search direction is only guaranteed to be well-defined in a narrow neighbourhood around the central path. The search direction corresponds to applying Newton's method to (24) in a scaled space and then scaling the resulting search direction back to the original space so that it is well-defined. A matrix

W

is a scaling matrix if it satisfies the conditions

W ≻ 0

and

W Q W = Q

where

W ≻ 0

means

W

is symmetric and positive definite and

Q

is a symmetric block diagonal matrix composed by so called reflection matrices

Q_{i}

with the following definition:

\begin{array}{l} Q_{i} ≔ (\begin{matrix} 1 & 0 & \dots & 0 \\ 0 & - 1 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & - 1 \end{matrix}) for quadratic cone, \\ Q_{i} ≔ (\begin{matrix} 0 & 1 & 0 & \dots & 0 \\ 1 & 0 & 0 & \dots & 0 \\ 0 & 0 & - 1 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & 0 & \dots & - 1 \end{matrix}) for rotated quadratic cone. \end{array}

It is easy to see that if we scale

x

W x

z

W^{- 1} z

A

A W^{- 1}

, and

c

W^{- 1} c

, the resulting primal and dual pair is equivalent to (19) and (20), see Alizadeh and Goldfarb (2003) for more details.

An important issue is the choice of the scaling matrix

W

. According to Andersen et al. (2003), the best results are obtained for the Nesterov–Todd (NT) scaling suggested by Nesterov and Todd (1997). In the NT scaling,

W

is chosen such that

W x = \bar{x} = \bar{z} = W^{- 1} z .

Then the resulting Newton system to be solved to get direction

(Δ x, Δ τ, Δ y, Δ z, Δ κ)

\begin{matrix} A Δ x - b Δ τ & = & (γ - 1) (A x^{0} - b τ^{0}), \\ A^{T} Δ y + Δ z - c Δ τ & = & (γ - 1) (A^{T} y^{0} + z^{0} - c τ^{0}), \\ - c^{T} Δ x + b^{T} Δ y - Δ κ & = & (γ - 1) (- c^{T} x^{0} + b^{T} y^{0} - κ^{0}), \\ W x^{0} \circ W^{- 1} Δ z + W^{- 1} z^{0} \circ W Δ x & = & - W x^{0} \circ W^{- 1} z^{0} + γ μ^{0} e, \\ τ^{0} Δ κ + κ^{0} Δ τ & = & - τ^{0} κ^{0} + γ μ^{0} . \end{matrix}

(25)

11.3 Mehrotra's Predictor-Corrector Method

When Newton's method is applied to the perturbed complementarity conditions in (24), the quadratic terms are neglected. Instead of neglecting the quadratic term Mehrotra (1992) suggested using a second-order correction of the search direction which increases the efficiency of the algorithm significantly in practice (Lustig et al. (1992)).

To implement this idea, we first solve (24) for

γ = 0

to get an affine scaling direction and a maximum step size

α_{n}^{max}

to the boundary. Then use these directions to estimate the quadratic terms

W Δ x \circ W^{- 1} Δ z and Δ τ Δ κ

from (24) and use

α_{n}^{max}

to choose

γ = \min (δ, {(1 - α_{n}^{max})}^{2}) (1 - α_{n}^{max}),

where

δ \in (0, 1)

is a constant. Therefore we can choose

γ

dynamically depending on how much progress can be made in the pure Newton (affine scaling) direction.

11.4 Solving the KKT System

The solution of the Newton system of equations (25) is the most computationally costly operation. To reduce the system, we need the following definition. Associated with each vector

x = (x_{0}; \bar{x}) \in ℝ^{n}

there is an arrow-shaped matrix

Arw (x)

defined as:

Arw (x) ≔ (\begin{matrix} x_{0} & {\bar{x}}^{T} \\ \bar{x} & x_{0} I \end{matrix}),

where

I

is the identity matrix of dimension

n - 1

. Together with the definition in (22), it is not hard to see that

x \circ z = Arw (x) z .

In practice, system (25) is reduced to the augmented system by eliminating

Δ z

and

Δ κ

from the system as follows:

(\begin{matrix} - W^{2} & A^{T} \\ A & 0 \end{matrix}) (\begin{matrix} g_{1} \\ g_{2} \end{matrix}) = (\begin{matrix} r_{2} - W {(Arw (W x^{0}))}^{- 1} r_{4} \\ r_{1} \end{matrix})

(26)

and

(\begin{matrix} - W^{2} & A^{T} \\ A & 0 \end{matrix}) (\begin{matrix} h_{1} \\ h_{2} \end{matrix}) = (\begin{matrix} c \\ b \end{matrix})'

(27)

where

r_{1}, \dots, r_{4}

(

r_{5}

eliminated) are the corresponding right-hand side in (25) and we have that

Δ τ = \frac{r_{3} - c^{T} g_{1} + b^{T} g_{2}}{{(τ^{0})}^{-^{1}} κ^{0} + c^{T} h_{1} - b^{T} h_{2}}

and

(\begin{matrix} Δ x \\ Δ y \end{matrix}) = (\begin{matrix} g_{1} \\ g_{2} \end{matrix}) + (\begin{matrix} h_{1} \\ h_{2} \end{matrix}) Δ τ .

Linear systems (26) and (27) are systems of

m + n

variables, symmetric and indefinite. Submatrix

W

is block diagonal and positive definite. Note that systems (26) and (27) have the same coefficient matrix so we only need to perform factorization once per iteration.

The system (27) can be further reduced by eliminating

g_{1}

and

h_{1}

, to a positive definite system usually called normal equations defined as

(A W^{- 2} A^{T}) h_{2} = b + A W^{- 2} c,

(28)

also system (26) can be reduced similarly.

Typically, formulation (28) is preferred for many problems as the system matrix can be factorized by a sparse Cholesky. However, this brings some well-known disadvantages: ill-conditioning of the system is often observed during the final stages of the algorithm. If matrix

A

contains dense columns (columns with relatively many nonzeros), then

A W^{- 2} A^{T}

has many nonzeros, which in turn makes the factorization expensive. On the other hand, solving the augmented system by Bunch–Parlett type factorization is usually slower, but it normally avoids the fill-in caused by dense columns.

e04ptf can detect and handle dense columns in the KKT system effectively. Since matrix

W^{- 2}

in (28) is block diagonal, so dense columns also come as a linear combination of some columns in

A

. Depending on the number and the density of the ‘dense’ columns, the solver may either choose to directly use an augmented system formulation or to treat these columns separately in a product-form Cholesky factorization as described by Goldfarb and Scheinberg (2005). It is also possible to manually override the automatic choice via the optional parameter SOCP System Formulation and let the solver use a normal equations or an augmented system formulation.

Badly scaled optimal solutions may present numerical challenges, therefore iterative refinement is employed for reducing the roundoff errors produced during the solution of the system. When the condition number of the system

A W^{- 2} A^{T}

prevents the satisfactory use of iterative refinement, e04ptf switches automatically to an augmented system formulation, reporting RS (Restart) in the last column of the iteration log (I). Furthermore, e04ptf provides several scaling techniques to adjust the numerical characteristics of the problem data, see optional parameter SOCP Scaling.

Finally, factorization of the system matrix can degrade sparsity, so the resulting fill-in can be large, therefore several ordering techniques are included to minimize it. e04ptf uses Harwell packages MA97 (see Hogg and Scott (2011) and HSL (2011)) for the underlying sparse linear algebra factorization and MC68 approximate minimum degree algorithm, and METIS (Karypis and Kumar (1998)) nested dissection algorithm for the ordering.

11.5 Stopping Criteria

11.5.1 Convergence-optimal termination

To measure the infeasibility, the following measures

\begin{array}{l} ρ_{P} ≔ \frac{{‖A x - b τ‖}_{\infty}}{\max (1, {‖A, b‖}_{\infty})}, & relative primal feasibility, \\ ρ_{D} ≔ \frac{{‖A^{T} y + z - c τ‖}_{\infty}}{\max (1, {‖A^{T}, I, - c‖}_{\infty})}, & relative dual feasibility, \\ ρ_{G} ≔ \frac{|- c^{T} x + b^{T} y - κ|}{\max (1, {‖- c^{T}, b^{T}, 1‖}_{\infty})}, & relative duality gap \end{array}

are defined to measure the relative reduction in the primal, dual and gap infeasibility, respectively. In addition, an extra measure is considered to quantify the accuracy in the objective function, which is given by

ρ_{A} ≔ \frac{|c^{T} x - b^{T} y|}{τ + |b^{T} y|} .

The iteration is considered nearly feasible and optimal, and the interior point algorithm is stopped when the following conditions

\max (ρ_{P}, ρ_{D}) \leq ε_{1} and ρ_{A} \leq ε_{2}

are satisfied. Here

ε_{1}

and

ε_{2}

may be set using SOCP Stop Tolerance and SOCP Stop Tolerance 2, respectively.

Premature termination is triggered and the returned solution is considered as an optimal solution if the current iteration exhibits fast convergence and the optimality measures lie within a small range of desired precision. In particular, the self-dual algorithm is stopped if the above termination conditions are met within a small factor and

τ > 1000 κ

. This measure is tracked after the first

10

iterations.

In addition, the solver stops prematurely and reports suboptimal solution when it predicts that the current estimate of the solution will not be improved in subsequent iterations. In most cases the returned solution should be acceptable.

11.5.2 Infeasibility/Unboundedness Detection

The problem is concluded to be primal or dual infeasible if one of the following conditions hold:

1. $\max (ρ_{P}, ρ_{D}, ρ_{G}) \leq ε_{1} and τ \leq ε_{2} \max (1, κ)$ .
2. $μ \leq ε_{2} μ_{0} and τ \leq ε_{2} \max (1, κ)$ .

Then the problem is declared dual infeasible if

c^{T} x < 0

or primal infeasible otherwise.

11.6 Further Details

e04ptf includes an advance preprocessing phase (called presolve) to reduce the dimensions of the problem before passing it to the solver. The reduction in problem size generally improves the behaviour of the solver, shortening the total computation time. In addition, infeasibility may also be detected during preprocessing. The default behaviour of the presolve can be modified by optional parameter SOCP Presolve.

12 Optional Parameters

Several optional parameters in e04ptf define choices in the problem specification or the algorithm logic. In order to reduce the number of formal arguments of e04ptf these optional parameters have associated default values that are appropriate for most problems. Therefore, you need only specify those optional parameters whose values are to be different from their default values.

The remainder of this section can be skipped if you wish to use the default values for all optional parameters.

The optional parameters can be changed by calling e04zmf anytime between the initialization of the handle and the call to the solver. Modification of the optional parameters during intermediate monitoring stops is not allowed. Once the solver finishes, the optional parameters can be altered again for the next solve.

The option values can be retrieved by calling e04znf.

The following is a list of the optional parameters available. A full description of each optional parameter is provided in Section 12.1.

Defaults
Infinite Bound Size
Monitoring File
Monitoring Level
Print File
Print Level
Print Options
Print Solution
SOCP Iteration Limit
SOCP Monitor Frequency
SOCP Presolve
SOCP Scaling
SOCP Stop Tolerance
SOCP Stop Tolerance 2
SOCP System Formulation
Stats Time
Task
Time Limit

12.1 Description of the Optional Parameters

For each option, we give a summary line, a description of the optional parameter and details of constraints.

The summary line contains:

the keywords;
a parameter value, where the letters $a$ , $i$ and $r$ denote options that take character, integer and real values respectively;
the default value, where the symbol $ε$ is a generic notation for machine precision (see x02ajf).

All options accept the value

DEFAULT

to return single options to their default states.

Keywords and character values are case and white space insensitive.

Defaults

This special keyword may be used to reset all optional parameters to their default values. Any argument value given with this keyword will be ignored.

Infinite Bound Size

r

Default

= 10^{20}

This defines the ‘infinite’ bound

bigbnd

in the definition of the problem constraints. Any upper bound greater than or equal to

bigbnd

will be regarded as

+ \infty

(and similarly any lower bound less than or equal to

- bigbnd

will be regarded as

- \infty

). Note that a modification of this optional parameter does not influence constraints which have already been defined; only the constraints formulated after the change will be affected.

Constraint:

Infinite Bound Size \geq 1000

Monitoring File

i

Default

= - 1

i \geq 0

, the unit number for the secondary (monitoring) output. If set to

- 1

, no secondary output is provided. The following information is output to the unit:

–a listing of the optional parameters;
–problem statistics, the iteration log, and the final status as set by Monitoring Level;
–the solution if set by Print Solution.

Constraint:

Monitoring File \geq - 1

Monitoring Level

i

Default

= 4

This parameter sets the amount of information detail that will be printed by the solver to the secondary output. The meaning of the levels is the same as with Print Level.

Constraint:

0 \leq Monitoring Level \leq 5

Print File

i

Default

= advisory message unit number

i \geq 0

, the unit number for the primary output of the solver. If

Print File = - 1

, the primary output is completely turned off independently of other settings. The default value is the advisory message unit number as defined by x04abf at the time of the optional parameters initialization, e.g., at the initialization of the handle. The following information is output to the unit:

–a listing of optional parameters if set by Print Options;
–problem statistics, the iteration log, and the final status from the solver as set by Print Level;
–the solution if set by Print Solution.

Constraint:

Print File \geq - 1

Print Level

i

Default

= 2

This parameter defines how detailed information should be printed by the solver to the primary output.

$i$	Output
$0$	No output from the solver
$1$	Only the final status and the primal and dual objective value
$2$	Problem statistics, one line per iteration showing the progress of the solution with respect to the convergence measures, final status and statistics
$3$	As level $2$ but each iteration line is longer, including step lengths and errors
$4, 5$	As level $3$ but further details of each iteration are presented

Constraint:

0 \leq Print Level \leq 5

Print Options

a

Default

= YES

Print Options = YES

, a listing of optional parameters will be printed to the primary output.

Constraint:

Print Options = YES

NO

Print Solution

a

Default

= NO

Print Solution = X

, the final values of the primal variables are printed on the primary and secondary outputs.

Print Solution = YES

ALL

, in addition to the primal variables, the final values of the dual variables are printed on the primary and secondary outputs.

Constraint:

Print Solution = YES

NO

X

ALL

SOCP Iteration Limit

i

Default

= 100

The maximum number of iterations to be performed by e04ptf. Setting the option too low might lead to

ifail = 22

Constraint:

SOCP Iteration Limit \geq 1

SOCP Monitor Frequency

i

Default

= 0

This parameter defines the frequency of how often subroutine monit is called. If

i > 0

, the solver calls monit at the end of every

i

th iteration. If it is set to

0

, the subroutine is not called at all.

Constraint:

SOCP Monitor Frequency \geq 0

SOCP Presolve

a

Default

= FULL

This parameter allows you to reduce the level of presolving of the problem or turn it off completely. If the presolver is turned off, the solver will try to handle the problem as given by you. In such a case, the presence of fixed variables or linear dependencies in the constraint matrix can cause numerical instabilities to occur. In normal circumstances, it is recommended to use the full presolve which is the default.

Constraint:

SOCP Presolve = FULL

BASIC

NO

SOCP Scaling

a

Default

= NONE

This parameter controls the type of scaling to be applied on the constraint matrix

A

before solving the problem. More precisely, the scaling procedure will try to find diagonal matrices

D_{1}

and

D_{2}

such that the values in

D_{1} A D_{2}

are of a similar order of magnitude. The solver is less likely to run into numerical difficulties when the constraint matrix is well scaled.

Constraint:

SOCP Scaling = ARITHMETIC

GEOMETRIC

NONE

SOCP Stop Tolerance

r

Default

= \sqrt{ε}

This parameter sets the value

ε_{1}

which is the tolerance for the convergence measures in the stopping criteria, see Section 11.5.

Constraint:

SOCP Stop Tolerance > ε

SOCP Stop Tolerance 2

r

Default

= \sqrt{ε}

This parameter sets the additional tolerance

ε_{2}

used in the stopping criteria, see Section 11.5.

Constraint:

SOCP Stop Tolerance 2 > ε

SOCP System Formulation

a

Default

= AUTO

As described in Section 11.4, e04ptf can internally work either with the normal equations formulation (28) or with the augmented system (26) and (27). A brief discussion of advantages and disadvantages is presented in (27). Setting the option value to

AUTO

leaves the decision to the solver based on the structure of the constraints and it is the recommended value. This will typically lead to the normal equations formulation unless there are many dense columns or the system is significantly cheaper to factorize as the augmented system. Note that in some cases even if

SOCP System Formulation = NORMAL EQUATIONS

the solver might switch the formulation through the computation to the augmented system due to numerical instabilities or computational cost.

Constraint:

SOCP System Formulation = AUTO

AUGMENTED SYSTEM

AS

NORMAL EQUATIONS

NE

Stats Time

a

Default

= NO

This parameter allows you to turn on timings of various parts of the algorithm to give a better overview of where most of the time is spent. This might be helpful for a choice of different solving approaches. It is possible to choose between CPU and wall clock time. Choice

YES

is equivalent to

WALL CLOCK

Constraint:

Stats Time = YES

NO

CPU

WALL CLOCK

Task

a

Default

= MINIMIZE

This parameter specifies the required direction of the optimization. If

Task = FEASIBLE POINT

, the objective function (if set) is ignored and the algorithm stops as soon as a feasible point is found with respect to the given tolerance. If no objective function is set, Task reverts to

FEASIBLE POINT

automatically.

Constraint:

Task = MINIMIZE

MAXIMIZE

FEASIBLE POINT

Time Limit

r

Default

= 10^{6}

A limit to the number of seconds that the solver can use to solve one problem. If during the convergence check this limit is exceeded, the solver will terminate with

ifail = 23

Constraint:

Time Limit > 0

NAG Library Manual, Mark 27

Interfaces: FL CL CPP AD

NAG FL Interface Introduction

E04 (Opt) Chapter Contents

E04 (Opt) Chapter Introduction

e04pt: FL CL CPP AD

NAG FL Interfacee04ptf (handle_​solve_​socp_​ipm)

▸▿ Contents

1 Purpose

2 Specification

3 Description

3.1 Structure of the Lagrangian Multipliers

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Formulating Problems as SOCPs

9.1.1 Quadratically Constrained Quadratic Programming

9.1.2 Norm Minimization Problems

9.1.3 Robust Linear Programming

9.2 Description of the Printed Output

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

11 Algorithmic Details

11.1 Homogeneous Self-Dual Algorithm

11.2 The Nesterov–Todd Search Direction

11.3 Mehrotra's Predictor-Corrector Method

11.4 Solving the KKT System

11.5 Stopping Criteria

11.5.1 Convergence-optimal termination

11.5.2 Infeasibility/Unboundedness Detection

11.6 Further Details

12 Optional Parameters

12.1 Description of the Optional Parameters

NAG FL Interface
e04ptf (handle_solve_socp_ipm)