Module factor 

Chapter 5.3 — Sparse LU factorization of the basis: Lu

Purpose

Every simplex iteration solves two or three linear systems with the basis matrix B (Chapter 5.4). We factorize B once into sparse triangular factors and keep them across iterations, modified by the Forrest–Tomlin update (Chapter 5.5) until a fresh factorization is cheaper or safer. Simplex bases are special: most of their columns are logical or nearly triangular, so most pivots can be found without any arithmetic by peeling off singletons; only a small bump (the kernel) needs genuine Gaussian elimination, done with Markowitz pivoting under a threshold.

Notation

B ∈ R^{m×m} has rows i and positions k (Chapter 5.1). The factorization is a sequence of m steps; step s pivots on the entry (p_s, q_s) (row p_s, position q_s). With the permutation matrices P (row s of PB is row p_s of B) and Q (column s of BQ is column q_s of B) we compute

  P B Q = L U,

L unit lower triangular and U upper triangular. We store both in the original indices, never permuting values:

• L as one column eta per step with a nonzero subdiagonal: the pivot row p_s and the multipliers l_is for rows i pivoted later;
• U as its diagonal u_{p_s q_s} (indexed by the row p_s) and its off-diagonal entries u_ij (row i, position j, with step(i) < step(j)), held twice: by rows and by columns. FTRAN needs U by columns, BTRAN by rows, and the Forrest–Tomlin update must delete a column and a row of U — so both copies are dynamic (Chapter 5.2) and are kept identical.

The active submatrix S during the elimination has row counts r_i and column counts c_j. A column singleton is a column with c_j = 1, a row singleton a row with r_i = 1. The Markowitz merit of a candidate pivot is (r_i − 1)(c_j − 1), an upper bound on the fill it creates. With the threshold u (LuOptions::rel_pivot_tol ∈ (0, 1]), an entry is acceptable iff |s_ij| ≥ u · max_k |s_kj| and |s_ij| > τ_abs (LuOptions::abs_pivot_tol). Entries of magnitude at most ε_drop (LuOptions::drop_tol) are negligible.

Algorithm L (LU factorization of a simplex basis)

Given the columns of B through a crate::ColumnSource (in the simplex method: the basis heading and A).

L1. [Load.] Validate the source (Chapter 5.1, B1). Build S = B in a column file with values and a row file with the pattern only; record nnz(B) and max |b_ij|. Bucket columns and rows by count. L2. [Empty columns.] Every active column with c_j = 0 is dead: it can never be pivoted. Remove it (its position joins the deficiency list). L3. [Column singleton.] If some column j has c_j = 1, with entry (i, v): if |v| > τ_abs pivot on (i, j) (step L6) — the L column of this step is empty and row i is copied to U — else declare j dead. Go to L2. L4. [Row singleton.] Otherwise, if some row i with r_i = 1 has its only entry (i, j) above τ_abs, pivot on it (step L6) — the U row is empty, the multipliers s_kj / s_ij form the L column, and no entry of S is updated. Go to L2. Steps L3–L4 constitute the triangular phase: no arithmetic beyond one division per L entry. Singleton pivots are forced (the only entry of their row or column in any completion of the elimination), so no threshold is applied to them; they cannot create growth in U. L5. [Markowitz.] Otherwise (the bump), remove the empty columns as in L2; if no active column remains, go to L7. Run Algorithm M to obtain a pivot and a list of dead columns; remove the dead columns, pivot (step L6) and repeat L5. L5′. [Pre-ordered bump; only with Lu::set_pre_order(Some(PreOrder::Colamd)), task L1.3.] Instead of L5: compute the COLAMD order (Chapter 5.10, Algorithm C) of the columns of the bump on its pattern; for each column c in that order that is still active (after removing the empty columns as in L2): if γ_c = max_i |s_ic| ≤ τ_abs declare c dead; otherwise pivot (step L6) on the acceptable entry (i, c) of least row count r_i (ties: the larger |s_ic|, then the first in the column) — threshold partial pivoting inside a fixed column order, as UMFPACK and SuperLU do. With pre_order = None step L5 runs exactly as before, so the factors are bit for bit those of the second edition. L6. [Pivot on (p, c).] Append the step (p, c) and the diagonal u_pc = s_pc. For every other (i, s_ic) of column c record l_ic = s_ic / s_pc; for every other (p, s_pj) of row p record u_pj = s_pj. Remove row p and column c from S. Then, for each j of the U row, find the positions in column j of the rows of the L column (one pass over column j against a map of those rows) and set s_ij ← s_ij − l_ic u_pj for every i of the L column, in its order, appending fill-in to both files; entries that cancel to |s_ij| ≤ ε_drop are dropped, in the order of column j. Only the entries just updated or filled can cancel — plus, the first time column j is updated, entries of B that were negligible from the start, so such a column is scanned whole once. Update the counts of the touched rows and columns. L7. [Finish.] If any column is dead, the basis is singular: report its positions and the rows that never received a pivot (LuError::Singular). Otherwise build the row and column copies of U, the inverse step maps, and the statistics, including the growth factor max |u_ij| / max |b_ij| and max |u_ij| itself, the reference of the update policy (Chapter 5.5). Record the pivot sequence (p_s, q_s) as the SymbolicLu of Chapter 5.12 (into the arrays of the previous one: 2m copies, no allocation once warm). Step L7 reads only the output of the elimination (Elim), so Algorithm R of Chapter 5.12 shares it.

Algorithm M (Markowitz search with threshold)

Let best be the acceptable entry of least merit found so far (ties go to the larger ratio |s_ij| / max_k |s_kj|), and e the number of lines examined that contained an acceptable entry.

M1. For k = 1, 2, …, m: M2. [Columns of count k.] For each active column j with c_j = k: compute γ_j = max_i |s_ij|; if γ_j ≤ τ_abs, mark j dead and go on. Otherwise evaluate every acceptable entry with merit (r_i − 1)(k − 1). Every unexamined candidate now lies in a row with r_i ≥ k (rows of smaller count were examined whole) and a column with c_j ≥ k, so its merit is at least (k − 1)²: stop if best reaches that bound, or if e ≥ search_limit and best exists. M3. [Rows of count k.] The same for each active row i with r_i = k, scanning its columns j (merit (k − 1)(c_j − 1)); the bound for unexamined candidates is now k (k − 1). M4. Return best (possibly none) and the dead columns.

The search limit (default 4, LuOptions::search_limit) is Zlatev’s restriction: the merit found after a few lines is almost always the minimum, and the full search would cost O(nnz(S)) per pivot. If M returns no pivot and no dead column, every active column was examined and found empty, so L5 terminates.

Invariants

Checked by LuFactor::check_invariants in the tests and after every factorization in debug builds. (F1) The non-vacant slots of (p_s) and (q_s) are permutations of 0..m, and the inverse maps step(p_s) = step(q_s) = s agree with them (right after a factorization there are no vacant slots; each update adds one, Chapter 5.5). (F2) Every L multiplier l_is of step s has step(i) > s, where step is the order of the last factorization (updates reorder U but leave L alone). (F3) Every off-diagonal u_ij has step(i) < step(j); the row and column copies of U hold exactly the same entries; every diagonal entry is nonzero. (F4) Right after a factorization, L U = P B Q up to rounding; tested by property tests against the dense product.

The two files of S (task P3.37)

The column file and the row file of S are cross-linked (Chapter 5.2, Algorithm T): every entry holds a hint of the position of its twin in the other file, exact in long lines. Every deletion of L2, L6 and of the cancellations of L6, and every lookup of s_ij from the row side (L4, M3), is then O(1) in a long line instead of a linear search along it (short lines are still searched, which is cheaper than keeping their hints exact). The deletions keep the lines in the same order as the searches did, so the elimination — pivots, fill, the order of L and U — is the same, bit for bit. Invariant (S3) of Chapter 5.2 is checked after the load and after the triangular phase in debug builds.

The objects (second edition)

Lu<T, U> is the factorization object of the contract: the options LuOptions, the stored factors Factors<T> (the steps, L, U by rows and by columns, the spike, the work arrays and the statistics LuStats), the basis-update policy U (Chapter 5.5, crate::ForrestTomlin by default, which owns the row etas) and the moving averages of the density policy (Chapter 5.4). The simplex method programs against the trait LuFactor. In the first edition all of this was one struct, Factor; splitting it moves no arithmetic: the update reaches the factors through &mut Factors<T>, and the type parameter U is resolved at compile time.

LuOptions merges the first edition’s FactorOptions with the two tolerances the factorization borrowed from the solver’s Tolerances: markowitz became LuOptions::rel_pivot_tol and zero became LuOptions::drop_tol (the solver copies them in when it creates a factor).

Complexity

The triangular phase costs O(nnz(B)) plus one division per L entry. A Markowitz pivot costs O(Σ_j c_j) over the columns of its U row, plus the bounded search. (Before task P3.37 every deletion also searched its line, O(Σ_i r_i²) per factorization: on hypothyroid-k1 half the time of a factorization.)

References

Markowitz, “The elimination form of the inverse and its application to linear programming”, Management Sci. 3, 1957; Suhl & Suhl, “Computing sparse LU factorizations for large-scale linear programming bases”, ORSA J. Comput. 2(4), 1990 (triangular phase, count lists, search); Zlatev, “On some pivotal strategies in Gaussian elimination by sparse technique”, SIAM J. Numer. Anal. 17, 1980 (restricted search); Maros (2003), ch. 5 (§5.4 LU of the basis) and ch. 8 (§8.2 bump); Davis, Direct Methods for Sparse Linear Systems (2006), ch. 6; Duff, Erisman & Reid (1986), ch. 7 and §9.4.

Structs

Factors

The stored factors P, L, U, Q of a basis (row and column files, MA28/Suhl style), the recorded spike, the work arrays of the solves and the statistics. Opaque: the basis-update policies reach it through BasisUpdate.

Lu

Sparse LU factors of a basis with a basis-update policy U (Forrest–Tomlin by default): the contract’s Lu<T, U>.

LuOptions

Options of the factorization, of the solves and of the update policy.

LuStats

Statistics of the current factor, for logs and the phase-1 benchmark. The update counters are in UpdateStats.

Enums

PreOrder

Column pre-orders available to the dynamic LU (Lu::set_pre_order; the contract’s PreOrder, which it placed in LuOptions::pre_order).

Traits

LuFactor

What every LU of the crate exposes; the simplex method programs against this trait (the contract’s LuFactor<T>).