Graph Diffusion Residuals for Control-Function Instrumental Variables
arXiv:2606.14636v1 [cs.LG] 12 Jun 2026
Rui Wu∗ Zongyuan Chen ∗ Hong Xie † Defu Lian Enhong Chen
[email protected] [email protected] [email protected] [email protected] [email protected]
School of Computer Science and Engineering University of Science and Technology of China
Abstract Control-function instrumental variable estimators need a first-stage residual, not merely a first-stage prediction. High-capacity first stages can interpolate treatment and leave too little residual information for the outcome equation. We study Adaptive Anisotropic Instrumental Heat Flow (A-IHF), a deterministic graph-diffusion residual extractor for flexible control functions. A-IHF treats treatment as a signal on a graph of first-stage features, uses pilot diffusion to detect large treatment jumps, attenuates conductance across those jumps, and computes the generated control with a sparse graph resolvent. Its observational selection rule uses only (Z, X), combining graph generalized cross-validation, roughness, residualized-treatment relevance, and graph-admissibility filtering. The analysis decomposes error into structural leakage, residual attenuation, and residualized treatment variation, yielding finite-sample bounds, graph-admissibility rates under latent piecewisesmooth geometry, and finite-path selection calibration. Across 54 synthetic benchmark cells with tuned graph, kernel, tree, boosting, series, and neural control-function baselines, guarded observational A-IHF has the lowest average structural-response MSE; the A-IHF family beats the best non-A-IHF baseline in 32 cells. Performance is strongest when the graph captures piecewise-smooth first-stage structure. Keywords: instrumental variables, control functions, graph diffusion, generated controls, nonparametric first stages
1 Introduction Instrumental variables are used when treatment is endogenous and unobserved confounding prevents direct regression from having a causal interpretation (Imbens and Angrist, 1994; Angrist et al., 1996; Angrist and Pischke, 2009; Wooldridge, 2010). In a control-function model, the first stage supplies a residual that enters the outcome equation as a generated control (Rivers and Vuong, 1988; Newey et al., 1999; Blundell and Powell, 2003). The central point is easy to miss: the quality of this residual is not the same as the predictive accuracy of the first-stage regression. A first stage can predict treatment well and, by doing so, remove variation that the second stage needs. The control-function first stage is therefore not only a prediction problem. It is a residual-extraction problem: the estimator must ∗. These authors contributed equally. †. Corresponding author.
©2026 Rui Wu and Zongyuan Chen and Hong Xie and Defu Lian and Enhong Chen. License: CC-BY 4.0, see https://creativecommons.org/licenses/by/4.0/.
Wu, Chen, Xie, Lian, and Chen
A control-function first stage is a residual-extraction problem, not only a prediction problem X and fitted first stage
Prediction target: fit treatment
Residual target: fit systematic component
observed X prediction fit true g(Z)
1.5 1.0
1.0
0.5
0.5
0.0
0.0
interpolation leaves no generated control 0.5
0.5
Residual collapse: corr(V, U) 0 residual / control
0.2
conductance is suppressed at the jump
Residual preserved for the second stage
true U generated residual
0.4
observed X anisotropic graph fit true g(Z)
1.5
true U generated residual, corr=0.88
0.4 0.2
0.0
0.0
0.2
0.2
0.4
0.4 0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
first-stage feature Z first-stage feature Z Toy fractured SCM: the black dotted line marks a jump in the systematic first-stage component.
0.8
1.0
Figure 1: Prediction fit and generated-control recovery are different first-stage targets. In a toy fractured first-stage model, an interpolation-oriented first stage can fit the treatment signal and collapse the generated residual. A graph-resolvent smoother with conductance suppressed near the jump fits the systematic component while preserving residual variation aligned with the outcome-relevant control. The figure is illustrative; Proposition 1 gives a stylized interpolation limit, and the finite-sample theory below decomposes the residual error. remove the systematic treatment component without erasing the endogenous variation that the structural equation uses. This distinction is particularly consequential for flexible IV estimation. Modern estimators, ranging from kernels and neural networks to adversarial and minimax objectives, can match nonlinear first-stage maps with high accuracy (Hartford et al., 2017; Singh et al., 2019; Muandet et al., 2020; Bennett et al., 2019; Dikkala et al., 2020; Syrgkanis et al., 2019). They also change the selection problem. Tuning by hidden residuals, confounders, or structural-response loss uses information unavailable in observational applications, while tuning by factual prediction loss targets a different object. Proposition 1 gives a concrete failure mode: in an interpolation limit for kernel gradient flow, every positive-kernel empirical mode of the generated residual is exponentially damped. This failure mode motivates a first stage that smooths the systematic treatment component without making zero prediction residual the goal. We study Adaptive Anisotropic Instrumental Heat Flow (A-IHF), a gradient-free residual extractor that treats the treatment vector as a signal on a graph built from first-stage features. Smooth regions conduct; large pilot-stage treatment jumps act as barriers. A-IHF estimates the first-stage systematic component by anisotropic graph diffusion, obtains the fitted component through a symmetric resolvent, and uses the remaining signal as the generated control. Stage 1 has no gradient training, while the second stage may use any supervised learner because Y is observed. 2
Graph Diffusion Residuals for Control-Function IV
The paper studies additive control-function designs in which the first-stage systematic component is piecewise smooth on a first-stage feature graph. The claim is graphconditional: A-IHF is useful when the graph representation supports small structural leakage, small residual attenuation, and treatment variation after conditioning on the generated control. A latent low-dimensional geometry is one sufficient route to these quantities. The contributions are as follows. • A deterministic graph residual extractor for control-function IV, with fixed, observational, and graph-admissibility-filtered first-stage selection rules. • A finite-sample theory based on leakage, attenuation, first-stage noise mismatch, and residualized treatment variation, plus graph-admissibility rates under piecewise smooth geometry. • Experiments that connect performance and failure modes to these quantities in synthetic benchmarks, empirical covariate geometries, and a real IV application.
2 Related Work Control functions and nonparametric IV. Control-function methods express endogeneity through a generated residual from a first-stage equation. The approach appears in limited-information models (Rivers and Vuong, 1988) and in nonparametric triangular systems (Newey et al., 1999; Blundell and Powell, 2003). Nonparametric IV estimation also has a conditional-moment formulation, where the structural function is recovered through an ill-posed inverse problem (Newey and Powell, 2003; Hall and Horowitz, 2005; Darolles et al., 2011). Those results study identification, regularization, and inference for nonparametric structural functions. A-IHF uses the control-function decomposition instead. Its contribution is narrower: it constructs a first-stage residual by graph diffusion when the first-stage systematic treatment component is piecewise smooth on a first-stage feature graph. Machine learning IV and orthogonal estimation. Recent IV estimators use kernels, neural networks, adversarial losses, and minimax conditional-moment formulations (Hartford et al., 2017; Singh et al., 2019; Muandet et al., 2020; Bennett et al., 2019; Dikkala et al., 2020). Orthogonal and doubly robust procedures reduce sensitivity to nuisance estimation in the second stage (Chernozhukov et al., 2018; Syrgkanis et al., 2019). These methods address broad IV estimation problems. A-IHF addresses the first-stage residual in a control-function design. It can be combined with orthogonal second-stage scores, as stated in Remark 26. The non-orthogonal control-function regressions reported below are point-estimation results unless stronger inference conditions hold. Graph regularization and diffusion. The graph part of A-IHF is related to manifold regularization, diffusion maps, spectral graph methods, and graph signal processing (Belkin and Niyogi, 2003; Belkin et al., 2006; Coifman and Lafon, 2006; von Luxburg, 2007; Shuman et al., 2013). Existing graph-Laplacian theory studies convergence of graph operators under low-dimensional geometric assumptions (Hein et al., 2005; Belkin and Niyogi, 2006; von Luxburg et al., 2008; Ting et al., 2010; Garcia Trillos and Slepčev, 2018). A-IHF adds a pilot-dependent anisotropic conductance inspired by anisotropic diffusion (Perona and Malik, 3
Wu, Chen, Xie, Lian, and Chen
1990; Weickert, 1998). The theory below treats baseline graph convergence through graphadmissibility conditions and isolates the perturbation induced by adaptive conductance. Positioning. A-IHF has a narrow role: it extracts a generated control in additive controlfunction designs. The graph must encode neighborhoods relevant to the first-stage systematic component. This premise is shared by spectral and manifold-based estimators. The representation experiments test its empirical consequences. The paper records how A-IHF residuals enter standard orthogonal-score arguments under product-rate conditions. The reported non-orthogonal regressions are point-estimation results.
3 Problem Formulation We consider data {(zi , xi , yi )}ni=1 generated by X = g(Z) + V ⋆ ,
V ⋆ = U + ηX ,
Y = f0 (X) + γU + ϵY ,
(1)
where Z ∈ RdZ denotes the first-stage feature vector, X ∈ R is the treatment, V ⋆ = X −g(Z) is the first-stage residual, U is the endogenous component entering the outcome, and ηX is first-stage noise. In the simplest case, Z is the excluded instrument. With included exogenous controls W , the first-stage feature vector may contain both excluded instruments and controls, and the outcome equation may be written as Y = f0 (X, W )+γU +ϵY . To keep notation light, we absorb these first-stage features into Z and write f0 (X) throughout. The usual exclusion and relevance assumptions apply to the excluded-instrument component. The noiseless control-function model is the special case ηX = 0, so V ⋆ = U . The synthetic designs use a small first-stage noise term. The nonlinear outcome is generated with U , while graph residual checks can be reported for both U and V ⋆ . A-IHF is a first-stage estimator under (1); the usual IV assumptions still govern causal interpretation. Weak instruments remain a separate source of instability (Staiger and Stock, 1997; Stock et al., 2002). Throughout, ∥ · ∥2 denotes the Euclidean norm on finite sample vectors. The scaled quantity n−1/2 ∥ · ∥2 is the empirical L2 norm. We write unscaled norms in finite-sample algebra and scaled norms in rate statements. Let x = (x1 , . . . , xn )⊤ , v⋆ = (v1⋆ , . . . , vn⋆ )⊤ , u = (u1 , . . . , un )⊤ , and g = (g(z1 ), . . . , g(zn ))⊤ . A control-function first stage estimates g and forms v̂ = x − ĝ.
(2)
The second stage estimates the structural response while conditioning on v̂. The graph problem is to recover the first-stage residual V ⋆ . The causal control problem is to recover U . Their difference is the first-stage noise term ηX . 3.1 Selection Information Sets We distinguish the information used by the first stage from the information available to the full observational procedure: I1 = σ{(zi , xi ) : 1 ≤ i ≤ n},
Iobs = σ{(zi , xi , yi ) : 1 ≤ i ≤ n}.
(3)
The A-IHF first-stage rule is I1 -measurable. The second stage may use Iobs . In a simulator one may also observe v⋆ , U , ηX , or f0 ; this larger information set is denoted by Ioracle . We use four labels throughout: 4
Graph Diffusion Residuals for Control-Function IV
• fixed : parameters are specified before evaluation; • observational : first-stage parameters are selected using I1 , with Y used only in the second stage or in evaluation; • graph-admissibility-filtered observational : observational selection after applying firststage graph diagnostics, still using only I1 ; • oracle: parameters are selected using hidden simulator quantities. These labels are part of the estimator description. Oracle rows are reference benchmarks. 3.2 Prediction and Residual Extraction A first-stage learner optimized only for prediction of X from Z can remove the residual that the control-function second stage needs. The following proposition is a formal failure mode for interpolation-oriented first-stage learning. It records the kernel-gradient-flow limit in the neural-tangent-kernel regime (Jacot et al., 2018); the spectral-bias perspective in neural networks gives complementary intuition for why low-frequency prediction can dominate learned structure (Rahaman et al., 2019). The statement is not a claim that every neural first stage interpolates in finite samples; it shows why a first-stage objective that drives prediction residuals to zero can destroy the generated control. Proposition 1 (Residual collapse under kernel gradient flow) Let ft denote the prediction of an over-parameterized neural network trained on squared loss for the first-stage regression X on Z. In the neural tangent kernel limit, suppose the empirical kernel matrix Σ is positive definite with eigenpairs (µj , ϕj ), µj > 0. If v̂N N (t) = x − ft (Z), then ⟨v̂N N (t), ϕj ⟩ = e−µj t ⟨v̂N N (0), ϕj ⟩.
(4)
Thus every positive-kernel empirical mode of the generated residual is exponentially damped, and v̂N N (t) → 0 as t → ∞. Appendix A.1 gives the proof. The implication is not that flexible first stages are unusable, but that their objective must be aligned with residual extraction. We therefore seek a first-stage procedure that smooths the systematic component of treatment while preserving a nondegenerate generated control for the second stage.
4 Adaptive Anisotropic Instrumental Heat Flow 4.1 Graph Construction and Pilot Diffusion We view treatment as heat on a first-stage feature graph. Within a smooth region of the first-stage mechanism, neighboring observations should conduct information; across a large jump, conductance should be attenuated so that the systematic component does not bleed across the boundary. Given feature vectors z1 , . . . , zn , we construct a symmetric K-nearestneighbor graph with nonnegative affinity matrix A. In the implementation, distances are converted to radial-basis weights using the median nonzero graph-edge distance, and the graph is symmetrized. The construction follows the graph-based view used in manifold 5
Wu, Chen, Xie, Lian, and Chen
learning, spectral methods, and graph signal processing (Belkin and Niyogi, 2003; Coifman and Lafon, 2006; von Luxburg, 2007; Shuman et al., 2013). For any symmetric weight matrix W , let DW be the diagonal degree matrix and let L(W ) =
DW − W , d¯W
1 d¯W = Tr(DW ). n
(5)
This scaled combinatorial Laplacian is symmetric positive semidefinite and has the constant vector in its nullspace on each connected component. The scaling keeps the diffusion parameter comparable across graphs with different mean degrees. The resolvent form is the graph analogue of Tikhonov-type smoothing (Tikhonov and Arsenin, 1977; Wahba, 1990). A pilot smoother is computed as x̃ = (I + τ L(A))−1 x.
(6)
The pilot is used only to estimate which graph edges cross large treatment jumps. 4.2 Anisotropic Conductance For each edge (i, j) in the graph, define (x̃i − x̃j )2 Cij = exp − γ
.
(7)
The scale γ is the p-th percentile of the nonzero squared pilot differences on graph edges. The final weights are Wij = Aij Cij , (8) with an optional small-conductance cutoff used for numerical sparsity. This is a graph version of anisotropic diffusion (Perona and Malik, 1990; Weickert, 1998): edges with small pilot differences retain conductance, while edges with large pilot differences are attenuated. 4.3 Resolvent Residual For a parameter tuple h = (K, τ, p, λ), let Wh be the anisotropic weight matrix and define Sh = (I + λL(Wh ))−1 .
(9)
The fitted first-stage systematic component and the generated control are ĝh = Sh x,
v̂h = (I − Sh )x.
(10)
No gradient descent is used in the first stage. 4.4 Observational First-Stage Selection A-IHF can be used with fixed defaults or with an observational rule over a finite candidate family Hn . The rule uses only (Z, X). For each h ∈ Hn , let Qh = [1, v̂h ],
† ⊤ Mh = I − Qh (Q⊤ h Qh ) Qh ,
6
(11)
Graph Diffusion Residuals for Control-Function IV
and define the residualized treatment variation κn (h) =
1 ⊤ x Mh x. n
(12)
The observational score is ⊤
2
b obs (h) = ∥(I − Sh )x∥2 /n + α (Sh x) L(Wh )(Sh x) . Q {1 − Tr(Sh )/n}2 ∥x∥22 /n + ε
(13)
The first term is graph generalized cross-validation for the first-stage smoother (Craven and Wahba, 1979; Golub et al., 1979). The fixed-smoother identity in Proposition 4 gives the idealized motivation for this term. The second term penalizes graph roughness of the fitted first-stage systematic component. The trace can be estimated by stochastic trace probes (Hutchinson, 1989). The score is an observational graph-admissibility proxy: it rewards residual size, controls effective degrees of freedom, and discourages rough fitted systematic components inside the chosen candidate family. The classical scalar-variance argument assumes a fixed linear smoother and homoskedastic noise. Here the smoother is data-adaptive and the residual may be heteroskedastic or graph-correlated. The correlatedresidual experiments test this limitation. For graph-admissibility filtering, let ℓn (h) be the largest post-cut connected-component fraction of the final graph induced by Wh , and let dmin,n (h) be its minimum final weighted degree after the numerical conductance cutoff. Given thresholds ωn ∈ [0, 1] and d0,n ≥ 0, define Hnκ = h ∈ Hn : κn (h) ≥ cκ ∥x − x̄1∥22 /n , (14) Hnadm (ωn , d0,n ) = {h ∈ Hnκ : ℓn (h) ≥ ωn ,
dmin,n (h) ≥ d0,n } .
(15)
The unfiltered observational rule selects over Hnκ , b obs (h). ĥobs ∈ arg minκ Q h∈Hn
(16)
The graph-admissibility-filtered observational rule selects ĥgobs ∈ arg
min
adm (ω ,d h∈Hn n 0,n )
b obs (h). Q
(17)
We take α ≥ 0, ε > 0, cκ > 0, ωn , and d0,n as fixed before evaluation. These constants have scale meanings rather than oracle meanings. The term ε is a numerical stabilizer. The relevance floor cκ requires a minimum fraction of centered treatment variation to remain after projection on the generated control. The graph-admissibility thresholds require the final graph to retain a macroscopic connected region and nondegenerate local degrees. The roughness weight α fixes the scale of the graph-energy penalty within the pre-specified family Hn ; it is not selected using Y , v⋆ , or the structural response. The candidate family and admissibility filters are part of the procedure. A largestcomponent screen is analogous to a giant-component diagnostic in random-graph and percolation analyses (Bollobás, 2001; Grimmett, 1999): the selected graph should retain a macroscopic connected region rather than shatter into many microscopic pieces. In Section 7
Wu, Chen, Xie, Lian, and Chen
7.3.1, we report that allowing sparse graphs and sharp cuts can over-fragment the graph. We treat Hn and Hnadm as admissible graph families, not unconstrained hyperparameter searches. The score in (13) is a finite-path rule. It balances residual size, effective degrees of freedom, and fitted-signal roughness inside the specified candidate family. A minimax analysis would require a continuum of graph scales and a uniform risk approximation for adaptive smoothers. Section 6.3 gives a split-sample calibration inequality for the finite family used by the procedure. Algorithm 1 Adaptive Anisotropic Instrumental Heat Flow 1: Input: first-stage features z1 , . . . , zn , treatment vector x, and parameters (K, τ, p, λ). 2: Construct a symmetric K-nearest-neighbor affinity matrix A. 3: Form L(A) and compute x̃ = (I + τ L(A))−1 x. 4: Set γ to the p-th percentile of nonzero (x̃i − x̃j )2 over graph edges. 5: For graph edges, define Cij = exp{−(x̃i − x̃j )2 /γ} and Wij = Aij Cij ; set Wij = 0 otherwise. 6: Form S = (I + λL(W ))−1 . 7: Output: ĝ = Sx and v̂ = (I − S)x. The observational versions apply Algorithm 1 to each h ∈ Hn and select h by either (16) or the graph-admissibility-filtered rule (17).
STAGE 1 Stage 1 uses first-stage information only. 𝑽∗ , 𝑼, 𝒇𝟎 , structural-response loss and counterfactual outcomes are not included.
input
Graph
Pilot
features 𝒁 treatment 𝑿
K-NN affinity
𝑨 = 𝑨(𝒁)
Conductance
Control
= 𝑺𝝉 𝑨 𝒙 𝒙
𝑾𝒊𝒋 = 𝑨𝒊𝒋 𝑪𝒊𝒋
edge jumps
cut large jumps
ෝ = 𝑺𝝀 𝑾 𝒙 𝒈 ෝ= 𝒙 − 𝒈 ෝ 𝒗
Observational path selection
STAGE 2
ෝ, 𝒀) fit response with (𝑿, 𝒗
GCV + graph roughness screen by 𝑲𝒏 (𝒉)
Figure 2: A-IHF mechanism. The first stage constructs a graph on first-stage features, computes a pilot diffusion, attenuates conductance across large pilot treatment jumps, and forms a generated control by a final graph resolvent. The observational tuning rules use only (Z, X), relevance, and graph-admissibility diagnostics; hidden simulator quantities are not part of Stage 1 selection.
5 Finite-Sample Algebraic Properties This section records algebraic properties that hold after a graph smoother has been realized. The results hold without a graph limit. They describe the residual tradeoff, the secondstage perturbation, and the empirical certificate used in the experiments. Section 6 gives large-sample conditions under which these quantities vanish. 8
Graph Diffusion Residuals for Control-Function IV
5.1 Standing Assumptions and Targets Assumption 2 (Additive first stage) On the observed sample, x = g + v⋆ ,
v⋆ = u + η X ,
(18)
where g = (g(z1 ), . . . , g(zn ))⊤ is the first-stage systematic treatment component, v⋆ is the first-stage residual, u is the outcome-relevant endogenous control, and η X is first-stage noise. The graph residual error is evaluated against v⋆ . The second-stage perturbation depends on the distance from v̂ to u. Assumption 3 (Linear control-function second stage) For the theoretical perturbation bound, yi = β0 xi + γ0 ui + ϵi , E[ϵi | xi , ui ] = 0, i = 1, . . . , n. (19) Equivalently, y = β0 x + γ0 u + ϵ in vector notation. This condition is used only for the finitesample coefficient bound. The nonlinear response experiments are empirical evaluations of the same generated controls. The finite-sample bounds below condition on the realized smoother Sh . Assumption 13, introduced in Section 6.2, gives one set of sufficient conditions for such smoothers to recover the target residual. 5.2 Fixed-Smoother Selection Identities Proposition 4 (Fixed-smoother risk identity) Assume x = g + η, where E[η | Z] = 0 and E[ηη ⊤ | Z] = σ 2 I. Let S be a symmetric smoother treated as fixed conditional on Z. Then E ∥(I − S)x∥22 | Z = ∥(I − S)g∥22 + σ 2 Tr{(I − S)2 }. (20) Proof Since (I − S)x = (I − S)g + (I − S)η, h i E ∥(I − S)x∥22 | Z = ∥(I − S)g∥22 + E η ⊤ (I − S)2 η | Z = ∥(I − S)g∥22 + σ 2 Tr{(I − S)2 }.
(21) (22)
Proposition 4 is an idealization. It motivates using residual size and effective degrees of freedom to rank smoothers. Equation (13) turns this identity into a finite-family admissibility score for data-adaptive graph smoothers, where IV residuals need not have scalar conditional variance. 5.3 Oracle Selection Is Not Observational Selection Proposition 5 (Information-set separation) Let P and P ′ be two mechanisms that induce the same law for the observed sample (Z, X, Y ). Suppose that, under a coupling with the same observed sample, their oracle selectors over Hn satisfy P{h⋆P ̸= h⋆P ′ } > 0, 9
(23)
Wu, Chen, Xie, Lian, and Chen
where h⋆P and h⋆P ′ may depend on hidden residual targets. Then no Iobs -measurable map can equal both oracle selectors with probability one. Proof Let φ be an Iobs -measurable selector. Under the coupling, φ takes the same value under P and P ′ because the observed sample is the same. On the event {h⋆P ̸= h⋆P ′ }, it cannot equal both oracle selectors. This event has positive probability. This proposition separates two selection targets. Oracle-tuned rows use hidden residual information. Observational first-stage rules are functions of I1 . 5.4 Residual Error Decomposition Proposition 6 (Exact graph residual decomposition) Under Assumption 2, let v̂h = (I − Sh )x. Then v̂h − v⋆ = (I − Sh )g − Sh v ⋆ . (24) | {z } | {z } structural leakage
residual attenuation
Proof Substitute x = g + v⋆ into v̂h = (I − Sh )x and subtract v⋆ . Equation (24) gives the algebraic tradeoff. A first-stage smoother should make (I − Sh )g small and keep Sh v⋆ small. Isotropic smoothing can fail when g has discontinuities on the graph. Anisotropic conductance reduces edges across large pilot jumps and can lower structural leakage across such discontinuities. If v⋆ has low-frequency graph structure, then Sh v⋆ need not be small; this is the spectral leakage failure mode. Proposition 7 (Spectral leakage for correlated residuals) Let S = (I + λL)−1 for a fixed graph Laplacian L with eigenpairs (µj , ϕj ). If E[v⋆ (v⋆ )⊤ ] = ΣV , then E∥Sv⋆ ∥22 =
n X ϕ⊤ j ΣV ϕj . (1 + λµj )2
(25)
j=1
When residual variation has energy in low-frequency graph modes, the denominator in (25) is close to one. Graph smoothing cannot remove that component without also smoothing low-frequency structure. 5.5 Generated-Control Coefficient Error Let β̂(v̂) be the least-squares coefficient on x when regressing y on (x, v̂) and an intercept. Proposition 8 (Generated-control coefficient decomposition) Under Assumption 3, let v̂ be any generated control. Let Mv̂ be the residual-maker after projecting onto [1, v̂], and suppose 1 κn (v̂) = x⊤ Mv̂ x > 0. (26) n Then x⊤ Mv̂ (u − v̂) x⊤ Mv̂ ϵ β̂(v̂) − β0 = γ0 + . (27) x⊤ Mv̂ x x⊤ Mv̂ x 10
Graph Diffusion Residuals for Control-Function IV
Consequently, β̂(v̂) − β0 −
x⊤ Mv̂ ϵ |γ0 | ∥v̂ − u∥2 . ≤p ⊤ x Mv̂ x nκn (v̂)
(28)
The relevance screen in (16) controls a collinearity failure. If Sh ≈ 0, then v̂h = (I − Sh )x ≈ x; after projecting on the generated control, little treatment variation remains and κn (h) is small. This is distinct from first-stage interpolation, where Sh ≈ I and v̂h ≈ 0. The GCV term and the relevance screen address different degeneracies. Corollary 9 (Graph-to-causal perturbation) Under Assumptions 2 and 3, let v̂h = (I − Sh )x and suppose κn (h) > 0. Define Eϵ,h = Then
x⊤ Mv̂h ϵ . x⊤ Mv̂h x
|γ0 | |β̂(v̂h ) − β0 − Eϵ,h | ≤ p (∥(I − Sh )g∥2 + ∥Sh v⋆ ∥2 + ∥v⋆ − u∥2 ) . nκn (h)
(29)
(30)
Equivalently, 1/2 |γ0 | |β̂(v̂h ) − β0 − Eϵ,h | ≤ p 3∥(I − Sh )g∥22 + 3∥Sh v⋆ ∥22 + 3∥v⋆ − u∥22 . nκn (h)
(31)
Proof By Proposition 6, ∥v̂h − u∥2 ≤ ∥v̂h − v⋆ ∥2 + ∥v⋆ − u∥2 ≤ ∥(I − Sh )g∥2 + ∥Sh v⋆ ∥2 + ∥v⋆ − u∥2 .
(32)
Substitute this inequality into Proposition 8. The second display follows from (a + b + c)2 ≤ 3(a2 + b2 + c2 ). Corollary 9 is the link between graph signal processing and the causal second stage. The first term is structural leakage from removing part of g(Z). The second term is residual attenuation from smoothing away part of v⋆ . The third term is first-stage noise in the generated residual target. It vanishes in the noiseless control-function model. The denominator records whether treatment still has variation after conditioning on the generated control. 5.6 Graph-Admissibility Certificate For a realized smoother Sh , define Leakn (h) = n−1/2 ∥(I − Sh )g∥2 ,
(33)
−1/2
(34)
Attenn (h) = n
⋆
∥Sh v ∥2 ,
Noisen = n−1/2 ∥v⋆ − u∥2 ,
(35)
Reln (h) = κn (h).
(36)
The first quantity measures structural leakage. The second measures residual attenuation. The third measures first-stage noise in the control target. The fourth measures residualized treatment variation after conditioning on the generated control. 11
Wu, Chen, Xie, Lian, and Chen
Proposition 10 (Graph-admissibility certificate) Under Assumptions 2 and 3, it follows that if Reln (h) > 0, then n−1/2 ∥v̂h − v⋆ ∥2 ≤ Leakn (h) + Attenn (h),
(37)
and, with Eϵ,h defined in (29), |γ0 | {Leakn (h) + Attenn (h) + Noisen } . |β̂(v̂h ) − β0 − Eϵ,h | ≤ p Reln (h) Proof Equation (37) is Proposition 6 after dividing by
√
(38)
n. Also,
n−1/2 ∥v̂h − u∥2 ≤ Leakn (h) + Attenn (h) + Noisen .
(39)
Equation (38) follows from Proposition 8. In a simulator, g, u, and v⋆ are known, so (33)–(36) can be reported directly. In an observational application, Leakn (h), Attenn (h), and Noisen are hidden. A-IHF uses graph roughness, GCV, graph connectivity, and Reln (h) as observable checks. The failure modes in the experiments correspond to the terms above: cross-boundary smoothing raises structural leakage, low-frequency residual structure raises attenuation, first-stage noise separates V ⋆ from U , and weak instruments reduce relevance.
6 Asymptotic Theory and Inference This section first reduces consistency to graph-admissibility conditions on leakage, attenuation, noise mismatch, and relevance. It then gives a sufficient route through latent piecewisesmooth geometry and a local graph scale. The geometric result also states how the algorithmic regularization parameter λn maps to the continuum smoothing scale. 6.1 Consistency Under Graph-Admissible Sequences The finite-sample bounds in Section 5 imply the abstract reduction below. The conditions are stated for a generic sequence of graph smoothers, so they also apply to other graph filters and to approximate solvers.
Assumption 11 (Graph-admissible smoother sequence) Consider a triangular array satisfying Assumptions 2 and 3. Let Sn = Shn be a sequence of A-IHF smoothers, where hn may be fixed, observationally selected, or otherwise measurable with respect to the first-stage sample. Let v̂n = (I − Sn )xn . 12
(40)
Graph Diffusion Residuals for Control-Function IV
There exist deterministic sequences ag,n → 0, av,n → 0, and ax,n ≥ 0, and a constant κ0 > 0, such that n−1/2 ∥(I − Sn )gn ∥2 = Op (ag,n ), −1/2
(41)
n ∥Sn vn⋆ ∥2 = Op (av,n ), n−1/2 ∥vn⋆ − un ∥2 = Op (ax,n ),
(42)
Pr{κn (v̂n ) ≥ κ0 } → 1,
(44)
x⊤ n Mv̂n ϵn ⊤ xn Mv̂n xn
(45)
= Op (n−1/2 ).
(43)
Condition (41) is structural preservation: the smoother must not erase the first-stage systematic treatment component. Condition (42) is residual removal: the part of V ⋆ retained by the smoother must vanish in root-mean-square norm. Condition (43) records the distance between the first-stage residual and the outcome-relevant control. Condition (44) is a relevance condition for the second stage. Condition (45) is the usual sampling requirement for the outcome noise term after partialling out the generated control. Proposition 12 (Consistency reduction to graph admissibility) Under Assumption 11, n−1/2 ∥v̂n − vn⋆ ∥2 = Op (ag,n + av,n ). (46) If ag,n + av,n → 0, then the generated control is root-mean-square consistent for V ⋆ . If also ax,n → 0, it is root-mean-square consistent for U . Under the linear second-stage model in Assumption 3, |β̂(v̂n ) − β0 | = Op (ag,n + av,n + ax,n ) + Op (n−1/2 ). (47) p
In particular, if ax,n → 0, then β̂(v̂n ) → − β0 . Proof By Proposition 6, v̂n − vn⋆ = (I − Sn )gn − Sn vn⋆ .
(48)
Taking norms and applying (41)–(42) gives (46). For the second stage, use the exact decomposition in the proof of Proposition 8: β̂(v̂n ) − β0 = γ0
x⊤ x⊤ n Mv̂n (un − v̂n ) n Mv̂n ϵn + . ⊤ ⊤ xn Mv̂n xn xn Mv̂n xn
(49)
On the event κn (v̂n ) ≥ κ0 , x⊤ ∥Mv̂n xn ∥2 ∥un − v̂n ∥2 n Mv̂n (un − v̂n ) ≤ ⊤ nκn (v̂n ) xn Mv̂n xn
(50)
n−1/2 ∥un − v̂n ∥2 p κn (v̂n )
(51)
= Op (ag,n + av,n + ax,n ).
(52)
=
13
Wu, Chen, Xie, Lian, and Chen
The sampling term is controlled by (45). This proves (47). Proposition 12 is a reduction result. Its role is to show that structural leakage, residual attenuation, first-stage noise mismatch, and second-stage relevance determine generatedcontrol consistency for the coefficient. Cross-boundary smoothing violates (41). Lowfrequency graph correlation in V ⋆ violates (42). Persistent first-stage noise violates (43). A generated control nearly collinear with X violates (44). The experiments below are organized around these quantities. 6.2 Concrete Graph-Admissibility Conditions We now state graph-geometric conditions that control the residual-recovery part of graph admissibility. Relevance, first-stage noise mismatch, and outcome-noise sampling remain separate second-stage conditions. Let Ti ∈ M be a latent coordinate, where M is a compact m-dimensional Riemannian manifold with volume density bounded above and below. The graph may be built from Zi directly or from a representation ψ(Zi ). The conditions require this representation to preserve local neighborhoods; representation learning itself is treated as a separate problem. Assumption 13 (Latent geometry and piecewise smoothness) The variables Ti are i.i.d. on M. There is a representation ψ and constants 0 < cψ < Cψ < ∞ such that, with probability tending to one, cψ dM (Ti , Tj ) ≤ ∥ψ(Zi ) − ψ(Zj )∥2 ≤ Cψ dM (Ti , Tj )
(53)
for all pairs with dM (Ti , Tj ) ≤ c0 . The manifold is partitioned into finitely many closed regions R [ M= Mℓ , (54) ℓ=1
with piecewise C 1 boundaries. On each region, g is s-Hölder with constant Lg , where 0 < s ≤ 2. If an edge connects observations from different regions and is not contained in an o(1) boundary band, then |g(Ti ) − g(Tj )| ≥ ∆ (55) for a constant ∆ > 0. Remark 14 (Geometry and included controls) Assumption 13 is a condition on the metric representation used by the graph. It is not a claim that raw excluded instruments, discrete controls, and fixed effects jointly form a smooth Riemannian manifold. Discrete included controls can be handled by stratification or localization before graph construction. In applications that concatenate standardized controls with excluded instruments, the graph is an empirical first-stage feature graph. The geometric rate theory should then be read conditionally on a representation that preserves the relevant local neighborhoods. The Card schooling check is interpreted in this empirical sense, not as evidence for the manifold assumption. 14
Graph Diffusion Residuals for Control-Function IV
Let K = Kn and define the intrinsic graph scale Kn 1/m εn = . n
(56)
We use the nearest-neighbor regime Kn → ∞, log n
Kn log n → 0. n
(57)
Under standard concentration conditions for random geometric graphs, graph edges have intrinsic length Op (εn ) away from boundary effects (Belkin and Niyogi, 2006; Hein et al., 2005; Singer, 2006; Ting et al., 2010; Garcia Trillos and Slepčev, 2018). The first condition is the usual connectivity-side condition for local graphs. The second keeps neighborhoods local. After anisotropic weighting, the final graph may disconnect across jump boundaries by construction; the connectivity requirement applies to the separated within-region graph. The scaled Laplacian in (5) is not divided by ε2n . On smooth functions supported away from region boundaries, the low-frequency scaling has the form L0n f = cL ε2n (−∆M f ) + en,f ,
∥en,f ∥L2 (Pn ) ≤ ζn ∥f ∥S ,
(58)
for a graph-dependent constant cL > 0 and a smoothness norm ∥ · ∥S . Hence the final resolvent (I + λn Ln )−1 corresponds to a continuum heat scale tn ≍ λn ε2n ,
tp,n ≍ τn ε2n ,
(59)
for the final and pilot diffusions. Constants such as cL are absorbed into tn . The rate statements below use tn as shorthand for the algorithmic scale λn ε2n . Assumption 15 (Pilot accuracy and threshold separation) Let x̃ be the pilot diffusion in (6). There is a sequence rp,n → 0 such that ∥x̃ − g∥∞ ≤ rp,n
(60)
with probability tending to one. The conductance scale γn satisfies (Lg εsn + 2rp,n )2 → 0, γn
γn → 0. (∆ − 2rp,n )2
(61)
This is a high-level sufficient condition used to analyze jump separation by the pilot conductance. It is not an observational diagnostic and is not claimed for arbitrary noisy first stages. Lemma 16 (Jump separation by anisotropic conductance) Suppose Assumptions 13 and 15 hold. For graph edges within a single smooth region and away from the boundary band, min Cij → 1 (62) (i,j)∈Ein
in probability. For graph edges that cross a jump satisfying (55), max (i,j)∈Ejump
Cij → 0
in probability. 15
(63)
Wu, Chen, Xie, Lian, and Chen
Lemma 17 (Percentile threshold) Let Dij = (x̃i − x̃j )2 on graph edges and let γ̂n be the empirical θ-quantile of {Dij : Aij > 0}, with θ = p/100. Suppose there are deterministic sequences an < bn and a constant cγ < ∞ such that, with probability tending to one, max Dij ≤ cγ γ̂n ,
(64)
γ̂n /bn → 0,
(65)
Dij ≥ bn .
(66)
(i,j)∈Ein
min (i,j)∈Ejump
Then the percentile rule yields min Cij ≥ e−cγ + op (1),
(i,j)∈Ein
max (i,j)∈Ejump
Cij → 0.
(67)
Lemma 16 gives uniform within-region conductance for a deterministic threshold. Lemma 17 states the weaker condition used by the percentile implementation: within-region weights remain comparable, while cross-jump weights vanish. This is sufficient for the separatedgraph rate below. Assumption 18 (Bridge condition for the separated graph) Let Bn be the set of observations within intrinsic distance cεn of a region boundary, and define |Bn | 1/2 . (68) ξn = n Let Wn0 be the oracle separated weight matrix obtained by retaining only within-region edges outside Bn , with weights comparable to the original affinities. Let L0n = L(Wn0 ),
Sn0 = (I + λn L0n )−1 .
(69)
Let tn ≍ λn ε2n be the continuum diffusion scale induced by the algorithmic parameter λn , and let ζn be the fixed-graph approximation error for the separated graph. The bias term below is the heat-semigroup truncation bias associated with the graph-to-manifold scaling in (58); tn is not an additional algorithmic input. Assume ζn → 0 and n−1/2 ∥(I − Sn0 )gn ∥2 ≤ C{ts/2 n + ζn + ξn },
(70)
Tr{(Sn0 )2 } ≤ Ct−m/2 . n
(71) s/2
Here Sn is the A-IHF smoother after the anisotropic conductance step. The term tn is the heat-smoothing bias on an m-dimensional geometry. The term ζn records fixed-graph approximation. The term ξn records the boundary band. For piecewise C 1 boundaries and 1/2 bounded density, ξn = Op (εn ). For fixed local graphs with nonadaptive weights, graph-Laplacian and graph-heat approximations are well studied under geometric regularity conditions (Belkin and Niyogi, 2006; von Luxburg et al., 2008; Garcia Trillos and Slepčev, 2018). The percentile conductance used by A-IHF is data-dependent. The theorem below does not prove a continuum limit for that adaptive thresholding rule. It reduces the adaptive-graph part to the perturbation quantity ηn . 16
Graph Diffusion Residuals for Control-Function IV
Lemma 19 (Resolvent perturbation by leaked graph mass) Let Sn = (I + λn Ln )−1 and Sn0 = (I + λn L0n )−1 , where Ln and L0n are symmetric positive semidefinite. Define ηn = λn Ln − L0n op .
(72)
n−1/2 ∥(Sn − Sn0 )xn ∥2 = Op (ηn ).
(73)
If n−1/2 ∥xn ∥2 = Op (1), then
If 0 < cd ≤ d¯0n , d¯n ≤ Cd < ∞, and
L0n op ∨ ∥Ln ∥op ≤ CL ,
|d¯n − d¯0n | ≤ Cℓn , d¯0
X 1 0 ℓn = ¯0 max |Wn,ij − Wn,ij |, dn i j
(74)
n
(75)
then Ln − L0n op ≤ Cℓn . Thus ηn ≤ Cλn ℓn . Assumption 20 (Residual spectral regularity) Conditional on the latent geometry and on the separated graph, vn⋆ is mean zero and sub-Gaussian with covariance σV2 In . More generally, the analysis uses the following operative condition: s −m/2 t n . n−1/2 ∥Sn0 vn⋆ ∥2 = Op (76) n The isotropic sub-Gaussian case implies (76): conditionally on the separated graph, E{∥Sn0 vn⋆ ∥22 | Z} = σV2 Tr{(Sn0 )2 } ≤ CσV2 t−m/2 . n
(77)
Markov’s inequality gives the stated order. Assumption 20 is the separation requirement between the graph smoothness of g and that of V ⋆ . Correlated residuals can violate it. Proposition 7 gives the finite-sample form of that violation. The next theorem is conditional on the adaptive graph satisfying the perturbation control encoded by ηn . Theorem 21 (Residual recovery rate under piecewise graph smoothness) Suppose Assumptions 2 and 13–20 hold, and suppose n−1/2 ∥xn ∥2 = Op (1). Let ηn be defined by (72), let tn ≍ λn ε2n , and let s −m/2
rn = ts/2 n +
tn
n
+ ζ n + ξn + η n .
(78)
Then the A-IHF generated control satisfies n−1/2 ∥v̂n − vn⋆ ∥2 = Op (rn ). 17
(79)
Wu, Chen, Xie, Lian, and Chen
m/2
If tn → 0, ntn → ∞, ζn → 0, ξn → 0, and ηn → 0, then v̂n is root-mean-square consistent for vn⋆ . If ζn , ξn , and ηn are of smaller order and tn ≍ n−2/(2s+m) , equivalently λn ≍ n−2/(2s+m) ε−2 n ,
(80)
n−1/2 ∥v̂n − vn⋆ ∥2 = Op n−s/(2s+m) .
(81)
then
The rate separates four terms: heat-smoothing bias, effective residual dimension, fixedgraph approximation and boundary error, and leaked cross-jump graph mass. Assumption 18 abstracts the baseline spectral convergence of the separated geometric graph, which is standard for local graphs under regularity conditions. The theorem is a reduction from adaptive graph construction to separated-graph heat approximation, not a continuum limit for the adaptive percentile graph. Once the adaptive graph is close to a separated local graph in the sense of ηn , the residual recovery rate follows. A continuum theory for percentileselected, data-adaptive anisotropic graph diffusion remains a separate graph-limit problem. Corollary 22 (Explicit second-stage rate) Under the conditions of Theorem 21 and Assumption 3, suppose also that (82)
Pr{κn (v̂n ) ≥ κ0 } → 1 for some κ0 > 0. Define Eϵ,n =
x⊤ n Mv̂n ϵn . x⊤ n Mv̂n xn
(83)
Then |β̂(v̂n ) − β0 − Eϵ,n | = Op (rn + δx,n ),
δx,n = n−1/2 ∥vn⋆ − un ∥2 .
(84)
If Eϵ,n = Op (n−1/2 ), then |β̂(v̂n ) − β0 | = Op (rn + δx,n ) + Op (n−1/2 ),
(85)
Corollary 23 (Root-n linear inference under negligible generated-control error) In the linear second-stage model of Assumption 3, suppose p
∥v̂n − un ∥2 = op (1),
κn (v̂n ) → − κ0 > 0.
(86)
1 ⊤ p x Mv̂n xn → − Q > 0. n n
(87)
Suppose further that 1 √ x⊤ Mv̂n ϵn ⇝ N (0, Ω), n n Then
√ n{β̂(v̂n ) − β0 } ⇝ N (0, Ω/Q2 ). 18
(88)
Graph Diffusion Residuals for Control-Function IV
Proof Use the decomposition in the proof of Proposition 8. The generated-control term is bounded by √ |γ0 | np ∥v̂n − un ∥2 = Op (∥v̂n − un ∥2 ) = op (1). (89) nκn (v̂n ) The sampling term converges by (87) and Slutsky’s theorem. Condition (86) is stronger than the residual rates used for point estimation. Theorem 21 controls n−1/2 ∥v̂n − vn⋆ ∥2 ; the condition above controls the unnormalized distance to un . It requires both accurate graph residual recovery and negligible first-stage noise mismatch. Corollary 23 is a benchmark statement for negligible generated-control error. The orthogonal construction below gives a standard route when first-stage rates are slower. 6.3 Selection Calibration The score in (13) is data-adaptive because the conductance uses the pilot-smoothed treatment. For a finite-sample calibration statement, define a split version. One part of the sample constructs the candidate graphs and conductances. Conditional on that split, the candidate smoothers {Sh : h ∈ Hn } are fixed on the evaluation split. The experiments use the corresponding transductive full-sample implementation. For fixed candidate smoothers, define the sample residual-recovery loss Rn (h) =
1 1 ∥Sh x − g∥22 = ∥v̂h − v⋆ ∥22 , n n
(90)
where the equality follows because Sh x − g = −(v̂h − v⋆ ). Define also the roughness penalty Jn (h) =
(Sh x)⊤ L(Wh )(Sh x) . ∥x∥22 /n + ε
(91)
Since L(Wh ) is positive semidefinite and ε > 0, Jn (h) ≥ 0. Let Hnκ be the candidates satisfying the relevance screen in (16). Proposition 24 (Calibration inequality for the split observational selector) Assume α ≥ 0, ε > 0, Hnκ is nonempty, Hn is finite, each candidate smoother satisfies 0 ⪯ Sh ⪯ I, and 1 1 − dh ≥ cdf > 0 (92) dh = Tr(Sh ), n for all h ∈ Hn . Suppose that, conditional on the graph-construction split, b obs (h) − {Rn (h) + αJn (h) + σV2 } ≤ ∆n sup Q
κ h∈Hn
(93)
with probability at least 1 − δn . Then the split selector satisfies Rn (ĥobs ) ≤ inf κ {Rn (h) + αJn (h)} + 2∆n h∈Hn
with probability at least 1 − δn . 19
(94)
Wu, Chen, Xie, Lian, and Chen
In addition, suppose x = g + η on the evaluation split, the entries of η are conditionally independent, mean zero, sub-Gaussian with variance proxy bounded by ση2 , E[ηη ⊤ | Z] = σV2 I, and n−1 ∥g∥22 ≤ Cg . Define h i b obs (h) − {Rn (h) + αJn (h) + σ 2 } Z . bsel,n = sup E Q (95) V κ h∈Hn
Then, for all 0 < δ < 1, there is a constant C depending only on cdf , ση , and Cg such that ) (r log(2|Hn |/δ) log(2|Hn |/δ) (96) + ∆n = bsel,n + C n n is valid with conditional probability at least 1 − δ. If Tr(Sh ) is replaced by an estimator c h ) satisfying suph |dbh − dh | ≤ etr,n < cdf /2, the same bound gains an additive term Tr(S Cetr,n . For A-IHF candidates, 0 ⪯ Sh ⪯ I holds because L(Wh ) is symmetric positive semidefinite and Sh = (I + λL(Wh ))−1 . The bias term bsel,n records the mismatch between the GCV proxy and residual-recovery risk. It may be non-negligible for adaptive graphs, heteroskedastic residuals, or graph-correlated residual variation. Thus Proposition 24 is a finite-family calibration result. It becomes an oracle-type statement only when bsel,n and the traceestimation error are small. Remark 25 (Asymptotic behavior of the selection-bias term) The idealized case in which bsel,n is expected to vanish is narrow but informative. If the candidate smoothers are fixed conditional on Z, the first-stage residual V ⋆ is conditionally homoskedastic and mean independent of the graph features, and its conditional covariance has no systematic lowfrequency alignment with the candidate graph eigenspaces, then the classical GCV identity approximates the residual-recovery risk up to sampling and trace-estimation errors. Under a finite or slowly growing candidate family and a uniform law of large numbers for the quadratic forms {n−1 ∥(I −Sh )x∥22 : h ∈ Hnκ }, this gives bsel,n = op (1). These conditions are deliberately stronger than the empirical setting. Their role is to identify what the observational score is trying to approximate; when heteroskedasticity or graph-correlated residual variation violates them, bsel,n becomes a selection-bias term rather than a vanishing remainder. The correlatedresidual design is included for exactly this failure mode. Establishing bsel,n → 0 for fully adaptive graph GCV remains a separate uniform risk problem. The proposition compares the split selector with the best candidate in the admissible finite family for the penalized residual-recovery target. The candidate family and the relevance screen are part of the estimator. 6.4 Orthogonal Score Implication The OLS control-function coefficient in Proposition 8 gives a point-estimation perturbation bound. Root-n inference requires stronger residual rates or an orthogonal score. A standard construction uses cross-fitting and solves n
1X ψ(Wi ; β, η̂−k(i) ) = 0, n i=1
20
(97)
Graph Diffusion Residuals for Control-Function IV
where Wi = (Yi , Xi , Zi ), the nuisance η contains the generated control and second-stage regressions, and k(i) denotes the fold not containing observation i. Remark 26 (Orthogonal-score normality under product-rate conditions) Let parameter β0 satisfy E[ψ(W ; β0 , η0 )] = 0. Suppose the score is Neyman orthogonal at η0 , has nonsingular Jacobian G = ∂β E[ψ(W ; β, η0 )]|β=β0 , and satisfies a central limit theorem with variance Ω = E[ψ(W ; β0 , η0 )ψ(W ; β0 , η0 )⊤ ]. (98) If the cross-fitted nuisance estimates obey the usual product-rate condition, in particular if the A-IHF control enters with L2 rate ρn = o(n−1/4 ) relative to the outcome-relevant control and the remaining nuisance errors have compatible rates, then the solution β̃ of (97) satisfies √
n(β̃ − β0 ) ⇝ N 0, G−1 ΩG−⊤ .
(99)
A plug-in sandwich estimator using cross-fitted scores is consistent under the same conditions. This remark records how the generated-control rate enters a standard orthogonal-score expansion (Chernozhukov et al., 2018). It is not a theorem for the nonlinear neural response learner used in the structural-response MSE benchmark. A comparable nonparametric structural-response guarantee would require additional stability conditions on the second-stage learner, such as Lipschitz or entropy control for the learned response class, or a genuinely orthogonal nonlinear score whose nuisance errors satisfy compatible productrate conditions. The benchmark therefore evaluates nonlinear second-stage performance empirically, while the formal causal error propagation in this paper is through the linear control-function and orthogonal-score statements above. The Card check below reports bootstrap intervals for estimate stability. Remark 27 (Intrinsic dimension and orthogonal inference) The rate condition ρn = o(n−1/4 ) is restrictive. Under the rate in Theorem 21, ignoring fixed-graph, boundary, and leakage terms gives ρn ≍ n−s/(2s+m) . (100) The product-rate condition for root-n orthogonal inference then requires s 1 > , 2s + m 4
equivalently
m < 2s.
(101)
Since the graph-Laplacian rate stated here uses a second-order smoother with s ≤ 2, this condition supports root-n orthogonal inference only below intrinsic dimension four, unless additional structure or stronger nuisance estimation is available. In higher intrinsic dimensions, generated-control error can dominate the score expansion. Conventional score intervals may then undercover. The finite-sample undercoverage in Section B.2 is consistent with this rate restriction. 21
Wu, Chen, Xie, Lian, and Chen
6.5 Approximate Resolvent Solves The preceding statements use the exact resolvent Sh = (I + λL(Wh ))−1 . Large graphs may use an iterative solver or a polynomial graph-filter approximation. The effect of this numerical approximation is additive. We use conjugate-gradient checks because the final linear systems are symmetric positive definite (Shewchuk, 1994; Saad, 2003). Proposition 28 (Perturbation from approximate graph solves) Let S̃h be a numerical approximation to Sh and define ṽh = (I − S̃h )x.
(102)
n−1/2 ∥(S̃h − Sh )x∥2 ≤ δn ,
(103)
n−1/2 ∥ṽh − v⋆ ∥2 ≤ n−1/2 ∥v̂h − v⋆ ∥2 + δn .
(104)
If then For a triangular array, suppose the exact smoother Sn satisfies n−1/2 ∥(I − Sn )gn ∥2 = Op (ag,n ),
(105)
n−1/2 ∥Sn vn⋆ ∥2 = Op (av,n ),
(106)
−1/2
∥vn⋆ − un ∥2 = Op (ax,n ).
(107)
n−1/2 ∥(S̃n − Sn )xn ∥2 = Op (δn ),
(108)
n Suppose also that
and that the relevance and sampling conditions in Assumption 11 hold for ṽn = (I − S̃n )xn . Then n−1/2 ∥ṽn − vn⋆ ∥2 = Op (ag,n + av,n + δn ), (109) and |β̂(ṽn ) − β0 | = Op (ag,n + av,n + ax,n + δn ) + Op (n−1/2 ).
(110)
If only the centered coefficient is needed, the same statement holds with the sampling term subtracted as in Proposition 8. Proof Since ṽh − v̂h = (Sh − S̃h )x,
(111)
the triangle inequality gives (104). The exact-smoother decomposition gives v̂n − vn⋆ = (I − Sn )gn − Sn vn⋆ .
(112)
Combining this identity with the approximate-solve error gives (109). The coefficient statement follows from Proposition 8 applied to ṽn , the relevance condition, and the sampling condition. Thus conjugate-gradient, preconditioned, or polynomial implementations preserve the reduction result when their empirical solve error δn is small. The condition is stated in the norm relevant for the generated control, not as an operator-norm requirement. 22
Graph Diffusion Residuals for Control-Function IV
7 Experiments 7.1 Protocol We evaluate first-stage residual recovery and second-stage structural response estimation. Hidden residuals are used only for evaluation or for explicitly labeled oracle references. In nonlinear synthetic tables, control correlation is computed against the outcome-relevant control U . Certificate tables report the graph residual target V ⋆ = X − g(Z) and, where needed, its distance from U . When the structural response f0 is known, we report MSEresp =
1 |Xtest |
X
{fˆ(x) − f0 (x)}2 .
(113)
x∈Xtest
We call this quantity structural-response MSE. In nonlinear synthetic experiments, fˆ is the structural-response component learned by the second-stage control-function regressor after fixing the generated control. The quantity is available in synthetic and semi-synthetic experiments, but not in the Card schooling check. These nonlinear response tables are boundarysensitive point-estimation benchmarks for the generated controls. The finite-sample coefficient bound in Section 5 is a linear control-function perturbation result; extending the same propagation uniformly to arbitrary neural second-stage learners would require additional stability conditions on the response class and training algorithm. The experiments use the quantities appearing in the theory. The benchmark asks whether the graph residual improves downstream response estimation. Mechanism experiments decompose performance into leakage, attenuation, noise mismatch, and relevance. Selection experiments examine the first-stage rule without oracle quantities. The remaining experiments replace synthetic first-stage feature geometry with empirical covariate clouds and a Card schooling stability check. The benchmark contains six synthetic designs: fractured, smooth, multi-fracture, weakinstrument, correlated-residual, and high-dimensional nuisance. Each synthetic design embeds a one-dimensional latent first-stage coordinate into an ambient observed feature vector of dimension dZ . Thus dZ = 50 is an ambient representation stress, not an intrinsicdimension claim; the latent geometry has m = 1, which lies inside the intrinsic-dimension regime discussed in Remark 27. For each design, we use n ∈ {800, 1500, 3000}, dZ ∈ {5, 20, 50}, and ten random seeds, giving 540 runs. The main fractured design uses n = 800 and dZ = 50. Additional experiments at n = 800, dZ = 50 cover first-stage cross-validation, split selection, mechanism ablations, graph certificates, and boundary stress. Appendix B reports orthogonal linear inference, meta-parameter sensitivity, component runtime, and approximate resolvent solves. The semi-synthetic real-Z experiment builds graphs from tabular covariates while simulating X, Y , U , and V ⋆ . The Card schooling check has no observed structural response and is reported as an estimate-stability study rather than as evidence for the latent-geometry theory. The fixed A-IHF defaults are K = 15, τ = 2, λ = 30, and p = 80. The observational rule uses the candidate family K ∈ {10, 15, 20},
τ ∈ {1, 2},
λ ∈ {10, 30, 50}, 23
p ∈ {70, 80, 90}.
(114)
Wu, Chen, Xie, Lian, and Chen
This gives 54 candidates per run. The graph-admissibility-filtered observational variant uses (17) with nondegenerate final degrees and ωn = 0.5, so the largest post-cut component must contain at least half of the sample. These diagnostics use only the first-stage graph. The asymptotic theory uses a local-graph sequence Kn ; the experiments use a finite candidate family of K values as part of the finite-path estimator. The trace in the GCV term is estimated by Hutchinson probes. In the main nonlinear benchmark, the second stage is the same additive neural regressor for all control-function methods. The gradient-free claim concerns only the construction of v̂ in Stage 1. The linear alignment experiment in Section 7.2.3 uses ordinary least squares. The full synthetic benchmark compares A-IHF with tuned graph ridge, graph spectral, series, kernel ridge, random forest, histogram gradient boosting, XGBoost, and neural control-function baselines. Tuning for these first-stage baselines uses only (Z, X), typically through first-stage cross-validation or graph GCV. Fixed deep and deep-ensemble controls are included as additional neural references. Appendix B.1 gives the baseline protocol and seed variation for the main fractured cell. Linear IV and DeepIV are kept as separate diagnostic references where explicitly stated, but they are not part of the tuned SOTA aggregate table. A neural control function selected by held-out first-stage prediction error for X given Z is included in the tuned baseline suite. Its selection target is predictive risk for X. 7.2 Synthetic Benchmark and First-Stage Selection 7.2.1 Main Fractured Design Table 1: Main fractured design. Selected methods from the tuned SOTA suite, averaged over 10 seeds for n = 800, dZ = 50. First-stage selection uses only (Z, X). Method
Selection
A-IHF A-IHF A-IHF Graph ridge CF Graph spectral CF Random forest CF XGBoost CF Kernel ridge CF CV-tuned Deep CF Deep ensemble CF
observational with guardrail observational fixed graph GCV graph GCV first-stage CV first-stage CV first-stage CV first-stage CV fixed ensemble
Control corr. ↑ Structural-response MSE ↓ 0.948 0.948 0.953 0.853 0.858 0.827 0.793 0.750 0.764 0.715
1.732 1.732 1.810 2.163 2.179 2.277 2.603 2.883 3.153 3.589
Table 1 reports gains in residual recovery and structural-response estimation for the fractured setting. The guarded and unguarded observational rules select the same effective candidate in this cell, have lower structural-response MSE than the fixed default, and use no hidden simulator quantities. The strongest non-A-IHF competitors are graph GCV and cross-fitted tree baselines. 24
Graph Diffusion Residuals for Control-Function IV
Main fractured design: n=800, d_Z=50
Structural-response MSE
3.5 3.0
Graph fixed KRR
HistGBDT XGBoost
2.5 RF
Graphridge spectral Graph
2.0 1.5
A-IHF fixed A-IHF obs. guardrail
0.70
0.75
0.80 0.85 0.90 Correlation with outcome-relevant control
0.95
Figure 3: Residual recovery and downstream error in the main fractured design. Guarded observational A-IHF uses only (Z, X) for first-stage selection and lies in the high-correlation, low-error region.
7.2.2 Observational Neural Cross-Validation This diagnostic compares A-IHF with a neural control function tuned by held-out first-stage prediction error. The full benchmark below includes CV-tuned Deep CF; this smaller table isolates the prediction-tuning comparison at n = 800, dZ = 50. Table 2: First-stage cross-validation. Entries report mean structural-response MSE, with mean correlation against U in parentheses, over 10 seeds for n = 800, dZ = 50. The CV-tuned neural control function selects architecture, weight decay, learning rate, and stopping epoch by held-out prediction error for X given Z. Design
A-IHF obs. A-IHF fixed CV-tuned Deep CF Unreg. Deep CF
Fractured 1.668 (0.948) High-dimensional nuisance 2.270 (0.767) Smooth 1.035 (0.971)
1.779 (0.953) 2.789 (0.637) 0.822 (0.977)
3.098 (0.756) 4.881 (0.647) 0.787 (0.971)
4.433 (0.652) 10.867 (0.093) 9.395 (0.528)
Table 2 separates predictive first-stage tuning from fixed regularization. Cross-validation improves the neural control function relative to unregularized training. Its target is held-out prediction error for X given Z. Corollary 9 instead involves ∥v̂−u∥2 and κn (v̂). In fractured and high-dimensional nuisance designs, predictive cross-validation does not close the gap to 25
Wu, Chen, Xie, Lian, and Chen
A-IHF. In the smooth design, the cross-validated neural first stage is competitive with fixed A-IHF. 7.2.3 Linear Second-Stage Alignment The perturbation bound in Corollary 9 is stated for a linear control-function second stage. To align the experiment with that setting, we reuse the same first-stage designs and generate a linear outcome Ylin = β0 X + γ0 U + ϵ, β0 = 1, γ0 = 2.5. (115) For each generated control V̂ , we estimate β0 by ordinary least squares on (1, X, V̂ ). Table 3 reports the mean absolute coefficient error over 10 seeds. Table 3: Linear second-stage alignment. Entries are mean absolute errors |β̂ − β0 | for Ylin = β0 X + γ0 U + ϵ, n = 800, dZ = 50. Design
2SLS CV Deep CF Regularized Deep CF A-IHF fixed A-IHF obs.
Fractured High-dimensional nuisance Smooth Correlated residual
0.048 0.096 0.262 0.447
0.085 0.172 0.146 0.421
0.061 0.056 0.065 0.485
0.054 0.419 0.162 0.398
0.031 0.253 0.085 0.409
The fractured design aligns with the mechanism in Corollary 9. Among the methods in Table 3, observational A-IHF has the smallest coefficient error, correlation 0.953 with U , and residualized treatment variation κn = 5.555. In the high-dimensional nuisance design, observational tuning improves fixed A-IHF, but the regularized neural baseline and 2SLS have smaller coefficient error. Control correlation alone is not a sufficient summary of secondstage behavior. A-IHF has higher correlation with U than the CV-tuned neural control function (0.771 versus 0.651), but leaves less residualized treatment variation (κn = 1.950 versus 2.864). This is the denominator in (30). The full 54-cell SOTA table reports a broader linear coefficient diagnostic separately. In the correlated-residual design, all methods have large coefficient error, as predicted by the spectral-leakage analysis. Appendix B.2 reports a finite-sample orthogonal-score check. 7.2.4 Full Benchmark Across all benchmark cells, guarded observational A-IHF has the lowest average structuralresponse MSE among the tuned first-stage baselines. It is the best method in 24 of 54 cells and appears in the top three in 36 cells. Fixed A-IHF is best in 8 cells, and the A-IHF family beats the best non-A-IHF method in 32 of 54 cells. The last column shows that this advantage is not a uniform statement about every diagnostic: graph spectral, ensemble deep, XGBoost, and series baselines have lower average error in the linear coefficient diagnostic. The two rightmost uses of the generated control should therefore be read separately. The structural-response columns evaluate a nonlinear second-stage point-estimation task. The linear column probes a low-capacity coefficient map for which the perturbation bound in Corollary 9 is directly relevant. Linear coefficient error is not a direct measurement of the global L2 distance between V̂ and U : it also depends on the direction of the residual 26
Graph Diffusion Residuals for Control-Function IV
Table 4: Aggregate results over 54 benchmark cells. Each cell is averaged over 10 seeds; the table reports summaries of those cell averages. Structural-response MSE is the main nonlinear benchmark metric. The last column is a linear coefficient diagnostic. Method A-IHF, guarded observational A-IHF, observational A-IHF, fixed Graph ridge CF Random forest CF Graph spectral CF Kernel ridge CF XGBoost CF Deep ensemble CF Series CF CV-tuned Deep CF Fixed Deep CF HistGBDT CF Graph ridge CF, fixed
Mean response MSE ↓ Median response MSE ↓ Median corr. ↑ Mean linear abs. error ↓ 2.676 2.688 2.772 2.912 2.973 3.049 3.148 3.219 3.385 3.417 3.444 3.537 3.606 3.822
2.533 2.533 2.710 2.756 2.825 3.168 3.034 3.429 3.545 3.622 3.472 3.719 3.566 3.187
0.904 0.904 0.872 0.839 0.826 0.850 0.781 0.798 0.740 0.741 0.762 0.728 0.779 0.729
0.203 0.204 0.266 0.142 0.148 0.115 0.138 0.120 0.120 0.121 0.134 0.136 0.186 0.390
error, the projection geometry, and the relevance denominator in (30). Table 4 therefore shows a diagnostic tradeoff rather than uniform dominance. Smoother spectral and series controls can be stable for OLS projections, while guarded anisotropic residuals are more effective in the nonlinear fractured-response benchmark. The theory supports the residualerror, relevance, and orthogonal-score mechanisms; the nonlinear response advantage is an empirical comparison of downstream generated controls. The nonlinear benchmark has a simple stability interpretation. If a second-stage response map is locally Lipschitz in the generated-control coordinate, then replacing U by V̂ changes the oracle response by an amount proportional to ∥v̂ − u∥2 on the evaluation distribution. In fractured designs, however, a global average norm can hide boundary-local errors: a smoother may have low OLS coefficient error while still mixing observations across a jump in g(Z). Such cross-boundary contamination is precisely the structural-leakage term in Proposition 10. A-IHF should therefore be read as trading some global smoothness for boundary-sensitive control recovery, not as proving a uniformly smaller residual norm. Under algorithmic stability of the second-stage learner, the same residual error enters the learned response up to the learner’s estimation error. The certificate terms are therefore relevant to nonlinear point estimation, while all nonlinear rows use the same second-stage architecture and training schedule to keep optimization and architecture effects aligned across methods. 7.2.5 Performance by Design A-IHF is best on average in fractured, multi-fracture, high-dimensional nuisance, and correlatedresidual designs. Its relative margin is clearest in fractured and high-dimensional nuisance designs. In smooth designs, standard smooth or neural baselines are competitive or better. In weak-instrument designs, graph ridge performs best on average. In correlated-residual designs, residual recovery remains difficult because the residual has low-frequency graph energy, as predicted by (25); fixed A-IHF is slightly better than the guarded observational variant there. 27
Wu, Chen, Xie, Lian, and Chen
Table 5: Mean structural-response MSE by data-generating design. Values average over n, dZ , and seeds. Design
A-IHF guardrail A-IHF fixed Best non-A-IHF Best non-A-IHF method
Fractured Multi-fracture High-dimensional nuisance Smooth Weak instrument Correlated residual
A-IHF guardrail
3.528 2.045 3.264 1.162 0.806 5.252
3.718 2.299 3.470 1.198 0.738 5.212
A-IHF fixed
Best non-A-IHF
4.228 2.107 3.868 1.000 0.619 5.450
Graph ridge/spectral CF Graph ridge CF RF CF Deep ensemble CF Graph ridge CF Graph spectral CF
Mean structural-response MSE
5 4 3 2 1 0
ured Fract
ure
-fract
Multi
dim High-
nce
th
Smoo
nuisa
idual ment d res e t a l e Corr
instru Weak
Figure 4: Mean structural-response MSE by data-generating design. The comparison separates guarded observational A-IHF, fixed A-IHF, and the best non-AIHF baseline in each design. A-IHF has its largest relative gains in fractured and high-dimensional nuisance designs. It is not uniformly best in smooth or weakinstrument settings.
7.3 Mechanism, Selection, and Failure Modes 7.3.1 Selection and Graph-Admissibility Checks In the 540-run benchmark, all observational selections satisfy the relevance screen in (16). The unguarded observational rule selects p = 70 in all runs and K = 10 in 529 of 540 runs. The guarded rule selects K = 10 in 487 runs, K = 15 in 37 runs, and K = 20 in 16 runs; it selects p = 70 in 480 runs and keeps the median largest post-cut component fraction at 1.0. The criterion favors sharp graphs inside the pre-specified admissible family, while the guardrail prevents the most fragmented choices. 28
Graph Diffusion Residuals for Control-Function IV
An extended-grid experiment uses K ∈ {5, 10, 15, 20},
p ∈ {50, 60, 70, 80, 90},
(116)
and a fixed response evaluation grid Xtest ∈ [−5, 5]. In this experiment, the rule selected p = 50 in all 540 runs and K = 5 in 527 runs. Fixed A-IHF had average structural-response MSE 3.791, while the observational rule over the extended grid had average structuralresponse MSE 5.878. This is an over-fragmentation failure of the graph family. The extended-grid experiment is repeated with graph-admissibility constraints: the final graph must have mean degree at least 4, and the number of post-cut connected components √ must be no larger than n. This guardrail moves the selected graph away from the sparsest candidate: K = 15 is selected in 488 of 540 runs, the median number of post-cut components is 1, and the median largest-component fraction is 1.0. The average structural-response MSE of observational A-IHF decreases from 5.878 to 4.521 under the same fixed-grid evaluation. Table 6 reports the comparison. The guardrail is not a reversal of the anisotropic mechanism. Anisotropy is meant to reduce conductance across large pilot-stage treatment jumps while preserving the macroscopic neighborhood graph on each smooth region. If a data-adaptive candidate cuts the graph into many small components, the procedure is no longer separating structural boundaries; it is isolating finite-sample noise. The degree and component checks are therefore graph-admissibility filters. They rule out shattering while still allowing local cross-boundary conductance to be downweighted. The main SOTA suite uses the graph-admissibility filter defined in Section 4.4, with largest-component threshold ωn = 0.5. This value encodes a benchmark prior that the tested fractured designs contain a dominant macroscopic geometry after legitimate boundary attenuation. It is not a statistical constant and is not used in the graph-limit theory. In applications with several comparable disconnected regimes, such as three equally sized subpopulations, this screen should be replaced by a domain-appropriate component-size rule or removed and reported as a graph-certificate diagnostic. The estimator is therefore the candidate family together with its admissibility filters, not an unrestricted search over arbitrarily shattered graphs. Table 6: Fixed-grid graph-admissibility check. Values average over the same 54 benchmark cells and use Xtest ∈ [−5, 5]. Procedure A-IHF, fixed A-IHF, observational broad grid A-IHF, observational with graph guardrail
Mean structural-response MSE ↓ Selected graph summary 3.791 5.878 4.521
fixed K = 15, p = 80 K = 5 in 527/540; p = 50 in 540/540 K = 15 in 488/540; median components 1
The guardrail improves the broad grid but leaves the fixed default better on this fixedgrid evaluation. Graph connectivity constraints reduce one failure mode of observational selection. The main observational results use a pre-specified admissible family. Broader grids are treated as a graph-admissibility problem. 29
Wu, Chen, Xie, Lian, and Chen
7.3.2 Split-Selection Check Proposition 24 concerns a split selector with fixed candidate smoothers on an evaluation split. The implemented split experiment has a narrower purpose. One split selects A-IHF parameters from (Z, X), and the selected parameters are then used to refit the transductive graph smoother on the full sample. The experiment measures selection stability and computational cost. Table 7: Split-selection check. Entries report mean nonlinear structural-response MSE over 10 seeds for n = 800, dZ = 50. The split-selected row selects parameters on a subsample and refits the graph smoother on the full sample. The last column is the mean runtime ratio relative to full observational selection. Design
A-IHF full obs. A-IHF split-selected Response-MSE difference Runtime ratio
Fractured High-dimensional nuisance Smooth Weak instrument Correlated residual
1.263 2.059 0.955 0.719 6.385
1.264 2.067 1.384 1.114 6.385
0.001 0.008 0.429 0.395 0.000
0.42 0.29 0.42 0.41 0.46
Table 7 reports similar results for split-selected and full observational parameters in fractured, high-dimensional nuisance, and correlated-residual designs. The split-selected rule is less reliable in smooth and weak-instrument designs. The selected percentile is p = 70 in all 50 split runs. The split selector chooses K = 10 in 44 of 50 runs, λ = 10 in 43 runs, and τ = 1 in 28 runs. 7.3.3 Mechanism Ablation We isolate the main graph operations at n = 800, dZ = 50. The ablations use the same seeds and downstream learner as the split-selection experiment. Table 8: A-IHF mechanism ablation. Entries report mean nonlinear structural-response MSE over 10 seeds. The full observational rule uses the 54-candidate family. Design Fractured High-dimensional nuisance Smooth Weak instrument Correlated residual
Fixed A-IHF Isotropic No pilot Obs. full Obs. no roughness 1.308 2.556 0.846 0.606 5.990
2.583 2.164 0.963 0.902 9.197
1.766 2.442 0.877 0.836 5.952
1.263 2.059 0.955 0.719 6.385
1.266 2.041 1.367 1.568 6.385
The ablation isolates the role of anisotropy and the pilot stage. In the fractured design, the isotropic variant has lower control correlation (0.727 versus 0.953) and higher structuralresponse MSE. Removing the pilot step raises structural-response MSE from 1.308 to 1.766. The hard conductance cutoff has negligible effect in this grid; its structural-response MSE is 1.308 in the fractured design and 2.556 in the high-dimensional nuisance design, matching fixed A-IHF to three decimals. It is a sparsification device rather than a source of statisti30
Graph Diffusion Residuals for Control-Function IV
cal performance. Removing the roughness penalty leaves fractured and correlated-residual selections unchanged in this grid, but worsens smooth and weak-instrument designs. 7.3.4 Representation Check The representation experiment varies the graph input and addresses Assumption 13. Table 9 reports observational A-IHF under an oracle latent representation, the observed firststage feature representation, principal-component representations, and observed features with additional irrelevant noise. Appendix B.3 reports the sensitivity of the observational score to α and cκ . Table 9: Representation stress test. Mean structural-response MSE of observational AIHF over 10 seeds for n = 800, dZ = 50. The oracle latent representation is a simulation reference and is not used in the main benchmark. Design Fractured High-dimensional nuisance
Latent oracle Observed Z PCA-5 PCA-10 Noisy Z Noisy PCA-5 1.689 1.143
1.767 2.469
2.070 2.637
3.507 2.805
3.190 2.903
3.301 3.213
The latent representation has the lowest structural-response MSE. The observed representation remains usable, while degraded representations can be much worse. The fixed A-IHF row is unstable for PCA-10 in this experiment, with mean structural-response MSE 22.862 in the fractured design and 48.631 in the high-dimensional nuisance design; observational selection mitigates this representation failure. A-IHF is a residual extractor conditional on a graph representation. 7.3.5 Graph Certificate and Boundary Stress The preceding experiments use downstream response error. We also report the certificate terms in Proposition 10. The quantities are normalized root-mean-square terms for the graph residual target V ⋆ = X − g(Z). They are observable only in simulation because they use g and V ⋆ . Table 10 separates graph leakage, residual attenuation, and relevance. In the correlatedresidual design, leakage is small but residual attenuation is close to one. The graph smoother removes little of the low-frequency residual component. In the high-dimensional nuisance design, observational selection lowers leakage and raises relevance relative to the fixed default. In weak-instrument designs, correlations remain high while κn is small. Residual fidelity and second-stage relevance are reported together. The next experiment stresses the graph representation and first-stage relevance. For each stress, Table 11 reports averages over five designs and ten seeds. The endpoint values are the low and high settings in each path. For instrument strength, the left endpoint is the weak-instrument case. The stress tests separate two failure modes. Irrelevant coordinates and isotropic representation noise raise the certificate bound and lower residualized treatment variation. Observational selection reduces the degradation in several cases, including RMSE and κn . Dependence on a geometry-preserving representation remains. Weak instruments show a 31
Wu, Chen, Xie, Lian, and Chen
Table 10: Graph-admissibility certificate. Mean over 10 seeds for n = 800, dZ = 50. Correlation and RMSE are computed against V ⋆ = X − g(Z). The bound column is leakage plus residual attenuation for the graph residual target; the coefficient certificate also includes first-stage noise mismatch. Design
Method
Fractured A-IHF obs. Fractured A-IHF fixed Smooth A-IHF obs. Smooth A-IHF fixed Weak instrument A-IHF obs. Weak instrument A-IHF fixed Correlated residual A-IHF obs. Correlated residual A-IHF fixed High-dimensional nuisance A-IHF obs. High-dimensional nuisance A-IHF fixed
Corr. RMSE Leakage Atten. Bound
κn
0.953 0.957 0.977 0.982 0.961 0.961 0.170 0.287 0.771 0.639
5.555 5.232 0.909 0.848 0.563 0.488 7.091 7.016 1.950 1.017
0.372 0.291 0.247 0.207 0.292 0.281 0.993 0.971 0.656 1.077
0.133 0.236 0.110 0.131 0.212 0.249 0.021 0.046 0.629 1.100
0.328 0.158 0.218 0.162 0.215 0.162 0.992 0.970 0.242 0.093
0.461 0.394 0.327 0.292 0.427 0.411 1.013 1.016 0.871 1.194
Table 11: Boundary stress tests. Values average over five designs and ten seeds at n = 800, dZ = 50. Correlation and RMSE are computed against V ⋆ = X − g(Z). Each cell reports the low-endpoint value followed by the high-endpoint value. For instrument strength, the endpoint path is 0.2 → 1.0; for noise dimensions, 0 → 200; for representation noise, 0 → 2. Stress path
Method
Corr.
RMSE
Bound
κn
Instrument strength A-IHF obs. 0.820 → 0.766 0.379 → 0.512 0.451 → 0.620 0.321 → 3.214 Instrument strength A-IHF fixed 0.830 → 0.765 0.367 → 0.566 0.437 → 0.661 0.310 → 2.920 Irrelevant dimensions A-IHF obs. 0.766 → 0.545 0.512 → 1.129 0.620 → 1.400 3.214 → 0.684 Irrelevant dimensions A-IHF fixed 0.765 → 0.564 0.566 → 1.383 0.661 → 1.543 2.920 → 0.174 Representation noise A-IHF obs. 0.766 → 0.538 0.512 → 1.251 0.620 → 1.514 3.214 → 0.215 Representation noise A-IHF fixed 0.765 → 0.556 0.566 → 1.493 0.661 → 1.628 2.920 → 0.028
different limitation. At instrument strength 0.2, residual correlations remain high, but κn is near zero. The denominator in Corollary 9 is then the limiting term. 7.4 Empirical Covariate Geometry and Stability Checks 7.4.1 Semi-Synthetic Real-Z Experiment We replace the synthetic first-stage feature cloud with real tabular covariates from the diabetes, breast_cancer, and digits data sets. Outcomes remain simulated. For each data set, we generate fractured, smooth, and weak-instrument first-stage designs with n = 400 and ten seeds. This keeps U , V ⋆ , and the structural response observable while testing the graph rule on empirical covariate geometry. The nonlinear response uses U as the outcome-relevant control. The coefficient column in Table 12 uses a linear experiment with 32
Graph Diffusion Residuals for Control-Function IV
V ⋆ as the control target, so it measures residual-target recovery rather than the same linear design reported in Section 7.2.3. Table 12: Semi-synthetic real-Z experiment. Mean over three data sets, three simulated designs, and ten seeds. The covariates are real; X, Y , U , and V ⋆ are simulated.
Method
Structural-response MSE ↓ |β̂ − β0 | ↓ Control corr. ↑
Kernel ridge CF A-IHF, observational Random forest CF CV-tuned Deep CF 2SLS Regularized Deep CF Gradient boosting CF A-IHF, fixed Unregularized Deep CF
6.094 6.760 7.469 9.866 10.598 10.886 12.295 17.920 30.043
0.248 0.360 0.236 0.455 0.446 0.497 0.622 0.977 1.394
0.884 0.860 0.860 0.853 0.805 0.839 0.823 0.821 0.140
Table 12 is less favorable to A-IHF than the synthetic benchmark. Observational AIHF ranks second in structural-response MSE and improves the fixed A-IHF default by 11.160 points on average. Kernel ridge has the lowest structural-response MSE and the highest control correlation. Random forest has the smallest mean coefficient error. The unregularized neural first stage remains unstable. Table 13: Semi-synthetic real-Z results by simulated design. Entries are mean nonlinear structural-response MSE.
Design Real-Z fractured Real-Z smooth Real-Z weak instrument
A-IHF obs. KRR CF RF CF CV Deep CF 1.968 4.863 13.448
2.198 5.420 10.665
2.601 8.087 11.719
3.293 8.051 18.253
Table 13 separates the three simulated designs. Observational A-IHF is best on the fractured and smooth real-Z designs, but not under weak instruments. In the weak-instrument design, its selected residual has high control correlation (0.913) and low residualized treatment variation (κn = 0.115). This matches the denominator in (30): a high-fidelity generated control is insufficient if little treatment variation remains after conditioning on it. The selection pattern resembles the synthetic benchmark. The observational rule selects p = 70 in all 90 semi-synthetic runs, K = 10 in 89 runs, λ = 10 in 55 runs, and λ = 30 in 35 runs. The admissible family remains part of the estimator. 33
Wu, Chen, Xie, Lian, and Chen
7.4.2 Card Schooling Tabular Stability Check We use the Card schooling data as a tabular real-IV stability check (Card, 1995; Wooldridge, 2010). The outcome is log wage, the treatment is years of education, and the excluded instrument is proximity to a four-year college. The controls are experience, experience squared, race, southern residence, metropolitan indicators, and 1966 region indicators. These variables are mostly tabular and partly discrete, so the Euclidean K-nearest-neighbor graph is not expected to resemble the connected local-geometry regime used in the rate theory. The graph for A-IHF and the nonparametric control-function baselines is built from the standardized excluded instrument and controls. We report the original-sample coefficient on education and a nonparametric bootstrap with 200 resamples. The structural response is unobserved, so this experiment reports estimate stability and graceful degradation on a difficult tabular geometry rather than structural-response MSE or graph-admissibility evidence.
Table 14: Card schooling tabular stability check. Coefficient on education in a logwage equation. Intervals are the 2.5% and 97.5% bootstrap quantiles over 200 resamples. Method
Selection
Original Bootstrap mean Bootstrap interval
OLS 2SLS Linear CF KRR CF RF CF GBDT CF A-IHF, fixed A-IHF, observational
observational observational observational observational observational observational fixed observational
0.0747 0.1315 0.1315 0.0691 0.0626 0.0682 0.0677 0.0650
0.0748 0.1394 0.1394 0.0716 0.0677 0.0685 0.0704 0.0731
[0.0678, 0.0828] [0.0409, 0.2810] [0.0409, 0.2810] [0.0370, 0.1039] [0.0422, 0.0900] [0.0414, 0.0890] [0.0551, 0.0837] [0.0596, 0.0843]
Table 14 reports a conservative estimate relative to 2SLS. Observational A-IHF gives an original-sample estimate of 0.0650 and a bootstrap interval [0.0596, 0.0843], close to OLS and the nonlinear control-function baselines. The 2SLS estimate is larger and has a wider bootstrap interval. The first-stage F statistic is 13.25 in the original sample; across bootstrap rows its median is 13.31 and its 2.5% and 97.5% quantiles are 3.60 and 31.16. The A-IHF observational rule selects p = 70, τ = 1, and λ = 10 in all 201 samples, and selects K = 10 in 200 of 201 samples. The mean bootstrap value of κn (v̂) is 1.500 for observational A-IHF and 1.155 for fixed A-IHF. The selected graphs are often fragmented, with typical largestcomponent fractions around 0.1–0.2. This is a diagnostic of stability under an unfavorable tabular graph, not a positive example of the geometric assumptions. These graphs do not resemble the connected local-geometry regime used in the rate theory. The Card results should not be read as evidence for graph admissibility. They report estimate stability under the implemented observational rule. 34
Graph Diffusion Residuals for Control-Function IV
Card schooling application OLS 2SLS Linear CF KRR CF RF CF GBDT CF A-IHF fixed Original sample Bootstrap mean and 95% interval
A-IHF obs. 0.05
0.10
0.15
0.20
0.25
Coefficient on education
Figure 5: Card schooling bootstrap intervals. Filled circles show bootstrap means and open circles show original-sample estimates. A-IHF gives an estimate close to other nonlinear control-function baselines on this tabular stability check, while the linear IV estimate has a wider interval.
7.5 Numerical and Evaluation Notes Appendix B.4 reports first-stage component timings. Appendix B.5 reports conjugategradient checks for the final graph resolvent. The empirical solve error is below 5 × 10−5 at n = 10000 in the tested designs. Appendix B.6 records the response-grid convention used in the benchmark.
8 Limitations A-IHF is representation-conditioned. This is both a limitation and a diagnostic statement. The graph must encode neighborhoods relevant to g(Z); when it does not, the graph certificate and downstream checks reveal that the generated control is not geometrically supported. Boundary and representation experiments show the mechanism: irrelevant coordinates and representation noise raise structural leakage and reduce residualized treatment variation. For text, images, genomic measurements, or other unstructured first-stage features, a representation-learning step or a domain metric is needed before graph construction. This is the same premise used by spectral and manifold estimators. A-IHF is a residualextraction module once a manifold-respecting representation is available. The method also requires separation between the graph spectrum of g(Z) and that of Low-frequency residual variation is retained by the smoother and appears as residual attenuation. The correlated-residual experiments illustrate this case. Weak instruments V ⋆.
35
Wu, Chen, Xie, Lian, and Chen
create a different limitation: the generated control can be close to the target while κn (v̂) remains small. The selection theory is finite-path. Proposition 24 compares a split selector with the best candidate in the chosen admissible family. Minimax optimality would require a continuum of graph scales, smoothness classes, and a risk approximation for adaptive graphs. Broad candidate families can over-fragment the graph, so the graph family is part of the estimator. The observational score is best understood as a geometrically principled admissibility filter. Its GCV component is a first-stage smoothness proxy, and Proposition 24 contains an explicit selection-bias term that delineates when this proxy approximates residual-recovery risk. The term is expected to be small under fixed or stable candidate graphs, approximately homoskedastic first-stage residuals, and weak alignment between residual variation and lowfrequency graph modes. It can be persistent under heteroskedasticity or graph-correlated residuals. The method therefore supplies a finite-path observational rule with diagnostic failure modes rather than an oracle selector. Theorem 21 gives rates under latent piecewise-smooth geometry and graph approximation conditions. It abstracts the baseline spectral convergence of the separated graph and isolates the perturbation from anisotropic conductance. A full continuum limit for percentileselected, pilot-dependent graph diffusion remains a separate graph-limit problem. The inference results have the usual nonparametric rate bottleneck. Root-n orthogonal inference requires o(n−1/4 ) control error, which restricts the intrinsic dimension under second-order graph smoothing. The nonlinear structural-response tables use a neural second stage and are interpreted as empirical point-estimation comparisons under a common response learner. The formal propagation result in Section 5 is for a linear control-function coefficient, with a standard orthogonal-score route under product-rate conditions. Extending the result to nonlinear neural second stages requires additional stability assumptions on the response learner or an orthogonal nonparametric score tailored to the structural target. The semi-synthetic real-Z experiment simulates treatment, residuals, and outcomes. It tests empirical covariate geometry. The Card schooling check has no observed structural response or true generated control, and its discrete tabular features often produce fragmented Euclidean graphs; its intervals measure estimate stability under an unfavorable representation rather than graph admissibility. This degradation should not be read only as an algorithmic failure. It is evidence for the main premise: control-function residuals are not extractable from first-stage prediction alone, and a representation with the wrong local topology can make graph-based residual extraction unreliable. Exact sparse solves are feasible at the sample sizes used here. Larger samples require monitoring the generated-control solve error δn together with relevance κn (v̂).
9 Conclusion This paper treats the control-function first stage as a residual-extraction problem. A-IHF uses anisotropic graph diffusion to estimate the first-stage systematic treatment component and uses the remaining signal as a generated control. The finite-sample theory identifies three quantities that govern the second stage: structural leakage, residual attenuation, and 36
Graph Diffusion Residuals for Control-Function IV
residualized treatment variation. The asymptotic theory gives rates when the graph family is admissible and the representation has piecewise-smooth geometry. The experiments are consistent with this interpretation. In the 54-cell tuned SOTA benchmark, guarded observational A-IHF has the lowest average structural-response MSE, and the A-IHF family beats the best non-A-IHF method in 32 cells. Its clearest margins occur in fractured and high-dimensional nuisance designs, where graph cuts reduce cross-boundary leakage. It remains best on average in multi-fracture and correlated-residual designs, although correlated-residual errors remain large in absolute terms. It is less distinctive in smooth and weak-instrument designs, and it is not the leading method for the linear coefficient-error diagnostic. Mechanism checks, graph certificates, representation stress tests, semi-synthetic real-Z experiments, and the Card tabular stability check point to the same limitations. A-IHF is a non-oracle first-stage baseline for settings where graph geometry is informative and the reported quantities are stable.
Appendix A. Proofs A.1 Proof of Proposition 1 Proof [Proof of Proposition 1] Under continuous-time gradient flow in the NTK regime, d ft (Z) = Σ(x − ft (Z)). dt
(117)
d v̂N N (t) = −Σv̂N N (t). dt
(118)
Let v̂N N (t) = x − ft (Z). Then
Expanding in the eigenbasis of Σ gives d ⟨v̂N N (t), ϕj ⟩ = −µj ⟨v̂N N (t), ϕj ⟩. dt
(119)
Solving this scalar differential equation yields the stated expression. Since µj > 0 for every j, each coefficient converges to zero.
A.2 Proof of Proposition 7 Proof Since L is symmetric positive semidefinite, write L=
n X
µj ϕ j ϕ ⊤ j .
(120)
j=1
Then 2
−2
S = (I + λL)
=
n X j=1
37
(1 + λµj )−2 ϕj ϕ⊤ j .
(121)
Wu, Chen, Xie, Lian, and Chen
Using E[v⋆ (v⋆ )⊤ ] = ΣV , E∥Sv⋆ ∥22 = E[(v⋆ )⊤ S 2 v⋆ ]
(122) (123)
2
= Tr(S ΣV ) n X ϕ⊤ j ΣV ϕj . = (1 + λµj )2
(124)
j=1
A.3 Proof of Proposition 8 Proof Let Mv̂ project off the intercept and v̂. The coefficient on x in the regression of y on (1, x, v̂) is x⊤ Mv̂ y β̂(v̂) = ⊤ . (125) x Mv̂ x Substitute y = β0 x + γ0 u + ϵ: β̂(v̂) − β0 = γ0 Because Mv̂ v̂ = 0,
x⊤ Mv̂ ϵ x⊤ Mv̂ u + . x⊤ Mv̂ x x⊤ Mv̂ x
x⊤ Mv̂ u = x⊤ Mv̂ (u − v̂).
(126) (127)
This proves (27). By Cauchy–Schwarz, |x⊤ Mv̂ (u − v̂)| ≤ ∥Mv̂ x∥2 ∥u − v̂∥2 .
(128)
Since ∥Mv̂ x∥22 = x⊤ Mv̂ x = nκn (v̂), the generated-control term is bounded by (28). The remaining term is the sampling term.
A.4 Proof of Lemma 16 Proof On the event in Assumption 15, for an edge (i, j) within one smooth region, |x̃i − x̃j | ≤ |g(Ti ) − g(Tj )| + 2rp,n ≤ Lg εsn + 2rp,n . Thus
(Lg εsn + 2rp,n )2 1 ≥ Cij ≥ exp − γn For a cross-jump edge satisfying (55),
→ 1.
|x̃i − x̃j | ≥ |g(Ti ) − g(Tj )| − 2rp,n ≥ ∆ − 2rp,n . Therefore
(129) (130)
(131)
(∆ − 2rp,n )2 Cij ≤ exp − → 0. (132) γn The bounds are uniform over the stated edge sets on the same high-probability event.
38
Graph Diffusion Residuals for Control-Function IV
A.5 Proof of Lemma 17 Proof For within-region edges, (64) gives Cij = exp{−Dij /γ̂n } ≥ e−cγ
(133)
on the stated high-probability event. For cross-jump edges, (65) and (66) give Cij = exp{−Dij /γ̂n } ≤ exp{−bn /γ̂n } → 0.
(134)
The two bounds are uniform over the corresponding edge sets.
A.6 Proof of Lemma 19 Proof The resolvent identity gives Sn − Sn0 = (I + λn Ln )−1 λn (L0n − Ln )(I + λn L0n )−1 .
(135)
Since Ln and L0n are symmetric positive semidefinite, both resolvents have operator norm at most one. Hence Sn − Sn0 op ≤ λn Ln − L0n op = ηn . (136) Multiplying by n−1/2 ∥xn ∥2 = Op (1) gives (73). 0 For the degree bound, write ∆W = Wn − Wn0 and P ∆D = Dn − Dn . The operator norm of ∆D − ∆W is bounded by P a constant times maxi j |∆Wij | by Gershgorin’s theorem and the bound ∥∆W ∥op ≤ maxi j |∆Wij | for symmetric matrices. The normalization gives Dn − Wn Dn0 − Wn0 − . (137) d¯n d¯0n P The first difference is controlled by maxi j |∆Wij |/d¯0n . The second is controlled by |d¯n − d¯0n |/d¯0n and the operator-norm bound in (74). Condition (74) gives Ln − L0n op ≤ Cℓn . Ln − L0n =
A.7 Proof of Theorem 21 Proof By Proposition 6,
v̂n − vn⋆ = gn − Sn xn .
(138)
Add and subtract Sn0 xn : v̂n − vn⋆ = (I − Sn0 )gn − Sn0 vn⋆ + (Sn0 − Sn )xn . √ Taking norms and dividing by n gives n−1/2 ∥v̂n − vn⋆ ∥2 ≤ n−1/2 ∥(I − Sn0 )gn ∥2 + n−1/2 ∥Sn0 vn⋆ ∥2 + n−1/2 ∥(Sn0 − Sn )xn ∥2 .
(139) (140) (141)
The first term is controlled by (70), the second by (76), and the third by Lemma 19. This s/2 proves (79). The consistency statement follows from the displayed rate. Balancing tn with −m/4 n−1/2 tn gives tn ≍ n−2/(2s+m) and rate n−s/(2s+m) .
39
Wu, Chen, Xie, Lian, and Chen
A.8 Proof of Corollary 22 Proof Theorem 21 gives n−1/2 ∥v̂n − vn⋆ ∥2 = Op (rn ).
(142)
Also, n−1/2 ∥v̂n − un ∥2 ≤ n−1/2 ∥v̂n − vn⋆ ∥2 + n−1/2 ∥vn⋆ − un ∥2 = Op (rn + δx,n ).
(143)
On the event κn (v̂n ) ≥ κ0 , Proposition 8 bounds the generated-control contribution by a constant multiple of this term. This gives (84). The last display follows when Eϵ,n = Op (n−1/2 ).
A.9 Proof of Proposition 24 Proof Let ĥ = ĥobs . On the event (93), for any h ∈ Hnκ , b obs (ĥ) + ∆n Rn (ĥ) + αJn (ĥ) + σV2 ≤ Q b obs (h) + ∆n ≤Q ≤ Rn (h) + αJn (h) + σV2 + 2∆n .
(144) (145) (146)
Taking the infimum over h ∈ Hnκ and dropping the nonnegative term αJn (ĥ) gives (94). It remains to justify (96). Work conditional on the graph-construction split and on Z, so the matrices Sh are fixed. Set (I − Sh )2 Ah = . (147) (1 − dh )2 The degrees-of-freedom condition (92) gives √ ∥Ah ∥F ≤ c−2 n. df
∥Ah ∥op ≤ c−2 df ,
(148)
Since 0 ⪯ Sh ⪯ I, the same type of bounds hold for Sh2 and for Ah − Sh2 . Write x = g + η. The part of b obs (h) − {Rn (h) + αJn (h) + σ 2 } Q V that is random after centering is o 2 1n ⊤ η Bh η − E[η ⊤ Bh η | Z] + a⊤ η, n n h
(149)
where Bh = Ah − Sh2 ,
ah = {Ah + (I − Sh )Sh } g.
(150)
The deterministic remainder is included in bsel,n . The roughness term Jn (h) cancels because it appears in the empirical criterion and in the comparison target with the same coefficient α. 40
Graph Diffusion Residuals for Control-Function IV
By the sub-Gaussian quadratic-form inequality, for every u > 0, √ P η ⊤ Bh η − E[η ⊤ Bh η | Z] > C{∥Bh ∥F u + ∥Bh ∥op u} Z ≤ 2e−u . Using (148) gives a contribution bounded by r u u . C + n n The linear term satisfies
√ ∥ah ∥2 ≤ C∥g∥2 ≤ C n,
so the sub-Gaussian linear-form bound gives r 2 ⊤ u P Z ≤ 2e−u . |a η| > C n h n
(151)
(152)
(153)
(154)
Take u = log(2|Hn |/δ) and apply a union bound over Hnκ . Adding bsel,n proves (96). If the trace is estimated, the map d 7→ (1 − d)−2 is Lipschitz on d ≤ 1 − cdf /2, so the perturbation from replacing dh by dbh is bounded by Cetr,n .
A.10 Derivation for Remark 26 Proof Use a first-order expansion of the cross-fitted sample moment around (β0 , η0 ): n
1X ψ(Wi ; β0 , η0 ) + G(β̃ − β0 ) + Rn . 0= n
(155)
i=1
Cross-fitting removes first-order empirical-process terms from estimating η0 . Neyman orthogonality makes the deterministic first-order nuisance derivative zero. The assumed product-rate condition gives Rn = op (n−1/2 ). Hence √
n
1 X n(β̃ − β0 ) = −G−1 √ ψ(Wi ; β0 , η0 ) + op (1). n
(156)
i=1
The central limit theorem gives (99). Consistency of the sandwich estimator follows from cross-fitted plug-in consistency and the law of large numbers.
Appendix B. Additional Experimental Details The synthetic first stage is generated from a one-dimensional latent instrument embedded in dZ observed dimensions with sinusoidal random projections and Gaussian noise. The intrinsic dimension of the synthetic geometry is therefore m = 1; the reported values dZ ∈ {5, 20, 50} are ambient feature dimensions. The fractured design uses a step-like structural component. The multi-fracture design contains several discontinuities. The smooth design 41
Wu, Chen, Xie, Lian, and Chen
removes the discontinuity. The weak-instrument design scales the first-stage systematic component. The correlated-residual design draws residuals from a Gaussian process on the latent instrument. The high-dimensional nuisance design adds irrelevant ambient variation. These designs test graph recovery when a low-dimensional geometry is embedded in higherdimensional observations. For A-IHF, the affinity graph is a symmetric K-nearest-neighbor graph with RBF weights. The default fixed configuration is K = 15, τ = 2, λ = 30, p = 80, and cutoff 10−6 . The observational rule uses α = 0.05, ε = 10−8 , cκ = 0.02, and Hutchinson trace estimation with 16 probes. These values are fixed across designs and seeds. The graph-admissibility-filtered observational variant, abbreviated as guarded observational A-IHF in tables, filters candidate graphs using first-stage graph diagnostics before applying the observational score; in the SOTA suite it requires non-degenerate final degrees and a largest post-cut component fraction of at least 0.5. B.1 Baseline Protocol and Seed Variation Table 15 records the main SOTA baseline protocol. For control-function rows in nonlinear synthetic tables, the second stage is the same additive neural regressor for all generated controls in a given experiment. It uses separate one-hidden-layer ELU components (Clevert et al., 2016) for X and the generated control, each with hidden width 64. The network is trained by full-batch Adam (Kingma and Ba, 2015) with learning rate 0.01 and weight decay 10−4 . The main nonlinear benchmark and semi-synthetic experiments use 500 stage-2 epochs; other diagnostics keep the epoch count fixed within each experiment. The epoch count is not selected by method. Table 16 reports seed variation for the main fractured cell in Table 1. Values are means with standard deviations across ten seeds. The linear second-stage alignment experiment uses the same first-stage designs, but replaces the nonlinear response with Ylin = X + 2.5U + ϵ. The second stage is ordinary least squares on (1, X, V̂ ). This experiment evaluates the coefficient perturbation bound in Corollary 9. The orthogonal linear-inference check uses the same outcome and fivefold cross-fitting. It residualizes Ylin and X on the generated control by a linear nuisance regression and reports empirical score intervals. The split-selection check uses five designs at n = 800, dZ = 50: fractured, smooth, weak instrument, correlated residual, and highdimensional nuisance. It selects parameters on one split and then refits the transductive graph smoother on the full sample. The ablation experiment uses the same five-design grid and compares fixed anisotropic A-IHF with isotropic diffusion, no pilot diffusion, no hard cutoff, full observational selection, no relevance screen, and no roughness penalty. The meta-parameter experiment varies α and cκ while holding the candidate graph family fixed. The representation experiment rebuilds the graph using the observed first-stage features, the latent instrument, principal components, and noisy variants. The graph-certificate experiment reports structural leakage, residual attenuation, firststage noise mismatch, and residualized treatment variation as defined in Proposition 10. These terms use simulator information and are hidden in observational applications. The boundary stress tests use five designs, n = 800, dZ = 50, and ten seeds. The instrumentstrength path uses {0.2, 0.4, 0.6, 0.8, 1.0}. The irrelevant-dimension path appends {0, 25, 50, 42
Graph Diffusion Residuals for Control-Function IV
Table 15: Baseline protocol. First-stage targets and tuning rules used in the synthetic SOTA benchmark. All control-function methods use the generated residual as the second-stage control, while A-IHF uses the generated graph residual. Method
First-stage or structural target
Graph ridge CF
Graph-resolvent estimate of the first-stage Fixed default or graph GCV over K and λ systematic component Truncated graph spectral smoother for X on Graph GCV over spectral rank first-stage features Polynomial series on principal-component Five-fold first-stage CV over PC dimension, scores degree, and ridge penalty RBF kernel ridge regression for X on first- Five-fold first-stage CV over bandwidth and stage features ridge penalty Random forest regression for X on first-stage Five-fold first-stage CV over leaf size and feafeatures ture subsampling Histogram gradient boosting regression for X Five-fold first-stage CV over learning rate, on first-stage features leaves, leaf size, and ℓ2 penalty Gradient-boosted trees for X on first-stage Five-fold first-stage CV over depth, learnfeatures ing rate, subsampling, column sampling, and regularization One-hidden-layer ELU network for X on Fixed architecture, learning rate, training first-stage features length, and weight decay ELU networks for X on first-stage features Held-out prediction loss over widths, depths, learning rates, and weight decay Three fixed deep first stages with different Fixed architecture and averaging across enseeds semble members Anisotropic graph-resolvent estimate of the Fixed defaults, observational finite-path rule, first-stage systematic component or guarded observational rule
Graph spectral CF Series CF Kernel ridge CF Random forest CF HistGBDT CF XGBoost CF
Fixed Deep CF CV-tuned Deep CF Deep ensemble CF A-IHF
Tuning rule
100, 200} independent Gaussian coordinates to the observed first-stage feature representation. The representation-noise path adds isotropic Gaussian noise with scale {0, 0.25, 0.5, 1, 2} and rescales the representation. All boundary tests use the same fixed and observational A-IHF configurations as the main n = 800, dZ = 50 experiment. The semi-synthetic real-Z experiment uses covariates from the diabetes, breast_cancer, and digits data sets distributed with scikit-learn (Pedregosa et al., 2011). Covariates are standardized before graph construction. Each run samples n = 400 observations with replacement when needed. The simulated first stage is built from low-dimensional scores of the real covariate matrix: the fractured design applies a discontinuous transformation, the smooth design removes the discontinuity, and the weak-instrument design scales down the first-stage systematic component. The nonlinear response uses U as the control component. The linear coefficient experiment in Table 12 uses Ylin = β0 X + γ0 V ⋆ + ϵ. This design preserves access to U , V ⋆ , and f0 for evaluation while replacing the synthetic first-stage feature cloud by empirical tabular geometry. The Card schooling tabular stability check uses the Wooldridge version of the Card data. The outcome is lwage, the treatment is educ, and the excluded instrument is nearc4. Controls are exper, expersq, black, smsa, south, smsa66, and 1966 region indicators. AIHF, KRR, RF, and GBDT control-function baselines build first-stage residuals from the standardized first-stage feature matrix. Bootstrap samples are drawn with replacement at 43
Wu, Chen, Xie, Lian, and Chen
Table 16: Seed variation in the main fractured design. Values are mean (standard deviation) over 10 seeds for n = 800, dZ = 50.
Method A-IHF, guarded observational A-IHF, observational A-IHF, fixed Graph ridge CF Graph spectral CF Random forest CF XGBoost CF HistGBDT CF Kernel ridge CF Graph ridge CF, fixed CV-tuned Deep CF Series CF Deep ensemble CF Fixed Deep CF
Control corr. ↑ Structural-response MSE ↓ 0.948 (0.011) 0.948 (0.011) 0.953 (0.014) 0.853 (0.023) 0.858 (0.018) 0.827 (0.025) 0.793 (0.025) 0.770 (0.027) 0.750 (0.026) 0.727 (0.037) 0.764 (0.031) 0.751 (0.019) 0.715 (0.031) 0.701 (0.032)
1.732 (0.844) 1.732 (0.844) 1.810 (0.733) 2.163 (1.085) 2.179 (1.117) 2.277 (1.151) 2.603 (1.296) 2.822 (1.330) 2.883 (1.389) 2.948 (1.266) 3.153 (1.625) 3.205 (1.420) 3.589 (1.668) 3.860 (1.711)
the observation level. The bootstrap uses the same pre-specified observational candidate family as the main A-IHF benchmark. The component-runtime experiment uses the fixed A-IHF configuration at n = 800, dZ = 50 and decomposes first-stage graph time into affinity construction, pilot solve, conductance weighting, final solve, and Hutchinson trace estimation. The approximate-solver experiment uses the same A-IHF graph construction as the fixed configuration and compares exact sparse solves with conjugate-gradient solves for the final resolvent. The reported δn is computed relative to the exact sparse solution for the same graph. Solve times are wall-clock times for the final linear system only. Graph construction, pilot smoothing, trace estimation, and downstream regression are not included in those timings. B.2 Orthogonal Linear-Inference Check Remark 26 gives a root-n expansion under product-rate nuisance conditions. We run a finitesample check with the same linear outcome as Section 7.2.3. For each generated control, we use five-fold cross-fitting and residualize Ylin and X on the generated control by a linear nuisance regression. The generated control itself is computed before the second-stage crossfitting. The interval uses the empirical score variance. The Card schooling check uses a separate bootstrap. The oracle row has small error and coverage close to the nominal level. Observational A-IHF reduces average absolute error relative to fixed A-IHF and to the linear control function. Coverage remains below nominal in several designs, including correlated-residual and high-dimensional nuisance designs. The gap between the oracle and generated-control 44
Graph Diffusion Residuals for Control-Function IV
Table 17: Orthogonal linear-inference check. Entries report mean absolute error and empirical 95% coverage over 10 seeds for n = 800, dZ = 50. Oracle U uses hidden simulator information and is included only as a reference. Oracle U
Design
A-IHF obs.
A-IHF fixed
Linear CF
Abs. error Cov. Abs. error Cov. Abs. error Cov. Abs. error Cov. Fractured Smooth Weak instrument Correlated residual High-dimensional nuisance
0.003 0.021 0.008 0.004 0.004
1.00 0.80 1.00 1.00 1.00
0.030 0.089 0.117 0.410 0.252
0.60 0.40 0.70 0.10 0.00
0.053 0.162 0.270 0.399 0.420
0.30 0.10 0.10 0.10 0.00
0.047 0.261 0.276 0.449 0.100
0.60 0.10 0.00 0.00 0.40
Average
0.008
0.96
0.180
0.36
0.261
0.12
0.227
0.22
rows measures the finite-sample cost of control estimation. The pattern is consistent with the nuisance-rate condition in Remark 26 and the intrinsic-dimension restriction in Remark 27. B.3 Meta-Parameter Sensitivity The observational score contains two fixed scale constants, α and cκ . We varied α ∈ {0, 0.01, 0.05, 0.10, 0.20},
cκ ∈ {0.005, 0.01, 0.02, 0.05}.
(157)
Table 18 reports the fixed A-IHF default, the default observational rule, and the range of mean structural-response MSEs across the 20 meta-parameter pairs. The candidate graph family is the same as in the main benchmark. Table 18: Meta-parameter sensitivity. Mean structural-response MSE over 10 seeds for n = 800, dZ = 50. The last column gives the range across all tested (α, cκ ) pairs. Design Fractured High-dimensional nuisance Smooth Correlated residual
A-IHF fixed A-IHF obs. default A-IHF obs. range 1.710 2.240 1.056 5.312
1.804 2.793 0.993 5.042
1.693–1.891 2.186–2.743 1.032–1.509 4.911–5.532
The selected graph changes little across this experiment. The rule selects p = 70 in every run. It selects K = 10 in all high-dimensional nuisance and correlated-residual runs, and in most fractured and smooth runs. The main behavior is driven more by the candidate graph geometry than by fine adjustment of α or cκ . 45
Wu, Chen, Xie, Lian, and Chen
B.4 Component Runtime Check We measure the main A-IHF graph components at n = 800, dZ = 50, using the fixed configuration. The timings exclude downstream outcome regression and are reported only for the first-stage graph pipeline. Table 19: Component runtime decomposition. Mean over 10 seeds for fixed A-IHF at n = 800, dZ = 50. Times are seconds. Percentages are shares of the listed component total. Design
Total
Graph
Pilot solve
Conductance
Final solve
Trace
Fractured High-dimensional nuisance Smooth Weak instrument Correlated residual
0.0789 0.1405 0.0746 0.0764 0.0679
0.0469 (59.4%) 0.0406 (28.9%) 0.0409 (54.8%) 0.0448 (58.7%) 0.0429 (63.2%)
0.0040 (5.1%) 0.0110 (7.8%) 0.0038 (5.1%) 0.0038 (5.0%) 0.0043 (6.3%)
0.0023 (2.9%) 0.0023 (1.6%) 0.0020 (2.7%) 0.0020 (2.7%) 0.0020 (2.9%)
0.0036 (4.5%) 0.0109 (7.8%) 0.0043 (5.8%) 0.0036 (4.8%) 0.0037 (5.5%)
0.0222 (28.1%) 0.0758 (53.9%) 0.0235 (31.6%) 0.0220 (28.8%) 0.0150 (22.1%)
Table 19 decomposes runtime in the tested regime. Graph construction and Hutchinson trace estimation account for most of the fixed A-IHF time. The pilot and final sparse solves are small at this sample size. The large-n solver experiment below isolates the final resolvent. B.5 Approximate Resolvent Check The exact A-IHF implementation solves sparse linear systems for the pilot and final graph resolvents. Proposition 28 shows that an approximate solve affects the generated control through the empirical error δn = n−1/2 ∥(S̃h − Sh )x∥2 .
(158)
We evaluate conjugate-gradient solves for the final resolvent on three designs with dZ = 50, n ∈ {800, 3000, 10000}, and five seeds. Table 20 reports the largest sample size and tolerance 10−4 . The times are solve-component times for the final resolvent, not end-to-end pipeline times. Table 20: Approximate resolvent solves. Mean over five seeds at n = 10000, dZ = 50. CG uses relative tolerance 10−4 . Design Fractured Smooth High-dimensional nuisance
Exact corr. CG corr. 0.96290 0.98487 0.69430
0.96290 0.98487 0.69430
δn 4.36 × 10−5 2.79 × 10−5 3.32 × 10−5
Exact time CG time 1.115s 1.139s 6.867s
0.0145s 0.0150s 0.0164s
The generated residual is unchanged at the reported precision. Jacobi-preconditioned CG gives the same control correlations at the reported precision and smaller solve times in this experiment. Approximate resolvent solves can be used in larger graph instances when δn and κn (v̂) are monitored. 46
Graph Diffusion Residuals for Control-Function IV
Solve time at tolerance 10−4
Solve accuracy CG Jacobi CG
10−4
Final solve time (s, log scale)
Generated-control perturbation δn
10−3
10−5 10−6 10−7 10−8
Exact sparse CG Jacobi CG 1
0.1
0.01
10−9 10−8
10−7
10−6
10−5
10−4
10−3
Fractured
Relative tolerance
Smooth
High-dim. nuisance
Figure 6: Approximate final-resolvent solves at n = 10000. The left panel reports the generated-control perturbation δn across iterative-solve tolerances. The right panel compares exact sparse, CG, and Jacobi-preconditioned CG solve times for the final resolvent. These are solve-component times, not end-to-end pipeline times.
B.6 Evaluation Grid Note The main benchmark evaluates the structural response on the empirical support of X in each run. As n changes, this support can widen. Cross-n MSE comparisons combine estimation effects with support-width effects and are not used as empirical evidence for Proposition 12. The consistency statement is a graph-sequence reduction, not a monotonicity claim about these finite benchmark tables.
References Joshua D. Angrist and Jörn-Steffen Pischke. Mostly Harmless Econometrics: An Empiricist’s Companion. Princeton University Press, 2009. Joshua D. Angrist, Guido W. Imbens, and Donald B. Rubin. Identification of causal effects using instrumental variables. Journal of the American Statistical Association, 91(434): 444–455, 1996. Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural Computation, 15(6):1373–1396, 2003. Mikhail Belkin and Partha Niyogi. Convergence of laplacian eigenmaps. In Advances in Neural Information Processing Systems 19, pages 129–136, 2006. Mikhail Belkin, Partha Niyogi, and Vikas Sindhwani. Manifold regularization: A geometric framework for learning from labeled and unlabeled examples. Journal of Machine Learning Research, 7(85):2399–2434, 2006. 47
Wu, Chen, Xie, Lian, and Chen
Andrew Bennett, Nathan Kallus, and Tobias Schnabel. Deep generalized method of moments for instrumental variable analysis. In Advances in Neural Information Processing Systems 32, pages 3564–3574, 2019. Richard Blundell and James L. Powell. Endogeneity in nonparametric and semiparametric regression models. In Mathias Dewatripont, Lars Peter Hansen, and Stephen J. Turnovsky, editors, Advances in Economics and Econometrics: Theory and Applications, Eighth World Congress, Volume II, pages 312–357. Cambridge University Press, 2003. Béla Bollobás. Random Graphs. Cambridge University Press, 2nd edition, 2001. David Card. Using geographic variation in college proximity to estimate the return to schooling. In Louis N. Christofides, E. Kenneth Grant, and Robert Swidinsky, editors, Aspects of Labour Market Behaviour: Essays in Honour of John Vanderkamp, pages 201– 222. University of Toronto Press, 1995. Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins. Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1):C1–C68, 2018. Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter. Fast and accurate deep network learning by exponential linear units (elus). In International Conference on Learning Representations, 2016. Ronald R. Coifman and Stéphane Lafon. Diffusion maps. Applied and Computational Harmonic Analysis, 21(1):5–30, 2006. Peter Craven and Grace Wahba. Smoothing noisy data with spline functions: Estimating the correct degree of smoothing by the method of generalized cross-validation. Numerische Mathematik, 31(4):377–403, 1979. doi: 10.1007/BF01404567. Serge Darolles, Yanqin Fan, Jean-Pierre Florens, and Eric Renault. Nonparametric instrumental regression. Econometrica, 79(5):1541–1565, 2011. Nishanth Dikkala, Greg Lewis, Lester Mackey, and Vasilis Syrgkanis. Minimax estimation of conditional moment models. In Advances in Neural Information Processing Systems 33, pages 12248–12262, 2020. Nicolas Garcia Trillos and Dejan Slepčev. A variational approach to the consistency of spectral clustering. Applied and Computational Harmonic Analysis, 45(2):239–281, 2018. Gene H. Golub, Michael Heath, and Grace Wahba. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics, 21(2):215–223, 1979. Geoffrey Grimmett. Percolation. Springer, 2nd edition, 1999. Peter Hall and Joel L. Horowitz. Nonparametric methods for inference in the presence of instrumental variables. The Annals of Statistics, 33(6):2904–2929, 2005. 48
Graph Diffusion Residuals for Control-Function IV
Jason Hartford, Greg Lewis, Kevin Leyton-Brown, and Matt Taddy. Deep IV: A flexible approach for counterfactual prediction. In Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 1414–1423, 2017. Matthias Hein, Jean-Yves Audibert, and Ulrike von Luxburg. From graphs to manifolds: Weak and strong pointwise consistency of graph laplacians. In Learning Theory, volume 3559 of Lecture Notes in Computer Science, pages 470–485. Springer, 2005. Michael F. Hutchinson. A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in Statistics - Simulation and Computation, 18(3):1059–1076, 1989. Guido W. Imbens and Joshua D. Angrist. Identification and estimation of local average treatment effects. Econometrica, 62(2):467–475, 1994. Arthur Jacot, Franck Gabriel, and Clément Hongler. Neural tangent kernel: Convergence and generalization in neural networks. In Advances in Neural Information Processing Systems 31, 2018. Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015. Krikamol Muandet, Arash Mehrjou, Si Kai Lee, and Anant Raj. Dual instrumental variable regression. In Advances in Neural Information Processing Systems 33, pages 2710–2721, 2020. Whitney K. Newey and James L. Powell. Instrumental variable estimation of nonparametric models. Econometrica, 71(5):1565–1578, 2003. Whitney K. Newey, James L. Powell, and Francis Vella. Nonparametric estimation of triangular simultaneous equations models. Econometrica, 67(3):565–603, 1999. Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Édouard Duchesnay. Scikit-learn: Machine learning in python. Journal of Machine Learning Research, 12(85):2825–2830, 2011. Pietro Perona and Jitendra Malik. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(7):629–639, 1990. Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred A. Hamprecht, Yoshua Bengio, and Aaron Courville. On the spectral bias of neural networks. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 5301–5310, 2019. Douglas Rivers and Quang H. Vuong. Limited information estimators and exogeneity tests for simultaneous probit models. Journal of Econometrics, 39(3):347–366, 1988. 49
Wu, Chen, Xie, Lian, and Chen
Yousef Saad. Iterative Methods for Sparse Linear Systems. SIAM, 2nd edition, 2003. Jonathan Richard Shewchuk. An introduction to the conjugate gradient method without the agonizing pain. Technical report, Carnegie Mellon University, 1994. David I. Shuman, Sunil K. Narang, Pascal Frossard, Antonio Ortega, and Pierre Vandergheynst. The emerging field of signal processing on graphs: Extending highdimensional data analysis to networks and other irregular domains. IEEE Signal Processing Magazine, 30(3):83–98, 2013. Amit Singer. From graph to manifold laplacian: The convergence rate. Applied and Computational Harmonic Analysis, 21(1):128–134, 2006. Rahul Singh, Maneesh Sahani, and Arthur Gretton. Kernel instrumental variable regression. In Advances in Neural Information Processing Systems 32, pages 4595–4607, 2019. Douglas Staiger and James H. Stock. Instrumental variables regression with weak instruments. Econometrica, 65(3):557–586, 1997. James H. Stock, Jonathan H. Wright, and Motohiro Yogo. A survey of weak instruments and weak identification in generalized method of moments. Journal of Business & Economic Statistics, 20(4):518–529, 2002. Vasilis Syrgkanis, Victor Lei, Miruna Oprescu, Maggie Hei, Keith Battocchi, and Greg Lewis. Machine learning estimation of heterogeneous treatment effects with instruments. In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d’Alché-Buc, Emily B. Fox, and Roman Garnett, editors, Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pages 15167–15176, 2019. Andrei N. Tikhonov and Vasiliy Y. Arsenin. Solutions of Ill-Posed Problems. V. H. Winston & Sons, 1977. Daniel Ting, Ling Huang, and Michael I. Jordan. An analysis of the convergence of graph laplacians. In Proceedings of the 27th International Conference on Machine Learning, pages 1079–1086, 2010. Ulrike von Luxburg. A tutorial on spectral clustering. Statistics and Computing, 17(4): 395–416, 2007. Ulrike von Luxburg, Mikhail Belkin, and Olivier Bousquet. Consistency of spectral clustering. The Annals of Statistics, 36(2):555–586, 2008. Grace Wahba. Spline Models for Observational Data, volume 59 of CBMS-NSF Regional Conference Series in Applied Mathematics. SIAM, 1990. Joachim Weickert. Anisotropic Diffusion in Image Processing. B. G. Teubner, Stuttgart, 1998. Jeffrey M. Wooldridge. Econometric Analysis of Cross Section and Panel Data. MIT Press, 2nd edition, 2010. 50