Minimum degree ordering in Python
Greedily eliminate the vertex of smallest degree to limit fill-in.
- File
minimum_degree_ordering.py, 167 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/minimum-degree-ordering/. The authors are listed in prov/3FDE64.jsonl.
"""Minimum degree ordering of a sparse symmetric matrix.
A graph is a list adj with adj[v] the neighbours of v (no self loops), the
pattern of a symmetric matrix A. An ordering is a list order; vertex order[k]
is eliminated k-th, so P A P^T has (P A P^T)[k][l] = a[order[k]][order[l]].
minimum_degree works on the elimination graph (explicit fill edges).
minimum_degree_quotient works on the quotient graph of George and Liu, whose
storage never exceeds that of the original graph, with exact degrees or with
the approximate degree of Amestoy, Davis and Duff (AMD).
Ties are broken by the smallest vertex index, so all results are reproducible.
"""
import heapq
def graph_from_edges(n: int, edges) -> list[list[int]]:
adj = [set() for _ in range(n)]
for v, w in edges:
if v != w:
adj[v].add(w)
adj[w].add(v)
return [sorted(s) for s in adj]
def graph_from_dense(a) -> list[list[int]]:
n = len(a)
return [[j for j in range(n) if j != i and (a[i][j] != 0 or a[j][i] != 0)] for i in range(n)]
def minimum_degree(adj) -> list[int]:
"""Minimum degree on the elimination graph; ties by smallest index."""
n = len(adj)
g = [set(nb) for nb in adj] # @line L1
alive = [True] * n
heap = [(len(g[v]), v) for v in range(n)]
heapq.heapify(heap)
order: list[int] = []
while len(order) < n: # @line L2
d, v = heapq.heappop(heap) # @line L3
if not alive[v] or d != len(g[v]):
continue # stale heap entry
order.append(v) # @line L4
alive[v] = False
nb = g[v] # @line L5
for w in nb: # @line L6
g[w].discard(v)
g[w] |= nb - {w}
g[v] = set()
for w in nb: # @line L7
heapq.heappush(heap, (len(g[w]), w))
return order # @line L8
def minimum_degree_quotient(adj, degree: str = "exact", stats=None) -> list[int]:
"""Minimum degree on the quotient graph.
Variables i keep A[i] (variable neighbours) and E[i] (adjacent elements);
an eliminated vertex p becomes an element with variable list L[p].
degree = "exact" uses |reach(i)|; degree = "approximate" uses the AMD bound.
If stats is a dict, stats["storage"] gets the largest value of
sum |A_i| + sum |L_e| over the run, stats["checks"] a list of
(approximate, exact, number of elements of i) for every degree update.
"""
n = len(adj)
a = [set(nb) for nb in adj]
e = [set() for _ in range(n)]
elems: dict[int, set[int]] = {}
alive = [True] * n
deg = [len(nb) for nb in adj]
heap = [(deg[v], v) for v in range(n)]
heapq.heapify(heap)
order: list[int] = []
storage = sum(deg)
while len(order) < n:
d, p = heapq.heappop(heap)
if not alive[p] or d != deg[p]:
continue
order.append(p)
alive[p] = False
lp = set(a[p])
for el in e[p]:
lp |= elems.pop(el) # absorb the elements adjacent to p
lp.discard(p)
absorbed = e[p]
elems[p] = lp
a[p] = set()
e[p] = set()
remaining = n - len(order)
for i in lp:
e[i] -= absorbed
e[i].add(p)
a[i] -= lp
a[i].discard(p)
for i in lp:
reach = set(a[i])
for el in e[i]:
reach |= elems[el]
reach.discard(i)
exact = len(reach)
if degree == "exact":
deg[i] = exact
else:
ext = sum(len(elems[el] - lp) for el in e[i] if el != p)
deg[i] = min(remaining - 1, deg[i] + len(lp) - 1, len(a[i]) + len(lp) - 1 + ext)
if stats is not None:
stats.setdefault("checks", []).append((deg[i], exact, len(e[i])))
heapq.heappush(heap, (deg[i], i))
size = sum(len(a[v]) for v in range(n) if alive[v]) + sum(len(s) for s in elems.values())
storage = max(storage, size)
if stats is not None:
stats["storage"] = storage
return order
def _numbers(order):
f = [0] * len(order)
for k, v in enumerate(order):
f[v] = k
return f
def cholesky_pattern(adj, order) -> list[set[int]]:
"""Exact symbolic Cholesky: cols[k] = rows i > k with l_ik != 0 (new numbers)."""
f = _numbers(order)
n = len(adj)
later = [set(f[w] for w in adj[v]) for v in order]
cols = []
for k in range(n):
s = {i for i in later[k] if i > k}
cols.append(s)
for i in s:
later[i] |= s - {i}
return cols
def fill_and_flops(adj, order) -> tuple[int, int]:
"""(nonzeros of L below the diagonal, sum_k c_k^2) for the ordering."""
cols = cholesky_pattern(adj, order)
return sum(len(c) for c in cols), sum(len(c) ** 2 for c in cols)
def fill_edges(adj, order) -> int:
"""Number of fill entries: nonzeros of L below the diagonal minus edges of G."""
return fill_and_flops(adj, order)[0] - sum(len(nb) for nb in adj) // 2
def ldlt_pattern_exact(a, order):
"""Exact LDL^T of P A P^T with Fractions; returns the set of (i, k), i > k, with l_ik != 0."""
from fractions import Fraction
n = len(a)
b = [[Fraction(a[i][j]) for j in order] for i in order]
nz = set()
for k in range(n):
piv = b[k][k]
for i in range(k + 1, n):
if b[i][k] != 0:
nz.add((i, k))
m = b[i][k] / piv
for j in range(k + 1, i + 1):
if b[k][j] != 0:
b[i][j] -= m * b[k][j]
b[j][i] = b[i][j]
return nz
test_minimum_degree_ordering.py 180 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/minimum-degree-ordering/. The authors are listed in prov/3FDE64.jsonl.
import itertools
import random
from minimum_degree_ordering import (
cholesky_pattern, fill_and_flops, fill_edges, graph_from_dense, graph_from_edges,
ldlt_pattern_exact, minimum_degree, minimum_degree_quotient,
)
def random_graph(rng, n, p):
return graph_from_edges(n, [(i, j) for i in range(n) for j in range(i) if rng.random() < p])
def grid(k):
e = []
for i in range(k):
for j in range(k):
v = i * k + j
if j + 1 < k:
e.append((v, v + 1))
if i + 1 < k:
e.append((v, v + k))
return graph_from_edges(k * k, e)
def rcm(adj):
"""A plain reverse Cuthill-McKee from a minimum-degree vertex, for comparison."""
n = len(adj)
seen, order = [False] * n, []
for s in sorted(range(n), key=lambda v: (len(adj[v]), v)):
if seen[s]:
continue
seen[s] = True
order.append(s)
h = len(order) - 1
while h < len(order):
v = order[h]
h += 1
for w in sorted((w for w in adj[v] if not seen[w]), key=lambda w: (len(adj[w]), w)):
seen[w] = True
order.append(w)
return order[::-1]
def test_edge_cases():
assert minimum_degree([]) == []
assert minimum_degree_quotient([]) == []
assert minimum_degree([[]]) == [0]
k5 = graph_from_edges(5, [(i, j) for i in range(5) for j in range(i)])
assert fill_edges(k5, minimum_degree(k5)) == 0
empty = [[] for _ in range(4)]
assert minimum_degree(empty) == [0, 1, 2, 3]
def test_permutation_and_greedy_choice():
rng = random.Random(11)
for _ in range(150):
n = rng.randint(1, 25)
g = random_graph(rng, n, rng.uniform(0.05, 0.5))
order = minimum_degree(g)
assert sorted(order) == list(range(n))
# replay: each chosen vertex has least degree in the current elimination graph
cur = [set(nb) for nb in g]
alive = set(range(n))
for v in order:
assert len(cur[v]) == min(len(cur[w]) for w in alive)
assert v == min(w for w in alive if len(cur[w]) == len(cur[v]))
nb = cur[v]
for w in nb:
cur[w].discard(v)
cur[w] |= nb - {w}
alive.discard(v)
def test_quotient_graph_matches_elimination_graph():
rng = random.Random(12)
for _ in range(150):
n = rng.randint(1, 30)
g = random_graph(rng, n, rng.uniform(0.05, 0.4))
st = {}
assert minimum_degree_quotient(g, "exact", st) == minimum_degree(g)
# storage of the quotient graph never exceeds that of the original graph
assert st["storage"] <= sum(len(nb) for nb in g)
def test_approximate_degree_is_an_upper_bound_and_exact_with_two_elements():
rng = random.Random(13)
exact_cases = 0
for _ in range(150):
n = rng.randint(2, 30)
g = random_graph(rng, n, rng.uniform(0.05, 0.4))
st = {}
order = minimum_degree_quotient(g, "approximate", st)
assert sorted(order) == list(range(n))
assert st["storage"] <= sum(len(nb) for nb in g)
for approx, exact, ne in st.get("checks", []):
assert approx >= exact
if ne <= 2:
assert approx == exact
exact_cases += 1
assert exact_cases > 1000
def test_tree_has_no_fill():
rng = random.Random(14)
for _ in range(100):
n = rng.randint(1, 60)
g = graph_from_edges(n, [(i, rng.randrange(i)) for i in range(1, n)])
assert fill_edges(g, minimum_degree(g)) == 0
assert fill_edges(g, minimum_degree_quotient(g, "approximate")) == 0
def test_not_optimal_on_two_cliques():
def clique(s):
return [(a, b) for a in s for b in s if a < b]
g = graph_from_edges(9, clique([0, 1, 2, 3]) + clique([5, 6, 7, 8]) + [(4, 0), (4, 5)])
order = minimum_degree(g)
assert order[0] == 4 # the unique vertex of degree 2
assert fill_edges(g, order) == 1
assert fill_edges(g, [1, 2, 3, 6, 7, 8, 0, 4, 5]) == 0 # a perfect elimination ordering
def test_against_brute_force_minimum_fill():
rng = random.Random(15)
worse = 0
for _ in range(25):
n = 6
g = random_graph(rng, n, rng.uniform(0.3, 0.6))
best = min(fill_edges(g, list(p)) for p in itertools.permutations(range(n)))
md = fill_edges(g, minimum_degree(g))
assert md >= best
worse += md > best
assert worse < 25 # usually optimal on tiny graphs, but not always
def test_symbolic_pattern_matches_exact_ldlt():
rng = random.Random(16)
equal = 0
for _ in range(40):
n = rng.randint(1, 12)
g = random_graph(rng, n, rng.uniform(0.1, 0.5))
a = [[0] * n for _ in range(n)]
for i in range(n):
for j in g[i]:
if j < i:
a[i][j] = a[j][i] = -rng.randint(1, 9)
for i in range(n):
a[i][i] = 1 + sum(abs(x) for x in a[i]) # strictly diagonally dominant, so SPD
assert graph_from_dense(a) == g
order = minimum_degree(g)
sym = {(i, k) for k, col in enumerate(cholesky_pattern(g, order)) for i in col}
num = ldlt_pattern_exact(a, order)
assert num <= sym
equal += num == sym
assert equal >= 36 # cancellation is rare with random values
def test_grid_fill_against_natural_and_rcm():
for k, ratio in ((8, 0.8), (16, 0.7), (24, 0.6)):
g = grid(k)
n = k * k
nat = fill_and_flops(g, list(range(n)))
r = fill_and_flops(g, rcm(g))
md = fill_and_flops(g, minimum_degree(g))
amd = fill_and_flops(g, minimum_degree_quotient(g, "approximate"))
assert md[0] < r[0] <= nat[0]
assert md[1] < ratio * r[1]
assert abs(amd[0] - md[0]) <= 0.05 * md[0]
def test_three_by_three_grid_example():
g = grid(3)
order = minimum_degree(g)
assert order == [0, 2, 6, 8, 1, 3, 4, 5, 7]
assert fill_edges(g, order) == 5
assert fill_edges(g, list(range(9))) == 8
assert fill_and_flops(g, order) == (17, 39)
What the next level needs
- the tests did not run in this build