imsl.optimize.sparse_lp

sparse_lp(a, b, c, xb=None, obj_constant=0.0, preorder=0, presolve=0, max_iter=200, opt_tol=1e-10, primal_infeas_tol=1e-08, dual_infeas_tol=1e-08)

Solve a sparse linear programming problem.

The linear programming problem is of the form

Minimize c^T * x

Subject to:

b_l <= A * x <= b_u,

x_l <= x <= x_u

Parameters:
  • a (sparse matrix) – The constraint matrix of shape (M,N) in SciPy’s COO or CSC format. If no constraint matrix is required, either set a to None or to a matrix of shape (0,N).
  • b (tuple) –

    A three-element tuple (bl, bu, constr_type) containing the constraint bounds and information about the constraint type. If a constraint is two-sided then bl contains the lower bound and bu the upper bound of the constraint. If a constraint is one-sided, then bl must contain the bound and the corresponding entry in bu is ignored.

    bl can be defined in one of two forms:

    • A scalar : The same bound bl will be applied to all constraints.
    • (M,) array_like : Bound bl[i] will be applied to the i -th constraint.

    bu can be defined in one of two forms:

    • A scalar : The same upper bound bu will be applied to all constraints.
    • (M,) array_like : Upper bound bu[i] will be applied to the i -th constraint.

    constr_type can be defined in one of two forms:

    • A scalar : All constraints are of type constr_type.
    • (M,) array_like : The i -th constraint is of type constr_type[i].

    Let \(r_i = a_{i1}x_1+\ldots+a_{in}x_n\) be the i -th constraint, with lower bound \(bl_i\) and upper bound \(bu_i\). Then, the value of constr_type[i] signifies the following:

    constr_type[i] Constraint
    0 \(r_i = bl_i\)
    1 \(r_i \le bl_i\)
    2 \(r_i \ge bl_i\)
    3 \(bl_i \le r_i \le bu_i\)
    4 Ignore this constraint

    Note that constr_type[i] = 3 should only be used for constraints i with both finite lower and finite upper bound. For one-sided constraints, use constr_type[i] = 1 or constr_type[i] = 2. For free constraints, use constr_type[i] = 4.

    If no constraint matrix a is available, then b is ignored.

  • c ((N,) array_like) – Array containing the coefficients of the objective function.
  • xb (tuple, optional) –

    A two-element tuple (xlb, xub) containing the lower and upper bounds on the variables.

    xlb can be defined in one of two forms:

    • A scalar : The same lower bound xlb will be applied to all variables.
    • (N,) array_like : Each variable \(x_i\) will be lower bounded by xlb[i].

    Use constant UNBOUNDED_BELOW from the constants module to specify negative infinite bounds.

    xub can be defined in one of two forms:
    • A scalar : The same upper bound xub will be applied to all variables.
    • (N,) array_like : Each variable \(x_i\) will be upper bounded by xub[i].

    Use constant UNBOUNDED_ABOVE from the constants module to specify positive infinite bounds.

    Default: All variables are non-negative.

  • obj_constant (float, optional) – Value of the constant term in the objective function.
  • preorder (int, optional) –

    The variant of the Minimum Degree Ordering (MDO) algorithm used in the preordering of the normal equations or augmented system matrix:

    preorder Method
    0 A variant of the MDO algorithm using pivotal cliques.
    1 A variant of George and Liu’s Quotient Minimum Degree algorithm.
  • presolve (int, optional) –

    Presolve the LP problem in order to reduce the problem size or to detect infeasibility or unboundedness of the problem. Depending on the number of presolve techniques used, different presolve levels can be chosen:

    presolve Description
    0 No presolving
    1 Eliminate singleton rows
    2 Additionally to 1, eliminate redundant (and forcing) rows
    3 Additionally to 2, eliminate dominated variables.
    4 Additionally to 3, eliminate singleton columns.
    5 Additionally to 4, eliminate doubleton rows.
    6 Additionally to 5, eliminate aggregate columns.
  • max_iter (int, optional) – The maximum number of iterations allowed for the primal-dual solver.
  • opt_tol (float, optional) – The relative optimality tolerance.
  • primal_infeas_tol (float, optional) – The primal infeasibility tolerance.
  • dual_infeas_tol (float, optional) – The dual infeasibility tolerance.
Returns:

  • A named tuple with the following fields
  • sol ((N,) ndarry) – Array containing the primal solution of the LP problem.
  • obj (float) – Optimal value of the objective function.
  • dual ((M,) ndarray) – Array containing the dual solution.
  • n_iter (int) – The number of iterations required by the primal-dual solver.
  • infeas (tuple) – A 3-element tuple containing the primal and dual infeasibilities. Element infeas[0] contains the primal infeasibility of the solution, infeas[1] the violation of the variable bounds, and infeas[2] the dual infeasibility of the solution.
  • cp_ratios (tuple) – A 2-element tuple containing the ratios of the smallest (in cp_ratios[0]) and largest (in cp_ratios[1]) complementarity product to the average.

Notes

The function sparse_lp uses an infeasible primal-dual interior-point method to solve linear programming problems, i.e., problems of the form

\[\begin{split}\min_{x \in {R}^n} c^Tx \quad \text{subject to} \quad b_l \le Ax \le b_u\\ x_l \le x \le x_u\end{split}\]

where c is the objective coefficient vector, A is the coefficient matrix, and the vectors \(b_l, b_u, x_l\) and \(x_u\) are the lower and upper bounds on the constraints and the variables, respectively.

Internally, sparse_lp transforms the problem given by the user into a simpler form that is computationally more tractable. After redefining the notation, the new form reads

\[\begin{split}\begin{array}{lrll} \min{c^Tx}&&&\\ \text{subject to} & Ax&= b&\\ & x_i+s_i&=u_i, \quad x_i, s_i \ge 0, &i \in I_u\\ & x_j&\ge 0, &j \in I_s \end{array}\end{split}\]

Here, \(I_u \cup I_s = \{1,\ldots,n\}\) is a partition of the index set \(\{1,\ldots,n\}\) into upper bounded and standard variables.

To simplify the description, the following assumes that the problem above contains only variables with upper bounds, i.e. is of the form

\[\begin{split}\begin{array}{clrl} &\min{c^Tx}&&\\ (P)&\text{subject to} & Ax&= b\\ && x+s&=u,\\ && x,s&\ge 0 \end{array}\end{split}\]

The corresponding dual problem is

\[\begin{split}\begin{array}{clrl} &\max{b^Ty-u^Tw}&&\\ (D)&\text{subject to} & A^Ty+z-w&=c,\\ && z,w&\ge 0 \end{array}\end{split}\]

The Karush-Kuhn-Tucker (KKT) optimality conditions for (P) and (D) are

\[\begin{split}\begin{array}{rlr} Ax&=b,&\quad (1.1)\\ x+s&=u,&\quad (1.2)\\ A^Ty+z-w&=c,&\quad (1.3)\\ XZe&=0,&\quad (1.4)\\ SWe&=0,&\quad (1.5)\\ x,z,s,w&\ge 0,&\quad (1.6) \end{array}\end{split}\]

where \(X=diag(x)\), \(Z=diag(z)\), \(S=diag(s)\) and \(W=diag(w)\) are diagonal matrices, and \(e=(1,\ldots,1)^T\) is a vector of ones.

Function sparse_lp, like all infeasible interior point methods, generates a sequence

\[(x_k,s_k,y_k,z_k,w_k), \quad k=0,1,\ldots\]

of iterates that satisfy \((x_k,s_k,y_k,z_k,w_k)>0\) for all k, but are in general not feasible, i.e. the linear constraints (1.1)-(1.3) are only satisfied in the limiting case \(k \to \infty\).

The barrier parameter \(\mu\), defined by

\[\mu = \frac{x^Tz+s^Tw}{2n}\]

measures how well the complementarity conditions (1.4), (1.5) are satisfied.

Mehrotra’s predictor-corrector algorithm is a variant of Newton’s method applied to the KKT conditions (1.1)-(1.5). Function sparse_lp uses a modified version of this algorithm to compute the iterates \((x_k,s_k,y_k,z_k,w_k)\). In every step of the algorithm, the search direction vector

\[\Delta := (\Delta x, \Delta s, \Delta y, \Delta z, \Delta w)\]

is decomposed into two parts, \(\Delta = \Delta_a + \Delta_c^{\omega}\), where \(\Delta_a\) and \(\Delta_c^{\omega}\) denote the affine-scaling and a weighted centering component, respectively. Here,

\[\Delta_c^{\omega}:=(\omega_P\Delta x_c, \omega_P\Delta s_c, \omega_D \Delta y_c, \omega_D \Delta z_c, \omega_D \Delta w_c)\]

where \(\omega_P\) and \(\omega_D\) denote the primal and dual corrector weights, respectively.

The vectors \(\Delta_a\) and \(\Delta_c := (\Delta x_c, \Delta s_c, \Delta y_c, \Delta z_c, \Delta w_c)\) are determined by solving the linear system

\[\begin{split}\begin{bmatrix} A & 0 & 0 & 0 & 0\\ I & 0 & I & 0 & 0\\ 0 & A^T & 0 & I & -I\\ Z & 0 & 0 & X & 0\\ 0 & 0 & W & 0 & S \end{bmatrix} \begin{bmatrix} \Delta x\\ \Delta y\\ \Delta s\\ \Delta z\\ \Delta w \end{bmatrix} = \begin{bmatrix} r_b\\ r_u\\ r_c\\ r_{xz}\\ r_{ws} \end{bmatrix} \quad (2)\end{split}\]

for two different right-hand sides.

For \(\Delta_a\), the right-hand side is defined as

\[(r_b,r_u,r_c,r_{xz},r_{ws})=(b-Ax,u-x-s,c-A^Ty-z+w,-XZe,-WSe).\]

Here, \(r_b\) and \(r_u\) are the violations of the primal constraints and \(r_c\) defines the violations of the dual constraints.

The resulting direction \(\Delta_a\) is the pure Newton step applied to the system (1.1)-(1.5).

In order to obtain the corrector direction \(\Delta_c\), the maximum stepsizes \(\alpha_{Pa}\) in the primal and \(\alpha_{Da}\) in the dual space preserving nonnegativity of \((x,s)\) and \((z,w)\) respectively, are determined, and the predicted complementarity gap

\[g_a = (x+\alpha_{Pa}\Delta x_a)^T(z+\alpha_{Da}\Delta z_a)+ (s+\alpha_{Pa}\Delta s_a)^T(w+\alpha_{Da}\Delta w_a)\]

is computed. It is then used to determine the barrier parameter

\[\hat{\mu} = \left( \frac{g_a}{g} \right)^2 \frac{g_a}{2n},\]

where \(g=x^Tz+s^Tw\) denotes the current complementarity gap.

The direction \(\Delta_c\) is then computed by choosing

\[(r_b,r_u,r_c,r_{xz},r_{ws})=(0,0,0,\hat{\mu}e- \Delta X_a \Delta Z_a e,\hat{\mu} e-\Delta W_a \Delta S_ae)\]

as the right-hand side in the linear system (2).

Function sparse_lp now uses a linesearch to find the optimal weight \(\hat{\omega}=\left( \hat{\omega_P}, \hat{\omega_D} \right)\) that maximizes the stepsizes \(\left( \alpha_P, \alpha_D \right)\) in the primal and dual directions of \(\Delta = \Delta_a + \Delta_c^{\omega}\), respectively.

A new iterate is then computed using a step reduction factor \(\alpha_0 = 0.99995\):

\[\begin{split}\begin{array}{cl} (x_{k+1},s_{k+1},y_{k+1},z_{k+1},w_{k+1}) \quad = & (x_k,s_k,y_k,z_k,w_k)+\\ & \quad \alpha_0 \left( \alpha_P \Delta x, \alpha_P \Delta s, \alpha_D \Delta y, \alpha_D \Delta z, \alpha_D \Delta w \right) \end{array}\end{split}\]

In addition to the weighted Mehrotra predictor-corrector, sparse_lp also uses multiple centrality correctors to enlarge the primal-dual stepsizes per iteration step and to reduce the overall number of iterations required to solve an LP problem. The maximum number of centrality corrections depends on the ratio of the factorization and solve efforts for system (2) and is therefore problem dependent. For a detailed description of multiple centrality correctors, refer to [1].

The linear system (2) can be reduced to more compact forms, the augmented system (AS)

\[\begin{split}\begin{bmatrix} -\Theta^{-1} & A^T\\ A & 0 \end{bmatrix} \begin{bmatrix} \Delta x\\ \Delta y \end{bmatrix} = \begin{bmatrix} r\\ h \end{bmatrix} \quad (3)\end{split}\]

or further by elimination of \(\Delta x\) to the normal equations (NE) system

\[A \Theta A^T \Delta y = A \Theta r + h, \quad (4)\]

where

\[\Theta = \left( X^{-1}Z+S^{-1}W \right)^{-1}, r = r_c-X^{-1}r_{xz} +S^{-1}r_{ws}-S^{-1}Wr_u, h=r_b.\]

The matrix on the left-hand side of (3), which is symmetric indefinite, can be transformed into a symmetric quasidefinite matrix by regularization. Since these types of matrices allow for a Cholesky-like factorization, the resulting linear system can be solved easily for \((\Delta x, \Delta y)\) by triangular substitutions. For more information on the regularization technique, see [3]. For the NE system, matrix \(A \Theta A^T\) is positive definite, and therefore a sparse Cholesky algorithm can be used to factor \(A \Theta A^T\) and solve the system for \(\Delta y\) by triangular substitutions with the Cholesky factor L.

In function sparse_lp, both approaches are implemented. The AS approach is chosen if A contains dense columns, if there is a considerable number of columns in A that are much denser than the remaining ones or if there are many more rows than columns in the structural part of A. Otherwise, the NE approach is selected.

Function sparse_lp stops with optimal termination status if the current iterate satisfies the following three conditions:

\[\begin{split}\begin{array}{rl} \frac{\mu}{1+0.5(|c^Tx|+|b^Ty-u^Tw|)} & \le \texttt{opt_tol},\\[1ex] \frac{\|(b-Ax,x+s-u)\|}{1+\|(b,u)\|} & \le \texttt{primal_infeas_tol}, \\[1ex] \frac{\|c-A^Ty-z+w\|}{1+\|c\|} & \le \texttt{dual_infeas_tol}, \end{array}\end{split}\]

where primal_infeas_tol, dual_infeas_tol and opt_tol are primal infeasibility, dual infeasibility and optimality tolerances, respectively. The default value is 1.0e-10 for opt_tol and 1.0e-8 for the other two tolerances.

Function sparse_lp is based on the code HOPDM developed by Jacek Gondzio et al., see [2].

References

[1]Gondzio, Jacek (1994), Multiple Centrality Corrections in a Primal-Dual Method for Linear Programming, Logilab Technical Report 1994.20, Logilab, HEC Geneva, Section of Management Studies, Geneva.
[2]Gondzio, Jacek (1995), HOPDM - Modular Solver for LP Problems, User’s Guide to version 2.12, WP-95-50, International Institute for Applied Systems Analysis, Laxenburg, Austria.
[3]Altman, Anna, and Jacek Gondzio (1998), Regularized Symmetric Indefinite Systems in Interior Point Methods for Linear and Quadratic Optimization, Logilab Technical Report 1998.6, Logilab, HEC Geneva, Section of Management Studies, Geneva.

Examples

The linear programming problem

\[\begin{split}\min f(x) = 2x_1-8x_2+3x_3 \\ \text{subject to} \quad x_1 + 3x_2 \le 3 \\ 2x_2+3x_3 \le 6 \\ x_1+x_2+x_3 \ge 2 \\ -1 \le x_1 \le 5\\ 0 \le x_2 \le 7\\ 0 \le x_3 \le 9\end{split}\]

is solved.

>>> import numpy as np
>>> from scipy.sparse import coo_matrix
>>> from imsl.optimize import sparse_lp
>>> row = np.array([0, 0, 1, 1, 2, 2, 2])
>>> col = np.array([0, 1, 1, 2, 0, 1, 2])
>>> val = np.array([1.0, 3.0, 2.0, 3.0, 1.0, 1.0, 1.0])
>>> bl = np.array([3.0, 6.0, 2.0])
>>> bu = None
>>> constr_type = np.array([1, 1, 2])
>>> b = (bl, bu, constr_type)
>>> c = np.array([2.0, -8.0, 3.0])
>>> xlb = np.array([-1.0, 0.0, 0.0])
>>> xub = np.array([5.0, 7.0, 9.0])
>>> xb = (xlb, xub)
>>> # Create matrix in SCS (COO) format
>>> a = coo_matrix((val, (row, col)))
>>> lp_result = sparse_lp(a, b, c, xb)
>>> np.set_printoptions(precision=3)
>>> print("Solution vector:\n")
... 
Solution vector:
>>> print(str(lp_result.sol)) 
[-0.375  1.125  1.25 ]
>>> print("\nOptimal objective value: {0:10.4f}".format
...       (lp_result.obj)) 
Optimal objective value:    -6.0000