Linear Programming solvers: the state of the art Julian Hall School - PowerPoint PPT Presentation

Linear Programming solvers: the state of the art Julian Hall School of Mathematics, University of Edinburgh 4th ISM-ZIB-IMI MODAL Workshop Mathematical Optimization and Data Analysis Tokyo 27 March 2019

Overview LP background Serial simplex Interior point methods Solvers Parallel simplex For structured LP problems For general LP problems A novel method Julian Hall Linear Programming solvers: the state of the art 2 / 59

Solution of linear programming (LP) problems f = c T x minimize subject to A x = b x ≥ 0 Background Example Fundamental model in optimal decision-making Solution techniques ◦ Simplex method (1947) ◦ Interior point methods (1984) ◦ Novel methods Large problems have ◦ 10 3 –10 8 variables ◦ 10 3 –10 8 constraints STAIR: 356 rows, 467 columns and 3856 nonzeros Matrix A is (usually) sparse Julian Hall Linear Programming solvers: the state of the art 3 / 59

Solving LP problems: Necessary and sufficient conditions for optimality f = c T x minimize subject to A x = b x ≥ 0 Karush-Kuhn-Tucker (KKT) conditions x ∗ is an optimal solution ⇐ ⇒ there exist y ∗ and s ∗ such that x T s A x = b (1) x ≥ 0 (3) = 0 (5) A T y + s = ≥ c (2) s 0 (4) For the simplex algorithm (1–2 and 5) always hold Primal simplex algorithm: (3) holds and the algorithm seeks to satisfy (4) Dual simplex algorithm: (4) holds and the algorithm seeks to satisfy (3) For interior point methods (1–4) hold and the method seeks to satisfy (5) Julian Hall Linear Programming solvers: the state of the art 4 / 59

Solving LP problems: Characterizing the feasible region x 3 f = c T x minimize subject to A x = b x ≥ 0 x 2 A ∈ R m × n is full rank Solution of A x = b is a n − m dim. hyperplane in R n K Intersection with x ≥ 0 is the feasible region K A vertex of K has m basic components, i ∈ B given by A x = b n − m zero nonbasic components, j ∈ N where B ∪ N partitions { 1 , . . . , n } A solution of the LP occurs at a vertex of K x 1 Julian Hall Linear Programming solvers: the state of the art 5 / 59

Solving LP problems: Optimality conditions at a vertex f = c T x minimize subject to A x = b x ≥ 0 Karush-Kuhn-Tucker (KKT) conditions x ∗ is an optimal solution ⇐ ⇒ there exist y ∗ and s ∗ such that x T s A x = b x ≥ = 0 (1) 0 (3) (5) A T y + s = c s ≥ 0 (2) (4) � x B � � c B � � s B � � � Given B ∪ N , partition A as B N , x as , c as and s as x N c N s N If x N = 0 and x B = � b ≡ B − 1 b then B x B + N x N = b so A x = b (1) For (2) If y = B − T c B and s B = 0 then B T y + s B = c B � B T � � s B � � c B � c N ≡ c N − N T y then (2) holds If s N = � y + = N T s N c N Finally, x T s = x T B s B + x T N s N = 0 (5) Need � b ≥ 0 for (3) and � c N ≥ 0 for (4) Julian Hall Linear Programming solvers: the state of the art 6 / 59

Solving LP problems: Simplex and interior point methods Simplex method (1947) Interior point method (1984) Replace x ≥ 0 by log barrier Given B ∪ N so (1–2 and 5) hold Primal simplex method Solve n � Assume ˆ b ≥ 0 f = c T x + µ (3) maximize ln( x j ) c N ≥ 0 Force ˆ (4) j =1 subject to A x = b Dual simplex method KKT (5) changes: Assume ˆ c N ≥ 0 (4) Replace x T s = 0 by XS = µ e Force ˆ b ≥ 0 (3) X and S have x and s on diagonal Modify B ∪ N KKT (1–4) hold Combinatorial approach Satisfy (5) by forcing XS = µ e as Cost O (2 n ) iterations µ → 0 Practically: O ( m + n ) iterations Iterative approach Practically: O ( √ n ) iterations Julian Hall Linear Programming solvers: the state of the art 7 / 59

Simplex method

The simplex algorithm: Definition x 3 � � � b x 2 corresponding to B ∪ N At a feasible vertex x = 0 If � c N ≥ 0 then stop: the solution is optimal 1 K Scan � c j < 0 for q to leave N x 2 � − � � x+ α d a q a q = B − 1 N e q and d = Let � 3 e q Scan � b i / � a iq > 0 for α and p to leave B 4 Exchange p and q between B and N 5 x 1 Go to 1 6 Julian Hall Linear Programming solvers: the state of the art 9 / 59

Primal simplex algorithm: Choose a column RHS N Assume � b ≥ 0 Seek � c N ≥ 0 Scan � c i < 0 for q to leave N B c T � c q � N Julian Hall Linear Programming solvers: the state of the art 10 / 59

Primal simplex algorithm: Choose a row RHS N Assume � b ≥ 0 Seek � c N ≥ 0 Scan � c j < 0 for q to leave N � � a q b Scan � b i / � a iq > 0 for p to leave B B � � a pq b p Julian Hall Linear Programming solvers: the state of the art 11 / 59

Primal simplex algorithm: Update cost and RHS RHS N Assume � b ≥ 0 Seek � c N ≥ 0 Scan � c j < 0 for q to leave N � � a q b Scan � b i / � a iq > 0 for p to leave B B a T � � a pq � b p p Update: Exchange p and q between B and N Update � b := � α P = � c T b − α P � a q b p / � a pq � � c q N c T c T a T Update � N := � N + α D � α D = − � c q / � a pq p Julian Hall Linear Programming solvers: the state of the art 12 / 59

Primal simplex algorithm: Data required RHS N Assume � b ≥ 0 Seek � c N ≥ 0 Scan � c j < 0 for q to leave N � � a q b Scan � b i / � a iq > 0 for p to leave B B a T � � a pq � b p p Update: Exchange p and q between B and N Update � b := � α P = � c T b − α P � a q b p / � a pq � � c q N c T c T a T Update � N := � N + α D � α D = − � c q / � a pq p Data required a T p = e T p B − 1 N Pivotal row � a q = B − 1 a q Pivotal column � Julian Hall Linear Programming solvers: the state of the art 13 / 59

Primal simplex algorithm RHS N Assume � b ≥ 0 Seek � c N ≥ 0 Scan � c j < 0 for q to leave N � � a q b Scan � b i / � a iq > 0 for p to leave B B a T � � a pq � b p p Update: Exchange p and q between B and N Update � b := � α P = � c T b − α P � a q b p / � a pq � � c q N c T c T a T Update � N := � N + α D � α D = − � c q / � a pq p Data required Why does it work? a T � p = e T p B − 1 N Pivotal row � b p × � c q Objective improves by − each iteration a q = B − 1 a q a pq � Pivotal column � Julian Hall Linear Programming solvers: the state of the art 14 / 59

Simplex method: Computation Standard simplex method (SSM): Major computational component N − 1 Update of tableau: � N := � a T RHS N a q � � p a pq � where � N = B − 1 N � � B N b Hopelessly inefficient for sparse LP problems c T � Prohibitively expensive for large LP problems N Revised simplex method (RSM): Major computational components a T B T π p = e p p = π T Pivotal row via and � p N BTRAN PRICE Represent B − 1 Pivotal column via B � a q = a q FTRAN INVERT Update B − 1 exploiting ¯ B = B + ( a q − B e p ) e T UPDATE p Julian Hall Linear Programming solvers: the state of the art 15 / 59

Serial simplex: Hyper-sparsity

Serial simplex: Solve B x = r for sparse r Given B = LU , solve L y = r ; U x = y In revised simplex method, r is sparse: consequences? If B is irreducible then x is full If B is highly reducible then x can be sparse Phenomenon of hyper-sparsity Exploit it when forming x Exploit it when using x Julian Hall Linear Programming solvers: the state of the art 17 / 59

Serial simplex: Hyper-sparsity Inverse of a sparse matrix and solution of B x = r B − 1 has density of 58%, so B − 1 r is Optimal B for LP problem stair typically dense Julian Hall Linear Programming solvers: the state of the art 18 / 59

Serial simplex: Hyper-sparsity Inverse of a sparse matrix and solution of B x = r B − 1 has density of 0.52%, so B − 1 r Optimal B for LP problem pds-02 is typically sparse—when r is sparse Julian Hall Linear Programming solvers: the state of the art 19 / 59

Serial simplex: Hyper-sparsity Use solution of L x = b To illustrate the phenomenon of hyper-sparsity To demonstrate how to exploit hyper-sparsity Apply principles to other triangular solves in the simplex method Julian Hall Linear Programming solvers: the state of the art 20 / 59

Serial simplex: Hyper-sparsity Recall: Solve L x = b using function ftranL ( L , b , x ) When b is sparse r = b for all j ∈ { 1 , . . . , m } do Inefficient until r fills in for all i : L ij � = 0 do r i = r i − L ij r j x = r Julian Hall Linear Programming solvers: the state of the art 21 / 59

Serial simplex: Hyper-sparsity Better: Check r j for zero When x is sparse function ftranL ( L , b , x ) Few values of r j are nonzero r = b Check for zero dominates for all j ∈ { 1 , . . . , m } do if r j � = 0 then Requires more efficient identification for all i : L ij � = 0 do of set X of indices j such that r j � = 0 r i = r i − L ij r j x = r Gilbert and Peierls (1988) H and McKinnon (1998–2005) Julian Hall Linear Programming solvers: the state of the art 22 / 59