Numeric kernel — one kernel, two surfaces (CONCEPT:AU-KG.compute.numeric-kernel)¶

P1 of the Analytics Program (plans/epistemic-graph-analytics_program.md, design reports/epistemic-graph-numeric-kernel-handoff.md). A slim, BLAS/LAPACK-free Rust numeric kernel that serves both Python-side array math (replacing numpy in agent-utilities) and — in later phases — in-database analytics over engine-resident data. This doc covers the shipped P1 kernel; P2–P5 are on the roadmap.

Thesis¶

epistemic-graph already does compute-near-data for vectors (batch_cosine_similarity, semantic_search run server-side in Rust, zero numpy). The Analytics Program generalizes that proven pattern into one numeric kernel (crates/eg-numeric) exposed on two surfaces.

flowchart TD
    subgraph K["crates/eg-numeric — pure Rust kernel (rlib)"]
        direction TB
        ND["ndarray 0.16<br/>arrays · reductions · element-wise"]
        FA["faer 0.20<br/>svd · eigh · solve · pinv · lstsq · qr · cholesky<br/>(NO system BLAS/LAPACK)"]
        RN["rand / rand_distr<br/>seedable RNG"]
        ERR["NumericError → LinAlgError<br/>(numpy parity)"]
    end

    K -->|"feature: python<br/>(pyo3 + rust-numpy,<br/>zero-copy + allow_threads)"| SA
    K -->|"rlib link<br/>(NO pyo3 — python feature OFF)"| SB

    subgraph SA["Surface A — in-process Python"]
        M1["epistemic_graph.numeric<br/>(extension module)"]
        M2["agent_utilities.numeric.xp<br/>(np-shim, numpy fallback) — AU-KG.compute.surface-analytics-program"]
        M1 --> M2
    end

    subgraph SB["Surface B — engine operators (P4, 🗺)"]
        D1["DataFusion SQL UDFs/UDAFs<br/>pca · covariance · zscore · svd"]
        D2["graph / vector / timeseries analytics"]
        D3["cross-modal join → PCA/cluster in-engine"]
    end

    TA["transient data<br/>(finance dataframes, ad-hoc KG math)"] --> SA
    TB["engine-resident data<br/>(embeddings, graph, columnar, timeseries)"] --> SB

    classDef done fill:#d5f5e3,stroke:#1e8449;
    classDef todo fill:#fdebd0,stroke:#b9770e;
    class K,SA done;
    class SB todo;

Decision rule (which surface): data already in the engine → Surface B (compute-near-data, no FFI). Data transient in Python (API dataframes, in-memory arrays) → Surface A (in-process). Never round-trip transient data into the DB just to compute (anti-pattern).

Crate layout & feature gating¶

crates/eg-numeric is a workspace member with crate-type = ["cdylib", "rlib"]:

rlib — the pure kernel (reductions, elementwise, linalg, random, error). No pyo3. This is what the engine links for Surface B.
cdylib — the Surface-A Python extension, built only with --features python (pulls pyo3 + numpy/rust-numpy).

Feature discipline (the Pi contract):

build	pulls eg-numeric?	faer/ndarray?	pyo3?
`--no-default-features --features server` (lib/minimal)	❌	❌	❌
default / `--features full` (→ `numeric`)	✅ (rlib)	✅	❌
`maturin ... -m crates/eg-numeric/Cargo.toml --features python`	✅ (cdylib)	✅	✅

The top-level numeric cargo feature (numeric = ["dep:eg-numeric"]) is part of the one main build (default/full), so a full-featured engine links the pure faer/ndarray kernel. A minimal --no-default-features --features server build links no eg-numeric/faer/ndarray.
pyo3 is behind eg-numeric's own python feature (off in every engine build), so a numeric engine build links the kernel rlib but no Python extension — the Plan-01 "no Python extension in the engine binary" contract holds (scripts/check_no_pyo3.sh stays green; the main wheel remains bindings = "bin").

Op-surface (P1)¶

Curated from the real audit (grep np\. over agent_utilities/) — not "all of numpy":

Reductions / stats: sum · prod · mean · var(ddof) · std(ddof) · min · max · argmin · argmax · argsort · cumsum · cumprod · percentile · quantile.
Element-wise: sqrt · log · exp · abs · tanh · clip · maximum · minimum · where · nan_to_num · isnan.
Linalg (faer): norm(+ord) · dot · matmul · solve · svd · svdvals · eigh · pinv · lstsq · qr · cholesky · det · inv · matrix_power, plus a LinAlgError exception that mirrors numpy.linalg.LinAlgError (raised on singular / non-PD).
Random: normal · uniform · integers (seedable ChaCha20; deterministic, and distribution-parity — bit-for-bit parity with numpy's PCG64 is a non-goal).

Parity¶

Every op is asserted np.allclose vs numpy on randomized inputs, with mandatory edge cases (nan/inf, singular matrices, empty arrays). P1 landed 847 parity checks, 0 failures. The corpus lives in agent-utilities (tests/test_numeric_parity.py, AU-KG.compute.surface-analytics-program) and runs against the compiled kernel when present, else the numpy fallback. A second, engine-side corpus (crates/eg-numeric/tests/test_kernel_parity.py, CONCEPT:AU-KG.compute.is-installed-kernel-discovery) tests the compiled kernel DIRECTLY (not the shim) and is what the CI gate runs — see below.

Packaging: ONE published package — `epistemic-graph[numeric]` (CONCEPT:AU-KG.compute.is-installed-kernel-discovery / AU-KG.compute.shim-goes-kernel-live)¶

The Surface-A kernel ships as package data folded into the single epistemic-graph wheel — there is NO separate eg-numeric package on PyPI. pip install epistemic-graph gets the server; pip install epistemic-graph[numeric] additionally provides the compiled kernel importable as epistemic_graph.numeric, so agent_utilities.numeric.xp runs kernel-LIVE (HAVE_KERNEL == True) instead of the numpy fallback (AU-KG.compute.shim-goes-kernel-live). The [numeric] extra only pulls numpy (for the kernel's zero-copy interop); the .so itself is self-contained (pure-Rust, BLAS/LAPACK-free).

How the fold works (release build). The kernel is still a distinct maturin target — a pyo3 cdylib built from crates/eg-numeric --features python — so the engine binary stays pyo3-free. The wheels job in .github/workflows/release-build.yml builds the node server wheel, then builds the kernel cdylib for the same target and injects its compiled .so into the server wheel as epistemic_graph/numeric.abi3.so (scripts/inject_numeric_kernel.py, which recomputes RECORD). The result is one epistemic_graph-<ver> wheel carrying both the server binary and the numeric kernel.

Dev / parity build it standalone (editable, non-publishing) with maturin — the kernel crate compiles fast and this is what the parity CI gate uses:

# On CPython ≤ 3.13 (pyo3 0.22's supported range) no flag is needed.
# On a newer interpreter (e.g. 3.14) pyo3 0.22 refuses to build unless you opt in:
export PYO3_USE_ABI3_FORWARD_COMPATIBILITY=1
maturin build --release -m crates/eg-numeric/Cargo.toml --features python --out target/wheels
pip install --no-index --find-links target/wheels eg-numeric   # dev-only local wheel

The #[pymodule] is named numeric with m.add("__kernel__", "eg-numeric"). The folded product build homes the raw extension at epistemic_graph.numeric (the xp shim's primary probe); a standalone crate build installs it as a top-level numeric package. The xp shim probes both names and checks the __kernel__ marker, so either layout activates the kernel. With the kernel present:

>>> from agent_utilities.numeric import xp, HAVE_KERNEL, KERNEL_SOURCE
>>> HAVE_KERNEL, KERNEL_SOURCE
(True, 'epistemic_graph.numeric')   # every routed 1-D/2-D float64 op now hits faer/ndarray

The kernel is a SEPARATE maturin build target from the engine facade (bindings = "bin" — the server binary, no pyo3) even though both ship in one wheel. scripts/check_no_pyo3.sh guards the facade (src/, epistemic_graph/, top-level Cargo.toml/pyproject.toml) and never scans crates/, and the Pi contract holds unchanged — cargo tree | grep -ci pyo3 = 0 (pyo3 is only ever pulled by eg-numeric's own python feature, which no engine build enables). The folded .so is named numeric.abi3.so, which the guard's epistemic_graph/epistemic_graph*.so / _epistemic_graph*.so patterns do not match, and it exists only in the built wheel, never in the source tree.

numpy-parity CI gate¶

The numeric-parity job in .github/workflows/rust-ci.yml maturin-builds the wheel, installs it, and runs crates/eg-numeric/tests/test_kernel_parity.py — a self-contained corpus that asserts every compiled-kernel op equals its numpy reference (np.allclose) across the full op-surface, including the mandatory edge cases (nan/inf, singular matrix, empty). The job FAILS CI if the Rust kernel ever diverges from numpy, which is what makes it safe to run the xp shim kernel-live in production. CI uses Python 3.12 (inside pyo3-0.22's native range, so no forward-compat flag is needed there).

Surfaces in later phases¶

Surface A migration (P2–P3): mechanical import numpy as np → from agent_utilities.numeric import xp as np across the 32 light-op files then the 6 linalg files.
Surface B (P4): the SAME rlib exposed as DataFusion UDFs/UDAFs + graph/vector/ timeseries operators, re-homing KG-resident numerics (spectral_navigator, world_model) to compute-near-data; cross-modal joins then PCA/cluster in-engine. First increment shipped — see below.
P5: drop numpy/scipy from agent-utilities; the kernel is the dep, shipped as epistemic-graph[numeric] (folded into the one engine wheel — no separate eg-numeric package). The xp shim stays so future backend swaps remain mechanical.

Surface B — in-database analytics operators (P4, CONCEPT:EG-KG.query.surface-b-numeric-operators/EG-KG.compute.l2-normalize-batch-vectors/EG-KG.query.concept-6/EG-KG.query.svd-eg-pca-column/EG-KG.query.kmeans-clustering-half-one/EG-KG.query.eg-3)¶

The pure kernel rlib is wired into the engine's query surface so analytics run where the data lives — no fetch-to-Python, no FFI. Two reach paths:

flowchart TD
    subgraph kernel["eg-numeric (rlib, feature numeric — faer + ndarray, NO pyo3)"]
      K1["linalg: dot · norm · svd · batch_l2_normalize"]
      K2["reductions: mean · std · var"]
    end
    subgraph sql["Surface B / SQL  (eg-query, feature numeric ⊃ sql)"]
      U1["cosine_sim(a,b) → Float64  (scalar)"]
      U2["l2_normalize(v) → List&lt;Float32&gt;  (scalar)"]
      U3["zscore(col) → Float64  (scalar-over-batch)"]
      U4["covariance(a,b) → Float64  (UDAF)"]
      U5["svd(vec_col) → List&lt;Float64&gt;  (UDAF, col→matrix)"]
      U6["pca(vec_col,k) → List&lt;List&lt;Float64&gt;&gt;  (UDAF, col→matrix)"]
      U7["kmeans(vec_col,k) → List&lt;Int64&gt;  (UDAF, col→matrix)"]
    end
    subgraph rpc["Surface B / Method  (src/server/handlers)"]
      M1["BatchL2Normalize { vectors }"]
    end
    SQL["SELECT zscore(price) …\nSELECT svd(emb) …\nSELECT pca(emb,3) …\nSELECT kmeans(emb,2) …"] --> U1 & U2 & U3 & U4 & U5 & U6 & U7
    U1 & U2 & U3 & U4 & U5 & U6 & U7 --> kernel
    client["client.batch_l2_normalize(vectors)"] --> M1 --> kernel
    kernel --> pi{{"Pi-contract: numeric ∉ pi\n→ no eg-numeric/faer in the pi tree"}}

SQL operators (registered on the graph-exec AND obs-tables SessionContext, gated #[cfg(feature = "numeric")] in crates/eg-query/src/sql/exec.rs::register_numeric, implemented in crates/eg-query/src/sql/numeric.rs):

Operator	Kind	Kernel	Semantics
`cosine_sim(a, b)`	scalar → `Float64`	`linalg::dot`/`norm`	`a·b/(‖a‖‖b‖)`; the raw-similarity complement to EG-115 `vector_cosine` (distance). Accepts a stored `List<Float{32,64}>` column or a `'[1,2,3]'` text literal; NULL on a dimension mismatch.
`l2_normalize(v)`	scalar → `List<Float32>`	`linalg::norm`	unit vector `v/‖v‖` (pgvector type — feeds `cosine_sim`/ANN in-query); zero-norm returned unchanged.
`zscore(col)`	scalar-over-batch → `Float64`	`reductions::mean`/`std`	standardize `(x-mean)/std` (population `ddof=0`) over the materialized batch. Exact for the engine's single-partition MemTable/`NodesTableProvider` (one batch/table); a global two-pass is the `(x-avg(x) OVER())/stddev(x) OVER()` window form.
`covariance(a, b)`	UDAF → `Float64`	`reductions::mean`	sample covariance `Σ(aᵢ-ā)(bᵢ-b̄)/(n-1)`; buffers aligned non-null pairs, merge state = two `List<Float64>` columns.
`svd(vec_col)` (EG-KG.query.svd-eg-pca-column)	UDAF → `List<Float64>`	`linalg::svdvals`	column→matrix: stacks the aggregated vector column into an `n×d` matrix (each row = one matrix row; same operand forms as `cosine_sim` — a `List<Float{32,64}>` column or `'[..]'` text) and returns its singular values (descending).
`pca(vec_col, k)` (EG-KG.query.concept-6)	UDAF → `List<List<Float64>>`	`reductions::mean` + `linalg::eigh`	column→matrix: mean-centers the `n×d` matrix, eigendecomposes the `d×d` sample covariance (`ddof=1`), and returns the top-`k` principal-component DIRECTIONS as `k` unit vectors of length `d`, descending by explained variance (sign arbitrary; `k` clamped to `d`; projected coords = `X_centered·componentsᵀ` downstream).
`kmeans(vec_col, k)` (EG-KG.query.kmeans-clustering-half-one)	UDAF → `List<Int64>`	`cluster::kmeans_labels`	column→matrix: stacks the aggregated vector column into an `n×d` matrix and returns one hard cluster label (`0..k`) per row, in ingestion order. Pure-Rust Lloyd + k-means++ (`eg-numeric`'s ChaCha20 RNG, seeded → deterministic; no linfa/BLAS); `k` clamped to `n`, empty clusters re-seeded to the farthest point.

Column→matrix marshalling (the deferred-in-P4 bridge, svd/pca/kmeans): all three are UDAFs whose Accumulator (MatrixAcc) decodes each ingested row via row_to_vector (the same operand forms cosine_sim accepts) and buffers it row-major into a flat Vec<f64> + fixed dim (ragged rows skipped, NULL-safe); evaluate() reshapes to a dense ndarray::Array2 and runs the kernel (faer for svd/pca, the cluster k-means for kmeans). Partial-aggregate state is the flat buffer as a List<Float64> plus dim (and k for pca/kmeans, via MatrixOp::takes_k) as Int64, so multi-phase grouping merges losslessly. The List<Float64> / List<List<Float64>> / List<Int64> results render as structural JSON arrays via cell_to_json (alongside the pgvector List<Float32> path).

Example — standardize a price column, rank rows by similarity, and reduce embeddings, all in-engine:

SELECT id, zscore(price) AS z FROM nodes ORDER BY z DESC;
SELECT a.id, cosine_sim(a.emb, b.emb) AS sim FROM nodes a JOIN nodes b ON a.id <> b.id;
SELECT svd(emb) AS singular_values FROM nodes;        -- singular values of the emb matrix
SELECT pca(emb, 3) AS top3_components FROM nodes;     -- top-3 principal-component directions
SELECT kmeans(emb, 4) AS cluster_labels FROM nodes;  -- 4-way clustering, one label per row

Batch Method — Method::BatchL2Normalize { vectors } (src/server/handlers/graph_ops.rs, handler #[cfg(feature = "numeric")]) L2-normalizes a batch in-engine via eg_numeric::linalg::batch_l2_normalize; the kernel-backed successor to the deprecated BatchCosineSimilarity on the same client path (client.batch_l2_normalize()).

Feature wiring & Pi-contract. eg-query gains a numeric feature (["sql", "dep:eg-numeric", "dep:ndarray"]); the engine's top-level numeric = ["dep:eg-numeric", "eg-query?/numeric"] turns the SQL operators on whenever the query surface is also built (i.e. the main build). eg-numeric's pyo3 (python) feature is off in every engine build, so a numeric engine links faer/ndarray but no Python extension. Verified: cargo tree | grep -ci pyo3 = 0.

Shipped this increment: kmeans(vec_col, k) (EG-KG.query.kmeans-clustering-half-one) — the clustering column→matrix UDAF, backed by a new pure-Rust eg-numeric::cluster k-means kernel (Lloyd + k-means++, NO linfa) — and cross-modal join→analytics (EG-KG.query.eg-3, below). Prior increment: svd/pca (EG-KG.query.svd-eg-pca-column/335). Still deferred (later P4): graph-algo/timeseries unification under the one kernel (native Method surfaces beyond SQL).

This is the capability that lets epistemic-graph surpass numpy: join graph + vector + timeseries, then run PCA / k-means / covariance over the JOINED result set IN-ENGINE. numpy has no data layer — to do this it must first fetch each modality into a separate array and align them by hand in Python. Here the join and the analytics are one SQL statement over resident data, computed where the data lives (compute-near-data, no FFI, no round-trip).

flowchart LR
    subgraph engine["epistemic-graph engine — one SQL surface"]
      direction TB
      G["graph / relational<br/>nodes(id, x, emb)"]
      V["vector<br/>emb = per-node embedding"]
      T["timeseries<br/>readings(nid, ts, reading)"]
      G -. "emb prop" .-> V
      T --> AGG["ts: AVG(reading) per node<br/>(timeseries reduction)"]
      G --> JOIN["JOIN nodes ⋈ ts  ON id = nid"]
      AGG --> JOIN
      JOIN --> AN["analytics over the joined rows<br/>kmeans(emb,k) · pca(emb,k) · covariance(x, avg_reading)"]
    end
    AN --> R["result (in-engine)<br/>clusters · components · cross-modal cov"]

One query, three modalities (from crates/eg-query/tests/cross_modal_analytics.rs, the proof test — synthetic data with hand-computed answers):

WITH readings(nid, ts, reading) AS (VALUES         -- ── timeseries modality ──
    ('n1',1,1.0),('n1',2,3.0), ('n2',1,3.0),('n2',2,5.0), ('n3',1,5.0),('n3',2,7.0),
    ('n4',1,7.0),('n4',2,9.0), ('n5',1,9.0),('n5',2,11.0),('n6',1,11.0),('n6',2,13.0)),
     ts AS (SELECT nid, avg(reading) AS avg_reading FROM readings GROUP BY nid)
SELECT kmeans(json_get(n.props, 'emb'), 2)                        AS clusters,   -- vector
       covariance(json_get_f64(n.props, 'x'), t.avg_reading)      AS xcov        -- graph×ts
FROM nodes n JOIN ts t ON n.id = t.nid;                           -- ── graph ⋈ timeseries ──

The six nodes form two communities (embeddings near [10,10] / [-10,-10], all on the line y = x); each avg_reading = 2·x by construction. The in-engine result:

column	value	why
`clusters`	`[0,0,0,1,1,1]` (2 balanced clusters of 3)	k-means over the joined `emb` vectors recovers the two communities
`xcov`	`7.0`	`cov(x, 2·x) = 2·var_sample(1..6) = 2·3.5` — a statistic spanning the graph scalar `x` and the timeseries aggregate `avg_reading`, aligned by the join
`pca(emb,1)`	`±[1/√2, 1/√2]`	all variance lies on the `y=x` diagonal, so PC1 is the diagonal direction

json_get(props,'emb') recovers each node's embedding (a JSON-array prop) as the pgvector text row_to_vector decodes — so the vector modality is a resident graph property, the graph modality is the nodes table, and the timeseries modality is the per-node AVG(reading) aggregate; the JOIN fuses them and the kernel UDAFs analyze the joined rows. The cross-modal covariance(x, avg_reading) is the sharpest "impossible in numpy" moment: two different modalities correlated in a single expression over the joined result set.