Skip to main content

Crate xingo_sparse

Crate xingo_sparse 

Source
Expand description

§Chapter 5 — xingo-sparse: sparse linear algebra for optimization

Second edition (PLAN.md trilha L, task L1.2; D-0080 and its two addenda). The first edition of this chapter was the crate xingo-linalg: the LU factorization of the simplex basis with Forrest–Tomlin updates and hypersparse solves. This edition moves it, together with the compressed containers of xingo-core (formerly Chapter 1.2), into a crate that does not depend on the solver and can be published on its own. The migration has no numerical effect: instantiated with f64, every operation below performs the same floating-point operations in the same order as the first edition, so the solver’s results are bit for bit those of before (checked on the MIPLIB anchors and the nl9 set, D-0080 addendum 2).

§5.A The layer inversion

Before: xingo-linalg → xingo-core (the factorization borrowed the containers, the typed indices, the error type and the tolerances of the solver’s core). Now: xingo-core → xingo-sparse, and xingo-sparse depends on nothing in the workspace. Each borrowed item received a home here:

first edition (in xingo-core)second edition (here)
CscMatrix, CsrMatrix, SparseVecCscMatrix<T>, CsrMatrix<T>, SparseVec<T> (§5.8), T = f64 by default
Col, Row (Chapter 1.1)Col, Row (§5.7), the same newtypes
CoreErrorSparseError (containers) and LuError (factorization), §5.0
Tolerances::{markowitz, zero}LuOptions::rel_pivot_tol, LuOptions::drop_tol
Tolerances::{lt, approx_eq, is_zero, …}the comparison helpers of §5.6 (crate-private)

xingo-core re-exports the containers at their old paths as the f64 instances (xingo_core::CscMatrix = xingo_sparse::CscMatrix<f64>, and so on), re-exports Col and Row, and converts SparseError and LuError into its CoreError; the solver maps its Tolerances into LuOptions when it creates a factor.

§5.B Scalars and genericity

Every container and the factorization are generic in the scalar T: Real (§5.6), with T = f64 the default type parameter, so CscMatrix written without parameters is CscMatrix<f64>. Real is the small trait a threshold-pivoting LU needs — the field operations, a total order on the values that occur (PartialOrd, NaN excluded by contract), abs, max, is_finite, conversion from and to f64 — and is implemented here for f64 only; xingo-extprec (trilha L2) will implement it for double-double and quad-double. Intervals are not Real (no total order, no single rounding; D-0080, Q2).

Genericity costs nothing (D-0080, addendum 2): the code is generic by monomorphization only — no trait objects in the kernels — and every method of Real for f64 is an #[inline] call of the inherent f64 method, so Lu<f64> compiles to the same machine operations as the first edition’s Factor. Options and statistics stay in f64 (they are thresholds and counters, not data); a threshold enters the kernels through Real::from_f64, the identity for f64.

§5.C Comparisons

The rule of the solver — “never compare f64 outside xingo_core::Tolerances” — cannot hold in a crate below xingo-core. Its equivalent here: every comparison of floating-point values goes through the helpers of §5.6 (lt, le, approx_eq, is_zero, is_exact_zero, the same definitions as Tolerances, generic in T), and every threshold of the factorization is a field of LuOptions — never a bare literal in a kernel. The two constants that remain are policy, documented where they are defined: the clearing fraction of SparseVec::clear (0.3, Chapter 5.8.3) and the moving-average weight HYPER_HISTORY_WEIGHT (Chapter 5.4).

§5.D Contract and implementation

The approved contract (contracts/xingo-sparse/src/lib.rs, L1.1) is the specification of the whole v0.1. This edition realizes the part listed for L1.2; its divergences from the contract are recorded in the cycle’s decisions and listed for ratification.

contracthere
RealReal (adds Sum, needed for bit-identical sums)
CscMatrix<T>, CsrMatrix<T>, SparseVec<T>the same layouts; first-edition methods kept
Triplets<T>, BlockPartitionTriplets (§5.8.4; to_csc keeps (C3): a zero sum is not stored, zero inputs are legal; + from_csc, len, is_empty, clear, reserve, check_invariants), BlockPartition (§5.8.5; + new, from_sizes, check_invariants and the block accessors) (L1.7a)
SparseError, LuErrorSparseError, LuError (variants of the first edition, see §5.0)
ColumnRef, Sign, ColumnSource<T>ColumnRef, Sign, ColumnSource (static dispatch)
LuOptions, LuStats, SolveStats, UpdateStatsLuOptions (+ drop_tol), LuStats, SolveStats, UpdateStats
Lu<T, U = ForrestTomlin<T>>, LuFactor<T>, Factors<T>Lu, LuFactor, Factors
BasisUpdate<T>, ForrestTomlin<T>BasisUpdate, ForrestTomlin
UpdateOutcome, RefactorReasonUpdateOutcome, RefactorReason (unchanged)
Pattern, Perm, CscMatrix::{pattern, symmetric_pattern}Pattern (+ transpose, symmetric, ata, cholesky_nnz), Perm (L1.3)
Ordering, OrderingTarget, Natural, Amd, Colamdthe same (L1.3; Amd::default() has aggressive absorption on)
PreOrder, LuOptions::pre_order, SymbolicLu, LuStats::symbolic_reusesPreOrder, Lu::{set_pre_order, pre_order} (not an LuOptions field: see Lu::set_pre_order), SymbolicLu (+ singletons), LuStats::{symbolic_reuses, symbolic_fallbacks} (L1.3)
LuFactor::{refactorize_same_pattern, symbolic}the same (L1.3); the replay threshold Lu::{set_refactor_pivot_tol, refactor_pivot_tol} (L1.5b, default = rel_pivot_tol)
Equilibration, Scaling<T>Equilibration, Scaling (+ check_invariants; §5.13, L1.5)
Inertia, Regularization<T>, Refinement<T>, RefineReport<T>, KktForm, KktSolver<T>, KktErrorthe same (§5.14, L1.5), plus Default for the option structs, Regularization::{none, relative_to, is_quasi_definite}, KktForm::expected_inertia, KktError::Structure, EQUILIBRATION_PASSES; the augmented inertia reads {positive: m, negative: n}; L1.5b: Regularization::scale with ShiftScale (shifts in the equilibrated matrix by default), KktError::NotFinite, RefineWorkspace and KktSolver::solve_refined_with, refinement by GMRES
FaerLdl, FaerSymbolic, FaerNumeric<T> (feature faer)the same (§5.15, L1.5), plus read-only accessors
DenseBlock<T>, DenseKernel<T>, NaiveDense, FaerDense, DenseErrorthe same (§5.16, L1.5); sytrf also returns its permutation and 2×2 marks

§5.E What the tasks add

  • L1.3 (done) — Perm and Pattern (§5.9), the orderings Natural, Amd and Colamd behind Ordering (§5.10, our own; their engine, Algorithm Q, in §5.11), the reusable pivot sequence SymbolicLu with LuFactor::refactorize_same_pattern and the COLAMD pre-order of the bump Lu::set_pre_order (§5.12, step L5′ of §5.3). With pre_order = None (the default) Algorithm L is unchanged, bit for bit; the simplex uses none of it yet.
  • L1.4 (after P3.9) — BartelsGolub<T>, the second BasisUpdate; the trait already separates the update policy from the stored factors (Factors), so the simplex will not change.
  • L1.5 (done) — row/column scaling Scaling by the policies of Equilibration (§5.13); the KKT interface KktSolver with the regularization policy and the iterative refinement shared by every backend (§5.14); its v0.1 backend FaerLdl over faer’s sparse LDLᵀ/LBLᵀ ordered by our Amd (§5.15, feature faer); and the dense kernels DenseKernel, NaiveDense and FaerDense (§5.16). The simplex uses none of it, and faer is not in the solver’s build: the workspace depends on this crate without default features.
  • L1.5b (done) — before 0.1.0: a default regularization that holds at the end of an interior-point method (shifts in the equilibrated matrix, §5.14), refinement by right-preconditioned GMRES with reusable work vectors, non-finite factors reported, solves without allocation and exact bump counts (§5.15), the elimination-tree postorder before faer’s supernodes (§5.9, Algorithm P), a faster Algorithm Q (§5.11), and the replay threshold of Algorithm R (§5.12).
  • L1.7a (done) — what 0.1.0 ships besides the code: the last two containers of the contract, the coordinate builder Triplets (§5.8.4, Algorithm C′) and BlockPartition (§5.8.5, Algorithm B); the crate’s README, CHANGELOG and three runnable examples (examples/: a basis LU with a Forrest–Tomlin update, the orderings against the natural order, one KKT step of an interior-point method); complete item documentation. No solver path changes (nothing uses the two new types yet).
  • L1.6 — the own SupernodalLdl<T, K>, after the interior-point method has measured faer.

§Purpose (the first edition, unchanged)

The revised simplex method (Chapter 6) works with the basis matrix B only through three operations: factorize B, solve with B or Bᵀ, and replace one column of B. This crate provides them (PLAN.md §5 items 1–2, tasks P1.1–P1.2):

§Notation

B ∈ R^{m×m}; rows i and basis positions k both range over 0..m but are distinct index spaces. P B Q = L U with P, Q permutations, L unit lower and U upper triangular; after t updates B = Pᵀ L R₁⁻¹ ⋯ R_t⁻¹ U Qᵀ with row etas R_e and updated P, Q, U. Identifiers follow the chapters: l_ic, u_pj, mu (μ), alpha_r (α_r), tau (τ), sigma (σ). Indices are u32 in the stored arrays, pointers usize (PLAN.md §3.2; D-0080, Q3).

§Contract with the simplex method

  1. lu.factorize(&(&heading[..], &a)) (a ColumnSource); on LuError::Singular replace each cols[k] by ColumnRef::Logical(rows[k], _) and factorize again.
  2. ftran(&mut v) maps a row-indexed v to the position-indexed B⁻¹ v; btran(&mut v) maps a position-indexed v to the row-indexed B⁻ᵀ v. Both take &mut self: the factor owns the work arrays of the solves and the density statistics that choose between the dense and the hypersparse algorithms. Give the right-hand side in sparse mode with a correct pattern to allow the hypersparse algorithms; the result comes back in sparse mode, its pattern a superset of its support (Chapter 5.4, invariant (V1)).
  3. To pivot position r out for column a_q: ftran_for_update(&mut a_q) (giving α_q), then update(r, α_q[r]). On UpdateOutcome::RefactorNeeded (or any error) refactorize with the new heading; the factor never does so on its own (D-0016).

§Safety

Per D-0006 this is the one crate where unsafe is permitted (audited, // SAFETY: comments, Miri where executable). This edition, like the first, uses none; optimized kernels wait for measurements (task P1.10). The faer feature (L1.5) brings faer’s own unsafe kernels in as a dependency; the code of this crate stays free of unsafe.

§References

Suhl & Suhl (1990); Forrest & Tomlin (1972); Maros, Computational Techniques of the Simplex Method (2003), ch. 5 and 8; Davis, Direct Methods for Sparse Linear Systems (2006); Duff, Erisman & Reid (1986); Koberstein (2005), ch. 5; Gilbert & Peierls (1988); Hall & McKinnon (2005); Gustavson (1978); Amestoy, Davis & Duff (1996); Davis, Gilbert, Larimore & Ng (2004); the contract of the crate (contracts/xingo-sparse, D-0080).

Re-exports§

pub use basis::ColumnRef;
pub use basis::ColumnSource;
pub use basis::Sign;
pub use csc::CscMatrix;
pub use csr::CsrMatrix;
pub use dense::FaerDense;
pub use dense::DenseBlock;
pub use dense::DenseKernel;
pub use dense::NaiveDense;
pub use error::DenseError;
pub use error::KktError;
pub use error::LuError;
pub use error::SparseError;
pub use factor::Factors;
pub use factor::Lu;
pub use factor::LuFactor;
pub use factor::LuOptions;
pub use factor::LuStats;
pub use factor::PreOrder;
pub use faer_ldl::FaerLdl;
pub use faer_ldl::FaerNumeric;
pub use faer_ldl::FaerSymbolic;
pub use index::Col;
pub use index::Row;
pub use kkt::EQUILIBRATION_PASSES;
pub use kkt::Inertia;
pub use kkt::KktForm;
pub use kkt::KktSolver;
pub use kkt::RefineReport;
pub use kkt::RefineWorkspace;
pub use kkt::Refinement;
pub use kkt::Regularization;
pub use kkt::ShiftScale;
pub use ordering::Amd;
pub use ordering::Colamd;
pub use ordering::Natural;
pub use ordering::Ordering;
pub use ordering::OrderingTarget;
pub use partition::BlockPartition;
pub use perm::Pattern;
pub use perm::Perm;
pub use real::Real;
pub use scaling::Equilibration;
pub use scaling::Scaling;
pub use solve::HYPER_HISTORY_WEIGHT;
pub use solve::SolveStats;
pub use symbolic::SymbolicLu;
pub use triplets::Triplets;
pub use update::BasisUpdate;
pub use update::ForrestTomlin;
pub use update::RefactorReason;
pub use update::UpdateOutcome;
pub use update::UpdateStats;
pub use vector::SparseVec;

Modules§

basis
Chapter 5.1 — Basis columns
compressed
Chapter 5.8 — Sparse matrices and vectors
csc
Chapter 5.8.1 — Compressed sparse column storage
csr
Chapter 5.8.2 — Compressed sparse row storage
dense
Chapter 5.16 — Dense kernels for supernodal blocks
error
Chapter 5.0 — Errors
factor
Chapter 5.3 — Sparse LU factorization of the basis: Lu
faer_ldl
Chapter 5.15 — FaerLdl: the v0.1 KKT backend
index
Chapter 5.7 — Typed indices
kkt
Chapter 5.14 — KKT systems: analyze once, factorize many
ordering
Chapter 5.10 — Fill-reducing orderings: Natural, Amd, Colamd
partition
Chapter 5.8.5 — Block partitions
perm
Chapter 5.9 — Permutations and sparsity patterns
real
Chapter 5.6 — Scalars and comparisons
scaling
Chapter 5.13 — Row and column scaling
solve
Chapter 5.4 — FTRAN and BTRAN, dense and hypersparse
symbolic
Chapter 5.12 — The reusable pivot sequence: SymbolicLu
triplets
Chapter 5.8.4 — The coordinate builder
update
Chapter 5.5 — The Forrest–Tomlin update and the refactorization policy
vector
Chapter 5.8.3 — The sparse work vector