naginterfaces.library.opt.handle_​solve_​socp_​ipm

naginterfaces.library.opt.handle_solve_socp_ipm(handle, x=None, u=None, uc=None, monit=None, data=None, io_manager=None)

handle_solve_socp_ipm 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).

Note: this function uses optional algorithmic parameters, see also: handle_opt_set(), handle_opt_get().

For full information please refer to the NAG Library document for e04pt

https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e04/e04ptf.html

Parameters

handleHandle

The handle to the problem. It needs to be initialized (e.g., by handle_init()) and to hold a problem formulation compatible with handle_solve_socp_ipm. It must not be changed between calls to the NAG optimization modelling suite.

xNone or float, array-like, shape $\left(\text{nvar}\right)$, optional

The input of $x$ is reserved for future releases of the NAG Library and it is ignored at the moment.

uNone or float, array-like, shape $\left(\text{nnzu}\right)$, optional

Note: if $\text{nnzu}>0$, $u$ holds Lagrange multipliers (dual variables) for the bound constraints and linear constraints. If $\text{nnzu}=0$, $u$ will not be referenced.

The input of $u$ is reserved for future releases of the NAG Library and it is ignored at the moment.

ucNone or float, array-like, shape $\left(\text{nnzuc}\right)$, optional

Note: if $\text{nnzuc}>0$, $uc$ holds Lagrange multipliers (dual variables) for second-order cones as defined by handle_set_group(). If $\text{nnzuc}=0$, $uc$ will not be referenced.

The input of $uc$ is reserved for future releases of the NAG Library and it is ignored at the moment.

monitNone or callable monit(handle, rinfo, stats, data=None), optional

Note: if this argument is None then a NAG-supplied facility will be used.

$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 option ‘SOCP Monitor Frequency’ (the default is $0$, $monit$ is not called).

Parameters

handleHandle

The handle to the problem as provided on entry to handle_solve_socp_ipm. It may be used to query the model during the solve, and extract the current approximation of the solution by handle_set_get_real().

rinfofloat, ndarray, shape $\left(100\right)$

Error measures and various indicators at the end of the current iteration as described in $rinfo$.

statsfloat, ndarray, shape $\left(100\right)$

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.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

dataarbitrary, optional

User-communication data for callback functions.

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

xfloat, ndarray, shape $\left(\text{nvar}\right)$

The final values of the variables $x$.

ufloat, ndarray, shape $\left(\text{nnzu}\right)$

The final values of the variables $u$.

ucfloat, ndarray, shape $\left(\text{nnzuc}\right)$

The final values of the variables $\text{uc}$.

rinfofloat, ndarray, shape $\left(100\right)$

Error measures and various indicators of the algorithm (see Algorithmic Details for details) as given in the table below:

$0$	Value of the primal objective.
$1$	Value of the dual objective.
$2$	Flag indicating the system formulation used by the solver, $0$: augmented system, $1$: normal equation.
$3$	Factorization type, $3$: Cholesky, $4$: Bunch–Parlett.
$4$–$13$	Not referenced in this solver.
$14$	Relative primal infeasibility, see Convergence-optimal termination.
$15$	Relative duality gap, see Convergence-optimal termination.
$16$	Relative dual infeasibility, see Convergence-optimal termination.
$17$	Accuracy, see Convergence-optimal termination.
$18$	$\tau$, see [equation].
$19$	$\kappa$, see [equation].
$20$	Step length.
$21$–$99$	Reserved for future use.

statsfloat, ndarray, shape $\left(100\right)$

Solver statistics as given in the table below. Note that times are measured in seconds, see option ‘Stats Time’.

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

Other Parameters

‘Defaults’valueless

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

‘Infinite Bound Size’float

Default $={10}^{20}$

This defines the ‘infinite’ bound $\text{bigbnd}$ in the definition of the problem constraints. Any upper bound greater than or equal to $\text{bigbnd}$ will be regarded as $+\infty$ (and similarly any lower bound less than or equal to $-\text{bigbnd}$ will be regarded as $-\infty$). Note that a modification of this option does not influence constraints which have already been defined; only the constraints formulated after the change will be affected.

Constraint: $\text{‘Infinite Bound Size'}\ge 1000$.

‘Monitoring File’int

Default $=-1$

If $i\ge 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 options;

• problem statistics, the iteration log, and the final status as set by ‘Monitoring Level’;

• the solution if set by ‘Print Solution’.

Constraint: $\text{‘Monitoring File'}\ge -1$.

‘Monitoring Level’int

Default $=4$

This argument 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\le \text{‘Monitoring Level'}\le 5$.

‘Print File’int

Default $=\text{advisory message unit number}$

If $i\ge 0$, the unit number for the primary output of the solver. If $\text{‘Print File'}=-1$, the primary output is completely turned off independently of other settings. The default value is the advisory message unit number at the time of the options initialization, e.g., at the initialization of the handle. The following information is output to the unit:

• a listing of options 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: $\text{‘Print File'}\ge -1$.

‘Print Level’int

Default $=2$

This argument 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\le \text{‘Print Level'}\le 5$.

‘Print Options’str

Default $=\text{'YES'}$

If $\text{‘Print Options'}=\text{'YES'}$, a listing of options will be printed to the primary output.

Constraint: $\text{‘Print Options'}=\text{'YES'}$ or $\text{'NO'}$.

‘Print Solution’str

Default $=\text{'NO'}$

If $\text{‘Print Solution'}=\text{'X'}$, the final values of the primal variables are printed on the primary and secondary outputs.

If $\text{‘Print Solution'}=\text{'YES'}$ or $\text{'ALL'}$, in addition to the primal variables, the final values of the dual variables are printed on the primary and secondary outputs.

Constraint: $\text{‘Print Solution'}=\text{'YES'}$, $\text{'NO'}$, $\text{'X'}$ or $\text{'ALL'}$.

‘SOCP Factorization Method’str

Default $=\text{'MA86'}$

If the value of ‘SOCP System Formulation’ is ‘AUGMENTED SYSTEM’, then this parameter controls whether Harwell packages ‘MA86’ or ‘MA97’ is used for the sparse linear algebra factorization. Note that if the option value ‘SOCP System Formulation’ is set to ‘AUTO’ or ‘NORMAL EQUATIONS’, then specifying ‘MA86’ with this option will allow the solver to use this package in the case that the solver switches to the augmented system formulation.

Constraint: $\text{‘SOCP Factorization Method'}=\text{'MA97'}$ or $\text{'MA86'}$.

‘SOCP Iteration Limit’int

Default $=100$

The maximum number of iterations to be performed by handle_solve_socp_ipm. Setting the option too low might lead to $errno$ = 22.

Constraint: $\text{‘SOCP Iteration Limit'}\ge 1$.

‘SOCP Monitor Frequency’int

Default $=0$

This argument defines the frequency of how often function $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 function is not called at all.

Constraint: $\text{‘SOCP Monitor Frequency'}\ge 0$.

‘SOCP Presolve’str

Default $=\text{'FULL'}$

This argument 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: $\text{‘SOCP Presolve'}=\text{'FULL'}$, $\text{'BASIC'}$ or $\text{'NO'}$.

‘SOCP Scaling’str

Default $=\text{'NONE'}$

This argument 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: $\text{‘SOCP Scaling'}=\text{'ARITHMETIC'}$, $\text{'GEOMETRIC'}$ or $\text{'NONE'}$.

‘SOCP Stop Tolerance’float

Default $=\sqrt{ϵ}$

This argument sets the value ${ϵ}_1$ which is the tolerance for the convergence measures in the stopping criteria, see Stopping Criteria.

Constraint: $\text{‘SOCP Stop Tolerance'}>ϵ$.

‘SOCP Stop Tolerance 2’float

Default $=\sqrt{ϵ}$

This argument sets the additional tolerance ${ϵ}_2$ used in the stopping criteria, see Stopping Criteria.

Constraint: $\text{‘SOCP Stop Tolerance 2'}>ϵ$.

‘SOCP System Formulation’str

Default $=\text{'AUGMENTED SYSTEM'}$

As described in Solving the KKT System, handle_solve_socp_ipm can internally work either with the normal equations formulation [equation] or with the augmented system [equation] and [equation]. A brief discussion of advantages and disadvantages is presented in [equation]. Setting the option value to ‘AUTO’ leaves the decision to the solver based on the structure of the constraints. 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 $\text{‘SOCP System Formulation'}=\text{'NORMAL EQUATIONS'}$ the solver might switch the formulation through the computation to the augmented system due to numerical instabilities or computational cost.

Constraint: $\text{‘SOCP System Formulation'}=\text{'AUTO'}$, $\text{'AUGMENTED SYSTEM'}$, $\text{'AS'}$, $\text{'NORMAL EQUATIONS'}$ or $\text{'NE'}$.

‘Stats Time’str

Default $=\text{'NO'}$

This argument 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: $\text{‘Stats Time'}=\text{'YES'}$, $\text{'NO'}$, $\text{'CPU'}$ or $\text{'WALL CLOCK'}$.

‘Task’str

Default $=\text{'MINIMIZE'}$

This argument specifies the required direction of the optimization. If $\text{‘Task'}=\text{'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: $\text{‘Task'}=\text{'MINIMIZE'}$, $\text{'MAXIMIZE'}$ or $\text{'FEASIBLE POINT'}$.

‘Time Limit’float

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 $errno$ = 23.

Constraint: $\text{‘Time Limit'}>0$.

Raises

NagValueError

(errno $1$)

$handle$ has not been initialized.

(errno $1$)

$handle$ does not belong to the NAG optimization modelling suite, has not been initialized properly or is corrupted.

(errno $1$)

$handle$ has not been initialized properly or is corrupted.

(errno $2$)

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

(errno $2$)

The problem is already being solved.

(errno $4$)

On entry, $\text{nvar}=⟨value⟩$, expected $value=⟨value⟩$.

Constraint: $\text{nvar}$ must match the current number of variables of the model in the $handle$.

(errno $5$)

On entry, $\text{nnzu}=⟨value⟩$.

$\text{nnzu}$ does not match the size of the Lagrangian multipliers for constraints.

The correct value is either $0$ or $⟨value⟩$.

(errno $5$)

On entry, $\text{nnzu}=⟨value⟩$.

$\text{nnzu}$ does not match the size of the Lagrangian multipliers for constraints.

The correct value is $0$ for no constraints.

(errno $5$)

On entry, $\text{nnzuc}=⟨value⟩$.

$\text{nnzuc}$ does not match the size of the Lagrangian multipliers for second-order cone constraints.

The correct value is either $0$ or $⟨value⟩$.

(errno $5$)

On entry, $\text{nnzuc}=⟨value⟩$.

$\text{nnzuc}$ does not match the size of the Lagrangian multipliers for second-order cone constraints.

$\text{nnzuc}=0$ when there are no second-order cone constraints.

(errno $51$)

The problem was found to be primal infeasible.

(errno $52$)

The problem was found to be dual infeasible.

Warns

NagAlgorithmicWarning

(errno $50$)

Suboptimal solution.

NagAlgorithmicMajorWarning

(errno $22$)

Maximum number of iterations exceeded.

(errno $23$)

The solver terminated after the maximum time allowed was exceeded.

(errno $24$)

No progress, stopping early.

NagCallbackTerminateWarning

(errno $20$)

User requested termination during a monitoring step.

Notes

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

$$\begin{array}{c}\begin{array}{cc}{minimize}_{x\in R^n}& c^{T}x\\ \text{subject to}& l_A\le Ax\le u_A\text{,}\\ & l_x\le x\le u_x\text{,}\\ & x\in K\text{,}\end{array}\end{array}$$

where $K=K^{n_1}\times \cdots \times K^{n_r}\times R^{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\times 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 R^{n_i}\text{:}{z}_1^2\ge \sum _{j=2}^{n_i}{z}_j^2\text{,}{z}_1\ge 0\}\text{.}$$

• Rotated quadratic cone:

  $$K_r^{n_i}:=\{z=(z_1,z_2,\dots ,z_{n_i})\in R^{n_i}\text{:}2z_1{z}_2\ge \sum _{j=3}^{n_i}{z}_j^2\text{,}{z}_1\ge 0\text{,}{z}_2\ge 0\}\text{.}$$

handle_solve_socp_ipm 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 functions in the NAG optimization modelling suite. First, the problem handle is initialized by calling handle_init(). Then some of the functions handle_set_group(), handle_set_linobj(), handle_set_qconstr(), handle_set_qconstr_fac(), handle_set_quadobj(), handle_set_simplebounds() or handle_set_linconstr() may be called to formulate the quadratic cones, linear objective function, quadratic objective function, quadratic constraints, bounds of the variables, and the block of linear constraints, respectively. Alternatively, the whole model can be loaded from a file by handle_read_file(). When the handle is no longer needed, handle_free() should be called to destroy it and deallocate the memory held within. See the E04 Introduction for more details about the NAG optimization modelling suite.

The solver method can be modified by various options (see Other Parameters) which can be set by handle_opt_set() and handle_opt_set_file() 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 options.

The option ‘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 Algorithmic Details and Other Parameters for further details.

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 $uc$ for quadratic cone constraints.

If the simple bounds have been defined (handle_set_simplebounds() was successfully called), the first $2n$ 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 $2m$ elements of $u$ belong to multipliers for the linear constraints (if handle_set_linconstr() 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.

If convex quadratic constraints have been defined successfully by handle_set_qconstr() or handle_set_qconstr_fac(), denote the number of such constraints as $\text{nq}$, then the following $2\text{nq}$ elements of $u$ belong to multipliers for the convex quadratic constraints. The organization is the same as linear 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 $uc$ are the corresponding Lagrangian multipliers for the variables in the quadratic cone constraints that have been defined by handle_set_group(). All multipliers are stored next to each other in array $uc$ in the same order as the cone constraints were defined by handle_set_group(). 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$.

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, 2010, An indefinite sparse direct solver for large problems on multicore machines, RAL Technical Report. RAL-TR-2010-011

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, 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

See also

naginterfaces.library.examples.opt.handle_solve_socp_ipm_ex.main()