Exact Selective Inference for Spectral Clustering via Piecewise-Constant Graph Perturbation

Introduction

Clustering is routinely followed by a confirmatory step: having obtained clusters, the analyst asks whether two of them are truly different, for example whether their feature means differ. The naive approach, namely a two-sample test between the estimated clusters, is invalid, because the clusters were chosen by the data precisely to be separated, so the test rejects far too often even when no structure is present (the “double-dipping” problem [2]). Selective inference resolves this by conditioning on the event that the algorithm returned the clusters of interest [11, 6].

For hierarchical clustering with squared-Euclidean Ward linkage, Gao et al. [6] derived an exact selective test: under a spherical Gaussian model, the difference-in-means statistic follows a $\chi$ distribution truncated to the set of perturbations that preserve the clusters, and that set is an analytically computable union of intervals. Chen and Witten [4] extended the program to $k$-means by conditioning on the iterates of Lloyd’s algorithm. Subsequent work has extended post-clustering selective inference to convex clustering [1], dependent observations and features (González-Delgado et al., 2026), and functional or longitudinal data [20]; see also the recent review of conditional selective inference by Neufeld et al. [14]. A parallel line based on parametric programming by Le Duy and Takeuchi [10] computes the same conditioning set by line search and has recently been assembled into a framework for multi-stage clustering pipelines (Miyata et al. [12]), whose component library covers distance-based outlier detection, variance/correlation feature selection, $k$-means and DBSCAN. A method-agnostic alternative, data thinning ([13, 5]), splits the data into independent parts and avoids conditioning by testing clusters selected on an independent noisy copy. This is valid, but it targets that randomized selection copy rather than the deterministic clustering returned on the observed data, which is our target.

Conspicuously absent from all of these is spectral clustering ([18] (Ng et al. 2001)), despite its ubiquity. Graph-based representations and decompositions are also central in genomic interaction analysis, including graph-cover formulations for Hi–C data and comparative evaluation of topological-domain callers on RNA-associated interactomes ([17, 16]). Spectral clustering composes three operations: an affinity map, an eigendecomposition of a graph Laplacian, and a rounding step ($k$-means on the spectral embedding); and the eigendecomposition is the step that does not admit the polyhedral / linear-or-quadratic characterization of the selection event on which existing exact post-clustering inference relies ([6, 1, 10]). It is therefore natural to suspect that the eigendecomposition is the obstruction. Our central observation is that it is not; the affinity construction is. We make this precise and instantiate an exact-in-real-arithmetic selective test. Many spectral clustering pipelines construct an unweighted $\varepsilon$- or $k$-nearest-neighbor graph before forming a graph Laplacian; our theorem targets exactly this graph-based regime, while smooth weighted kernels require separate analysis.

Contributions.

1. Finite-state theorem. We prove that any polynomial finite-state graph selector yields a finite-union selection set for deterministic spectral clustering (§ 3).

2. Neighborhood graphs. Unweighted $\varepsilon$-graphs and OR-symmetrized $k$-NN graphs satisfy the theorem with quadratic predicates and quadratic-root breakpoints.

3. Selective tests. Under the stated known-scale spherical Gaussian model, we derive finite-sample exact truncated-$\chi$ tests in real arithmetic for selected spectral-cluster pairs, with deterministic gap, rounding, and target-selection conventions.

4. Feature-wise inference. We derive truncated-normal tests for marker or driver features after spectral clustering.

5. Direct extensions and scope. We record short extensions for known feature covariance, eigengap-selected $K$, and auxiliary-coordinate affinities, while keeping the main theorem focused on the fixed-$K$ neighborhood-graph case.

6. Empirical validation and scope. We report global- and partial-null calibration, operational power and usability, numerical-sensitivity and misspecification diagnostics, scaling behavior, and a working-model PBMC illustration (§ 5). To our

knowledge, this is the first finite-sample exact selective-inference procedure, at the theorem level, for spectral clustering with unweighted neighborhood-graph affinities.

Setup and selective test

Model. Let $\mathbf{X} \in\mathbb{R}^{n \times d}$ have independent rows $\mathbf{x}_i \sim N(\boldsymbol{\mu}_i, \sigma^2 I_d)$ with $\sigma$ known. A clustering algorithm is a map $\mathcal{C} : \mathbb{R}^{n \times d} \to\Pi_n$ to partitions of $\{1,\ldots,n\}$; write $\mathcal{C}_{\mathrm{spec}}$ for spectral clustering. For two clusters $C_a, C_b \in\mathcal{C}(\mathbf{X})$ we test

Contrast and statistic. Define $\boldsymbol{\nu} \in\mathbb{R}^n$ by $\nu_i = \mathbb{1}(i \in C_a)/|C_a| - \mathbb{1}(i \in C_b)/|C_b|$, so that $\boldsymbol{\nu}^{\mathsf{T}}\mathbf{X} = \bar{\mathbf{x}}_{C_a} - \bar{\mathbf{x}}_{C_b} \in\mathbb{R}^d$ and $\mathbb{E}[\boldsymbol{\nu}^{\mathsf{T}}\mathbf{X}] = \bar{\boldsymbol{\mu}}_{C_a} - \bar{\boldsymbol{\mu}}_{C_b}$. The test statistic is $T(\mathbf{X}) = \|\boldsymbol{\nu}^{\mathsf{T}}\mathbf{X}\|_2$. Under $H_0$, $\boldsymbol{\nu}^{\mathsf{T}}\mathbf{X} \sim N(0, \sigma^2\|\boldsymbol{\nu}\|^2 I_d)$ and hence

Here $\chi_d$ denotes the distribution of $\|G\|_2$ for $G \sim N(0,I_d)$; it is a chi, not a chi-squared, distribution. Note $\|\boldsymbol{\nu}\|^2 = 1/|C_a| + 1/|C_b|$, and $\boldsymbol{\nu}$ is supported only on $C_a \cup C_b$: points outside the tested pair (“spectators”) do not enter the contrast, though they still shape the graph.

Perturbation path. Let $\boldsymbol{\pi} := \operatorname{dir}(\boldsymbol{\nu}^{\mathsf{T}}\mathbf{X}) = \boldsymbol{\nu}^{\mathsf{T}}\mathbf{X}/T(\mathbf{X})$. Define, for $\phi\ge0$,

so that $\boldsymbol{\nu}^{\mathsf{T}}\mathbf{X}'(\phi) = \phi\boldsymbol{\pi}$, the orthogonal nuisance $\Pi_{\boldsymbol{\nu}}^\perp\mathbf{X}$ and direction $\boldsymbol{\pi}$ are preserved, and $\mathbf{X}'(T(\mathbf{X})) = \mathbf{X}$. Crucially, (3) is affine in $\phi$: row $i$ is $\mathbf{p}_i + \phi\mathbf{q}_i$ with $\mathbf{p}_i = \mathbf{x}_i - (\nu_i/\|\boldsymbol{\nu}\|^2)\boldsymbol{\nu}^{\mathsf{T}}\mathbf{X}$ and $\mathbf{q}_i = (\nu_i/\|\boldsymbol{\nu}\|^2)\boldsymbol{\pi}$; spectator rows have $\mathbf{q}_i = 0$.

Target-selection rule. The tested pair may itself be chosen after clustering, so we make that rule part of the selection event. Let $\mathcal{A}$ be a deterministic target-selection rule that maps the clustering output (and any deterministic diagnostics such as cluster sizes) to a set of unordered cluster pairs to be tested. For example, $\mathcal{A}_{\mathrm{app}}(\mathcal{C})$ returns all unordered pairs of clusters in $\mathcal{C}$, whereas $\mathcal{A}_{2\mathrm{lg}}(\mathcal{C})$ returns only the pair consisting of the two largest clusters, with a fixed tie-breaking rule. The observed target is the unordered pair $\{C_a,C_b\}$ of index sets. Once the deterministic validity conventions of § 4.1 are fixed, the matched selection event is that the spectral pipeline is valid and $\{C_a,C_b\} \in\mathcal{A}(\mathcal{C}_{\mathrm{spec}}(\mathbf{X}))$. If the reported hypothesis is simply this observed pair, $\mathcal{A}_{\mathrm{app}}$ is the Gao–Bien–Witten “pair appears” event; if the analysis protocol says to test the two largest clusters, $\mathcal{A}_{2\mathrm{lg}}$ is the correct event to condition on. More generally, any pre-specified rule that is a deterministic function of the finite graph state and pipeline output is covered without new graph roots. A rule using additional raw-data summaries requires its own selection decision to be added, and the finite-state conclusion persists only when that added decision is itself piecewise constant with finitely many tractable boundaries.

Level	Event fixed along the path	Role
Graph-conditioned	realized graph state	cheapest, most conservative
Full partition	exact spectral partition	intermediate diagnostic
Matched rule	$\{C_a,C_b\} \in\mathcal{A}(C(\cdot))$ for the reporting rule	primary test
Pair appearance	$\{C_a,C_b\}$ appears among cluster pairs	valid when the pair is user chosen

Selective $p$-value. Conditioning on the matched selection event together with the nuisances $(\Pi_{\nu}^{\perp}\mathbf{X},\pi)$ reduces the test to a one-dimensional truncated $\chi$ problem [6]. With

the selective $p$-value is

If $S_A$ is a union of intervals $\bigcup_m[\ell_m,u_m]$, (5) is a finite sum of $\chi_d$-CDF differences. The entire problem thus reduces to computing $S_A$ for the actual target-selection rule.

The spherical model $\mathbf{x}_i\sim N(\boldsymbol{\mu}_i,\sigma^2 I_d)$ can be relaxed when the feature covariance is known: one applies the entire pipeline after whitening, as summarized with the direct extensions in § 4.3.

The tractability dichotomy

Because (3) is affine, every pairwise squared distance along the path is a scalar quadratic in $\phi$:

with $\alpha_{ij}=\lVert q_i-q_j\rVert^2$, $\beta_{ij}=2(p_i-p_j)^{\mathsf{T}}(q_i-q_j)$, $\gamma_{ij}=\lVert p_i-p_j\rVert^2$. What happens next depends sharply on the affinity.

are constant on $I$, where the indicator is defined to be zero whenever the pipeline declares the graph state invalid.

Proofs of Lemma 2, Theorem 3, and Corollary 4 are given in the supplement.

Remark 5 (The eigendecomposition is not the obstruction). Theorem 3 never differentiates eigenvectors in $\phi$: on each finite-state interval the Laplacian is a fixed matrix, so its spectrum is fixed. All $\phi$-dependence is combinatorial, living in the graph selector and in the deterministic target-selection rule. This is why unweighted graphs are the tractable case.

Remark 6 (RBF affinities). For a Gaussian/RBF affinity $w_{ij}(\phi)=\exp(-d_{ij}(\phi)^2/2h^2)$, $W(\phi)$ is real-analytic rather than finite-state. Away from eigenvalue crossings and rounding ties, eigenvectors vary analytically by Rellich's theorem. The endpoints of the selection set are then zeros of analytic, generally transcendental functions, not roots of finitely many graph-state polynomials. A line-search method may localize such endpoints numerically under suitable bracketing and tie-free conditions, but the finite quadratic-root interval representation of Corollary 4 is no longer supplied by our finite-state argument.

A finite-state spectral selective test

We instantiate the construction with a single code path parameterized by the graph type (unweighted $\varepsilon$- or $k$-NN) and the Laplacian (unnormalized $L=D-W$, symmetric-normalized $L_{\mathrm{sym}}=I-D^{-1/2}WD^{-1/2}$ with row-normalized embedding à la Ng–Jordan–Weiss, or random-walk $L_{\mathrm{rw}}=I-D^{-1}W$ à la Shi–Malik), with $k$-means rounding for general $K$, testing a chosen pair of clusters. All variants share the same breakpoint set (Corollary 4) and the same gap condition below, since $L_{\mathrm{rw}}$ and $L_{\mathrm{sym}}$ have identical eigenvalues.

Deterministic spectral pipeline and degeneracy conventions

For the exact selection event, spectral clustering must be a single-valued map from a graph state to either a partition or an invalid symbol. We use the following convention throughout.

Graph construction. The $\varepsilon$-graph has $W_{ij} = 1\{d_{ij}^{2} \le\varepsilon^{2}\}$ for $i \ne j$ and zero diagonal, so the endpoint convention is closed. For the $k$-NN graph, self-neighbors are excluded and node $i$ orders candidates by the lexicographic key $(d_{ij}^{2},j)$; $N_k(i)$ is the first $k$ indices in that order. The undirected graph is the OR-symmetrization $W_{ij} = 1\{j \in N_k(i)\ \text{or}\ i \in N_k(j)\}$. Mutual $k$-NN is not used in the reported experiments.

Degrees, Laplacians, and zero rows. Let $d_i = \sum_j W_{ij}$. For normalized variants we set $d_i^{-1} = d_i^{-1/2} = 0$ when $d_i = 0$; thus isolated vertices contribute zero rows to the normalized adjacency. We then form $L = D - W$, $L_{\mathrm{sym}} = I - D^{-1/2} W D^{-1/2}$, or, for the random-walk embedding, the equivalent symmetric generalized problem obtained from $L_{\mathrm{sym}}$ and mapped back by $D^{-1/2}$. Zero-degree states are not random choices; they are handled by this convention and may still be excluded by the gap check below.

Gap and bottom-$K$ representation. A graph state is declared invalid unless $\lambda_K(L) < \lambda_{K+1}(L)$ for the chosen Laplacian eigenvalues in nondecreasing order. Repeated eigenvalues inside the bottom-$K$ subspace are allowed. To remove basis ambiguity, the mathematical embedding is defined from the bottom-$K$ spectral projector by the following deterministic procedure:

1. Scan the coordinate vectors $e_1,\ldots,e_n$ in index order.

2. Project each vector onto the bottom-$K$ invariant subspace.

3. Orthogonalize the projection against all previously accepted vectors.

4. Append the normalized residual whenever it is nonzero.

5. Choose the sign so the first nonzero coordinate is positive.

The resulting canonical orthonormal basis is a deterministic representation of the selected invariant subspace; away from degeneracies it agrees with the usual spectral-clustering target up to rounding. In floating-point experiments we use the same convention up to fixed numerical tolerances, declaring a graph invalid when the observed gap is below tolerance.

Embedding, rounding, and labels. For $L_{\mathrm{sym}}$ we row-normalize the $n \times K$ embedding; a zero spectral row remains zero. For $L_{\mathrm{rw}}$ and $L$ we use the unnormalized rows. Rounding is deterministic Lloyd $k$-means: a fixed pseudorandom seed gives the initialization, center-distance ties are assigned to the smallest center index, center updates are performed in label order, empty clusters retain their previous center, and the iteration stops when labels no longer change or after a fixed maximum number of iterations. The inferential target is expressed as an unordered pair of cluster index sets, so numeric label names have no statistical meaning; for reporting only, labels are ordered by the smallest observation index in each cluster.

Target-selection rules. $\mathcal{A}_{\mathrm{app}}$ returns every unordered pair of nonempty clusters in the partition. $\mathcal{A}_{2\mathrm{lg}}$ returns the unordered pair of the two largest nonempty clusters, breaking size ties by the smallest observation index in the cluster and then lexicographically. Other graph-output-based rules are allowed when specified before inference and made deterministic. For example, a fixed external annotation may be used to choose the pair having the largest pre-specified annotation contrast. A rule that inspects additional raw-data summaries is not automatically covered: its selection event must also be characterized and included.

Conditioning beyond the cluster pair. Two further data-defined events must be included in the conditioning set to make the procedure well-defined. They are deterministic, hence measurable, functions of the graph state; intersecting them with the selection event only removes additional values of $\phi$ from $S_A$, so the truncated-pivot validity argument is unchanged.

Under Assumptions 7–8 and the conventions in § 4.1, the conditioning event is $S_A$ with invalid graph states removed. The rule $\mathcal{A}$ must match the way the reported hypothesis was selected: the Gao–Bien–Witten pair-appearance event uses $\mathcal{A}_{\mathrm{app}}$, whereas experiments that test the two largest clusters use $\mathcal{A}_{2\mathrm{lg}}$.

Proof. [Proof sketch] On the line (3), fixing $U = (P_{\nu}^{\perp}\mathbf{X},\pi)$ leaves only the radius $T(\mathbf{X})$ free. Under $H_0$, $\nu^{\mathsf{T}}\mathbf{X} \sim N(0,\sigma^2\lVert\nu\rVert_2^2 I_d)$, and Gaussian orthogonality plus isotropy give $T(\mathbf{X}) \sim\sigma\lVert\nu\rVert\chi_d$ independent of $U$. For fixed $U$, the deterministic pipeline makes the matched selection event exactly $\{T(\mathbf{X}) \in S_A(U)\}$, so the conditional law is the same $\chi_d$ law truncated to $S_A(U)$. The ratio (5) is therefore the truncated upper-tail probability and is uniform by the probability integral transform; zero-mass truncation sets are assigned $p = 1$. Full details are in the supplement.

Feature-wise selective inference

The global test asks whether two selected clusters differ in their overall $d$-dimensional mean. In high-dimensional applications it is often more useful to ask which coordinates explain the selected separation. For feature $j \in\{1,\ldots,d\}$, consider

Under $H_{0j}$, $Z_j \sim N(0,\sigma^2\lVert\boldsymbol{\nu}\rVert^2)$. Condition on the other features $\mathbf{X}_{\cdot,-j}$ and on the nuisance projection $P_{\boldsymbol{\nu}}^{\perp}\mathbf{X}_{\cdot j}$. The one-dimensional path that varies only the tested feature is

where $e_j$ is the $j$th coordinate vector. Then $\boldsymbol{\nu}^{\top}\mathbf{X}'_j(\phi)_{\cdot j}=\phi$, all other coordinates are unchanged, and the nuisance $P_{\boldsymbol{\nu}}^{\perp}\mathbf{X}_{\cdot j}$ is preserved. Define the feature-wise selection set

Because (8) is affine in $\phi$, every pairwise squared distance along the path is again a quadratic polynomial in $\phi$. Hence Theorem 3 and Corollary 4 apply verbatim: for unweighted $\varepsilon$- and $k$-NN graphs, $S_{A,j}$ is an exact finite union of intervals whose finite endpoints are among the corresponding data-dependent quadratic roots.

Proof [Proof sketch] Fix the selected pair and condition on $U_j=(\mathbf{X}_{\cdot,-j},P_{\boldsymbol{\nu}}^{\perp}\mathbf{X}_{\cdot j})$. Under $H_{0j}$, $Z_j=\boldsymbol{\nu}^{\top}\mathbf{X}_{\cdot j}\sim N(0,\sigma^2\lVert\boldsymbol{\nu}\rVert^2)$ and is independent of $U_j$ by Gaussian orthogonality and the spherical feature model. For fixed $U_j$, the matched selection event is $\{Z_j\in S_{A,j}(U_j)\}$, and Corollary 4 makes this set a finite union of intervals for the graph classes considered here. The truncated-normal probability integral transform gives (10); one-sided tails are analogous. Full details are in the supplement.

Remark 11 (Multiple features) If all features are tested for a fixed selected pair, Proposition 10 gives valid marginal selective $p$-values for marker or driver features. The subsequent multiplicity correction is a separate layer: Bonferroni or Benjamini–Yekutieli control is valid under arbitrary dependence, while more powerful false-discovery-rate procedures require assumptions on the dependence among the selective $p$-values. If the analyst first screens features and then tests only the screened subset, that screening rule should also be included in the selection event or handled by sample splitting.

Direct extensions

The finite-state proof is stable under deterministic operations that are fixed on each interval where the neighborhood graph is fixed. We therefore keep the main construction focused on the fixed-$K$ spherical case above, and record only the application-facing extensions here. The supplement gives the full whitening and auxiliary-coordinate statements and proofs.

Data-driven number of clusters. Fix a maximum $K_{\max}$ and use the deterministic eigengap rule

with ties broken by the smallest such $k$. Because the Laplacian is fixed on each finite-state interval, the whole spectrum and hence $\widehat{K}$ are constant there.

The proof, given in the supplement, is immediate from Lemma 2 and Proposition 9: the event $\{\widehat{K}(\phi) = K_{\mathrm{obs}}\}$ is just one more measurable, path-piecewise-constant selection event. The reported experiments fix $K = 3$; (11) is the validity-preserving way to let the data set $K$.

Graphs built from auxiliary coordinates. If the graph is a deterministic function of auxiliary coordinates independent of the tested features—for example, spatial coordinates used to cluster cells before testing expression—then the graph is held fixed along the expression perturbation path. The selection event is independent of the contrast, the truncation set is the whole ray (or all of $\mathbb{R}$ for feature-wise tests), and the ordinary test is already exact. This boundary case is stated and proved in the supplement.

Algorithm 1: Finite-state computation of the selective $p$-value

Input: $\mathbf{X}$, $\sigma$, $K$, neighborhood size; observed labels and pair $(C_a,C_b)$; deterministic target-selection rule $\mathcal{A}$.

1. Form $\boldsymbol{\nu}$ and the affine path coefficients $\mathbf{p}_i,\mathbf{q}_i$.

2. Compute all graph-change roots from $\alpha_{ij},\beta_{ij},\gamma_{ij}$ in Eq. (6).

3. Set $B \leftarrow\{0\} \cup\{\text{positive finite roots returned by the solver}\} \cup\{\infty\}$; deduplicate roots by floating-point equality and sort $B$.

4. Initialize $S_{\mathcal{A}} \leftarrow\varnothing$ and $g_{\mathrm{prev}} \leftarrow\bot$.

5. For each consecutive interval $(\ell,u)$ in $B$: (a) Choose $\phi^\star$ as the midpoint of $(\ell,u)$, or $\ell+1$ when $u=\infty$.

   (b) Build $G$ from $\{d_{ij}(\phi^\star)^2\}$ and set $g \leftarrow$ the byte representation of the binary adjacency matrix.

   (c) If $g=g_{\mathrm{prev}}$, verify matrix equality and reuse the previous decision.

   (d) Otherwise, check Assumption 7. If it fails, set $\mathrm{ok} \leftarrow\mathrm{false}$; if it passes, compute the canonical embedding, run deterministic $k$-means, and set $\mathrm{ok} \leftarrow(\{C_a,C_b\} \in\mathcal{A}(\mathcal{C}(\mathbf{X}'(\phi^\star))))$.

   (e) Set $g_{\mathrm{prev}} \leftarrow g$. If ok, add $[\ell,u]$ to $S_{\mathcal{A}}$.

6. Merge adjacent accepted intervals using the fixed boundary tolerance.

7. Return the truncated-$\chi_d$ $p$-value over $S_{\mathcal{A}}$.

Algorithm. Algorithm 1 constructs the finite-interval representation of $S_{\mathcal{A}}$. Two optimizations reduce computation without changing the output: (a) spectator sparsity—under $\sigma^2 I$, $q_i \ne0$ only for $i \in C_a \cup C_b$, so only breakpoints from candidate pairs with at least one endpoint in $C_a \cup C_b$ need be enumerated; spectator–spectator distances are constant along the path; (b) graph deduplication—consecutive intervals frequently share the same graph (the swapped pair lies deep in the ordering), so the eigendecomposition and rounding are recomputed only when $G$ actually changes.

Homotopy line search. We implement Algorithm 1 as a homotopy that walks $\phi$ from $0$ to $\infty$ and recomputes the embedding and rounding only when the graph’s binary adjacency representation changes, reusing the decision across the (often many) intervals that share a graph. This is consequential for $k$-NN graphs, where most rank-swap breakpoints occur deep in the distance ordering and leave the top-$k$ graph unchanged: in our experiments the number of spectral solves was roughly $15\times$ smaller than the number of breakpoint intervals (§ 5). For $\varepsilon$-graphs every breakpoint is an edge flip, so solves equal intervals.

Floating-point conventions and numerical audit. The theorem is an exact real-arithmetic statement; the implementation is a floating-point proof of concept. In the initially submitted solver, coefficients with $|a| < 10^{-12}$ are treated as linear, negative computed discriminants are discarded rather than clipped, only positive finite roots above $10^{-12}$ are retained, and candidate roots are deduplicated by exact floating-point equality. The eigengap tolerance is $10^{-7}$, and adjacent accepted intervals are merged when their displayed boundary difference is below $10^{-9}$. Infinite intervals use $\ell+1$ as an interior evaluation point.

Graph reuse keys on the binary adjacency byte string and verifies matrix equality before reusing a decision. The revised camera-ready code also makes the theoretical lexicographic $k$-NN tie rule explicit; the audit below found no target or validity changes from that update. Endpoint membership is immaterial for the continuous reference laws.

We audited 400 fixed datasets against six alternative numerical policies. Of 2,400 attempted policy–dataset comparisons, 2,394 had an evaluable submitted target. No alternative changed overall evaluation status; among evaluable cases, neither the selected target nor eigengap validity changed. In 484 comparisons the two computed selection sets had nonzero symmetric-difference probability mass, but every difference was confined to extremely narrow intervals adjacent to selection-set boundaries; the largest such mass was $5.35 \times10^{-16}$, and the largest absolute $p$-value difference was $1.10 \times10^{-13}$. A 100-digit root/graph-state audit matched the submitted state sequence in 48/48 high-risk cases. These checks are empirical sensitivity evidence, not interval-arithmetic certification: eigenspaces and rounding remain float64, so the implementation does not provide certified upper/lower truncation probabilities. Full policies and results are in the supplement.

Conditioning level: over-conditioning vs. parametric programming. A single pass can yield valid $p$-values at several conditioning levels. The over-conditioned test fixes the realized graph, $S_{\mathrm{oc}} = \{\phi: G(\phi) = G(x^{\mathrm{obs}})\}$; the matched-rule test uses $S_A$ for the target-selection rule actually used to report $\{C_a,C_b\}$; and the pair-appearance special case is $S_{\mathrm{app}} = \{\phi: \{C_a,C_b\} \in\mathcal{A}_{\mathrm{app}}(C(X'(\phi)))\}$. Since the graph determines the clustering and the target-selection rule is a deterministic function of the clustering, all three sets are produced by the same homotopy. The less one conditions, the more power is available, but validity requires that the conditioning event match the reported hypothesis: if the protocol tests the two largest clusters, the correct power-preserving set is $S_{2\mathrm{lg}}$, not merely $S_{\mathrm{app}}$.

Experiments

Analysis protocol. Unless a robustness cell explicitly changes one component, graph type, neighborhood size, Laplacian, $K$, rounding seed, target rule, preprocessing, and tested features are fixed before inspecting the tested outcome. The primary and power experiments use OR-symmetrized $k$-NN with $k = 5$, $L_{\mathrm{sym}}$, $K = 3$, deterministic rounding, and the two-largest rule. The $\varepsilon$ robustness threshold is fixed at $1.3\sqrt{2d\sigma^2}$. Other data-dependent tuning choices must be fixed in advance, selected on independent data, or incorporated into the selection event.

Primary null calibration

We first verify type-I control under the global null (all $\mu_i$ equal, so no true clusters), where valid selective $p$-values must be uniform on $[0,1]$. The primary configuration uses the full Algorithm 1: OR-symmetrized $k$-NN graph with $k = 5$, $L_{\mathrm{sym}}$, $K = 3$, deterministic $k$-means rounding, and the matched two-largest-cluster target rule $\mathcal{A}_{2\mathrm{lg}}$. The primary setting is $n = 32,d = 5$; we retain $n = 18,d = 5$ as a small supporting configuration, with 5000 usable null replicates in each setting. A replicate is “usable” if the deterministic pipeline returns the specified target and passes the validity checks; failures are reported separately.

Setting	$N$	mean $p$	$\mathbb{P}(p \le.01)$	$\mathbb{P}(p \le.05)$	$\mathbb{P}(p \le.10)$	KS $p$
current_n18_d5	5000	0.497	0.010 [0.008, 0.014]	0.051 [0.045, 0.057]	0.096 [0.088, 0.105]	0.222
larger_n32_d5	5000	0.498	0.010 [0.008, 0.013]	0.053 [0.047, 0.059]	0.104 [0.096, 0.113]	0.911

The selective $p$-values are indistinguishable from uniform: KS $p$-values are 0.911 for $n = 32$ and 0.222 for $n = 18$, and the empirical rejection rates at levels .01, .05, and .10 all fall within their binomial confidence intervals. No observed replicate failed because of the spectral gap, zero-degree vertices, empty-cluster output, or numerical truncation-mass issues. The naive test rejected in 99.98% of the primary $n = 32$ null replicates and 95.2% of the supporting $n = 18$ replicates. For the primary $n = 32$ setting, the average number of candidate intervals was 7317, but graph deduplication reduced this to 119 spectral solves; the supporting $n = 18$ setting had 1277 intervals and 69 solves. Thus the number of interval blocks is not the number of eigendecompositions, so candidate intervals and eigendecompositions should be distinguished, although the explicit interval walk remains costly.

Laplacian variants, graph types, and conditioning level

We next exercise the unified procedure across graph type ($\varepsilon$, $k$-NN), Laplacian ($L_{\mathrm{sym}}$, $L_{\mathrm{rw}}$), and conditioning level under the global null. The supplement reports 2000 usable replicates per cell for both $n = 16$ and $n = 32$ ($d = 5$, $K = 3$). The matched rows use the same two-largest target-selection rule as the primary calibration; the graph-conditioned rows fix the realized graph for the main $k$-NN, $L_{\mathrm{sym}}$ configuration. All level .05 conditional tail-rate confidence intervals cover the nominal level. Under the target-usable stopping rule, the $n = 16$ $\varepsilon$ random-walk and symmetric cells required 2,099 and 2,101 attempts to obtain 2,000 usable replicates: usability was 0.953 and 0.952, conditional rejection was 0.050 and 0.0555, and joint usable-and-reject probability was 0.0476 and 0.0528. The extra attempts were observed-gap failures. All $k$-NN cells and all $n = 32$ cells were usable in 2,000/2,000 attempts. Because these runs stopped at 2,000 usable datasets, the attempt-based rates are diagnostics rather than pre-specified fixed-attempt binomial estimates.

Power: conditioning level matters. We next vary signal strength in three-cluster Gaussian alternatives ($k$-NN, $L_{\mathrm{sym}}$, $K = 3$, $d = 5$). The cluster centers are $0$, $\delta e_1$, and $\delta e_2$, with $\delta\in\{0,1,2,3,4,5\}$. We use 500 fixed attempts per separation for a balanced design with cluster sizes $(10,10,10)$ and a mildly unbalanced design with sizes $(6,10,14)$. Figure 2 compares the naive test, realized-graph conditioning, and matched $A_{2\mathrm{lg}}$ conditioning. All 500/500 attempts were usable in every cell, so conditional rejection and operational joint discovery coincide in these runs; the finite-sample lower 95% confidence bound on usability in each cell is 0.993. The plotted uncertainty bars are exact binomial 95% confidence intervals for rejection. At $\delta= 0$, the selective tests are near nominal while the naive test rejects almost always because the clusters were selected from the same noise. As separation grows, matched conditioning is materially more powerful than fixing the realized graph: at $\delta= 4$, matched power is 0.832 (balanced) and 0.842 (unbalanced), versus 0.418 and 0.512 for graph conditioning; at $\delta= 5$ the matched powers are 0.952 and 0.958. The supplement also reports a valid Gaussian data-thinning baseline, but it targets clusters selected from an independent noisy copy $X^{\mathrm{sel}}$ rather than the spectral clustering actually returned on $X$. In some power cells this different target is easier; for example, in Supplementary Table S7 the balanced $\delta= 3$ thinning rejection rate is 0.728, compared with 0.448 for the matched selective test.

Feature-wise partial null in structured data. We generated $n = 30$, $d = 5$ spherical Gaussian data with three balanced groups, signal in the first two coordinates, and a pre-specified fifth feature having mean zero in every true group. For separations $\delta= 0, 2, 4, 5$, all 1,000/1,000 attempts were usable. Selective level-.05 rejection was 0.050, 0.048, 0.043, 0.054, respectively, versus naive rates 0.447, 0.275, 0.059, 0.051. Median truncation masses were 0.0301, 0.101, 0.917, 1.000. Thus calibration persisted both when the selection set was restrictive and when genuine cluster structure was present; full exact confidence intervals and truncation-mass quantiles are in the supplement.

Model-scope diagnostics. In a fixed-attempt grid of 500 datasets per setting, rejection was 0.052 under the theorem-backed spherical Gaussian model. Outside the theorem, row heteroskedasticity also gave 0.052; row correlation $\rho= 0.2$ gave a conservative .05 tail rate of 0.028 but a non-uniform full distribution. For feature correlations $\rho= 0.1, 0.25, 0.5$, the rates were 0.066, 0.068, and 0.144. The submitted standardized-$t_5$ stress test rejected at 0.013 in 300 datasets. These rows are diagnostics and not theorem extensions; feature correlation can be materially anti-conservative. Pilot direct evaluations found no selection-set mismatch at 746 checked points in the new partial-null runs.

Computational scaling

Let $M$ denote candidate intervals and $U$ the distinct graph states that require spectral solves. Coefficient construction costs $O(n^{2}d)$ time and $O(n^{2})$ memory; the explicit root list has $M = O(n^{2})$ for $\varepsilon$ graphs and $M = O(n^{3})$ for $k$-NN, with $O(M \log M)$ sorting and $O(n^{2} + M)$ storage. The current code rebuilds each graph, costing $O(Mn^{2})$ for $\varepsilon$ and $O(Mn^{2} \log n)$ for $k$-NN, then performs $U$ dense $O(n^{3})$ spectral solves; $U \le M$. Testing $F$ features repeats the feature-wise homotopy approximately $F$ times.

A component diagnostic (three fixed attempts per graph/$n$ cell) confirms the current bottleneck. For $k$-NN at $n = 32, 64, 128$, mean candidate intervals were 6,890, 59,083, and 497,839, versus 107, 240, and 482 spectral solves; mean times were 0.18, 5.89, and 244.59 seconds, with the graph walk accounting for 83.9%, 97.3%, and 99.6%. For $\varepsilon$ graphs, the corresponding $n = 32, 64, 128$ means were 273, 1,186, and 4,989 intervals/solves and 0.06, 0.56, and 6.54 seconds, with eigendecomposition/rounding accounting for 87.8–94.5%. Note that these measurements characterize the recorded proof-of-concept engine, not hardware-independent constants or an optimized-complexity claim. Larger problems require skipping irrelevant rank swaps, incremental graph maintenance, and warm-started or partial eigensolvers.

Working-model PBMC algorithmic illustration

We use the PBMC 3k dataset ([19]) to illustrate execution of the feature-wise homotopy on a fixed preprocessed matrix; this is not to support confirmatory biological inference. The protocol fixes a pre-specified 16-gene panel, log-normalized library-size-adjusted counts, marginal gene standardization, a fixed-seed subsample of 60 cells, $k = 5$, $K = 3$, and the matched two-largest rule. The observed cluster sizes are (31, 21, 8), and the 16 homotopies finish in 55.5 seconds.

The calculation sets the working scale to $\sigma= 1$ after standardization; this is not a known sampling variance, and marginal standardization does not remove gene–gene correlation. The misspecification grid above shows that feature correlation can be strongly anti-conservative. Accordingly, the displayed quantities are working-model selective values and not valid confirmatory scientific $p$-values under a defensible PBMC sampling model. We make no marker significance or biological conclusion. The fixed-panel table and diagnostic plots are moved to the supplement, where they are labeled with the same caveat.

Discussion

Conditioning and power. The same graph-state walk yields several valid conditioning levels: fixing the realized graph is cheapest but most conservative, while conditioning only on the matched target-selection rule $\{C_a, C_b\} \in\mathcal{A}(C(\cdot))$ preserves the reported hypothesis with more power. Figure 2 illustrates the gap: at separation $\delta= 4$, matched-rule conditioning gains 33–41 percentage points over graph conditioning without extra graph enumeration. This is the same conditioning/power trade-off emphasized in parametric-programming selective inference and related post-selection problems.

Different targets and pipeline use. Data thinning is a valid baseline, but it clusters a randomized copy $X^{\mathrm{sel}}$ and tests on an independent copy $X^{\mathrm{inf}}$; our selective test instead makes inference for the deterministic spectral clustering returned on the observed matrix $X$. Theorem 3 and Corollary 4 also make the unweighted-graph spectral step a natural component for multi-stage clustering-pipeline selective inference, where linear/quadratic selection events are already handled for steps such as feature selection, $k$-means, and DBSCAN.

Assumptions, numerical status, and limitations. The finite-state graph reduction does not require known $\sigma$, but the stated truncated-$\chi$ and truncated-normal reference laws do. An independent or same-data scale estimate cannot simply be plugged in and called exact; an exact studentized construction would require a new conditional-law derivation. When a known covariance $\Sigma$ is available, validity applies only if the analysis itself is defined by whitening first and clustering/testing in the resulting Mahalanobis metric; it does not validate the Euclidean clusters that would have been obtained before whitening, and per-coordinate pivots under general or unknown covariance remain open. Fixed preprocessing and deterministic rounding are required. Data-driven $K$ under Proposition 13 is covered, but feature selection, graph type, neighborhood size, normalization, and other tuning choices must be fixed, independently selected, or added to the selection event. Smooth kernels break the finite quadratic-root representation. Finally, the numerical audit supports stability on the audited datasets but is not end-to-end certification, and the current implementation is a proof of concept for the studied problem sizes rather than a practical large-scale system.

Conclusion

Exact selective inference for spectral clustering is tractable at the statistical-theorem level, and the eigendecomposition is not the obstruction. For unweighted neighborhood graphs, clustering is piecewise constant along the selective perturbation, with finite breakpoints among roots of data-dependent quadratics. This yields finite-sample exact truncated pivots in real arithmetic under the stated known-scale spherical Gaussian model. The same graph-state argument gives feature-wise tests and supplies the previously missing spectral component for clustering-pipeline selective inference. Smooth kernels, unknown/correlated covariance, adaptive preprocessing, scalable implementation, and end-to-end numerical certification remain outside the present guarantee.

Supplementary Material

Expanded Proofs and Technical Details

Proof of Theorem 3 and Corollary 4

Proof [Proof of Lemma 2] On an interval where $G(\phi)=G_0$, every downstream operation is evaluated at the same graph state. The Laplacian construction, gap/validity check, canonical embedding map, rounding rule, label-set representation, and target-selection rule are fixed measurable functions of that state. Their composition is therefore a single constant value on the interval.

Proof [Proof of Theorem 3] Let

Because there are finitely many nonzero polynomials and each has finite degree, $B$ is finite. On any connected component $I$ of $[0,\infty)\setminus B$, every sign $\operatorname{sign}h_m(\phi)$ is constant, so the graph selector returns a fixed graph $G_I$ throughout $I$. By the deterministic spectral pipeline and the strict-gap convention, the graph-to-label map either returns a single partition or the invalid symbol. Lemma 2 therefore makes the target decision $\mathbb{1}\{\{C_a,C_b\}\in\mathcal{A}(C(X'(\phi)))\}$ constant on $I$, with invalid graph states contributing decision zero. Hence $S_A$ is the union of those components for which the constant decision is one, up to finitely many endpoints in $B$. These endpoints have no effect on the truncated $\chi$ or truncated normal pivots because the reference laws are continuous.

If the coefficients of the $h_m$'s are algebraic, each finite endpoint is a root of a polynomial with algebraic coefficients. Thus every finite endpoint is an algebraic number.

Proof [Proof of Corollary 4] Along the selective path (3), row $i$ can be written as

For $i\ne j$,

where

Thus every pairwise squared distance is a polynomial of degree at most two.

For an $\varepsilon$-graph, the edge $\{i,j\}$ is determined by the predicate $d_{ij}(\phi)^2\le\varepsilon^2$, whose boundary is the quadratic equation $d_{ij}(\phi)^2-\varepsilon^2=0$. Hence the graph can change only at roots of these $O(n^2)$ quadratic equations.

For an OR-symmetrized $k$-NN graph, node $i$ orders candidates by $(d_{ij}(\phi)^2,j)$, with the index tie-breaker fixed in advance. Away from ties, the top-$k$ set is determined by the order predicates $d_{ij}(\phi)^2<d_{im}(\phi)^2$ for $j,m\ne i$. Each boundary is

again a polynomial equation of degree at most two. There are $O(n^3)$ such candidate order-swap equations. Between consecutive candidate roots all edge predicates, order predicates, and hence the graph are constant. Theorem 3 then gives the finite union of intervals. The finite endpoints are roots of degree-$\le2$ polynomials with coefficients determined by the data.

Remark S1 (Degenerate quadratic coefficients) If $\alpha_{ij}=0$, the corresponding distance or order boundary is linear or constant. Constant nonzero predicates contribute no breakpoint. Identically zero boundaries correspond to exact distance ties on the path; the deterministic endpoint and lexicographic tie conventions assign these cases to a fixed graph state, and the boundary itself is null for the continuous pivot laws.

Proof of Proposition 9

Proof Let $T(\mathbf{X})=\|\nu^\top\mathbf{X}\|$ and $\pi=\nu^\top\mathbf{X}/T(\mathbf{X})$ on $\{T>0\}$, and let

The data decompose as

Conversely, replacing $T(\mathbf{X})$ by any $\phi\ge0$ gives exactly the path

so $(U,\phi)$ reconstructs the data matrix on the selective line. Under the null hypothesis,

For each feature, $Z$ is uncorrelated with $P_\nu^\perp\mathbf{X}$, and joint Gaussianity implies independence. Since $Z$ is isotropic under $H_0$, its radius and direction are independent, with

and $T(\mathbf{X})$ independent of $U=(P_\nu^\perp\mathbf{X},\pi)$. By construction of the deterministic pipeline, for fixed $U$ the matched selection event is exactly

where invalid graph states are excluded from $S_A(U)$. Therefore, conditional on $U$ and on selection, $T(\mathbf{X})$ follows a $\sigma\lVert\nu\rVert\chi_d$ law truncated to $S_A(U)$. The $p$-value (5) is the upper-tail probability under that truncated law. By the probability integral transform it is uniform conditional on $U$ and $E_A$ whenever the truncation set has positive measure. If the truncation mass is zero we set the $p$-value to one, which preserves super-uniformity. Averaging over $U$ conditional on selection yields

Proof of Proposition 10

Proof. Fix the selected pair $\{C_a,C_b\}$ and its contrast vector $\nu$. For feature $j$, write

The tested column has the orthogonal decomposition

Replacing $Z_j$ by a candidate scalar $\phi\in\mathbb{R}$ while holding $U_j$ fixed gives the feature-wise path (8). Thus $(U_j,\phi)$ is the one-dimensional reparametrization used to evaluate the selection event for feature $j$.

Under $H_{0j}$,

Moreover,

so Gaussianity implies $Z_j \perp P^\perp_\nu\mathbf{X}_{\cdot j}$. The spherical model also makes feature $j$ independent of $\mathbf{X}_{\cdot,-j}$, hence $Z_j \perp U_j$.

For fixed $U_j$, define

Because $\mathbf{X}=\mathbf{X}'(Z_j)$, the matched feature-wise selection event is $\{Z_j\in S_{A,j}(U_j)\}$. The finite-state graph argument applies to this path as well because every pairwise squared distance remains quadratic in $\phi$; hence $S_{A,j}$ is a finite union of intervals for the graph classes covered by Corollary 4.

Conditional on $U_j$ and on selection, $Z_j$ therefore follows $N(0,s_j^2)$ truncated to $S_{A,j}(U_j)$. For the two-sided statistic, the null distribution of $\lvert Z_j\rvert$ induced by this truncated normal law gives the selective tail probability (10). Applying the probability integral transform to that induced distribution gives a uniform conditional pivot whenever the truncation mass is positive; zero-mass cases are assigned $p_j=1$. Averaging over the selected nuisance distribution proves super-uniformity. One-sided pivots follow by replacing the two-sided absolute tail with the corresponding upper or lower truncated-normal tail.

Direct extension statements and proofs

Proof. Whitening is a fixed invertible linear map on features. If $\mathbf{x}_i \sim N(\boldsymbol{\mu}_i,\boldsymbol{\Sigma})$ independently, then $\widetilde{\mathbf{x}}_i = \boldsymbol{\Sigma}^{-1/2}\mathbf{x}_i \sim N(\boldsymbol{\Sigma}^{-1/2}\boldsymbol{\mu}_i,\boldsymbol{\Sigma}^{-1/2}\boldsymbol{\Sigma}\boldsymbol{\Sigma}^{-1/2}) = N(\boldsymbol{\Sigma}^{-1/2}\boldsymbol{\mu}_i,\mathbf{I}_d)$, so the rows of $\widetilde{\mathbf{X}} = \mathbf{X}\boldsymbol{\Sigma}^{-1/2}$ obey the unit-scale spherical model of the main text. Since $\boldsymbol{\Sigma}^{-1/2}$ is invertible, $\bar{\boldsymbol{\mu}}_{C_a}-\bar{\boldsymbol{\mu}}_{C_b}=0$ iff $(\bar{\boldsymbol{\mu}}_{C_a}-\bar{\boldsymbol{\mu}}_{C_b})\boldsymbol{\Sigma}^{-1/2}=\bar{\widetilde{\boldsymbol{\mu}}}_{C_a}-\bar{\widetilde{\boldsymbol{\mu}}}_{C_b}=0$, so $H_0$ is unchanged. Running the pipeline on $\widetilde{\mathbf{X}}$ means the graph is built from $\lVert\widetilde{\mathbf{x}}_i-\widetilde{\mathbf{x}}_j\rVert^2=(\mathbf{x}_i-\mathbf{x}_j)^{\mathsf{T}}\boldsymbol{\Sigma}^{-1}(\mathbf{x}_i-\mathbf{x}_j)$ and the statistic is $T(\widetilde{\mathbf{X}})=\lVert\boldsymbol{\nu}^{\mathsf{T}}\widetilde{\mathbf{X}}\rVert_2=\lVert\boldsymbol{\nu}^{\mathsf{T}}\mathbf{X}\rVert_{\boldsymbol{\Sigma}^{-1}}$. The selective path for $\widetilde{\mathbf{X}}$, $\widetilde{\mathbf{X}}'(\phi)=\widetilde{\mathbf{X}}+\frac{\phi-T(\widetilde{\mathbf{X}})}{\lVert\boldsymbol{\nu}\rVert^2}\boldsymbol{\nu}\widetilde{\boldsymbol{\pi}}^{\mathsf{T}}$ with $\widetilde{\boldsymbol{\pi}}=\boldsymbol{\nu}^{\mathsf{T}}\widetilde{\mathbf{X}}/T(\widetilde{\mathbf{X}})$, is affine in $\phi$, so every pairwise squared distance is quadratic in $\phi$ exactly as in (6). Corollary 4 therefore holds verbatim for $\widetilde{\mathbf{X}}$, and Proposition 9 applied to $\widetilde{\mathbf{X}}$ with $\sigma=1$ gives the exact selective $p$-value with reference law $\lVert\boldsymbol{\nu}\rVert\chi_d$ truncated to $S_A$.

Remark S2 (Why graph whitening is required). Proposition 9 requires the truncation set and the test statistic to be measurable functions of the same one-dimensional reduction $(U,\phi)$. Whitening only the statistic, while building a Euclidean graph on the unwhitened $\mathbf{X}$, breaks this: the graph breakpoints are then quadratics in the Euclidean distances of $\mathbf{X}$, whereas the reference law lives in the Mahalanobis metric, and the two no longer share the affine line. Clustering the whitened data restores a single metric and hence the exact reduction.

Proof [Proof of Proposition 13] On any connected component $I$ of $[0,\infty)\setminus B$ (with $B$ the finite breakpoint set of Theorem 3) the graph, hence the Laplacian $L$, is a fixed matrix; therefore its ordered spectrum $\lambda_1(L)\le\cdots\le\lambda_n(L)$ is a fixed tuple and the eigengap maximiser $\hat{K}$ of (11), with its fixed tie-break, is a constant on $I$. Augment the target decision of Lemma 2 by the factor $\mathbf{1}\{\hat{K}(\phi)=K_{\mathrm{obs}}\}$; being a deterministic function of the frozen spectrum, it is constant on $I$, so the augmented decision is constant on $I$. Theorem 3 then makes $S_A^K$ a finite union of intervals whose finite endpoints lie among roots of the data-dependent polynomial predicates. Validity is Proposition 9 verbatim: the event $\{\hat{K}(\phi)=K_{\mathrm{obs}}\}$ is one more measurable selection event that is piecewise constant along the path, so intersecting it into $E_A$ leaves $T(\mathbf{X})\sim\sigma\lVert\boldsymbol{\nu}\rVert\chi_d$ conditionally on $U$ truncated to $S_A^K(U)$, and the probability integral transform gives super-uniformity (with $p=1$ on zero truncation mass). The strict-gap requirement $\lambda_{\hat{K}+1}>\lambda_{\hat{K}}$ holds automatically on any interval where the maximal gap in (11) is positive.

perturbation path, or equivalently that the graph metric ignores the perturbed directions. Then the graph is constant along the path. If the auxiliary coordinates are independent of the tested feature matrix $\mathbf{X}$, the selected contrast is independent of the test statistic, the truncation set is vacuous, and the ordinary global and feature-wise tests are exactly valid.

Proof If the graph metric ignores the perturbed directions then, in the notation of Corollary 4, every distance or order predicate entering the selector has $q_i-q_j=0$, so its defining quadratic (6) is the constant $\gamma_{ij}$ and contributes no breakpoint. Hence $G(\phi)\equiv G$ and the finite-state selector has no roots on $[0,\infty)$: the single component $[0,\infty)$ carries one graph state. If that state is valid and $\{C_a,C_b\}\in\mathcal{A}(C_{\mathrm{spec}}(G))$ then $S_A=[0,\infty)$, so the truncation in (5) is vacuous and the pivot is the untruncated $\sigma\lVert\nu\rVert\chi_d$ tail, i.e. the ordinary test; the identical argument on the feature-wise path (8) gives $S_{A,j}=\mathbb{R}$ and the ordinary $N(0,\sigma^2\lVert\nu\rVert^2)$ test. When the auxiliary coordinates defining the graph are independent of $\mathbf{X}$, the induced $\nu$ (a function of those coordinates) is independent of $\nu^{\mathsf{T}}\mathbf{X}$, so $T(\mathbf{X})\sim\sigma\lVert\nu\rVert\chi_d$ holds conditionally on the selection and the naive test is exactly valid. ▪

Primary Null Calibration Diagnostics

Figure S1 shows the smaller primary global-null design that was displayed in the submitted main paper and is moved here to make room for the camera-ready additions. The experiment used 5,000 fixed attempts; all 5,000 were usable. The matched and graph-conditioned selective $p$-values are close to uniform, whereas the ordinary post-clustering test is severely anti-conservative.

Floating-Point Numerical-Sensitivity Audit

The finite-state characterization and pivots in the main paper are exact in real arithmetic under the stated model. The implementation is a floating-point proof of concept and is not an interval-arithmetic certification. The originally submitted numerical policy treats $|a| < 10^{-12}$ as linear, discards negative computed discriminants, retains positive finite roots above $10^{-12}$, deduplicates candidate roots by exact floating-point equality, uses an eigengap tolerance of $10^{-7}$, and merges adjacent accepted intervals at a boundary tolerance of $10^{-9}$.

We audited 400 fixed datasets, 50 in each of eight strata: $k$-NN global-null settings at $n = 18$ and 32, $\varepsilon$-graph global-null settings at $n = 16$ and 32, and balanced or mildly unbalanced power settings at separations 3 and 5. Each dataset was compared with six alternative policies: a stable quadratic formula, long-double root arithmetic, relaxed root handling, eigengap tolerances $10^{-10}$ and $10^{-6}$, and explicit deterministic $k$-NN tie handling. The camera-ready code adopts the explicit theoretical lexicographic tie rule. One $n = 16$ $\varepsilon$-graph dataset failed the observed eigengap check under the submitted policy and all alternatives; this leaves 399 evaluable datasets per policy and 2,394 evaluable comparisons overall.

Alternative policy	Attempted	Evaluable	Set differences	Discrete changes	Max. $|\Delta p|$	Max. symmetric mass
Stable quadratic formula	400	399	121	0	1.10e-13	5.35e-16
Long-double roots	400	399	121	0	1.10e-13	5.35e-16
Relaxed root merging	400	399	121	0	1.10e-13	5.35e-16
Eigengap tol. $10^{-10}$	400	399	0	0	0.00e+00	0.00e+00
Eigengap tol. $10^{-6}$	400	399	0	0	0.00e+00	0.00e+00
Explicit deterministic ties	400	399	121	0	1.10e-13	5.35e-16

Across the 2,394 evaluable comparisons, no alternative policy changed the overall evaluation status, selected target, or eigengap validity. There were 484 comparisons with numerically nonzero symmetric-difference mass; every such difference was confined to extremely narrow intervals adjacent to selection-set boundaries. The maximum probability mass under the reference law was $5.35 \times10^{-16}$, and the maximum absolute $p$-value difference was $1.10 \times10^{-13}$. The only root-count change occurred for one balanced separation-5 dataset under relaxed root merging (5,040 to 5,039 roots); it merged a redundant root and changed neither the selection set’s probability mass nor the $p$-value.

We additionally selected 48 risk cases using small root spacing, discriminants, graph margins, eigengaps, or truncation masses and recomputed roots and compressed graph-state sequences at 100 decimal digits. The submitted and long-double state sequences both matched the high-precision sequence in 48/48 cases. This high-precision check covers root and graph-state enumeration only: spectral eigendecomposition and rounding remain float64 and are not interval-enclosed. Consequently, the audit supplies empirical sensitivity evidence rather than certified upper and lower truncation probabilities.

Operational Usability and Power

The submitted power files contain fixed attempts rather than only successful replicates. Table S2 therefore reports the requested end-to-end quantities. Every cell had 500/500 usable attempts, so the conditional rejection probability and the operational joint probability agree in this grid. The exact finite-sample lower 95% confidence bound on usability in each cell is 0.993; this finite Monte Carlo result is not a universal usability guarantee.

Design	$\delta$	Attempts	Usable	Pr(reject | usable)	Pr(usable $\cap$ reject)
Balanced	0	500	500/500	0.046 [0.029, 0.068]	0.046 [0.029, 0.068]
	1	500	500/500	0.052 [0.034, 0.075]	0.052 [0.034, 0.075]
	2	500	500/500	0.158 [0.127, 0.193]	0.158 [0.127, 0.193]
	3	500	500/500	0.448 [0.404, 0.493]	0.448 [0.404, 0.493]
	4	500	500/500	0.832 [0.796, 0.864]	0.832 [0.796, 0.864]
	5	500	500/500	0.952 [0.929, 0.969]	0.952 [0.929, 0.969]
Mildly unbalanced	0	500	500/500	0.054 [0.036, 0.078]	0.054 [0.036, 0.078]
	1	500	500/500	0.058 [0.039, 0.082]	0.058 [0.039, 0.082]
	2	500	500/500	0.172 [0.140, 0.208]	0.172 [0.140, 0.208]
	3	500	500/500	0.486 [0.441, 0.531]	0.486 [0.441, 0.531]
	4	500	500/500	0.842 [0.807, 0.873]	0.842 [0.807, 0.873]
	5	500	500/500	0.958 [0.937, 0.974]	0.958 [0.937, 0.974]

Table S3 gives the analogous attempt-based quantities for the submitted graph/Laplacian robustness file. These runs used a target-usable stopping rule, so the attempt totals are diagnostic denominators rather than a prospectively fixed sample size. The only observed failures were in the $n = 16$ $\varepsilon$-graph cells: 99 eigengap failures before 2,000 usable random-walk cases and 101 before 2,000 usable symmetric-normalized cases. All other displayed cells were 2,000/2,000 usable.

$n$	Graph	Laplacian	Cond.	Attempts	Pr(usable)	Pr(reject $|$ usable)	Pr(usable $\cap$ reject)
16	$\varepsilon$	$L_{\mathrm{rw}}$	matched	2099	0.953 [0.943, 0.962]	0.050 [0.041, 0.060]	0.048 [0.039, 0.058]
16	$\varepsilon$	$L_{\mathrm{sym}}$	matched	2101	0.952 [0.942, 0.961]	0.056 [0.046, 0.066]	0.053 [0.044, 0.063]
16	$k$-NN	$L_{\mathrm{rw}}$	matched	2000	1.000 [0.998, 1.000]	0.048 [0.039, 0.058]	0.048 [0.039, 0.058]
16	$k$-NN	$L_{\mathrm{sym}}$	graph	2000	1.000 [0.998, 1.000]	0.044 [0.036, 0.054]	0.044 [0.036, 0.054]
16	$k$-NN	$L_{\mathrm{sym}}$	matched	2000	1.000 [0.998, 1.000]	0.048 [0.039, 0.058]	0.048 [0.039, 0.058]
32	$\varepsilon$	$L_{\mathrm{rw}}$	matched	2000	1.000 [0.998, 1.000]	0.050 [0.041, 0.060]	0.050 [0.041, 0.060]
32	$\varepsilon$	$L_{\mathrm{sym}}$	matched	2000	1.000 [0.998, 1.000]	0.054 [0.045, 0.065]	0.054 [0.045, 0.065]
32	$k$-NN	$L_{\mathrm{rw}}$	matched	2000	1.000 [0.998, 1.000]	0.053 [0.044, 0.064]	0.053 [0.044, 0.064]
32	$k$-NN	$L_{\mathrm{sym}}$	graph	2000	1.000 [0.998, 1.000]	0.045 [0.037, 0.056]	0.045 [0.037, 0.056]
32	$k$-NN	$L_{\mathrm{sym}}$	matched	2000	1.000 [0.998, 1.000]	0.043 [0.035, 0.053]	0.043 [0.035, 0.053]

Feature-Wise Partial-Null Calibration with True Structure

To complement the global-null experiments, we generated spherical Gaussian data with $n = 30$, $d = 5$ and three balanced true groups. The first two coordinates carried group separation $\delta$, while the fifth coordinate was pre-specified for testing and had mean zero in every true group. We applied the feature-wise matched-rule test to that null coordinate after spectral clustering. All 4,000 attempts were usable. Table S4 shows calibration across regimes ranging from restrictive truncation at $\delta= 0$ to nearly vacuous truncation at $\delta= 5$.

$\delta$	Usable	Selective rejection	Operational joint	Naive rejection	Median truncation mass
0	1000/1000	0.050 [0.037, 0.065]	0.050 [0.037, 0.065]	0.447 [0.416, 0.478]	0.03007
2	1000/1000	0.048 [0.036, 0.063]	0.048 [0.036, 0.063]	0.275 [0.248, 0.304]	0.101
4	1000/1000	0.043 [0.031, 0.057]	0.043 [0.031, 0.057]	0.059 [0.045, 0.075]	0.9171
5	1000/1000	0.054 [0.041, 0.070]	0.054 [0.041, 0.070]	0.051 [0.038, 0.067]	1

Selective rejection rates were 0.050, 0.048, 0.043, and 0.054 at $\delta= 0, 2, 4, 5$, respectively. The corresponding naive rates were 0.447, 0.275, 0.059, and 0.051. A stable interval-mass helper and the submitted direct formula differed by at most $3.70 \times10^{-10}$ over the 4,000 attempts. A pilot direct evaluation of the selection rule at 746 non-boundary points in 10 successful datasets found no mismatch with the enumerated intervals.

Target-Rule Equivalence Sanity Check for the Reported $K = 3$ Design

Table S5 reports a small null diagnostic for the main $k$-NN, $L_{\mathrm{sym}}$, $K = 3$ configuration. The reported target is the pair of the two largest clusters. This table is not intended to show a power difference among the first three conditioning events. In this $K = 3$ setting, once the two exact selected cluster index sets are fixed, the pair-appearance and full-partition events coincide with the matched two-largest event: the remaining cluster is the complement, and the selected pair retains the same sizes. Thus the first three rows are expected to agree. The graph-conditioned row is the stricter over-conditioned event used in the main paper’s power comparison. Larger-$K$ ablations can be generated with run_target_rule_comparison.py.

Conditioning event	$N$	mean $p$	$\mathbb{P}(p \le.01)$	$\mathbb{P}(p \le.05)$	$\mathbb{P}(p \le.10)$
Pair appears	500	0.526	0.012	0.046	0.104
Two largest	500	0.526	0.012	0.046	0.104
Full partition	500	0.526	0.012	0.046	0.104
Realized graph	500	0.518	0.008	0.042	0.096

Robustness Matrix

Table S6 reports the graph-type, Laplacian, and conditioning-level robustness matrix summarized in the main paper. The matched rows condition on the two-largest-cluster target-selection rule; graph rows fix the realized graph for the main $k$-NN, $L_{\mathrm{sym}}$ configuration.

$n$	Graph	Laplacian	Cond.	$N$	mean $p$	$\mathbb{P}(p \le.05)$	KS $p$
16	$k$-NN	$L_{\mathrm{sym}}$	Matched	2000	0.499	0.048 [0.039, 0.058]	0.896
16	$k$-NN	$L_{\mathrm{sym}}$	Graph	2000	0.499	0.044 [0.036, 0.054]	0.802
16	$k$-NN	$L_{\mathrm{rw}}$	Matched	2000	0.498	0.048 [0.039, 0.058]	0.977
16	$\varepsilon$	$L_{\mathrm{sym}}$	Matched	2000	0.499	0.056 [0.046, 0.066]	0.992
16	$\varepsilon$	$L_{\mathrm{rw}}$	Matched	2000	0.498	0.050 [0.041, 0.060]	0.571
32	$k$-NN	$L_{\mathrm{sym}}$	Matched	2000	0.502	0.043 [0.035, 0.053]	0.817
32	$k$-NN	$L_{\mathrm{sym}}$	Graph	2000	0.502	0.045 [0.037, 0.056]	0.735
32	$k$-NN	$L_{\mathrm{rw}}$	Matched	2000	0.499	0.053 [0.044, 0.064]	0.776
32	$\varepsilon$	$L_{\mathrm{sym}}$	Matched	2000	0.504	0.054 [0.045, 0.065]	0.724
32	$\varepsilon$	$L_{\mathrm{rw}}$	Matched	2000	0.491	0.050 [0.041, 0.060]	0.080

Gaussian Data-Thinning Baseline

Table S7 adds a simple independent-inference baseline for the power experiment. For a thinning fraction $\rho= 1/2$, draw $E \sim N(0, \sigma^{2}I)$ and set

Then $X^{\mathrm{sel}}$ and $X^{\mathrm{inf}}$ are independent, each with mean $\mu$ and variance $2\sigma^2 I$. We run the same k-NN spectral clustering pipeline on $X^{\mathrm{sel}}$, select the two largest clusters, and test their mean difference on $X^{\mathrm{inf}}$ with the ordinary $\chi_d$ test using the inflated scale $\sqrt{2}\sigma$. This baseline is valid because inference is performed on data independent of the clustering step, but it answers a different question: the selected clusters come from the randomized noisy copy, not from the deterministic spectral clustering of the observed matrix $X$. Power entries where thinning is larger should therefore be read as performance for this randomized target, not as a direct power comparison for the same reported clustering.

Design	Method	$\delta= 0$	1	2	3	4	5
balanced	Matched selective	0.046	0.052	0.158	0.448	0.832	0.952
balanced	Graph-conditioned	0.050	0.042	0.068	0.134	0.418	0.660
balanced	Gaussian thinning	0.046	0.086	0.234	0.728	0.964	0.998
mildly unbalanced	Matched selective	0.054	0.058	0.172	0.486	0.842	0.958
mildly unbalanced	Graph-conditioned	0.060	0.038	0.076	0.216	0.512	0.770
mildly unbalanced	Gaussian thinning	0.054	0.072	0.292	0.808	0.976	0.992

Scale-Sensitivity Check

Table S8 evaluates the same matched selection event under the global null while changing only the scale used in the truncated-$\chi_d$ pivot. The first row uses the true simulation scale and is the theorem-backed setting. The plug-in and independent-pilot rows ignore scale-estimation uncertainty and are included only as sensitivity checks; they should not be read as an unknown-variance theorem. The last two rows show the expected directional effect of a fixed 10% scale misspecification.

Scale used in pivot	Theorem-backed?	$N$	mean $p$	$\mathbb{P}(p \le.05)$	KS $p$
True $\sigma$	Yes	500	0.497	0.056 [0.038, 0.080]	0.765
Same-data plug-in $\hat{\sigma}$	No	500	0.496	0.052 [0.034, 0.075]	0.701
Independent-pilot $\hat{\sigma}$	No	500	0.496	0.056 [0.038, 0.080]	0.744
$0.9\sigma$	No	500	0.457	0.078 [0.056, 0.105]	0.003
$1.1\sigma$	No	500	0.533	0.042 [0.026, 0.063]	0.016

Distributional Stress Tests

Table S9 changes the null distribution while keeping the same selective calculation and nominal $\sigma= 1$ pivot. Only the spherical Gaussian row is covered by the theorem. The standardized $t_5$ row is a heavy-tailed non-Gaussian stress test; the correlated-Gaussian row keeps unit marginal variances but violates the spherical feature covariance assumption. These rows are included to show empirical sensitivity, not to claim validity outside the model. The correlated-feature result is deliberately cautionary: the $\chi_d$ pivot can be distorted when the feature covariance is not proportional to the identity. This is the same caveat that applies to the PBMC illustration: marginal standardization sets unit scale but does not make correlated genes spherical.

Table S9: Distributional stress tests under the global null ($n = 32$, $d = 5$, $K = 3$, $k$-NN, $L_{\mathrm{sym}}$). Brackets give exact binomial 95% confidence intervals for the level .05 tail probability.

Scenario	Theorem-backed?	$N$	mean $p$	$\mathbb{P}(p \le.05)$	KS $p$
Gaussian, spherical	Yes	300	0.486	0.073 [0.047, 0.109]	0.501
Standardized $t_5$	No	300	0.530	0.013 [0.004, 0.034]	0.081
Gaussian, correlated features	No	300	0.408	0.143 [0.106, 0.188]	0.000

Fixed-Attempt Model-Scope Diagnostics

Table S10 extends the distributional stress analysis with 500 fixed attempts per scenario. Only the first row is covered by the spherical-Gaussian theorem. The other rows deliberately retain the same selective calculation under row heteroskedasticity, row correlation, or feature correlation and should be read only as scope diagnostics. All attempts were usable. Row dependence was conservative in this particular design, whereas feature correlation became increasingly anti-conservative; this does not contradict the theorem because these rows violate its covariance assumptions.

Table S10: Fixed-attempt model-scope diagnostics at nominal level .05. Brackets are exact binomial 95% confidence intervals. The KS column tests the full $p$-value distribution against uniformity.

Scenario	Theorem	Usable	Rejection	Operational joint	KS $p$	Median truncation mass
Spherical Gaussian	Yes	500/500	0.052 [0.034, 0.075]	0.052	0.809	6.57e-05
Row heteroskedastic	No	500/500	0.052 [0.034, 0.075]	0.052	0.253	9.23e-05
Row-correlated ($\rho= .2$)	No	500/500	0.028 [0.015, 0.047]	0.028	0.00112	0.000604
Feature-correlated ($\rho= .1$)	No	500/500	0.066 [0.046, 0.091]	0.066	0.472	6.86e-05
Feature-correlated ($\rho= .25$)	No	500/500	0.068 [0.048, 0.094]	0.068	0.0214	2.71e-05
Feature-correlated ($\rho= .5$)	No	500/500	0.144 [0.114, 0.178]	0.144	4e-15	1.92e-05

Reproducibility Notes

The code and computational artifacts described in this section are distributed with the camera-ready PMLR supplementary package; the arXiv source package contains only the paper and this textual supplement. The camera-ready Supplementary Material contains the core implementation, environment specification, fixed seeds, unit tests, quick and full reproduction runners, scripts and raw and summary outputs for the originally submitted experiments, and the manuscript sources. The result-to-script map in RESULTS_MAP.md identifies the files used for each original table and figure.

A dedicated rebuttal-results directory adds the compact raw and summary outputs used for the camera-ready numerical audit, operational usability, partial-null, misspecification, and component-runtime results. It also contains a machine-readable reproducibility manifest recording the original execution environment, seeds, commands, numerical policies, source hash, and output checksums. The accompanying table-generation script validates the checked-in summaries while regenerating the new LaTeX tables. The detailed interval serialization file from the sensitivity audit is not required to rebuild the paper tables and is omitted from the compact camera-ready package because of its size; the included per-policy raw file, risk-case file, high-precision raw file, summaries, and manifest contain the quantities reported here.

All theorem-level claims concern real arithmetic. The included outputs and implementation use floating-point arithmetic; the numerical audit is empirical sensitivity evidence and not interval certification. The 100-digit audit covers root and graph-state sequences only, while eigendecomposition and rounding remain float64. Likewise, the partial-null experiment is theorem-backed only for its independent spherical-Gaussian design, and the non-spherical rows in Table S10 are scope diagnostics.

The principal original output map is:

run_primary_calibration.py      → si_spectral_primary_calibration_*,
run_robustness_matrix.py        → si_spectral_robustness_*,
run_power_curves.py             → si_spectral_power_curve_*,
run_target_rule_comparison.py  → si_spectral_target_rule_*,
run_thinning_baseline.py       → si_spectral_thinning_baseline_*,
run_scale_sensitivity.py       → si_spectral_scale_sensitivity_*,
run_distribution_stress.py     → si_spectral_distribution_stress_*,
run_scaling_experiments.py     → si_spectral_scaling_*,
run_pbmc_featurewise_demo.py   → pbmc_featurewise_*.

Additional Scaling Diagnostics

Table S12: Component timing for the current proof-of-concept implementation (three fixed attempts per graph/$n$ cell). Percentages are fractions of total wall-clock time in the recorded environment and are not hardware-independent complexity constants.

Graph	$n$	Candidate intervals	Spectral solves	Total sec.	Graph walk	Eigen/rounding
$\varepsilon$	32	273	273.3	0.06	9.4%	87.8%
$\varepsilon$	64	1,186	1,186.0	0.56	6.8%	92.8%
$\varepsilon$	128	4,989	4,989.0	6.54	5.5%	94.5%
$k$-NN	32	6,890	107.3	0.18	83.9%	13.9%
$k$-NN	64	59,083	240.3	5.89	97.3%	2.4%
$k$-NN	128	497,839	482.0	244.59	99.6%	0.3%

Table S11: Scaling diagnostics for the finite-state homotopy computation under the global null. Entries are averages over successful replicates; Reps is the number of successful replicates used in the averages, and Fail is the fraction of attempted replicates that failed the deterministic validity checks. Peak memory is Python-side tracemalloc peak memory.

Graph	$d$	$n$	Reps	Breakpoints	Intervals	Solves	Skipped	Time (s)	Peak MB	Fail
eps	5	16	3	72	73	73	0.0%	0.04	0.1	0.0%
eps	5	32	3	297	298	298	0.0%	0.21	0.2	0.0%
eps	5	64	3	1176	1177	1177	0.0%	1.22	0.7	0.0%
eps	5	128	3	4877	4878	4878	0.0%	25.53	2.5	0.0%
eps	20	16	2	39	40	40	0.0%	0.03	0.1	33.3%
eps	20	32	3	237	238	238	0.0%	0.16	0.3	0.0%
eps	20	64	3	1179	1180	1180	0.0%	1.23	1.2	0.0%
eps	20	128	3	4646	4647	4647	0.0%	25.56	4.5	0.0%
knn	5	16	3	829	830	57	93.1%	0.12	0.1	0.0%
knn	5	32	3	7242	7243	122	98.3%	0.58	0.7	0.0%
knn	5	64	3	60892	60893	242	99.6%	8.19	4.6	0.0%
knn	5	128	3	482753	482754	463	99.9%	299.73	24.9	0.0%
knn	20	16	3	894	895	55	93.8%	0.13	0.1	0.0%
knn	20	32	3	7809	7810	121	98.4%	0.62	0.6	0.0%
knn	20	64	3	59744	59745	205	99.7%	8.18	4.2	0.0%
knn	20	128	3	488983	488984	405	99.9%	302.92	26.5	0.0%

Working-Model PBMC Algorithmic Illustration

The PBMC calculation is retained only to demonstrate the feature-wise homotopy on a fixed preprocessed matrix. It is not a confirmatory biological analysis. The protocol fixes a 16-gene panel, log-normalized library-size-adjusted counts, marginal gene standardization, a fixed-seed subsample of 60 cells, $k=5$, $K=3$, and the matched two-largest rule. It then sets a working scale $\sigma=1$. This value is not a known sampling variance, and marginal standardization does not remove gene–gene correlation. In light of Table S10, the values below must not be interpreted as scientifically valid PBMC $p$-values or marker discoveries.

Gene	Mean diff.	Naive value	Selective value	Adjusted selective value
LTB	1.27	6.61e-06	9.72e-06	1.56e-04
CD3E	1.40	6.98e-07	2.44e-04	3.90e-03
IL7R	1.05	2.15e-04	4.29e-04	6.87e-03
CCR7	0.93	9.99e-04	4.44e-03	7.10e-02
CD74	-0.87	2.20e-03	4.95e-03	7.93e-02
CD3D	1.54	5.41e-08	7.14e-03	1.14e-01
CD3G	0.56	4.57e-02	3.36e-02	5.38e-01
CD8B	0.68	1.59e-02	4.18e-02	6.69e-01
NKG7	-0.57	4.56e-02	7.57e-02	1.00e+00
CD8A	0.43	1.31e-01	1.03e-01	1.00e+00
GZMB	-0.52	6.76e-02	1.25e-01	1.00e+00
GNLY	-0.41	1.50e-01	3.11e-01	1.00e+00