Solve Large Linear Optimization Problems

Yifan Hu
Wolfram Research, Inc.

The LinearProgramming function in Mathematica is unique in its ability to solve linear optimization problems of integer/rational coefficients and giving rational solutions. These algorithms however are not designed for solving large machine precision problems. In this talk we discuss new development in LinearProgramming, including a new interior point algorithm that is capable of solving large scale linear optimization problems with hundreds of thousands of variables and equations.

Outline of the talk

Motivation

Introduction to interior point method

Import data files (MPS format)

Numerical examples

DualLinearProgramming

Conclusions

\(\(\(\(\(Off[General::spell1]\) \); \(\(Off[
        General::spell]\) \); \(\(Off[
        ParametricPlot3D::"ppcom"]\) \);\)\(
\)\(
\)\(
\)\(
\)\(
\)\(
\)\(
\)\(
\)\(
\)\(
\)\(
\)
    \)\)

Motivation

The Linear Programming problem

LinearProgramming  (pre-5.0) implements Simplex/Revised Simplex algorithm using dense linear algebra.

Unique in solving  exact/extended precision problems

But slow for medium/large sized problem

For medium/large machine number problems, we could improve the performance with algorithms based on sparse linear algebra

Sparse Simplex or Interior Point algorithms.

Requires efficient sparse arithmetics and solvers, one of the consumer and driver of our sparse technology, detailed tomorrow.

Choose interior point -- the technology useful for other class of problems in the future

Introduction to interior point method

Simplex: exponential complexity in the worst case. Interior point: polynomial.
First proposed 1984 by Karmarkar (then aged 26, born in India), with substantial hype at the time

But substantially modified since

Wright found a formulation that is essentially the application of the classical barrier function method (1960's), which make interior point method much easier to understand/analyse

Actually there was a paper in the 1960s' which applied the barrier function method to Linear Programming. So interior point was not new!

It was the hype in 1984 which attracted sufficient interests to make it a serious alternative to Simplex: in 1994, 10 years after the Karmarkar's invention, over 1300 papers had been published on interior point methods

Introduction to interior point: barrier function method

Min x,  
    x  ≥   0

Plot[x, {x, 0, 10}, PlotRange -> {{0, 1}, {0, 1}}];

Min x  -  t Log[x],  
    x  ≥ 0

Plot[{x - 0.1  Log[x], x - 0.05  Log[x], x - 0.025  Log[x]}, {x, 
          0.000000001, 1}, PlotRange -> {{0, 1}, {0, 1}}];

Introduction to interior point: the ill-conditioning

A simple linear programming problem

The barrier function plot

t = 0.1;  epsilon = 1.*^-10;
surface = 
        Plot3D[x - t (Log[x]  +  Log[y]), {x, epsilon, 1}, {y, epsilon, 
            1.2}, Mesh -> False, 
          ViewPoint -> {1, 1, 1}];

The barrier function and the feasible region

alpha  =  0.9; 
    f[s_] :=  alpha {s, 1 - s, 
            s - t (Log[s]  +  Log[1 - s])} + (1 - alpha) {1, 1, 
            1};
line = 
        ParametricPlot3D[f[s], {s, epsilon, 1}, 
          ViewPoint -> {1, 1, 1}];

Show[{surface, line}]

The Jacobian or the Hessian of the barrier function

is very ill conditioned e.g., the Hessian is

Introduction to interior point: solving the barrier system

The problem

Barrier formulation

First order necessary condition:

Let z = t (X^(-1)) e, then

This is a set of 2m + n linear/nonlinear equations with constraints. Solving using Newton's method

and

Introduction to interior point method: numerical issues

ill-Conditioning

Jacobi matrix

Newton's method => solve

Example the following system, with solution (Subscript[x, 1])  =  0 and (Subscript[x, 2])  =  1

E.g., for the first problem in the netlib suite "25fv47.mps",

Solution: sparse symmetric solver, modified to overcome such singularity.

Interior point can only be expected to solve the problem to about tolerance = Sqrt[$MachineEpsilon]

Introduction to interior point method: numerical issues

dense columns

A (Z^(-1))X (A^T)is dense if A has a dense column

Solution: iterative refinement, conjugate gradient, Sherman-Morrison formula, ...

others issues

steplength: to maintain positive definiteness
the penalty parameters: Mehrotra's predictor-corrector scheme
choice of initial starting point
preprocessing

Import data files (MPS format)

MPS format is the most commonly used format for documenting linear programming problems

MPS originated from IBM Mathematical Programming System

A lot of collections, such as the netlib collection

Import["afiro.mps", "Text"]

problem = Import["afiro.mps"]

Apply[NMinimize, problem]

Import data files (MPS format)

For large problem, it is best to solve using LinearProgramming interface for efficiency.

Load the utility package

<< Optimization`MPSData`

We can import a problem as an MPS object, without convert to/from equation form. More efficient

Example: 60 MB of file

x = 
      Import["osa-60.mps", "MPSData"] // AbsoluteTiming

Not solving it here, but take about 2 minutes to solve on a 1.4 GHZ Pentium 4.

This now gives an MPSData object

afiro = Import["afiro.mps", "MPSData"]

Data suitable for LinearProgramming can be extracted using ToLinearProgrammingData method

{c, A, b, bound} = afiro@ToLinearProgrammingData[];

This solve it

xs = LinearProgramming[c, A, b, bound]

c . xs

Numerical examples: comparing with dense simplex

bandm = Import["bandm.mps", "MPSData"]

{c, A, b, bound} = ToLinearProgrammingData[bandm];

x = 
          LinearProgramming[c, A, b, bound, 
            Method -> RevisedSimplex]; // Timing

c . x // InputForm

xi = LinearProgramming[c, A, b, bound]; // Timing

c . xi // InputForm

Numerical examples: a larger example

Import as MPSData object

<< Optimization`MPSData`

(pds02 = Import["pds-02.mps", "MPSData"]) // 
      AbsoluteTiming

Convert to a form suitable for LinearProgramming

{c, A, b, bound} = ToLinearProgrammingData[pds02];

Solve

x = LinearProgramming[c, A, b, bound]; // Timing

c . x // InputForm

Numerical examples: infeasible/unbounded problem

The primal-dual interior point is design for feasible problem, for infeasible/unbound problem, it diverges badly.

We implemented a heuristics to catch infeasible/unbounded problems

LinearProgramming[{1. , 1}, {{1, 1}, {1, 1}}, {{1, -1}, {2, 1}}, 
      Method -> InteriorPoint]

Sometimes the heuristic can not tell if a problem is infeasible with certainty

forest6 = Import["forest6.mps", "MPSData"]

{c, A, b, bound} = ToLinearProgrammingData[forest6];

LinearProgramming[c, A, b, bound, 
        Method -> InteriorPoint];

LinearProgramming[c, A, b, bound];

DualLinearProgramming

This is a new functionality

Given general LP problem

its dual is:

Useful to know solutions for both for some applications

The relationship between primal and dual:

When both are feasible, the optimal value of (P) and (D) are the same, and the complementary conditions hold

DualInearProgramming gives the primal solution x, and the dual solution y,z,w, in the form

{x,y,z,w} = DualLinearProgramming[...]

DualLinearProgramming

Interior point process involved both the prime and dual variables, hence the dual solution comes naturally when the algorithm converges.

However interior point only kicks in after a preprocessing stage

The preprocessing can reduce the problem size by as much as 30%.

Each preprocessing move is recorded on a stack

During post processing the stack is unwounded to reverse these moves and recover the duals of the deleted variables/equations.

DualLinearProgramming: a simple example

{x, y, z, w} = 
      DualLinearProgramming[{3, 
          4}, {{1, 2}}, {5}, {{1, 4}, {1, 4}}]

{3, 4} . x

The dual of the constrain is y={2} ==> for one unit of increase in the RHS of constraint, we will get 2 unit of increase in the objective.

{x, y, z, w} = 
      DualLinearProgramming[{3, 
          4}, {{1, 2}}, {5 + 0.001}, {{1, 4}, {1, 4}}]

{3, 4} . x

Conclusions

Primal-dual interior point implemented

Large machine number linear programming can now be solved efficiently

Linear programming problem can now be imported in MPS format

A new function DualLinearProgramming is useful for parametric analysis of LP problem

Bibliography

I. I. Dikin, 1967, Iterative solution of problems of Linear and Quadratic Programming, Doklady Akademii Nauk SSSR 174, 747-748, Translated in: Soviet Mathematics Doklady 8, 674-675, 1967.