Module update 

Chapter 5.5 — The Forrest–Tomlin update and the refactorization policy

Purpose

A simplex iteration replaces one basic column: B' = B + (a_q − B e_r) e_rᵀ, where r is the leaving position. Refactorizing at every iteration would cost far more than the iteration itself, so the factors are updated. Forrest and Tomlin (1972) keep L untouched, replace one column of U by the transformed entering column (the spike) and restore triangularity by eliminating one row, recorded as a row eta. Their update keeps U as sparse as a fresh factorization in practice, which the product form of the inverse does not.

Notation

As in Chapter 5.4, B = Pᵀ L R₁⁻¹ ⋯ R_t⁻¹ U Qᵀ. Let s be the spike R_t ⋯ R_1 L⁻¹ a_q (indexed by rows, recorded by crate::LuFactor::ftran_for_update), α_q = B⁻¹ a_q and α_r its entry at position r. Let τ = step(r), p = p_τ the pivot row of that step and u_old = u_{p r} the diagonal entry of U there.

Algorithm U (Forrest–Tomlin update, classical variant)

U1. [Locate.] Find τ = step(r) and p = p_τ. Fail with LuError::NoSpike if no spike was recorded since the factor last changed. U2. [Spike row.] Copy row p of U (entries u_pj, step(j) > τ) into a dense work vector w indexed by positions, and let σ ← s_p. U3. [Eliminate.] Think of column r of U as replaced by s and of step τ as moved to the end of the pivot sequence (rows and columns shifted cyclically). Row p is then the only row below the diagonal and must be eliminated with the rows of the later steps, in step order: for the steps k > τ with j = q_k, if w_j is not negligible, set μ_k ← w_j / u_{p_k j}, w_j ← 0, w_l ← w_l − μ_k u_{p_k l} for the entries of row p_k, and σ ← σ − μ_k s_{p_k}. Since row p_k only has entries in later steps, the fill lands where it is eliminated later. U4. [Test.] The new diagonal is d = σ. Since det B' = α_r det B and the cyclic shifts of rows and columns cancel in sign, exact arithmetic gives d = α_r u_old. Let ζ = |s_p| + Σ_k |μ_k s_{p_k}|, the size of the terms summed into σ (its rounding error is O(ε ζ)). If |d| ≤ τ_abs max(1, ζ), the new basis is singular: return LuError::Singular with position r and row p. If |d − α_r u_old| > ε_FT max(|d|, |α_r u_old|), the update is numerically unreliable: return UpdateOutcome::RefactorNeeded(RefactorReason::Unstable). In both cases nothing has been modified: the factor still represents B. U5. [Commit.] Delete the old column r and the off-diagonal row p from both copies of U; insert the spike entries s_i (i ≠ p, |s_i| > ε_drop) as the new column r; set u_{pr} ← d; move step τ to the end of the sequence (its slot τ becomes vacant and the step is appended, so the step numbers grow beyond m − 1 until the next factorization; only their order matters); store the row eta R = I − e_p μᵀ if μ ≠ 0. Now R L⁻¹ P B' Q' = U' with U' upper triangular in the new step order. U6. [Policy.] Count the update and decide whether to ask for a refactorization (below).

The variant is the classical one: the spike is stored as a column of U and the elimination of row p uses a dense position-indexed work vector. The steps of U3 with a nonzero w_j are visited through a min-priority queue of steps (task P1.4) rather than by scanning τ+1..m, U5 inserts the spike from its recorded pattern and moves step τ to the end by leaving a vacant slot (no renumbering), so the cost of an update is proportional to the entries involved, as with the linked permutation lists of Suhl’s variant.

Refactorization policy

The factor never refactorizes by itself — the columns of B belong to the caller. After a successful update, crate::LuFactor::update returns UpdateOutcome::RefactorNeeded with the first applicable reason, the update already applied and the factor valid:

• RefactorReason::UpdateLimit — refactor_interval updates (default 100) since the last factorization;
• RefactorReason::FillGrowth — nnz(L) + nnz(U) + nnz(etas) exceeds fill_factor (default 3) times its value right after the factorization: solves have become more expensive than a refactorization would make them;
• RefactorReason::Growth — the updates have made U grow by more than growth_limit (default 10¹⁰): with g_0 = max |u_ij| right after the factorization and β = max |b_ij|, max |u_ij| > growth_limit · max(β, g_0).

The growth trigger is measured from the fresh factorization, not from B: a refactorization can only undo the growth the updates added, so a basis whose fresh factors already grow past the limit (MAROS-R7 of Netlib with the Markowitz threshold 0.01: max |u| / max |b| > 10¹⁰ straight out of the Markowitz phase) must not ask for one at every update — the old test, max |u_ij| / β > growth_limit, did exactly that and turned the solve into a loop of refactorizations (278 instead of 854 iterations per second). The absolute growth max |u_ij| / β stays available as a diagnostic (crate::LuStats::growth_exceeded); taking the larger of β and g_0 as the reference makes the trigger identical to the old one when the fresh factor did not grow (g_0 ≤ β). The cure of a grown fresh factor is a tighter threshold, a decision of the caller (crate::LuFactor::set_options with a larger crate::LuOptions::rel_pivot_tol; the simplex method of Chapter 6 does it).

Only RefactorReason::Unstable comes with the update not applied. The caller should refactorize before the next solve in every case; with the first three it may postpone (the factor is exact up to rounding).

The update as a policy (second edition)

The contract makes the update algorithm a type parameter of the factorization, Lu<T, U = ForrestTomlin<T>>, through the trait BasisUpdate (D-0080, addendum): an update policy owns its own data — here the row etas R_e and the lists of the etas that read or target each row — and reaches the stored factors through &mut crate::Factors. Besides Algorithm U it provides the steps of the solves that apply its data (F2, B2, HF4 and HB4 of Chapter 5.4), called by the solves at the point where the first edition ran them inline. Bartels–Golub (task L1.4) will be the second implementation. The trait is implementable inside the crate only (its arguments are opaque outside), which keeps the factors’ invariants in one place.

Invariants

(F1)–(F3) of Chapter 5.3 hold after every successful update (checked by the property tests through crate::LuFactor::check_invariants); FTRAN and BTRAN after k updates agree with a fresh factorization of the updated basis (property test, relative error ≤ 1e-9).

References

Forrest & Tomlin, “Updated triangular factors of the basis to maintain sparsity in the product form simplex method”, Math. Programming 2(1), 1972; Suhl & Suhl, “A fast LU update for linear programming”, Annals of OR 43, 1993; Maros (2003), §8.3 (Forrest–Tomlin) and §8.5 (stability checks); Koberstein (2005), §5.4; Huangfu & Hall (2018) for the comparison d versus α_r u_old.

Structs

ForrestTomlin

Forrest & Tomlin (1972): row etas and a permuted, modified U (Algorithm U).

UpdateStats

Counters of the update policy since the last factorization.

Enums

RefactorReason

Why crate::LuFactor::update asks for a refactorization.

UpdateOutcome

Result of a successful update.

Traits

BasisUpdate

A basis-update policy over the stored factors (the contract’s BasisUpdate<T>; see the module documentation). Every implementation is transactional — on Err and on RefactorReason::Unstable the factors are exactly as before — and never modifies L.