The Quadrature-TreeSHAP Kernel

01 — OrientationWhat the file is

Four functions, no dependencies beyond the STL. It computes all Shapley interaction values of every order in [min_order, max_order] for a batch of rows against a whole forest, in a single depth-first pass per row.

The paper's Algorithm 1 is one recursive DFS over one tree for one sample at one fixed interaction order s. This file expands that in four directions at once — batch of rows, forest of trees, range of orders, iterative traversal — while shrinking the per-node arithmetic. Everything else in this document is a consequence of those four expansions plus one algebraic substitution.

The quantity being computed

The paper's chain is short. For a leaf v, the empty prediction and marginal multiplier are

Everything the kernel does at runtime is an efficient evaluation of the telescoping update, Eq. 11, which the paper states for a split e on feature i with head h(e):

02 — Notation bridgeEvery paper symbol, and where it lives

The single most useful table for reading the file. The kernel does not store q, α, or p_e anywhere — it stores a rescaled form of each. Section 03 explains why.

03 — The core substitutionHot bit and cold factor

If you understand one thing about this file, make it this. Algorithm 1 computes with reciprocal cover ratios. The kernel never forms a reciprocal at all — and this is what the field names act ("hot") and c_acc ("cold factor") are about.

Fix a feature f and a node whose incoming edge splits on f. Let

• h = act[node] ∈ {0,1} — does x satisfy every split on f seen so far on this path;
• c = c_acc[node] — the product of cover ratios w_e over the same-feature chain up to and including this edge.

Then the paper's p_e = h / c and α_e = p_e − 1 = (h − c)/c. Substituting:

The same substitution collapses the Banzhaf ratio:

Why u is the safer denominator

Worth being precise, because the intuitive reading — "reciprocals of small numbers are dangerous" — is not the mechanism. Dividing by w_e is well-conditioned: IEEE division carries a relative error of at most half an ulp and cannot suffer cancellation, so 1/w_e is exactly as accurate as w_e. Three other properties do the work.

a · The bound on live_g, derived

Fix one edge. Three quantities go in, and their ranges are all we need:

• h = act[node] ∈ {0, 1} — exactly two values;
• c = c_acc[node] ∈ (0, 1] — a product of cover ratios, each ≤ 1;
• t ∈ (0, 1) — a Gauss–Legendre node, strictly interior, and fixed before the traversal starts.

Because h has only two values, the whole analysis is two cases.

Gauss–Legendre nodes are symmetric about ½, so 1 − t_n = t₁ and the two cases give the same magnitude. Combining them:

b · The contrast with the paper's intermediates

The paper computes the same number, α/(1 + αt), and it is equally bounded — that is the whole point of §03's algebra. The difference is what the intermediates do on the way. With α = (h − c)/c, a shrinking cover product sends α and 1 + αt both to infinity while their quotient stays near 1/t:

Read the last two rows together. The paper's form is accurate right up to the edge of the exponent range and then fails discontinuously — α overflows, 1 + αt overflows, and ∞/∞ is NaN. The kernel's form never approaches that edge, because t₁ holds the denominator up. In isolation this is a narrow win, confined to the last twenty decades. It matters because it composes: the same rescaling is what keeps the path accumulator in range, which is where the failures actually occur.

c · Two accumulators whose ranges cancel

This is the one that bites in practice. Algorithm 1 carries H_v(t) ∼ 1/∏w_e, which grows, and w_prod = ∏w_e, which shrinks, and multiplies them only at the leaf. Each can leave double's range while their product sits comfortably inside it — and ∞ × 0 is NaN. The rescaling multiplies them at every edge instead, so there is one accumulator, it is bounded by 1, and it is the answer. Evaluating both forms in double at the 8-node rule against a log-space reference:

Each row is one root-to-leaf path, evaluated at the top quadrature node t₈ = 0.980. "N distinct features" means a path of depth N where every edge splits on a different feature — so D = d = N in the paper's Table 1 notation, and each feature's chain is a single edge. The fourth row is the opposite extreme: one feature split ten times. The two middle columns are the intermediates that Algorithm 1 keeps apart until the leaf.

In every failing row the true value is an ordinary double — only the intermediates leave the range, and they leave it in opposite directions, which is why the product would have been fine. Relative error of the kernel form against a log-space reference stays at ≤ 4e−12 across all five rows.

d · Degrading to the limit instead of to NaN

The three properties combine into one useful behaviour. When a hot chain's cover product underflows to exactly 0, the kernel evaluates (1 − 0)/(t + 0) = 1/t — and 1/t is precisely lim_W→0 α/(1 + αt), the mathematically correct value at that limit. The formula does not merely avoid failing; it returns the right answer at the point where the alternative returns NaN. That is a consequence of h·t surviving when c does not, which is the same structural fact as the floor in (a).

The path invariant that follows

Group the path's edges by feature. Telescoping (§07) makes the net contribution of feature f to A equal to u at its last edge, which is W_f(1 + α_ft) with W_f the full chain cover product. Multiplying over features:

04 — InputsQuadTree and split semantics

QuadTree L20–56 is a flat, forest-wide struct of arrays. Two of its fields are the kernel's whole preprocessing story; the rest are ordinary tree plumbing.

goes_left L38–55

bool goes_left(double feature_value, int node) const{    if (std::isnan(feature_value))                      // missing → model's default branch        return children_left_default[node] != 0;    if (cat_size[node] > 0)                             // categorical → set membership    {        const int64_t category = static_cast<int64_t>(feature_value);        const int64_t *begin = cat_values + cat_start[node];        return std::binary_search(begin, begin + cat_size[node], category);    }    if (decision_type == Q_LESS_THAN)        return feature_value < thresholds[node];    return feature_value <= thresholds[node];}

Algorithm 1 compresses all of this into s_e ← [x satisfies edge v→u]. The kernel spends fifteen lines on it because real converted models carry NaN policy, categorical sets, and two incompatible comparison conventions. Note that the categorical branch requires cat_values to be sorted per node — binary_search is silently wrong otherwise.

05 — StateThe workspace

QuadWorkspace L65–161 is allocated once per kernel invocation and reused across every tree and every row. Nothing in the hot loop allocates.

06 — OutputFlat blocks and lexicographic rank

Algorithm 1 writes into an abstract map Φ[S]. The kernel writes into one flat double* per row, and the layout is a hard contract with the Python side.

The constructor L104–134 builds three tables:

1. binom — Pascal's triangle up to max_order, (n_feats+1) × (max_order+1). Entries with k > n are left at 0, which is the correct value.
2. order_offsets[s] — where order s's block starts. Blocks are laid out consecutively for s = min_order … max_order, each of length C(n_feats, s). Orders outside the range keep offset 0, which is why the first-order write is gated on min_order == 1 L240 — without the guard it would alias into whatever block starts at 0.
3. cum[r][v] = ∑_{u < v} C(n_feats − 1 − u, r) — the prefix sums that turn combination ranking into O(s) work instead of O(s·n_feats).

merged_rank L139–160

Given the sorted subset prefix chosen[0..size) and a pivot feature not in it, this merges the pivot into place and returns the position of the resulting tuple in itertools.combinations(range(n_feats), s) — lexicographic order.

    int s = size + 1;    int64_t rank = 0;    int prev = -1;    for (int i = 0; i < s; ++i)    {        const int64_t *cum_r = cum.data() + (size_t)(s - 1 - i) * (n_feats + 2);        rank += cum_r[merged[i]] - cum_r[prev + 1];        prev = merged[i];    }    return rank;

Each iteration counts the tuples that share the prefix merged[0..i) but have a smaller i-th element — a difference of two prefix sums rather than a loop.

Order 1 skips merged_rank entirely — the rank of the singleton {i} is i, so the write is just out[order_offsets[1] + feature] L242.

07 — Traversalquadrature_inference and the stage machine

Algorithm 1's DFS is recursive and returns a vector. quadrature_inference L257–393 is an explicit stack machine over four stages, with A and E living in depth-indexed rows instead of being passed by value.

Frame scheduling L324–336

int right = tree.children_right[node];if (left >= 0){    bool go_left = tree.goes_left(x[tree.features[node]], node);    ws.act[left]  = go_left;          // hot bit for the child edges    ws.act[right] = !go_left;    ws.stack.push_back({node,  depth,     3});   // leave      (popped 5th)    ws.stack.push_back({node,  depth,     2});   // after right (popped 4th)    ws.stack.push_back({right, depth + 1, 0});   //            (popped 3rd)    ws.stack.push_back({node,  depth,     1});   // after left  (popped 2nd)    ws.stack.push_back({left,  depth + 1, 0});   //            (popped 1st)}

The scheduling is what makes the depth-indexed rows safe: row d+1 is read at stage 1 before the right subtree is allowed to touch it, and row d−1 (the parent's A) is never written by any descendant, because descendants only write rows ≥ d.

Stage 0 — the edge update L287–323

int ancestor = tree.ancestors[node];if (node != root){    if (ancestor >= 0)        ws.act[node] = ws.act[node] && ws.act[ancestor];   // hot chain is transitive    int feature = tree.features[tree.parents[node]];    double h = ws.act[node] ? 1.0 : 0.0;    double c = tree.c_acc[node];    const double *A_prev = ws.A.data() + (size_t)(depth - 1) * n_quad;    double *g_row = ws.live_g.data() + (size_t)feature * n_quad;    if (ancestor >= 0)    {        double h0 = ws.act[ancestor] ? 1.0 : 0.0;        double c0 = tree.c_acc[ancestor];        for (int m = 0; m < n_quad; ++m)        {            double u_new = h  * t[m] + c  * (1.0 - t[m]);            double u_old = h0 * t[m] + c0 * (1.0 - t[m]);            A_row[m] = A_prev[m] * (u_new / u_old);      // telescope: divide the old factor out            g_row[m] = (h - c) / u_new;        }    }    else    {        for (int m = 0; m < n_quad; ++m)        {            double u_new = h * t[m] + c * (1.0 - t[m]);            A_row[m] = A_prev[m] * u_new;                // first edge on this feature            g_row[m] = (h - c) / u_new;        }        ws.path_feats.insert(/* sorted position */, feature);    }}

if p[f] = UNSEEN then
    p_e ← 1/w_e if s_e, else 0;   p↑ ← 1;
else if p[f] = 0 then
    p_e ← 0;   p↑ ← 0;
else
    p_e ← p[f]/w_e if s_e, else 0;   p↑ ← p[f];
α_e ← p_e − 1;
c′ ← c ⊙ (1 + α_e t);
if p[f] ≠ UNSEEN and |p[f] − 1| > ε then
    c′ ← c′ ⊘ (1 + (p[f] − 1) t);
Save p[f] and whether f ∈ A; set p[f] ← p_e; add f to A;

ancestor = ancestors[node]      // precomputed
act[node] &= act[ancestor]      // hot chain
h = act[node];  c = c_acc[node]

if ancestor >= 0:
    A[d] = A[d-1] * (u_new / u_old)
else:
    A[d] = A[d-1] * u_new
    path_feats.insert(feature)

live_g[feature] = (h - c) / u_new

Leaf and unwind L337–361, L372–390

A leaf writes E[d] = A[d]·value, calls quad_extract, then restores. An internal node does the same at stage 3, after both children have folded into E[d]. The restore is identical in both places:

• Repeated feature (ancestor >= 0): recompute live_g[feature] from the ancestor's h₀, c₀. The feature stays in path_feats — it is still live above.
• First occurrence: erase the feature from path_feats. live_g[feature] is left stale, which is harmless because membership is gated by path_feats, and any sibling subtree that splits on f again will overwrite the row at its own stage 0.

This is Algorithm 1's Restore p[f]; if f was not previously active, remove f from A — with the recomputation replacing the save.

Forest and batch L271–276, L395–418

roots[] indexes into one flat node array covering the entire ensemble, so the workspace is sized once from tree.num_nodes and tree.max_depth and every tree accumulates into the same output row. act and path_feats are cleared once per row L268–269, not per tree — safe because the DFS restores them symmetrically, and because act for any node is written by its parent before it is ever read. quadrature_tree_shap then loops rows with out_stride spacing, reusing the single workspace throughout.

08 — Extractionquad_extract

L204–255. Called once per edge, with E[depth] already holding the subtree sum. It computes Eq. 11's bracket, emits the first-order value, and hands the rest to the subset enumerator.

The cold-chain early return L214–216

int ancestor = tree.ancestors[node];if (ancestor >= 0 && !ws.act[ancestor])    return;  // both telescoping terms cancel exactly once the hot chain is broken

This has no counterpart in Algorithm 1, and it is exact rather than approximate:

δ and the first-order write L217–243

if (ancestor >= 0){    double h0 = ws.act[ancestor] ? 1.0 : 0.0;    double c0 = tree.c_acc[ancestor];    for (int m = 0; m < n_quad; ++m)    {        double u0 = h0 * t[m] + c0 * (1.0 - t[m]);        ws.delta[m] = w[m] * E_row[m] * (g_new[m] - (h0 - c0) / u0);        sum1 += ws.delta[m];    }}else{    for (int m = 0; m < n_quad; ++m)    {        ws.delta[m] = w[m] * E_row[m] * g_new[m];   // p↑ = 1 ⟹ the correction term is 0        sum1 += ws.delta[m];    }}if (ws.min_order == 1)    out[ws.order_offsets[1] + feature] += sum1;

Three things are folded into delta that Algorithm 1 keeps separate until the innermost loop: the quadrature weight w_m, the subtree sum H_u, and the telescoping difference. As a result the per-subset cost in §09 is a bare dot product of delta against gamma, with nothing else to multiply.

sum1 — accumulated in the same pass — is the order-1 value: it is exactly Eq. 11 with S = ∅, where the empty product is 1. The kernel gets it for free rather than routing the empty subset through the enumerator.

The ancestor < 0 branch drops the correction term entirely, matching Algorithm 1's p_↑ ← 1 on a feature's first appearance: (1−1)/(1+(1−1)t) = 0.

Candidate set L244–254

candidates = path_feats \ {feature}, rebuilt per edge. This is Algorithm 1's A \ {f}. Because path_feats is kept sorted and the pivot is removed by value, candidates is sorted too — which is what lets quad_enumerate generate subsets in lexicographic order and lets merged_rank assume a sorted chosen.

09 — Enumerationquad_enumerate

L166–202. One recursion over the subset lattice of candidates that emits every order from min_order to max_order, sharing the γ product between a subset and its prefixes.

for (size_t i = start; i < candidates.size(); ++i){    int j = candidates[i];    const double *g_j = ws.live_g.data() + (size_t)j * n_quad;    double *gamma_level = ws.gamma.data() + (size_t)level * n_quad;    const double *gamma_prev = (level == 0) ? nullptr                              : ws.gamma.data() + (size_t)(level - 1) * n_quad;    for (int m = 0; m < n_quad; ++m)        gamma_level[m] = (level == 0 ? 1.0 : gamma_prev[m]) * g_j[m];   // extend the prefix    chosen[level] = j;    int order = level + 2;  // subset of size level+1 plus the pivot feature    if (order >= ws.min_order && order >= 2)    {        double contribution = 0.0;        for (int m = 0; m < n_quad; ++m)            contribution += weighted[m] * gamma_level[m];   // dot product, nothing else        int64_t rank = ws.merged_rank(chosen, level + 1, feature);        out[ws.order_offsets[order] + rank] += contribution;    }    if (order < ws.max_order)        quad_enumerate(ws, candidates, i + 1, level + 1, chosen, weighted, feature, out);}

Two details worth pinning down:

• order = level + 2 — the subset has level+1 members and the pivot feature makes it level+2. The order >= 2 guard is therefore always true; it is there to keep the first-order block, which is written elsewhere, out of reach.
• Emit and descend are independent. A node emits if order >= min_order and descends if order < max_order. With min_order = 2, max_order = 4 the same descent writes into three different output blocks.

The recursion is a genuine function call per level, but the levels are bounded by max_order, not by tree depth — so the recursion depth is a small constant even on very deep trees.

10 — The comparisonFourteen departures from Algorithm 1

Everything below is a deliberate deviation, not a divergence in what is computed. Given the same tree, sample, and quadrature rule, the kernel and Algorithm 1 produce the same numbers — up to the floating-point differences the reparametrisation is designed to create.

The two that change the numbers

Departures 1 and 5 are not just faster paths to the same floating-point result — they produce different floating-point results, and that is the point. Algorithm 1 accumulates H (which grows like 1/∏w_e) and w_prod (which shrinks like ∏w_e) separately, and multiplies them only at the leaf; each can leave double's range while their product sits well inside it. The kernel multiplies them at every edge, so the running value never leaves [0,1]. The mechanism is not that a division by a small cover was avoided — that division is well-conditioned. See §03 for the measurements.

Departures 8, 9, and 10 change the complexity. Algorithm 1 as written costs O(C(|A|, s−1) · s · n) per edge for order s; the kernel costs O(∑_k C(|A|, k) · n) for all orders up to max_order together, with the cold-chain return removing edges entirely.

What is unchanged

The recursion structure (enter, descend, extract on the way out), the telescoping idea itself, the placement of the extraction at the child edge rather than the node, and the mathematical content of Eq. 11 are all carried over intact. If you read Algorithm 1 with §03's substitution in hand, the kernel is line-for-line recognisable.

11 — While readingInvariants worth trusting

Choosing n_quad

Theorem 2 of the paper bounds the required node count for exactness at order s: n ≥ ⌈(d − s + 1)/2⌉, where d is the maximum number of unique features on a root-to-leaf path. Higher orders need fewer points, so a kernel run configured for a range of orders is governed by its smallest order — min_order, not max_order. The paper's practical finding is that 8 points reach the float32 noise floor regardless; this kernel works in double throughout and takes n_quad, t, and w from the caller, so the choice is the caller's to make.

12 — ScopeWhat the paper has that this file doesn't

• The converter. ancestors and c_acc are inputs; nothing in this file computes them, and nothing in this local repository does either. They come from upstream shapiq: tree/quadrature/computer.py builds ancestors through create_edge_tree(), derives the Gauss–Legendre nodes and weights with _gauss_legendre_unit(), and calls this kernel from _explain_cpp(). Reading that file alongside this one is the fastest way to see what the kernel expects.
• Explicit SIMD. §7 of the paper describes vectorising 4 points per SSE register on CPU and 32 per warp on GPU. This file writes plain for (m = 0; m < n_quad; ++m) loops over contiguous double arrays with __restrict__ on the pointers, and leaves the vectorisation to the compiler. The layout is deliberately vector-friendly — every buffer is n_quad-contiguous — but the intrinsics are not here.
• Parallelism. No OpenMP. The row loop in quadrature_tree_shap L412 is embarrassingly parallel apart from the shared QuadWorkspace, which would need to be per-thread.
• float32. The paper's precision analysis is about the float32 output format; this kernel is double end to end.
• The GPU path. CPU only.