# Systems Optimization Laboratory

## LUSOL: Sparse LU for Ax = b

• AUTHOR: Michael Saunders
• CONTRIBUTORS: Philip Gill, Walter Murray, Margaret Wright, Michael O'Sullivan, Kjell Eikland, Yin Zhang, Nick Henderson, Ding Ma
• CONTENTS: A sparse LU factorization for square and rectangular matrices A, with Bartels-Golub-Reid updates for column replacement and other rank-1 modifications. Typically used for a sequence of linear equations as in the simplex method:

```        Solve Ax = b and/or A'y = c
Replace a column of A
Repeat with different b, c
```
The matrix A may have any shape and rank. Rectangular LU factors may be used to form a sparse null-space matrix operator.

Special feature 1: Three sparse pivoting options in the Factor routine:

```        Threshold partial pivoting (TPP)
Threshold rook pivoting    (TRP)
Threshold complete pivoting (TCP)
```
All options choose row and column permutations as they go, balancing sparsity and stability according to different rules. TPP is normally most efficient for solving Ax = b. TRP and TCP are rank-revealing factorizations. In practice, TRP is an effective method for estimating rank(A). TCP tends to be too dense and expensive to be useful, although MINOS and SNOPT switch from TPP to TRP and even TCP if necessary in case of persistent numerical difficulty.

Special feature 2: Multiple update routines:

```        Add, delete, or replace a column of A
Add, delete, or replace a row    of A
Add a general (sparse) rank-1 matrix to A
```

Numerical stability: LUSOL maintains LU factors with row and column permutations P, Q such that A = LU with PLP' lower triangular (with unit diagonals) and PUQ upper triangular. The condition of L is controlled throughout by maintaining |Lij| <= factol (= 10 or 5 or 2 or 1.1, ...), so that U tends to reflect the condition of A. This is essential for subsequent Bartels-Golub-type updates (which are implemented in a manner similar to John Reid's LA05 and LA15 packages in the HSL library).

If a fresh factorization is thought of as A = LDU (with unit diagonals on PLP' and PUQ), then TRP and TCP control the condition of both L and U by maintaining |Lij| <= factol and |Uij| <= factol, so that D reflects the condition of A. This is why TRP and TCP have rank-revealing properties.

Proven applications: LUSOL is the basis factorization package (BFP) for MINOS, SQOPT, SNOPT, lp_solve, AMPL/PATH, GAMS/PATH.

Shortcomings:
Factor: No special handling of dense columns.
Solve: No special treatment of sparse right-hand sides.
Documentation: No user's manual. Primary documentation is in-line comments within the f77 source code (and the more recent f90 version).

• REFERENCES:

J. K. Reid (1982). A sparsity-exploiting variant of the Bartels-Golub decomposition for linear programming bases, Mathematical Programming 24, 55-69.

P. E. Gill, W. Murray, M. A. Saunders and M. H. Wright (1987). Maintaining LU factors of a general sparse matrix, Linear Algebra and its Applications 88/89, 239-270.

P. E. Gill, W. Murray and M. A. Saunders (2005). SNOPT: An SQP algorithm for large-scale constrained optimization, SIGEST article, SIAM Review 47(1), 99-131. (See sections 4 and 5.)