Package com.imsl.math

Class SparseLP

java.lang.Object
com.imsl.math.SparseLP
All Implemented Interfaces:
Serializable, Cloneable
public class SparseLP extends Object implements Serializable, Cloneable
Solves a sparse linear programming problem by an infeasible primal-dual interior-point method.

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

$$\begin{array}{cll} \min{c^Tx} & \mbox{subject to} & b_l\le Ax\le b_u\\ {x\in {R}^n} && \\ & & x_l \le x \le x_u\\ \end{array}$$

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, SparseLP 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{array}{llrlrr} \min{c^Tx} & \mbox{subject to} & Ax&= b & & \\ & & x_i+s_i&=u_i, & x_i, s_i \ge 0, & i\in I_u \\ & & x_j&\ge 0, & & j\in I_s \end{array}$$

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.

In order to simplify the description it is assumed in the following that the problem above contains only variables with upper bounds, i.e. is of the form

$$\begin{array}{cllc} (P) &\min{c^Tx}& \mbox{subject to} & Ax= b\\ &&& x+s=u,\\ &&& x,s\ge 0 \end{array}$$

The corresponding dual problem is

$$\begin{array}{cllr} (D)&\max{b^Ty-u^Tw} & \mbox{subject to} & A^Ty+z-w=c,\\ &&& z,w\ge 0 \end{array}$$

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

$$\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}$$

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

Class SparseLP, 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 good 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). Class SparseLP 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 the weighted centering components, 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{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)$$

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

Class SparseLP now uses a line search 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{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}$$

In addition to the weighted Mehrotra predictor-corrector, SparseLP 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 Gondzio (1994).

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

$$\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)$$

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 Altman and Gondzio (1998). 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 class SparseLP, both approaches are implemented. The AS approach is chosen if A contains dense columns, if there are 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.

Class SparseLP stops with optimal termination status if the current iterate satisfies the following three conditions:

$$\frac{\mu}{1+0.5\left(\left|c^Tx\right|+\left|b^Ty-u^Tw\right|\right)} \le \mbox{optimalityTolerance}$$

$$\frac{\parallel \left(b-Ax,x+s-u\right)\parallel}{1+\parallel (b,u) \parallel} \le \mbox{primalTolerance}$$

and

$$\frac{\parallel c-A^Ty-z+w \parallel }{1+\parallel c \parallel} \le \mbox{dualTolerance}$$

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

Class SparseLP is based on the code HOPDM developed by Jacek Gondzio et al., see the HOPDM User's Guide (1995).

See Also:

  • Constructor Details

    • SparseLP

      public SparseLP(SparseMatrix a, double[] b, double[] c)
      Constructs a SparseLP object.
      Parameters:
      a - a SparseMatrix object containing the location and value of each nonzero coefficient in the constraint matrix A. If there is no constraint matrix, set a = null.
      b - a double array of length m, the number of constraints, containing the right-hand side of the constraints. If there are limits on both sides of the constraints, then b contains the lower limit of the constraints.
      c - a double array of length n, the number of variables, containing the coefficients of the objective function

    • SparseLP

      public SparseLP(MPSReader mps)
      Constructs a SparseLP object using an MPSReader object.
      Parameters:
      mps - an MPSReader object specifying the Linear Programming problem

    • SparseLP

      public SparseLP(int[] colPtr, int[] rowInd, double[] values, double[] b, double[] c)
      Constructs a SparseLP object using Compressed Sparse Column (CSC), or Harwell-Boeing format. See Compressed Sparse Column (CSC) Format, Chapter 1.
      Parameters:
      colPtr - an int array containing the location in values in which to place the first nonzero value of each succeeding column of the constraint matrix A. colPtr, rowInd and values specify the location and value of each nonzero coefficient in the constraint matrix A in CSC format.
      rowInd - an int array containing a list of the row indices of each column of the constraint matrix A. colPtr, rowInd and values specify the location and value of each nonzero coefficient in the constraint matrix A in CSC format.
      values - a double array containing the value of each nonzero coefficient in the constraint matrix A. colPtr, rowInd and values specify the location and value of each nonzero coefficient in the constraint matrix A in CSC format.
      b - a double array of length m, number of constraints, containing the right-hand side of the constraints; if there are limits on both sides of the constraints, then b contains the lower limit of the constraints
      c - a double array of length n, the number of variables, containing the coefficients of the objective function

  • Method Details

    • solve

      public void solve() throws SparseLP.DiagonalWeightMatrixException, SparseLP.CholeskyFactorizationAccuracyException, SparseLP.PrimalUnboundedException, SparseLP.PrimalInfeasibleException, SparseLP.DualInfeasibleException, SparseLP.InitialSolutionInfeasibleException, SparseLP.TooManyIterationsException, SparseLP.ProblemUnboundedException, SparseLP.ZeroColumnException, SparseLP.ZeroRowException, SparseLP.IncorrectlyEliminatedException, SparseLP.IncorrectlyActiveException, SparseLP.IllegalBoundsException
      Solves the sparse linear programming problem by an infeasible primal-dual interior-point method.
      Throws:
      SparseLP.DiagonalWeightMatrixException - is thrown if a diagonal element of the diagonal weight matrix is too small
      SparseLP.CholeskyFactorizationAccuracyException - is thrown if the Cholesky factorization failed because of accuracy problems
      SparseLP.PrimalUnboundedException - is thrown if the primal problem is unbounded
      SparseLP.PrimalInfeasibleException - is thrown if the primal problem is infeasible
      SparseLP.DualInfeasibleException - is thrown if the dual problem is infeasible
      SparseLP.InitialSolutionInfeasibleException - is thrown if the initial solution for the one-row linear program is infeasible
      SparseLP.TooManyIterationsException - is thrown if the maximum number of iterations has been exceeded
      SparseLP.ProblemUnboundedException - is thrown if the problem is unbounded
      SparseLP.ZeroColumnException - is thrown if a column of the constraint matrix has no entries
      SparseLP.ZeroRowException - is thrown if a row of the constraint matrix has no entries
      SparseLP.IncorrectlyEliminatedException - is thrown if one or more LP variables are falsely characterized by the internal presolver
      SparseLP.IncorrectlyActiveException - is thrown if one or more LP variables are falsely characterized by the internal presolver
      SparseLP.IllegalBoundsException - is thrown if the lower bound is greater than the upper bound

    • setUpperBound

      public void setUpperBound(double[] upperBound)
      Sets the upper bound on the variables.
      Parameters:
      upperBound - a double array of length n containing the upper bound \(x_u\) on the variables. If there is no upper bound on a variable, then 1.0e30 should be set as the upper bound.

      Default: None of the variables has an upper bound

    • getUpperBound

      public double[] getUpperBound()
      Returns the upper bound on the variables.
      Returns:
      a double array containing the upper bound on the variables

    • setLowerBound

      public void setLowerBound(double[] lowerBound)
      Sets the lower bound on the variables.
      Parameters:
      lowerBound - a double array of length n containing the lower bound \(x_l\) on the variables. If there is no lower bound on a variable, then -1.0e30 should be set as the lower bound.

      Default: lowerBound[i] = 0

    • getLowerBound

      public double[] getLowerBound()
      Returns the lower bound on the variables.
      Returns:
      a double array containing the lower bound on the variables

    • setConstraintType

      public void setConstraintType(int[] constraintType)
      Sets the types of general constraints in the matrix A.
      Parameters:
      constraintType - an int array of length m containing the types of general constraints in the matrix A. Let \(r_i = a_{i1}x_1 + \ldots + a_{in}x_n\). Then, the value of constraintType[i] signifies the following:

      constraintType

      Constraint
      0 \({\rm {r}}_i = {\rm {b}}_i\)
      1 \({\rm {r}}_i \le {\rm {bu}}_i\)
      2 \({\rm {r}}_i \ge {\rm {b}}_i\)
      3 \({\rm {b}}_i \le {\rm {r}}_i \le {\rm {bu}}_i\)
      4 Ignore this constraint

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

      Default: constraintType[i] = 0

    • getConstraintType

      public int[] getConstraintType()
      Returns the types of general constraints in the matrix A. See setConstraintType.
      Returns:
      an int array containing the types of general constraints in the matrix A

    • getOptimalValue

      public double getOptimalValue()
      Returns the optimal value of the objective function.
      Returns:
      a double scalar containing the optimal value of the objective function

    • getSolution

      public double[] getSolution()
      Returns the solution x of the linear programming problem.
      Returns:
      a double array containing the solution x of the linear programming problem

    • setConstant

      public void setConstant(double c0)
      Sets the value of the constant term in the objective function.
      Parameters:
      c0 - a double scalar containing the value of the constant term in the objective function

      Default: c0 = 0

    • getConstant

      public double getConstant()
      Returns the value of the constant term in the objective function.
      Returns:
      a double scalar containing the value of the constant term in the objective function

    • setMaxIterations

      public void setMaxIterations(int maxIterations)
      Sets the maximum number of iterations allowed for the primal-dual solver.
      Parameters:
      maxIterations - an int scalar containing the maximum number of iterations allowed for the primal-dual solver

      Default: maxIterations = 200

    • getMaxIterations

      public int getMaxIterations()
      Returns the maximum number of iterations allowed for the primal-dual solver.
      Returns:
      an int scalar containing the maximum number of iterations allowed for the primal-dual solver

    • setRelativeOptimalityTolerance

      public void setRelativeOptimalityTolerance(double optimalityTolerance)
      Sets the relative optimality tolerance.
      Parameters:
      optimalityTolerance - a double scalar containing the relative optimality tolerance

      Default: optimalityTolerance = 1.0e-10

    • getRelativeOptimalityTolerance

      public double getRelativeOptimalityTolerance()
      Returns the relative optimality tolerance.
      Returns:
      a double scalar containing the relative optimality tolerance

    • setPrimalInfeasibilityTolerance

      public void setPrimalInfeasibilityTolerance(double primalTolerance)
      Sets the primal infeasibility tolerance.
      Parameters:
      primalTolerance - a double scalar containing the primal infeasibility tolerance

      Default: primalTolerance = 1.0e-8

    • getPrimalInfeasibilityTolerance

      public double getPrimalInfeasibilityTolerance()
      Returns the primal infeasibility tolerance.
      Returns:
      a double scalar containing the primal infeasibility tolerance

    • setDualInfeasibilityTolerance

      public void setDualInfeasibilityTolerance(double dualTolerance)
      Sets the dual infeasibility tolerance.
      Parameters:
      dualTolerance - a double scalar containing the dual infeasibility tolerance

      Default: dualTolerance = 1.0e-8

    • getDualInfeasibilityTolerance

      public double getDualInfeasibilityTolerance()
      Returns the dual infeasibility tolerance.
      Returns:
      a double scalar containing the dual infeasibility tolerance

    • setPreordering

      public void setPreordering(int preorder)
      Sets the variant of the Minimum Degree Ordering (MDO) algorithm used in the preordering of the normal equations or augmented system matrix.
      Parameters:
      preorder - an int scalar containing the variant of the Minimum Degree Ordering (MDO) algorithm used in the preordering of the normal equations or augmented system matrix
      preorderMethod
      0 A variant of the MDO algorithm using pivotal cliques.
      1 A variant of George and Liu's Quotient Minimum Degree algorithm.

      Default: preorder = 0

    • getPreordering

      public int getPreordering()
      Returns the variant of the Minimum Degree Ordering (MDO) algorithm used in the preordering of the normal equations or augmented system matrix. See setPreordering.
      Returns:
      an int scalar containing the variant of the Minimum Degree Ordering (MDO) algorithm used in the preordering of the normal equations or augmented system matrix

    • setPrintLevel

      public void setPrintLevel(int printLevel)
      Sets the print level.
      Parameters:
      printLevel - an int containing the print level
      printLevel Action
      0 No printing.
      1 Prints statistics on the LP problem and the solution process.

      Default: printLevel = 0

    • getPrintLevel

      public int getPrintLevel()
      Returns the print level. See setPrintLevel.
      Returns:
      an int scalar containing the print level

    • setPresolve

      public void setPresolve(int presolve)
      Sets the presolve option.
      Parameters:
      presolve - an int containing the the presolve option to resolve 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 In addition to 1, eliminate redundant (and forcing) rows.
      3 In addition to 1 and 2, eliminate dominated variables.
      4 In addition to 1, 2, and 3, eliminate singleton columns.
      5 In addition to 1, 2, 3, and 4, eliminate doubleton rows.
      6 In addition to 1, 2, 3, 4, and 5, eliminate aggregate columns.

      Default: presolve = 0

    • getPresolve

      public int getPresolve()
      Returns the presolve option. See setPresolve.
      Returns:
      an int scalar containing the presolve option

    • getIterationCount

      public int getIterationCount()
      Returns the number of iterations used by the primal-dual solver.
      Returns:
      an int scalar containing the number of iterations used by the primal-dual solver

    • getTerminationStatus

      public int getTerminationStatus()
      Returns the termination status for the problem.
      Returns:
      an int scalar containing the termination status for the problem

      status

      Description
      0 Optimal solution found.
      1 The problem is primal infeasible (or dual unbounded).
      2 The problem is primal unbounded (or dual infeasible).
      3 Suboptimal solution found (accuracy problems).
      4 Iterations limit maxIterations exceeded.
      5 An error outside of the solution phase of the algorithm, e.g. a user input or a memory allocation error.

    • getDualSolution

      public double[] getDualSolution()
      Returns the dual solution.
      Returns:
      a double array containing the dual solution

    • getSmallestCPRatio

      public double getSmallestCPRatio()
      Returns the ratio of the smallest complementarity product to the average.
      Returns:
      a double scalar containing the ratio of the smallest complementarity product to the average

    • getLargestCPRatio

      public double getLargestCPRatio()
      Returns the ratio of the largest complementarity product to the average.
      Returns:
      a double scalar containing the ratio of the largest complementarity product to the average

    • getDualInfeasibility

      public double getDualInfeasibility()
      Returns the dual infeasibility of the solution.
      Returns:
      a double scalar containing the dual infeasibility of the solution, \(\parallel c - A^Ty - z + w \parallel\)

    • setUpperLimit

      public void setUpperLimit(double[] bu)
      Sets the upper limit of the constraints that have both a lower and an upper bound.
      Parameters:
      bu - a double array of length m containing the upper limit \(b_u\) of the constraints that have both a lower and an upper bound. If such a constraint exists, then method setConstraintType must be used to define the type of the constraints. If constraintType[i] != 3, i.e. if constraint i is not two-sided, then the corresponding entry in bu, bu[i], is ignored.

      Default: None of the constraints has an upper limit

    • getUpperLimit

      public double[] getUpperLimit()
      Returns the upper limit of the constraints that have both a lower and an upper bound.
      Returns:
      a double array containing the upper limit of the constraints that have both a lower and an upper bound. Returns null if the upper limit has not been set.

    • getPrimalInfeasibility

      public double getPrimalInfeasibility()
      Returns the primal infeasibility of the solution.
      Returns:
      a double scalar containing the primal infeasibility of the solution, \(\parallel x + s - u \parallel\)

    • getViolation

      public double getViolation()
      Returns the violation of the variable bounds.
      Returns:
      a double scalar containing the violation of the variable bounds, \(\parallel b - Ax \parallel\)