GPU Reasoning: Batching the Fixpoint

📍 Where we are: Part IV · Scale and Systems — Chapter 10. Materialization versus Rewriting weighed storing the closure against rewriting the query; this chapter takes the materialization lane and asks how fast the closure itself can be built when the fixpoint becomes matrix algebra.

Volume 4 had a moment worth remembering: in the Neural LP chapter, one Horn chain rule became a matrix product over relation adjacencies, and suddenly a piece of logic ran on the same hardware as a neural network. That demonstration hand-picked the friendliest possible rule, one binary relation composed with another. A real reasoning calculus is messier: Volume 2's EL completion keeps two different data structures (concept labels and role edges), fires six rule shapes that read and write both, and must run to an exact fixpoint, not a fixed number of hops. This chapter shows that the tensor move survives the upgrade. The companion module gpu_fixpoint.py re-expresses the entire completion calculus as boolean matrix algebra, iterates it to saturation, and then does something this volume insists on: it asserts cell-for-cell equality against Volume 2's own reasoner, the very oracle it vectorizes, on the real academic TBox and on three seeded synthetic ones. Nothing is learned and nothing is approximated; the payoff is purely architectural. A state that lives in matrices can be updated by whole-matrix operations, delta-filtered like a production Datalog engine, and stacked one hundred ontologies deep so a single pass closes them all. The chapter ends where the field keeps ending up: the measured, cited contrast between exact batched reasoning and a neural network trained to imitate a reasoner faster, at the price of soundness.

The simple version

Imagine knitting a scarf stitch by stitch versus weaving it on a loom. The knitter's hands visit one stitch at a time, deciding what to do next by looking at the neighboring stitches: that is Volume 2's reasoner, chasing pointers from fact to fact. A loom works differently: the pattern is set up once in the harness, and then one pull of the beam applies it to every thread across the whole row simultaneously. Same fabric, same pattern, radically different motion. This chapter re-rigs the completion rules as loom pulls: each rule becomes one matrix operation that fires everywhere at once, and a full pass of all the rules is one row of weaving. Two refinements complete the picture. First, a careful weaver does not re-weave rows that have not changed, only the fresh frontier (that is the delta trick). Second, if you have a hundred scarves to make with nearly the same pattern, you can thread all hundred into the loom at once and weave them with the same pulls (that is batching, and it is the entire reason graphics hardware is fast). The one thing the loom never does is change the fabric: every scarf is checked thread for thread against the hand-knitted original.

What this chapter covers

• The state as matrices: how the completion algorithm's subsumer sets and role edges become one boolean label matrix $S$ plus one adjacency matrix per role, every index decoded before any operation, and why the float32 trick underneath the boolean product is exact rather than approximate.
• Every rule a product: each completion rule derived into its matrix form and then quoted from the committed code, CR1 as a label-closure product, CR2 as a gather-and-join, CR3 as an outer-product emission, CR4 as a two-step join, the bottom rule as a matrix-vector product, and the role chain as graph composition; why the algebra is the OR-AND boolean semiring and where idempotence does the quiet work.
• The oracle standard: exact subsumption-set equality with el_completion.classify on the real 14-axiom TBox and on three seeded synthetic TBoxes, the asserts quoted, and why synthetic cross-checks matter when a single test ontology could hide rule-coverage gaps.
• Semi-naive deltas: the production-engine discipline restated and derived in matrix terms, with the committed operation counts, 71 derivation events naive against 28 with deltas, and the wavefront picture that connects straight back to Volume 1's fixpoint chapter.
• Batching as the GPU thesis: one hundred TBox variants closed in a single stacked tensor pass with no Python loop over the batch, the committed loop-versus-stacked timing table labelled honestly as a CPU proxy, and what real systems add on top: worklists, sparse kernels, and 128-core materialization.
• The honest ledger: what dense boolean matrices cannot see (memory that squares, provenance that vanishes, rule shapes that need new operators), and the neural-emulation alternative whose order-of-magnitude speedup buys throughput with soundness, the trade this volume keeps refusing in the exact lane.

From one chain rule to a whole calculus

Start by recalling what must be represented. Volume 2's completion algorithm saturates two structures over a finite universe of concept names: $S(A)$, the set of concepts that label $A$ (reading "$B$ labels $A$" as the derived subsumption $A \sqsubseteq B$), and $R(r)$, the set of derived role edges $(A, B)$ meaning every instance of $A$ has an $r$-successor in $B$. The rules are the classic consequence-based calculus of the EL family, the lineage that first classified SNOMED CT's 379,691 concepts and by 2006 did so in 1,782 seconds, roughly twice as fast as the one tableau reasoner that finished at all (FaCT++, at 3,859 seconds) while the other two ran out of memory [1], later modernized with concurrent worklist queues that brought the same saturation down to seconds [2]. Volume 2's el_completion.py implements that calculus in miniature the naive way: Python sets, membership tests, a changed flag (el_completion.py lines 177–261). It is the oracle for everything below, and its rule set is worth having in front of us. Decode the symbols it uses before reading it: ⊓ is concept conjunction ($A_1 \sqcap A_2$ is the concept of things that are both $A_1$ and $A_2$), $\exists r.B$ is the concept of things with at least one $r$-successor in $B$, $r \circ s$ is role composition (an $r$-edge followed by an $s$-edge), ⊥ reads "bottom", the impossible concept nothing can satisfy, and ⟹ reads "then add" (el_completion.py lines 25–30):

CR1  A' ∈ S(A),  A' ⊑ B          ⟹  add B to S(A)
CR2  A₁,A₂ ∈ S(A),  A₁ ⊓ A₂ ⊑ B  ⟹  add B to S(A)
CR3  A' ∈ S(A),  A' ⊑ ∃r.B       ⟹  add (A, B) to R(r)
CR4  (A,B) ∈ R(r),  B' ∈ S(B),  ∃r.B' ⊑ C  ⟹  add C to S(A)
CR⊥  (A,B) ∈ R(r),  ⊥ ∈ S(B)     ⟹  add ⊥ to S(A)   (the bottom rule)
CRχ  (A,B) ∈ R(r),  (B,C) ∈ R(s),  r ∘ s ⊑ t  ⟹  add (A,C) to R(t)

The matrixization is a change of data structure, nothing more. Fix the universe of concept names exactly as the oracle fixes it (every name occurring in the normalized axioms, plus ⊤, read "top", the concept that contains everything, and the bottom concept ⊥), sort it for determinism, and let $n$ denote its size, the number of concept names. Assign each name an integer index from $0$ to $n-1$. Then:

• the label state becomes one boolean matrix $S \in \lbrace 0,1 \rbrace^{n \times n}$ (an $n$-by-$n$ grid of true/false cells), with $S[A, B] = 1$ exactly when $B \in S(A)$, that is, when $A \sqsubseteq B$ has been derived;
• each role $r$ becomes an adjacency matrix $R_r \in \lbrace 0,1 \rbrace^{n \times n}$, with $R_r[A, B] = 1$ exactly when $(A, B) \in R(r)$;
• the axioms themselves become boolean tables: $H_1[A', B] = 1$ when the TBox states $A' \sqsubseteq B$, one existential-emission table $E^3_r[A', B] = 1$ per role when the TBox states $A' \sqsubseteq \exists r.B$, one trigger table $E^4_r[B', C] = 1$ per role when the TBox states $\exists r.B' \sqsubseteq C$, and an index array plus projection matrix for the conjunction axioms, described with CR2 below (gpu_fixpoint.py lines 86–137).

Initialization mirrors the oracle's $S(A) = \lbrace A, \top \rbrace$: the diagonal of $S$ is set true (every concept labels itself) and the ⊤ column is set true (everything is subsumed by ⊤), while every $R_r$ starts all-false (gpu_fixpoint.py lines 140–148). On the running example the committed run reports the whole setup:

[1] matrixization of the normalized academic TBox
    14 axioms → 16 normal forms; universe n = 12 (8 named + 2 fresh + ⊤ + ⊥); roles: advises, authored, grandAdvisor
    S: 12×12 bool (S[A,B] = "A ⊑ B"); one 12×12 adjacency R_r per role
    axiom tables: H1 nnz=6, conjunctions m=2, E3 nnz=3, E4 nnz=4, chains=1

Decode the counts. Normalization (Volume 2's own elc.normalize, imported, never retyped) turns the 14 authored axioms into 16 normal forms and invents two fresh names for the compound subconcepts; of the ten declared concepts only eight occur in the TBox, so the universe is $n = 12$: eight named concepts, two fresh names, ⊤ and ⊥. Three roles appear: advises and authored from the axioms, grandAdvisor only as a chain target. The abbreviation nnz is the standard sparse-matrix term for the number of nonzero entries: $H_1$ carries the six told subsumptions, the two conjunction axioms fill the tables for CR2, and one chain axiom (advises ∘ advises ⊑ grandAdvisor) feeds CRχ.

One operation powers everything: the boolean matrix product, written $\odot$ here. Where the ordinary matrix product sums products of numbers, the boolean product ORs together ANDs of bits. For matrices $X$ and $Y$ and any row index $i$ and column index $j$, with $k$ ranging over the $n$ inner positions,

$$(X \odot Y)[i, j] \;=\; \bigvee_{k=0}^{n-1} \; X[i, k] \,\wedge\, Y[k, j],$$

where $\bigvee_k$ is the big OR over all values of $k$ ("does at least one $k$ work?") and $\wedge$ is AND. This is the matrix product over the boolean semiring, the algebra whose addition is OR and whose multiplication is AND. The companion computes it by a deliberate trick (gpu_fixpoint.py lines 73–81):

def bmm(a: np.ndarray, b: np.ndarray) -> np.ndarray:
    """Boolean matrix product over the (∨, ∧) semiring,

        (A ⊙ B)[i, j] = ⋁_k  A[i, k] ∧ B[k, j],

    computed as a float32 BLAS matmul of 0/1 matrices followed by "> 0".
    Exact: every partial sum is an integer ≤ n ≪ 2²⁴, representable in
    float32.  Broadcasts over leading batch dimensions like np.matmul."""
    return (a.astype(np.float32) @ b.astype(np.float32)) > 0.0

The exactness claim deserves its two-line proof, because "float" and "exact" do not usually share a sentence. Cast to float32, every entry of the two matrices is $0.0$ or $1.0$, so each term $X[i,k] \cdot Y[k,j]$ is $0$ or $1$ and the ordinary matrix product computes the integer $c = \sum_k X[i,k]\, Y[k,j]$, the count of witnesses $k$ for which both bits are set. That count is at most $n$. A float32 number carries a 24-bit significand, so every integer of magnitude at most $2^{24} = 16{,}777{,}216$ is represented exactly, and adding $1$ to an integer below that bound is again exact; with $n = 12$, or indeed with $n$ in the millions, no rounding ever occurs. The final comparison $c > 0$ is true exactly when at least one witness exists, which is precisely the OR of the ANDs. The float path is not an approximation of the boolean product; it is the boolean product, routed through the one kernel (dense matrix multiply) that every accelerator on earth has spent twenty years optimizing. That routing is the entire thesis of the chapter.

The completion calculus re-rigged as linear algebra: the saturation state as boolean matrices, each rule as one masked or composed matrix product iterated to a fixpoint, and the same pass closing one hundred stacked ontology variants at once, every slice oracle-exact. Original diagram by the authors, created with AI assistance.

Every rule becomes a product

Take the rules one at a time; each gets a derivation first and its committed code second.

CR1, the label closure. The rule says: if $A'$ labels $A$ and the TBox states $A' \sqsubseteq B$, then $B$ labels $A$. Fix a row $A$ and a target column $B$ and ask when the rule fires for some witness $A'$. Renaming the witness to the inner index $k$: the rule fires exactly when there exists a $k$ with $S[A, k] = 1$ (concept $k$ labels $A$) and $H_1[k, B] = 1$ (the TBox says $k \sqsubseteq B$). Existence over a finite index set is a disjunction, so

$$\text{fires}(A, B) \;=\; \bigvee_{k} \; S[A, k] \wedge H_1[k, B] \;=\; (S \odot H_1)[A, B],$$

where the second equality is nothing but the definition of $\odot$ read backwards. One product therefore evaluates the rule for every pair $(A, B)$ and every witness $k$ simultaneously: all $n^3$ potential rule instances in one kernel call. The state update is the saturating OR-assignment $S \leftarrow S \vee (S \odot H_1)$, and the code is one line plus bookkeeping (gpu_fixpoint.py lines 179–181):

        # CR1: A' ∈ S(A), A' ⊑ B ⟹ B ∈ S(A)   —   S ← S ∨ S⊙H1.
        out = bmm(S, H1)
        matmuls += 1; derivations += int(out.sum()); S |= out

CR2, the conjunction join. The rule needs both conjuncts of an axiom $A_1 \sqcap A_2 \sqsubseteq B$ to label $A$. The axioms are held as an integer array $C_2$ of shape $(m, 2)$, where $m$ counts the conjunction axioms and row $m'$ holds the column indices of that axiom's two conjuncts, plus a projection matrix $P_2 \in \lbrace 0,1 \rbrace^{m \times n}$ with $P_2[m', B] = 1$ exactly when axiom $m'$'s head is $B$. The NumPy gather S[..., :, C2] builds a three-index array of shape $(n, m, 2)$ whose entry at $[A, m', j]$ is $S[A, \text{conjunct}_j(m')]$; reducing with .all over the last axis ANDs the two slots, giving the firing mask $\text{fires}[A, m'] = S[A, A_1^{m'}] \wedge S[A, A_2^{m'}]$, which is word for word CR2's premise "both conjuncts label $A$". A final boolean product routes each firing to its head column: $(\text{fires} \odot P_2)[A, B] = \bigvee_{m'} \text{fires}[A, m'] \wedge P_2[m', B]$. In database terms this is an AND-join of the state against the axiom table, executed as gather, reduce, multiply (gpu_fixpoint.py lines 183–187):

        # CR2: A1, A2 ∈ S(A), A1 ⊓ A2 ⊑ B ⟹ B ∈ S(A) — conjunction mask
        # (both conjunct columns true) projected onto the target columns.
        fires = S[..., :, mats["C2"]].all(axis=-1)       # (..., n, m)
        out = bmm(fires, mats["P2"])
        matmuls += 1; derivations += int(out.sum()); S |= out

CR3, the existential emission. If $A'$ labels $A$ and the TBox states $A' \sqsubseteq \exists r.B$, then the edge $(A, B)$ joins $R(r)$. With the axiom table $E^3_r$, the same disjunction argument as CR1 gives $R_r \leftarrow R_r \vee (S \odot E^3_r)$. There is a second, more physical way to read this product. Any boolean matrix is the OR of its set cells, so $E^3_r = \bigvee_{(A', B)} \mathbf{e}_{A'} \mathbf{e}_B^{\top}$, where $\mathbf{e}_j$ is the $j$-th standard basis column (all zeros except a single one in position $j$), the transpose $\top$ lays a column on its side as a row (a superscript operation on vectors that happens to share its glyph with the top concept ⊤, which it is not), and the outer product $\mathbf{e}_{A'} \mathbf{e}_B^{\top}$ is the matrix with a single one at cell $(A', B)$. Because $\odot$ distributes over $\vee$ (that distributivity is one of the semiring axioms, and it holds cell by cell because AND distributes over OR),

$$S \odot E^3_r \;=\; \bigvee_{(A', B)\,:\, A' \sqsubseteq \exists r.B} \; (S\, \mathbf{e}_{A'})\, \mathbf{e}_B^{\top},$$

and $S \mathbf{e}_{A'}$ is just column $A'$ of $S$: the set of all concepts currently labelled by $A'$. Each existential axiom stamps that whole column into column $B$ of the role matrix, an outer-product update crossing the current labels with the axiom template. The code is the CR1 pattern with the other table (gpu_fixpoint.py lines 189–192).

CR4, the two-step join. The premise chains three facts: an edge $(A, B) \in R(r)$, a label $B' \in S(B)$, and an axiom $\exists r.B' \sqsubseteq C$. Two nested existentials mean two products. First eliminate the witness $B'$: define $T = S \odot E^4_r$, so that

$$T[B, C] \;=\; \bigvee_{B'} \; S[B, B'] \wedge E^4_r[B', C],$$

which reads "some label of $B$ triggers an $r$-existential axiom landing on $C$". Then eliminate the witness $B$: $(R_r \odot T)[A, C] = \bigvee_B R_r[A, B] \wedge T[B, C]$, which reads "follow the $r$-edge from $A$ to some $B$ whose labels reach $C$". Substituting the first display into the second recovers CR4's premise exactly, with both existentials realized as the two inner disjunctions. The code charges itself two matmuls accordingly (gpu_fixpoint.py lines 194–198):

        # CR4: (A,B) ∈ R(r), B' ∈ S(B), ∃r.B' ⊑ C ⟹ C ∈ S(A)
        #   — the two-step join  S ← S ∨ R_r ⊙ (S ⊙ E4_r).
        for r, E in mats["e4"].items():
            out = bmm(R[r], bmm(S, E))
            matmuls += 2; derivations += int(out.sum()); S |= out

CR⊥, the bottom rule. An $r$-edge into an unsatisfiable concept poisons its source: if $(A, B) \in R(r)$ and $\bot \in S(B)$, then $\bot \in S(A)$. Slice the ⊥ column of $S$ as an $n \times 1$ matrix and multiply: $(R_r \odot S[:, \bot])[A] = \bigvee_B R_r[A, B] \wedge S[B, \bot]$, a matrix-vector product whose output ORs into the ⊥ column (gpu_fixpoint.py lines 200–205). This is the rule that will carry ⊥ backwards along the advises edge in the committed run below.

CRχ, the role chain. For an axiom $r \circ s \sqsubseteq t$, the rule composes edges: $(R_r \odot R_s)[A, C] = \bigvee_B R_r[A, B] \wedge R_s[B, C]$, true exactly when some $B$ completes the two-hop path. The chain axiom is literally graph composition, the same operation Volume 4 used for its chain rule, now ORed into a third relation: $R_t \leftarrow R_t \vee (R_r \odot R_s)$ (gpu_fixpoint.py lines 207–211). One table collects the whole calculus:

rule	premise (oracle form)	matrix form	code
CR1	$A' \in S(A)$, $A' \sqsubseteq B$	$S \leftarrow S \vee (S \odot H_1)$	lines 179–181
CR2	$A_1, A_2 \in S(A)$, $A_1 \sqcap A_2 \sqsubseteq B$	gather + AND, then $\vee\,(\text{fires} \odot P_2)$	lines 183–187
CR3	$A' \in S(A)$, $A' \sqsubseteq \exists r.B$	$R_r \leftarrow R_r \vee (S \odot E^3_r)$	lines 189–192
CR4	$(A,B) \in R(r)$, $B' \in S(B)$, $\exists r.B' \sqsubseteq C$	$S \leftarrow S \vee \big(R_r \odot (S \odot E^4_r)\big)$	lines 194–198
CR⊥	$(A,B) \in R(r)$, $\bot \in S(B)$	$S[:, \bot] \leftarrow S[:, \bot] \vee (R_r \odot S[:, \bot])$	lines 200–205
CRχ	$(A,B) \in R(r)$, $(B,C) \in R(s)$, $r \circ s \sqsubseteq t$	$R_t \leftarrow R_t \vee (R_r \odot R_s)$	lines 207–211

Why the boolean semiring, and why it terminates. Two properties of OR do the quiet work. The first is idempotence, $x \vee x = x$: re-deriving a known fact changes nothing, so the OR-assignments above can fire the same rule instance every wave without corrupting the state. Run the same products over the counting semiring (ordinary $+$ and $\times$) and each cell would instead count derivation paths; iterated, those counts grow without bound whenever the TBox has a cycle, and the iteration has no fixpoint of values, only of supports. Saturation wants an algebra where "known" is absorbing, and OR is exactly that. The second property is monotonicity: every update only turns cells on, never off, so the sequence of states is an ascending chain in the finite lattice of boolean matrices, $\lbrace 0,1 \rbrace^{n \times n}$ for $S$ and one more per role. A strictly ascending chain there has length at most the total number of cells, $n^2 (1 + \lvert \text{roles} \rvert)$, where $\lvert \text{roles} \rvert$ is the number of role names (the vertical bars read "the size of the set"), so some wave must change nothing, and the loop detects exactly that by comparing total set-cell counts before and after the wave (gpu_fixpoint.py lines 213–216). This is Volume 1's Kleene iteration verbatim, a monotone operator on a finite lattice climbing to its least fixpoint; only the algebra carrying each wave has changed, from set union to whole-matrix OR.

The oracle standard: exact, twice over

The committed run closes the academic TBox and reports:

[2] naive dense fixpoint: every rule × the FULL state, every wave
    waves=3  matmuls=33  derivations=71  →  closure |S|=39 cells, |R|=7 edges
    oracle equality: S and R match el_completion cell-for-cell; 8 named subsumptions == classify()
    unsatisfiable == oracle: ['TenuredStudent', 'TenuredStudentAdvisor']

Decode the numbers before trusting them. Three waves suffice, and each wave issues 11 boolean products: one for CR1, one for CR2, one CR3 product for the single role with existential-emission axioms (advises), two-product CR4 joins for each of the two roles carrying trigger axioms (advises and authored), three bottom-rule products (one per role), and one chain composition; $3 \times 11 = 33$ matmuls in all. The initial state has 23 set cells in $S$ (the 12-cell diagonal plus the 12-cell ⊤ column, overlapping at the single cell for ⊤ itself); the closure has 39, so saturation derived 16 new label cells plus the 7 role edges. The unsatisfiable pair is the calculus working end to end: TenuredStudent collapses to ⊥ through the Professor/Student disjointness via CR2, and CR⊥ then drags ⊥ backwards across TenuredStudentAdvisor's advises-edge into the advisor.

The spine of the chapter is what happens next in run(): not a spot check but set equality on everything. The matrices are read back as name-pair sets and compared against a fresh run of Volume 2's reasoner on the identical normalized axioms (gpu_fixpoint.py lines 405–408):

    pairs, roles = matrix_sets(S_nv, R_nv, mats["names"])
    o_pairs, o_roles = oracle_sets(normalized)
    assert pairs == o_pairs, "vectorized S differs from el_completion's S"
    assert roles == o_roles, "vectorized R differs from el_completion's R"

Equality here means every cell of $S$ and of every $R_r$, trivial facts, fresh names and all, not a sample of famous entailments. A second assert then rebuilds a dictionary state from the matrix, pushes it through el_completion's own reporting filters, and demands that the reader-facing classification match classify() exactly, eight named subsumptions and the two-element unsatisfiable set (gpu_fixpoint.py lines 412–422). The dense engine and the pointer-chasing engine are two implementations of one mathematical object, and the assert treats any daylight between them as a build failure.

One test ontology is still one test ontology. The real TBox was written for pedagogy: its tables are tiny ($H_1$ has six entries, there are two conjunctions and one chain), and a vectorization bug that only bites in a configuration the academic world never produces, say a conjunction head that feeds a chain, or ⊥ arriving through an existential rather than a disjointness, would sail through. So the module carries a small seeded generator, synth_tbox, which emits random 14-axiom TBoxes over 8 concepts and 2 roles, mixing told subsumptions, ⊥-headed conjunctions, existentials on both sides, and role chains, deterministically per seed (gpu_fixpoint.py lines 333–360). Three seeds are normalized, matrixized, closed by both dense engines, and held to the identical set-equality standard against the oracle (gpu_fixpoint.py lines 428–445):

[4] seeded synthetic TBoxes — exact oracle equality on each
    seed  axioms   n   |S|  |R|  waves  vs oracle
       1      14  10    28    0      3  exact
       2      14  10    21    2      2  exact
       3      14  10    21    3      2  exact

The three closures are usefully unlike the academic one (seed 1 derives a 28-cell label closure but no role edges at all; seeds 2 and 3 derive edges from different rule mixes), which is the point: coverage the showcase TBox cannot provide. Honesty about what this buys: three seeds are a smoke battery, not a proof. The correctness argument is the derivation in the previous section, rule by rule; the seeds guard the code that implements it against the gap between the mathematics and the NumPy.

Semi-naive: the delta discipline in matrix terms

The naive loop above fires every rule against the full state, every wave. Look at its committed per-wave derivation counts (a derivation event is one set cell in one rule's raw output, re-derivations of already-known facts included): 19, then 26, then 26. Wave 3 derives 26 events and adds nothing; every one is a re-derivation of a fact the state already holds, and wave 2 already re-derived most of wave 1. This is the naive tax: once a fact becomes derivable it is re-derived in every subsequent wave until the fixpoint. The classic remedy, the engine discipline that production Datalog systems parallelize [3], is semi-naive evaluation: fire a rule only through premises derived in the previous wave, joining that delta against the full state for the rule's other premise.

The matrix derivation takes four lines. Write the state entering a wave as old-plus-delta, $X = X_{\text{old}} \vee \Delta X$ and $Y = Y_{\text{old}} \vee \Delta Y$, where $\Delta X$ holds exactly the cells that turned on last wave. Because $\odot$ distributes over $\vee$ on both sides,

$$X \odot Y \;=\; (X_{\text{old}} \odot Y_{\text{old}}) \;\vee\; (\Delta X \odot Y_{\text{old}}) \;\vee\; (X_{\text{old}} \odot \Delta Y) \;\vee\; (\Delta X \odot \Delta Y).$$

The first term uses only old premises, so it was already computed, and its output already ORed in, on an earlier wave: recomputing it is pure waste. The remaining three terms each contain at least one delta operand, and they are jointly covered by the two products the semi-naive wave actually runs, because expanding each against the full state gives

$$\Delta X \odot Y \;=\; (\Delta X \odot Y_{\text{old}}) \vee (\Delta X \odot \Delta Y), \qquad X \odot \Delta Y \;=\; (X_{\text{old}} \odot \Delta Y) \vee (\Delta X \odot \Delta Y),$$

whose union is exactly the three delta terms. So firing only $\Delta X \odot Y$ and $X \odot \Delta Y$ loses no derivation with at least one new premise and skips precisely the all-old term. The companion's close_seminaive instantiates this split for every rule, and its docstring is the formula sheet (gpu_fixpoint.py lines 224–240):

        CR1   ΔS ⊙ H1                     CR3   ΔS ⊙ E3_r
        CR2   (all conjuncts ∈ S) ∧ (some conjunct ∈ ΔS), projected
        CR4   ΔR_r ⊙ (S ⊙ E4_r)  ∨  R_r ⊙ (ΔS ⊙ E4_r)
        CR⊥   ΔR_r ⊙ S[:,⊥]      ∨  R_r ⊙ ΔS[:,⊥]
        CRχ   ΔR_r ⊙ R_s         ∨  R_r ⊙ ΔR_s

CR1 and CR3 have only one state operand, so their delta form is a single smaller product; CR2 keeps the full conjunction mask but additionally requires the delta to have touched some conjunct (gpu_fixpoint.py lines 264–268); the three two-operand rules split in two. At wave's end, the next delta is what was derived minus what was known, $\Delta^{i+1} = \text{new} \wedge \neg S^{i}$, where the superscript $i$ counts waves ($S^{i}$ is the state after wave $i$, not a matrix power), and the state absorbs it (gpu_fixpoint.py lines 299–305). A rule whose delta operand is entirely false is skipped outright. The committed comparison:

[3] semi-naive (delta) iteration — same closure, fewer derivations
    engine       waves  matmuls  derivations  per-wave derivations
    naive            3       33           71   19  26  26
    semi-naive       4       49           28    9  12   7   0
    (the delta split makes MORE, smaller products, but never fires
     on all-old premises: derivations 71 → 28, the semi-naive economy RDFox parallelizes)

Read this honestly, both columns. The matmul count went up, 33 to 49, because each two-operand rule now issues two products and the delta framing takes one extra round to drain; on a 12-name toy, where every product costs microseconds, semi-naive is not a wall-clock win. What drops is the derivation-event count, 71 to 28, and that is the number that scales. The closure needed 23 genuinely new facts (16 label cells plus 7 edges); semi-naive spent 28 events reaching them, an overhead of five, while naive spent 71, re-deriving the same consequences wave after wave. In a real engine, work is proportional to derivation events, not to how many kernel launches bracket them: a worklist engine touches memory once per event, and the module asserts the economy as a semantic fact, st_sn["derivations"] < st_nv["derivations"], right next to the assert that both engines produce identical closures (gpu_fixpoint.py lines 400–402 and 424–426). The per-wave columns are also a picture worth pausing on: 9, 12, 7, 0 is a wavefront, a frontier of new consequences that expands, crests, and empties. It is the same wavefront Volume 1 drew for the closure waves of the $T_P$ operator (Volume 1's immediate-consequence operator, which maps a set of facts to everything the program $P$ derives from them in one step) on the 23-fact knowledge base, and the same one Volume 2's oracle traverses with Python sets; the delta discipline is just the observation that only the front of the wave is worth visiting.

Batching: the actual GPU thesis

Everything so far would be a curiosity if it only closed one ontology at a time; a worklist engine already does that well [2]. The move that hardware rewards is closing many at once. The companion builds $B = 100$ variants of the academic TBox, each adding one extra told subsumption between distinct named concepts, so only the $H_1$ table differs across variants; the variants become one stacked boolean tensor of shape $(100, 12, 12)$, one $H_1$ slab per variant, sharing every other axiom table (gpu_fixpoint.py lines 365–383). Here the earlier design choices pay off at once. bmm broadcasts over leading batch dimensions exactly as np.matmul does, the initial state constructor takes batch dimensions, and close_naive accepts the stacked tensor as its $H_1$ override, so every rule of every wave becomes one whole-batch operation and the entire batch closes with no Python loop over the hundred variants (gpu_fixpoint.py lines 153–169). The batch reaches its fixpoint when its slowest member does; slices that saturate early simply stop changing, and idempotence makes the extra waves harmless no-ops for them. The verification is as unforgiving as before: every one of the 100 slices is asserted equal to its own per-variant closure and to the oracle replayed on the base axioms plus that variant's added subsumption (gpu_fixpoint.py lines 447–460). The committed run:

[5] batching: B = 100 TBox variants (one extra told subsumption each) closed in
    ONE stacked (100, 12, 12) tensor pass, no Python loop over B; every slice ==
    its per-variant closure == the replayed oracle
    batched pass: waves=5, matmuls=55 (each a whole-batch op), derivations=16152

Five waves, because the hardest variant needs them; 55 matmuls, each touching all hundred worlds; 16,152 derivation events across the batch. Wall-clock is the one number this module refuses to assert, because milliseconds are machine-specific; the suite instead commits one canonical run as an artifact so the quoted table traces to something real. That committed canonical timing (best of 5, this suite's README.md):

    timing, best of 5 (measured wall-clock: varies run to run, asserted NOWHERE):
      Python loop over 100 closures :   14.75 ms
      one stacked tensor pass       :    1.79 ms   (≈8.2× faster)

Label this honestly, as the module does. It is a CPU-NumPy proxy, not a GPU measurement: an eight-fold gain from removing Python-loop overhead and letting one BLAS (Basic Linear Algebra Subprograms, the standard optimized matrix-multiply library every numeric stack calls into) call amortize across the batch, on matrices so small that a real accelerator would barely wake up. What it demonstrates is the shape of the win, not its magnitude: a fixpoint expressed as dense regular algebra batches for free, which is the property GPU hardware pays for. The distance between this proxy and a production system is a list of named engineering, all of it living downstream of the same calculus contract [1]: sparse kernels instead of dense slabs (a real ontology's $S$ is overwhelmingly zero), hash joins and columnar layouts for the rule tables, and worklist scheduling so that idle slices stop costing anything at all. The multicore anchor for what that engineering achieves is the parallel-materialization line behind RDFox, which showed how to partition a Datalog fixpoint's rule instances across cores with mostly lock-free indexes [3]; its later system report [4] measured materialization speedups up to 87-fold on a 128-core machine, sustaining millions of triples per second. Those numbers belong to that engine and that hardware, and they are quoted here as the ceiling this toy points toward, never as anything the toy achieved.

The other lane: buying speed with soundness

There is a second way to make a reasoner fast on batch hardware, and this volume keeps meeting it: stop reasoning and train a network to imitate the reasoner. The mature version of that line trains a recursive architecture on an ontology's materialized closure and then predicts entailments in a forward pass, retrained per ontology and reporting inferable-relation F1 between 0.916 and 0.999 across its benchmark ontologies [5]; the earlier report in the same line measured up to two orders of magnitude faster than the RDFox materialization engine it emulated, at F1 near 0.95 [6]. Decode the trade in systems currency. The throughput is real: a forward pass is exactly the dense batched tensor program this chapter just built, minus the fixpoint loop, so it inherits every hardware advantage at issue here. The price has two lines. First, soundness: F1 near 0.95 means roughly one derived fact in twenty is wrong or missing, with no proof attached and no signal distinguishing the good nineteen; a missed unsatisfiability is a modelling bug silently shipped, the exact failure Volume 2's ⊥-machinery exists to catch. Second, generality: the emulator is trained per ontology and answers only for distributions it saw, where this chapter's engine ran unchanged on three synthetic TBoxes it had never met, because exactness by construction transfers for free. The comparison is not a strawman; it is the field's recurring bargain, and the honest statement is that both lanes now run on the same hardware. This volume refuses the trade in the exact lane, where the whole point of the closure is that every cell is a theorem. Part VII buys the speed back on different terms: symbolic attention will put the reasoning inside an attention-shaped kernel and then pay for its guarantees with proofs and probes rather than with F1.

What the matrices cannot see

The dense representation has three blind spots, and each has a name and a fix.

Memory squares with the universe. $S$ and every $R_r$ cost $n^2$ cells. At this chapter's $n = 12$ that is 144 booleans; at SNOMED CT's scale, $n \approx 3.8 \times 10^5$ concept names, a single dense matrix is $n^2 \approx 1.4 \times 10^{11}$ cells, dead on arrival, while the actual closure is sparse (even our toy's saturated $S$ fills only 39 of 144 cells). Sparse formats and worklists are therefore not an optimization at scale; they are a feasibility condition, which is why the production engines cited above are built around them [2][3].

Boolean OR erases provenance. The naive run's 71 derivation events collapse into 39 cells; which axioms produced which cell, and along which rule firings, is gone the moment idempotence absorbs the re-derivation. The fix is Volume 4's move: evaluate the same products over a richer semiring whose elements carry explanations, provenance polynomials over axiom indicators, as Scallop's Datalog does and as Volume 2's provenance chapter set up. The cost is real: cells are no longer one bit, addition is no longer idempotent, and termination needs care again; this volume's own justifications chapter paid that cost deliberately to get minimal axiom sets.

Rules beyond these shapes need their own operators. The module vectorizes exactly CR1–CR4, CR⊥, and the chain rule, and it refuses what it does not implement: an assert rejects simple role inclusions $r \sqsubseteq s$ loudly rather than shipping a rule no committed run exercises (gpu_fixpoint.py lines 96–99). Full EL++ adds nominals, concrete domains, and range restrictions, each demanding its own matrix operator and its own correctness argument. Every new rule shape is a new kernel; that maintenance surface is precisely what a mature reasoner is. The sibling instance of this whole discipline sits one layer down in Volume 4: KLay layerizes irregular arithmetic circuits into dense per-level tensor ops for accelerators [7], the same move made here for a saturation calculus, and the GPU-native chapter told that story; this chapter is its systems kin at the level of logic rather than circuits.

The unsolved part

Three gaps stay open, and the third is open in this chapter's own numbers. First, dynamic worklists on accelerators: the semi-naive frontier is the right thing to compute (9, 12, 7, 0), but it is sparse and it shrinks, while SIMD (single instruction, multiple data) hardware, which applies one operation to many values at once, wants wide, regular, predictable work; the batched pass here papers over that mismatch by recomputing densely, which is only sane because the matrices are tiny. Nobody has a settled design for a saturation frontier that stays GPU-resident as it sparsifies. Second, incremental materialization on device: the previous chapter showed deletion is the hard direction (DRed overdeletes and re-derives), and that entire machinery, overdeletion sets, re-derivation passes, provenance tests, has no accepted batched-tensor formulation at all; retraction on a GPU is essentially unstudied. Third, the missing cost model. This chapter exhibits three regimes with committed numbers, pointer-chasing (the oracle), dense naive, and dense batched, and a fourth cited one (parallel worklists), but nothing here or in the literature predicts, from TBox statistics like $n$, table densities, and batch width, which regime wins a given workload. Our 8.2× is one machine, one shape, one batch size; a system that could choose its execution strategy per workload would need the model this chapter, like the field, does not have.

Why it matters

The series has been converging on this picture from three directions. Volume 1 defined reasoning as a monotone operator climbing a finite lattice; Volume 2 instantiated it as EL completion and proved the calculus sound and complete; Volume 4 showed that logic evaluated over semirings runs on learning hardware. This chapter is where the three meet with the receipts attached: the same closure, computed by pointer-chasing sets and by batched boolean algebra, asserted identical cell for cell, so "which hardware should reasoning use" becomes an engineering question about sparsity and scheduling rather than a question about correctness. For your own research the transferable asset is the oracle standard itself: when you vectorize, delta-filter, or otherwise re-platform an exact method, hold the fast version to set-equality against an independent slow implementation, on the real input and on adversarial synthetic ones, and let asserts, not eyeballs, own the claim. And keep the two lanes distinct in your write-ups as this chapter had to: measured toy timings, labelled and never asserted; cited system numbers, attributed and never conflated; learned emulators, priced in the soundness they spend.

Key terms

• Matrixization: compiling a normalized TBox into boolean tables ($H_1$, $E^3_r$, $E^4_r$, conjunction indices and projections) so completion rules become matrix operations over an indexed concept universe.
• Boolean matrix product ($\odot$): the matrix product over the OR-AND semiring, $(X \odot Y)[i,j] = \bigvee_k X[i,k] \wedge Y[k,j]$; computed exactly by a float matmul of 0/1 matrices followed by a greater-than-zero test.
• Saturation / fixpoint: iterating monotone, idempotent updates until a wave changes nothing; guaranteed to terminate on the finite lattice of boolean matrices.
• Naive evaluation: firing every rule against the full state every wave; correct, and maximally wasteful, since known facts are re-derived each wave.
• Semi-naive (delta) evaluation: firing each rule only through premises new in the previous wave, joined against the full state for the other premise; skips exactly the all-old rule instances by distributivity of $\odot$ over $\vee$.
• Derivation event: one set cell in one rule application's raw output, re-derivations included; the work measure that semi-naive evaluation shrinks (71 to 28 here) and that real engines' costs track.
• Batching / broadcasting: stacking $B$ problem instances into one leading tensor dimension so every rule application processes the whole batch in a single kernel call, with no per-instance loop.
• Materialization throughput: the systems currency of this Part, derived facts per second at saturation; the cited multicore engine's reported scaling (up to 87-fold at 128 cores) is its reference point.
• Neural reasoner emulation: training a network to predict a reasoner's closure; fast and batched by construction, approximate by construction, and retrained per ontology.

Where this leads

This chapter counted operations per wave and made the waves cheap; it never asked how many waves there must be. That question has a shape of its own: work (total operations) against depth (the longest chain of dependent steps), and for reasoning it is a genuine trilemma, because the tricks that shrink depth, like squaring a relation instead of extending chains link by link, explode work, and the tricks that bound work stretch depth back out. The Work-Depth Trilemma makes this exact on the running example and on synthetic chains: a 512-node relation chain closed in 510 waves of cheap steps or nine squarings (plus one confirming pass) of enormous ones, a committed 10,270-fold work ratio for a 51-fold depth cut, and the honest statement of why no schedule gets both. It is the bridge from this Part's systems ledger to Part V, where depth stops being a scheduling choice and becomes a hard ceiling baked into an architecture.

---

Companion code: examples/frontier/gpu_fixpoint.py implements the matrixization, both fixpoint engines, the seeded synthetic generator, the batched tensor pass, and every oracle assert, importing Volume 2's examples/symbolic/el_completion.py unchanged as the ground truth. Run python3 examples/frontier/gpu_fixpoint.py to reproduce every number in this chapter; apart from the labelled wall-clock rows, the output is byte-identical run to run and every claim is guarded by an assert.