Module solve 

Chapter 5.4 — FTRAN and BTRAN, dense and hypersparse

Purpose

The simplex method needs α_q = B⁻¹ a_q (FTRAN, forward transformation) and ρ_r = B⁻ᵀ e_r, y = B⁻ᵀ c_B (BTRAN, backward transformation). This chapter applies the factors of Chapter 5.3 and the row etas of Chapter 5.5 to a right-hand side, in one of two ways:

• the dense Algorithms F and B run over every elimination step, skipping exact zeros: O(m + nnz(L) + nnz(U) + nnz(etas)) whatever the right-hand side;
• the hypersparse Algorithms HF and HB first compute, by depth-first search in the graphs of the factors, the set of entries that can become nonzero (the reach of the right-hand side, Gilbert & Peierls 1988), and then perform only the operations of the reach, in a topological order. Their cost is proportional to the number of entries and factor nonzeros actually involved — on the large sparse LPs (Kennington ken-*, pds-*) a handful out of m ≈ 10⁵, which is what Hall & McKinnon (2005) call hypersparsity.

The choice between the two is automatic, per call (policy below).

Notation

After t updates the factor represents

  B = Pᵀ L R₁⁻¹ ⋯ R_t⁻¹ U Qᵀ,

where R_e = I − e_{p} μ_eᵀ is the row eta of update e (target row p, multipliers μ_e on rows pivoted after p). FTRAN maps a vector indexed by rows to one indexed by positions; BTRAN the reverse. For a row p let q(p) = q_{step(p)} be the position pivoted with it, and for a position q let p(q) = p_{step(q)}.

The four graphs of the hypersparse solves have an arc a → b when the value at b is modified by an operation that reads the value at a:

• G_L on rows: p → i for every multiplier l_ip of the L eta with pivot row p (column-wise L, as stored);
• G_U on rows: p → i for every off-diagonal u_{i q(p)} of column q(p) of U (U by columns);
• G_Uᵀ on positions: q → j for every off-diagonal u_{p(q) j} of row p(q) of U (U by rows);
• G_Lᵀ on rows: i → p for every multiplier l_ip (L by rows, a transposed copy built once per factorization).

The row etas are applied in creation order (FTRAN) or in reverse creation order (BTRAN); their dependencies are held by two lists per row: the etas that read row k (k among their p_k) and the etas that target row p.

Algorithm F (dense FTRAN: solve B x = b)

F1. [L.] For each L eta in pivot order (pivot row p, multipliers l_ip): if y_p ≠ 0, set y_i ← y_i − l_ip y_p. F2. [R.] For each row eta in creation order (target p): set y_p ← y_p − Σ_k μ_k y_{p_k}. The vector y now holds the spike R L⁻¹ b — recorded by crate::LuFactor::ftran_for_update. F3. [U.] For the steps s in decreasing order, with pivot (p_s, q_s): x_{q_s} ← y_{p_s} / u_{p_s q_s}; if nonzero, y_i ← y_i − u_{i q_s} x_{q_s} for the entries of column q_s of U (all in rows of earlier steps).

Algorithm B (dense BTRAN: solve Bᵀ y = c)

B1. [Uᵀ.] For the steps s in increasing order: w_{p_s} ← c_{q_s} / u_{p_s q_s}; if nonzero, c_j ← c_j − u_{p_s j} w_{p_s} for the entries of row p_s of U (all in positions of later steps). B2. [Rᵀ.] For each row eta in reverse creation order (target p): if w_p ≠ 0, w_{p_k} ← w_{p_k} − μ_k w_p. B3. [Lᵀ.] For each L eta in reverse pivot order (pivot row p): w_p ← w_p − Σ_i l_ip w_i. Then y = w.

Correctness: F1–F3 apply L⁻¹, R_t ⋯ R_1 and U⁻¹ (with the permutations implicit in the indexing) to b; B1–B3 apply the transposes in the opposite order.

Algorithm D (reach: depth-first search, Gilbert & Peierls)

Given a graph G and a set S of start nodes, produce the nodes reachable from S in a topological order (every arc a → b between reached nodes has a before b).

D1. [Start.] Take a fresh mark (a stamp: no clearing). Empty the output list post. For each s ∈ S not yet marked, mark it and push (s, 0) on the stack; then run D2 until the stack is empty. D2. [Step.] Let (v, k) be the top of the stack. If v has a k-th successor w, advance k; if w is unmarked, mark it and push (w, 0). Otherwise pop v and append it to post (v is finished). D3. [Order.] A node finishes after every node reachable from it, so post read backwards is a topological order.

The stack is explicit (no recursion: the paths in G_L can be 10⁵ long). Cost: O(|reach| + arcs out of the reach).

Algorithm HF (hypersparse FTRAN)

The right-hand side b is given with its pattern. HF1. [Load.] Copy the nonzeros of b into the row work array y; empty b (it receives the result). HF2. [Reach of L.] Run Algorithm D on G_L from the rows of HF1. HF3. [Apply L.] For each reached row p in topological order: if p is the pivot of an L eta and y_p ≠ 0, set y_i ← y_i − l_ip y_p. The reached rows form the list Λ (marked). HF4. [Row etas.] Put in a min-priority queue every eta that reads a row of Λ. While the queue is not empty, take the smallest eta e (target t): y_t ← y_t − Σ_k μ_k y_{p_k}; if t ∉ Λ, add it to Λ and queue the etas after e that read t. (An eta before e that reads t was applied — or skipped — while y_t was still zero, exactly as in F2.) For crate::LuFactor::ftran_for_update, record the spike y on Λ. HF5. [Reach of U.] Run Algorithm D on G_U from Λ. HF6. [Apply U.] For each reached row p in topological order: x_{q(p)} ← y_p / u_{p q(p)}, y_p ← 0; if x_{q(p)} ≠ 0, set y_i ← y_i − u_{i q(p)} x_{q(p)} along column q(p). The pattern of the result is the set of positions q(p) of the reach.

Algorithm HB (hypersparse BTRAN)

HB1. [Load.] Copy the nonzeros of c into the position work array; empty c. HB2. [Reach of Uᵀ.] Run Algorithm D on G_Uᵀ from the positions of HB1. HB3. [Apply Uᵀ.] For each reached position q in topological order: w_{p(q)} ← c_q / u_{p(q) q}, c_q ← 0; if nonzero, c_j ← c_j − u_{p(q) j} w_{p(q)} along row p(q). The rows p(q) form the list Λ (marked). HB4. [Row etas, transposed.] Put in a max-priority queue every eta whose target is in Λ. While the queue is not empty, take the largest eta e (target t): if w_t ≠ 0, w_{p_k} ← w_{p_k} − μ_k w_t; for every p_k ∉ Λ, add it to Λ and queue the etas before e that target p_k. (The etas after e that target p_k already ran, with w_{p_k} = 0: no-ops in B2 too.) HB5. [Reach of Lᵀ.] Run Algorithm D on G_Lᵀ from Λ. HB6. [Apply Lᵀ.] For each reached row i in topological order: if w_i ≠ 0, set w_p ← w_p − l_ip w_i for every multiplier of row i of L. The result y = w has the reach as pattern.

The direction of the arcs is the delicate point. In FTRAN L and U are used by columns (scatter: a known value pushes into the rows it modifies), so the graphs are those of the stored columns. In BTRAN the dense algorithm uses U by rows (B1 is already a scatter, along row p(q): arcs q → j) but L by columns as a gather (B3 sums over the column of eta p): a gather cannot be driven by a reach, so HB6 turns it into a scatter along the rows of L (i → p), which needs the transposed copy of L. The row etas are the opposite: F2 is a gather (driven by the “readers” lists) and B2 a scatter (driven by the “targets” lists).

Choice of the algorithm (alternation policy)

Following Hall & McKinnon (2005, §4) and HiGHS, the hypersparse algorithm is used for a call iff (a) the right-hand side is in sparse mode (it has a pattern); (b) its density count/m is at most LuOptions::hyper_rhs_density (0.10); (c) the expected density of the result — an exponential moving average, with weight HYPER_HISTORY_WEIGHT for the newest result, of the result densities of the previous calls of the same kind (FTRAN, BTRAN, FTRAN for update; kept across factorizations, starting at 0) — is at most LuOptions::hyper_result_density (0.10).

When the result fills in, the reach costs more than the dense sweep (the search visits every arc once, then the numerical pass visits them again), so (c) switches a kind of operation to the dense algorithms after about three dense results, and back once results are sparse again (the dense algorithms measure their result density exactly). The per-kind statistics (SolveStats) are in crate::LuStats.

Division of labour (second edition)

Steps F2, B2, HF4 and HB4 — the row etas — belong to the basis-update policy (Chapter 5.5): they are the methods BasisUpdate::apply_ftran, BasisUpdate::apply_btran, BasisUpdate::apply_ftran_hyper and BasisUpdate::apply_btran_hyper of crate::ForrestTomlin, called at exactly the place where the first edition ran them inline, with the same work arrays. Everything else is here, on the stored factors.

Invariants

(V1) On return the SparseVec is in sparse mode. After a dense solve its pattern is the exact support of the result (in increasing order); after a hypersparse solve it is the symbolic reach, a superset of the support that may contain explicit zeros (allowed by the invariant (V2) of SparseVec; call SparseVec::tidy to drop zeros and tiny entries). Computing the exact support would cost one more pass for no benefit to the simplex method. (V2) Both work arrays are zero between calls. (V3) The dense and the hypersparse algorithms perform the same operations except for the order of independent updates (and HB6, which sums the terms of B3 one at a time): their results agree to rounding (property test: relative error ≤ 1e-12).

Complexity

Dense: O(m + nnz(L) + nnz(U) + nnz(etas)). Hypersparse: O(r + f + e log e) with r the reach sizes, f the factor nonzeros in the columns (or rows) of the reach and e the number of row etas involved.

References

J. R. Gilbert and T. Peierls, “Sparse partial pivoting in time proportional to arithmetic operations”, SIAM J. Sci. Stat. Comput. 9(5), 1988 (the reach by depth-first search); J. A. J. Hall and K. I. M. McKinnon, “Hyper-sparsity in the revised simplex method and how to exploit it”, Comput. Optim. Appl. 32, 2005 (hypersparse FTRAN/BTRAN, the density thresholds); T. A. Davis, Direct Methods for Sparse Linear Systems (2006), §3.2 (cs_reach, cs_dfs); Maros (2003), §5.4–5.5; Forrest & Tomlin (1972); Suhl & Suhl (1990), §4; Koberstein (2005), §5.3; HiGHS HFactor (studied, not copied).

Structs

SolveStats

Counters of one kind of solve (FTRAN, BTRAN or FTRAN for update), accumulated over the life of the Lu (not reset by a factorization). The contract’s mean densities are the methods SolveStats::mean_rhs_density and SolveStats::mean_result_density over exact sums (no division per call).

Constants

HYPER_HISTORY_WEIGHT

Weight of the newest result in the moving average of the result density of each kind of solve (policy (c) of the module documentation).