Fill-reducing ordering

Deep dive

How it works

Parts:

• Elimination graph — adjacency sets of A's graph.
• Ordering heuristic — picks the next vertex: minimum degree/AMD, nested dissection, RCM.
• Symbolic factorization — replays elimination on the pattern; the elimination tree (parent of j = first later vertex adjacent to j in the filled graph) gives column counts and nnz(L).
• Permutation P — the output; the solver factors PAPᵀ (sparse direct solver).

1. Graph. Vertex per unknown, edge per off-diagonal nonzero (A+Aᵀ if unsymmetric).
2. Fill rule. Eliminate v: for every pair of remaining neighbours (a, b), add edge a-b if absent. Each added edge is one fill nonzero.
3. Minimum degree. Eliminate a vertex of least current degree, repeatedly. AMD (Amestoy, Davis, Duff 1996) uses cheap upper bounds on degrees in a quotient graph (eliminated vertices merged into elements), so no clique is stored.
4. Nested dissection. Find a small separator S splitting the graph in halves; order both halves first, S last; recurse. No fill crosses S and the halves are independent. METIS finds separators by multilevel coarsening.
5. RCM. BFS from a peripheral vertex, visit neighbours by increasing degree, reverse: small bandwidth, fill stays in a narrow band.
6. Unsymmetric/QR. COLAMD orders the columns of A so LU of A(:,P) (or Cholesky of AᵀA) stays sparse, without forming AᵀA.

 natural order           nested dissection (k=7)
 x x x x x x x           A | S | B     S = middle line, ordered last
 x x x x x x x           A | S | B     A, B: independent subtrees
 x x x x x x x           ----S----     of the separator tree,
 ...                     recurse in A and B   factor in parallel

Code

use std::collections::BTreeSet;

type Graph = Vec<BTreeSet<usize>>;

/// 5-point Laplacian graph on a k x k grid (vertex = r*k + c).
fn grid(k: usize) -> Graph {
    let mut g = vec![BTreeSet::new(); k * k];
    for r in 0..k {
        for c in 0..k {
            let v = r * k + c;
            if c + 1 < k { g[v].insert(v + 1); g[v + 1].insert(v); }
            if r + 1 < k { g[v].insert(v + k); g[v + k].insert(v); }
        }
    }
    g
}

/// Eliminate v: its neighbours become a clique; returns the new edges.
fn eliminate(g: &mut Graph, v: usize) -> usize {
    let nb: Vec<usize> = std::mem::take(&mut g[v]).into_iter().collect();
    for &u in &nb { g[u].remove(&v); }
    let mut fill = 0;
    for (i, &a) in nb.iter().enumerate() {
        for &b in &nb[i + 1..] {
            if g[a].insert(b) { g[b].insert(a); fill += 1; }
        }
    }
    fill
}

/// Symbolic elimination: number of fill edges.
fn fill(g: &Graph, order: &[usize]) -> usize {
    assert!(is_perm(order, g.len()));
    let mut g = g.clone();
    order.iter().map(|&v| eliminate(&mut g, v)).sum()
}

fn is_perm(o: &[usize], n: usize) -> bool {
    let mut seen = vec![false; n];
    o.len() == n && o.iter().all(|&v| v < n && !std::mem::replace(&mut seen[v], true))
}

/// Exact minimum degree (ties: lowest index); O(n^2) scan, small n only.
fn min_degree(g: &Graph) -> Vec<usize> {
    let mut g = g.clone();
    let mut done = vec![false; g.len()];
    let mut order = Vec::with_capacity(g.len());
    for _ in 0..g.len() {
        let v = (0..g.len()).filter(|&v| !done[v]).min_by_key(|&v| (g[v].len(), v)).unwrap();
        done[v] = true;
        order.push(v);
        eliminate(&mut g, v);
    }
    order
}

/// Nested dissection: middle line of the longer side last, halves first.
fn nested_dissection(k: usize) -> Vec<usize> {
    fn rec(k: usize, r: (usize, usize), c: (usize, usize), out: &mut Vec<usize>) {
        let (h, w) = (r.1 - r.0, c.1 - c.0);
        if h == 0 || w == 0 { return; }
        if h * w <= 4 {
            for i in r.0..r.1 { for j in c.0..c.1 { out.push(i * k + j); } }
        } else if h >= w {
            let m = r.0 + h / 2; // separator row
            rec(k, (r.0, m), c, out);
            rec(k, (m + 1, r.1), c, out);
            for j in c.0..c.1 { out.push(m * k + j); }
        } else {
            let m = c.0 + w / 2; // separator column
            rec(k, r, (c.0, m), out);
            rec(k, r, (m + 1, c.1), out);
            for i in r.0..r.1 { out.push(i * k + m); }
        }
    }
    let mut out = Vec::with_capacity(k * k);
    rec(k, (0, k), (0, k), &mut out);
    out
}

fn main() {
    // Path 0-1-2-3: no fill in order; eliminating 1 first adds edge 0-2.
    let mut path = vec![BTreeSet::new(); 4];
    for v in 0..3 { path[v].insert(v + 1); path[v + 1].insert(v); }
    assert_eq!(fill(&path, &[0, 1, 2, 3]), 0);
    assert_eq!(fill(&path, &[1, 0, 2, 3]), 1);
    // Arrow, hub 0 and 5 leaves: hub first fills C(5,2) = 10, hub last 0.
    let mut arrow = vec![BTreeSet::new(); 6];
    for v in 1..6 { arrow[0].insert(v); arrow[v].insert(0); }
    assert_eq!(fill(&arrow, &[0, 1, 2, 3, 4, 5]), 10);
    assert_eq!(fill(&arrow, &[1, 2, 3, 4, 5, 0]), 0);

    for k in [15, 31] {
        let g = grid(k);
        let nat: Vec<usize> = (0..k * k).collect();
        let (md, nd) = (min_degree(&g), nested_dissection(k));
        assert!(is_perm(&md, k * k) && is_perm(&nd, k * k));
        let (fn_, fm, fd) = (fill(&g, &nat), fill(&g, &md), fill(&g, &nd));
        println!("k={k:2} n={:4} edges={} fill: natural={fn_} min-degree={fm} nested-dissection={fd}",
                 k * k, g.iter().map(|s| s.len()).sum::<usize>() / 2);
        assert!(fm < fn_ && fd < fn_);
    }
}

Compiled and run (release); all asserts pass, output below.

Complexity and constants

Measured with the code above (fill = edges added during elimination; the 5-point grid has 2k(k-1) edges, so nnz of L below the diagonal is edges + fill):

• Natural row order is the banded case (fill grows about n^1.5). Exact minimum degree won at these sizes; nested dissection wins asymptotically (2D: nnz(L) O(n log n), flops O(n^1.5), George 1973). The crude middle-line separators here skip the leaf refinement real codes use.
• This code orders in O(n·(n+fill)); production AMD is near linear in nnz(A) in practice (order of).
• Cost drivers: nnz(L) bytes (8 B per value plus indices), flops ~ Σ colcount².

Efficiency

Cost: a one-off per sparsity pattern, small next to the factorization it shrinks; pointer-chasing graph work, memory-latency bound, while the factorization is compute-bound on dense blocks.

Levers:

• Any fill-reducing order — 3.1× less fill (min-degree), 2.2× (dissection) at k=31 (measured above) — always.
• Approximate degrees, quotient graph (AMD) — no stored cliques (Amestoy et al. 1996) — large matrices.
• Reuse permutation and symbolic result — no re-ordering — pattern fixed, values change.
• Dense rows last — a hub first costs a full clique (arrow test: 10 fill vs 0) — arrow-like rows.
• Separator tree — independent subtrees factor on separate threads or nodes — parallel direct solvers.

Naive → standard → state of the art:

• Natural order · fill ~ n^1.5 in 2D
• AMD (COLAMD if unsymmetric) · near-linear ordering time, good fill
• Multilevel nested dissection (METIS), minimum degree at the leaves · best fill on big meshes, tree-parallel

Implementation notes

• Work on the pattern only (CSR/CSC index arrays, no values); symmetrize A+Aᵀ first. AMD, METIS and COLAMD take these arrays directly.
• After ordering, read nnz(L) and flops off the elimination tree; check against RAM.

Pitfalls

• Using RCM to cut fill: it targets bandwidth/profile and is usually worse than AMD for fill.
• Unsymmetric A: order with COLAMD or the pattern of A+Aᵀ, not one triangle.

Variants and alternatives

In production

• Structural engineering: FEM codes order stiffness matrices with AMD or METIS before CHOLMOD/MUMPS solves.
• Circuit simulation: KLU (SuiteSparse) is designed for SPICE-like circuit matrices and uses AMD or COLAMD after a block-triangular permutation.
• Robotics SLAM: GTSAM's Ordering supports COLAMD, METIS, natural and custom orders for factor graphs.
• Libraries: Eigen's simplicial Cholesky defaults to AMDOrdering; SciPy spsolve to COLAMD (permc_spec); MUMPS via ICNTL(7).

Implementations

Sources

• Patrick R. Amestoy, Timothy A. Davis, Iain S. Duff, An approximate minimum degree ordering algorithm, SIAM J. Matrix Anal. Appl. 17(4), 1996. paper · code
• Alan George, Nested dissection of a regular finite element mesh, SIAM J. Numer. Anal. 10(2), 1973. paper
• George Karypis, Vipin Kumar, A fast and high quality multilevel scheme for partitioning irregular graphs, SIAM J. Sci. Comput. 20(1), 1998. paper · code
• Elizabeth Cuthill, James McKee, Reducing the bandwidth of sparse symmetric matrices, Proc. 24th ACM National Conference, 1969. paper
• Timothy A. Davis, John R. Gilbert, Stefan I. Larimore, Esmond G. Ng, A column approximate minimum degree ordering algorithm, ACM TOMS 30(3), 2004. paper
• Mihalis Yannakakis, Computing the minimum fill-in is NP-complete, SIAM J. Algebraic Discrete Methods 2(1), 1981. paper