Perturbative methods for non-parametric instrumental variable

# Perturbative methods for non-parametric instrumental variable
Source: [https://arxiv.org/html/2606.00322](https://arxiv.org/html/2606.00322)
###### Abstract

We introduce a perturbative approach for nonparametric instrumental variable \(NPIV\) estimation\. By drawing inspiration from perturbation theory in physics, we extend standard kernel ridge methods with systematic higher perturbation order corrections that significantly improve estimation accuracy\. Spectrally, the perturbation introduces mixing between different eigenmodes of the expectation integral operator, which becomes especially useful when the integral equation is ill\-defined\. One source for such ill\-definedness can be the curse of dimensionality\. Our method performs across various dimensionality regimes, particularly when the dimensionality parameterβ\\betawhich is defined through the number of samplesnnand dimensionddasnβ=dn^\{\\beta\}=d, becomes large\. Experimental results show that our first\-order perturbative corrections can reduce prediction error by up to 99% in high\-dimensional ill\-defined cases \(β\>0\.7\\beta\>0\.7\) compared to standard ridge regression approaches\. The performance improvement is maintained across a wide range of dimensions, with the advantage becoming more pronounced as dimensionality increases\.

Machine Learning, ICML

## 1Introduction

Nonparametric instrumental variable \(NPIV\) estimation has emerged as a fundamental tool for causal inference in the presence of unmeasured confounding\. The method leverages instrumental variables–variables that affect the treatment but not the outcome directly–to identify causal effects without imposing restrictive parametric assumptions on the functional form of the causal relationship\. Despite its theoretical appeal, NPIV estimation remains a challenging problem especially when condition number of the expectation integral operator is high\. This can be the result of multiple different factors: the presence of weak instrumental, misalignment between target and kernel spectrum and curse of dimensionality\. Informatively speaking, we consider a setting to be ”high dimensional” when sample size is small compared to the dimensionality of the variables\. As the dimension of the feature space increases relative to sample size, estimation accuracy typically deteriorates\.

Classical approaches to NPIV estimation typically formulate the problem as an ill\-posed integral equation, where the causal function is estimated by solving a conditional regression problem\(Newey and Powell,[2003](https://arxiv.org/html/2606.00322#bib.bib21); Hall and Horowitz,[2005](https://arxiv.org/html/2606.00322#bib.bib16)\)\. Kernel\-based methods, particularly those employing reproducing kernel Hilbert spaces \(RKHS\), have been applied to nonparameteric estimation due to their flexibility and theoretical guarantees\(Darolleset al\.,[2011](https://arxiv.org/html/2606.00322#bib.bib33); Singhet al\.,[2020](https://arxiv.org/html/2606.00322#bib.bib19)\)\. However, these methods encounter fundamental limitations whenβ≫1\\beta\\gg 1\(high treatment dimension\) even in the plain kernel ridge regression case with rotationally invariant kernels, leading to poorly conditioned kernel matrices and unstable estimations\(Donhauseret al\.,[2021](https://arxiv.org/html/2606.00322#bib.bib15)\)\.

The curse of dimensionality in kernel methods manifests in several ways: kernel matrices become increasingly ill\-conditioned, with exponentially decaying eigenvalue spectrum \(high condition number\), which leads to worse performance\(Steinwart and Christmann,[2008](https://arxiv.org/html/2606.00322#bib.bib34); Donhauseret al\.,[2021](https://arxiv.org/html/2606.00322#bib.bib15)\)in the kernel ridge solution\. This leads to ridge solutions that primarily retain high eigenvalue modes while discarding low eigenvalue directions, similar limitations have been noted in low\-rank kernel approximations\(Bach,[2013](https://arxiv.org/html/2606.00322#bib.bib14)\)\. In the context of NPIV estimation, the layer of conditional expectation adds a further layer of complication\. The inherently ill\-posedness of the integral equation, where small perturbations in the data can lead to large changes in the estimated causal function\(Carrascoet al\.,[2007](https://arxiv.org/html/2606.00322#bib.bib35)\)\. The ill\-posedness of integral equations can be simply understood as attempting to recover information of the original function from the convolved function, which for a generic convolution kernel is not fully invertible\. Existing regularization techniques such as Tikhonov regularization provide general suppression over all unstable eigendirections of the kernel function, can result in over\-smoothed estimates that fail to capture complex nonlinear relations\. This impedes the performance especially when those eigendirections align with the target causal function111Rather than kernel spectrum which we are discussing here, misaligned eigenspectrum in the conditional integral operator in IV was studied in\(Meunieret al\.,[2025](https://arxiv.org/html/2606.00322#bib.bib8)\)\.\.

Recent advances in understanding the fundamental limits of NPIV estimation have highlighted the polynomial approximation barrier where standard kernel methods require sample sizes that grow polynomially with dimension to maintain estimation accuracy\(Donhauseret al\.,[2021](https://arxiv.org/html/2606.00322#bib.bib15)\)\. This limit has motivated the development of alternative approaches, including deep learning methods\(Xuet al\.,[2023](https://arxiv.org/html/2606.00322#bib.bib36); Kimet al\.,[2025](https://arxiv.org/html/2606.00322#bib.bib9)\)and other regularization schemes\(Muandetet al\.,[2020](https://arxiv.org/html/2606.00322#bib.bib20)\)\.

In this paper, we introduce a novel perturbative renormalization approach to NPIV estimation that addresses the situations where the kernel matrix is low ranked or with high condition number \(caused mainly by high dimensionality\) while preserving the kernel framework’s desirable properties\. Our method draws inspiration from quantum physics, where perturbative expansions and renormalization techniques are used to handle divergent series and extract finite meaningful results from seemingly intractable calculations\(Peskin and Schroeder,[1995](https://arxiv.org/html/2606.00322#bib.bib37); de Faria and de Melo,[2010](https://arxiv.org/html/2606.00322#bib.bib44)\)\. The key insight is that the kernel coefficients in NPIV estimation can be expanded as a power series in a coupling parameter, allowing for systematic correction of the base kernel ridge regression solution\. Our approach consists of three main components:\(1\)\(1\)a perturbative expansion that introduces higher\-order moment interactions to capture complex dependencies in high dimensional spaces,\(2\)\(2\)a renormalization procedure that tames the expansion by adaptively rescaling the higher order terms in the power series, and\(3\)\(3\)a resurgence technique that handles cases where the naive perturbative series becomes factorially divergent\. From a spectral perspective, we demonstrate how these higher moment interactions induce perturbative solution effectively boost the contribution of eigenmodes with small eigenvalues and nonlinearly mix contributions from different eigenmodes222A different approach attempting to regularizing the kernel eigenmodes were discussed in the context of kernel MMD flow\(Chenet al\.,[2024](https://arxiv.org/html/2606.00322#bib.bib13); Hagrasset al\.,[2024](https://arxiv.org/html/2606.00322#bib.bib11)\)\.\.

We demonstrate the effectiveness of our approach through extensive experiments on challenging high\-dimensional NPIV problems\. Our results show that Gaussian RBF kernels see minimal improvement due to their rotational invariance which almost zero them out in high dimensions, the family of fractional Brownian kernels\(Sejdinovicet al\.,[2013](https://arxiv.org/html/2606.00322#bib.bib31)\)achieve substantial performance gains, with improvements of up to99%99\\%in mean square error compared to base kernel ridge regression\. The method is particularly effective in regimes where the dimensionality grows rapidly with sample size, precisely where traditional NPIV methods struggle most\. Our approach opens new possibilities for tackling ill\-conditioned problems in causal inference and suggests broader applications of perturbative methods in machine learning\.

This paper is structured in the following way\. In section[2](https://arxiv.org/html/2606.00322#S2)and[3](https://arxiv.org/html/2606.00322#S3), we introduce the basic approach to NPIV using kernel ridge regression\. In section[4](https://arxiv.org/html/2606.00322#S4), we describe the novel perturbative corrections we add to standard ridge regression and how we regularize the result\. In section[5](https://arxiv.org/html/2606.00322#S5), we discuss the implementations of the perturbative method and in section[6](https://arxiv.org/html/2606.00322#S6)we present and discuss the experiments\.

### 1\.1Related works

Nonparametric instrumental variable estimation has a rich history in econometrics and statistics\. Early work by\(Newey and Powell,[2003](https://arxiv.org/html/2606.00322#bib.bib21)\)and\(Hall and Horowitz,[2005](https://arxiv.org/html/2606.00322#bib.bib16)\)established the theoretical foundations, while subsequent research has focused on addressing the ill\-posedness of the inverse problem inherent in NPIV\. Kernel\-based approaches to NPIV have been developed by\(Singhet al\.,[2020](https://arxiv.org/html/2606.00322#bib.bib19)\)and\(Muandetet al\.,[2020](https://arxiv.org/html/2606.00322#bib.bib20)\), leveraging reproducing kernel Hilbert spaces \(RKHS\) to represent the unknown structural function\. These methods typically rely on Tikhonov regularization to ensure stability, but their performance degrades rapidly in high dimensions\. The challenges of high\-dimensional kernel ridge regression were systematically analyzed by\(Donhauseret al\.,[2021](https://arxiv.org/html/2606.00322#bib.bib15)\), who established minimax optimal rates and identified fundamental limits on the performance of polynomial approximation methods when the dimensionality parameterβ\\betaexceeds certain thresholds\. Their analysis showed that standard methods face a ”polynomial approximation barrier” as dimensionality increases\. Various approaches have been proposed to address high\-dimensional challenges in related contexts\.\(Belloniet al\.,[2015](https://arxiv.org/html/2606.00322#bib.bib22)\)developed LASSO\-based methods for high\-dimensional instrumental variables, focusing on sparse linear models\. Neural network approaches have been explored by\(Hartfordet al\.,[2017](https://arxiv.org/html/2606.00322#bib.bib17)\)and\(Bennettet al\.,[2020](https://arxiv.org/html/2606.00322#bib.bib23); Xuet al\.,[2023](https://arxiv.org/html/2606.00322#bib.bib36); Kimet al\.,[2025](https://arxiv.org/html/2606.00322#bib.bib9)\), leveraging the representation power of deep learning for instrumental variable estimation\.

Our perturbative renormalization approach differs from existing methods by systematically incorporating higher\-order corrections while ensuring stability through quantum field theory\-motivated renormalization\. This allows us to achieve performance improvements in dimensionality regimes where traditional methods struggle, without sacrificing interpretability or theoretical foundations\.

## 2Causal Framework

We consider a causal model with the following structure:

ZZXXYYUUFigure 1:Causal diagram:ZZis an instrument forXX, whileUUrepresents unobserved confounding betweenXXandYY\.In this framework,ZZis an instrumental variable that affectsXXbut notYYdirectly\.XXis the treatment/exposure variable that causally affectsYY, which is the outcome of interest\. Finally,UUrepresents unobserved confounding that affects bothXXandYY\.

The key assumptions are:ZZaffectsYYonly throughXX\(exclusion restriction\);ZZhas a non\-zero effect onXX;ZZis independent of the unobserved confoundersUU:U⟂⟂ZU\\perp\\\!\\\!\\\!\\perp Zand𝔼​\[U\|Z\]=0\\mathbb\{E\}\[U\|Z\]=0and the confounding effect ofUUonYYis additive333If this is not the case, one could use proxy learning\(Miaoet al\.,[2018](https://arxiv.org/html/2606.00322#bib.bib12)\)which also leads to an integral equation with similar issues, the algorithm described in this paper would also improve the performance\.\. The structural equation for the outcomeYYcan be written as:

Y=g​\(X\)\+U,Y=g\(X\)\+U\\,,\(1\)whereg​\(X\)g\(X\)is the causal effect of interest,UUis the unobserved confounder\. Our goal is to estimate the functiong​\(X\)g\(X\)that represents the causal effect ofXXonYY\.

## 3Nonparametric Instrumental Variable Estimation \(NPIV\)

Due to the confounding byUU, standard regression ofYYonXXwill not identify the causal effectg​\(X\)g\(X\)\. Instead, we use the instrumentZZto form moment conditions\. The key moment equation is:

𝔼​\[Y\|Z\]=𝔼​\[g​\(X\)\|Z\],\\mathbb\{E\}\[Y\|Z\]=\\mathbb\{E\}\[g\(X\)\|Z\]\\,,\(2\)after using the fact thatU⟂⟂ZU\\perp\\\!\\\!\\\!\\perp Zand𝔼​\[U\|Z\]=0\\mathbb\{E\}\[U\|Z\]=0\.

We represent the unknown functiong​\(X\)g\(X\)using the reproducing kernel Hilbert space \(RKHS\) approach, with finitennsamples, the representer theory for kernel ridge regression dictates that:

g​\(x\)=∑i=1nαi​K​\(x,xi\),g\(x\)=\\sum\_\{i=1\}^\{n\}\\alpha\_\{i\}K\(x,x\_\{i\}\)\\,,\(3\)whereK​\(x,x′\)K\(x,x^\{\\prime\}\)is a positive definite kernel,xix\_\{i\}are stage 2 samples conditioned onZZin the standard two\-stage approach andαi\\alpha\_\{i\}are coefficients to be determined from the standard representer theorem\(Schölkopfet al\.,[2001](https://arxiv.org/html/2606.00322#bib.bib10)\)\.

To estimateg​\(X\)g\(X\), we formulate the following objective function:

S0=𝔼​\[\(𝔼​\[Y\|Z\]−∫g​\(x\)​f​\(x\|Z\)​𝑑x\)2\],S\_\{0\}=\\mathbb\{E\}\\left\[\\left\(\\mathbb\{E\}\[Y\|Z\]\-\\int g\(x\)f\(x\|Z\)dx\\right\)^\{2\}\\right\]\\,,\(4\)wheref​\(x\|Z\)f\(x\|Z\)is the conditional density ofXXgivenZZ, the outer expectation is taken with respect toZZ\.

Substituting the kernel representation, we get:

S0=𝔼​\[\(𝔼​\[Y\|Z\]−∫∑i=1nαi​K​\(x,xi\)​f​\(x\|Z\)​d​x\)2\]\.S\_\{0\}=\\mathbb\{E\}\\left\[\\left\(\\mathbb\{E\}\[Y\|Z\]\-\\int\\sum\_\{i=1\}^\{n\}\\alpha\_\{i\}K\(x,x\_\{i\}\)f\(x\|Z\)dx\\right\)^\{2\}\\right\]\\,\.\(5\)To ensure stability of the solution, we add a quadratic regularization term:

Sridge=S0\+λ​‖g​\(x\)‖ℋ2,S\_\{\\text\{ridge\}\}=S\_\{0\}\+\\lambda\\\|g\(x\)\\\|^\{2\}\_\{\\mathcal\{H\}\}\\,,\(6\)where‖g​\(x\)‖ℋ2\\\|g\(x\)\\\|^\{2\}\_\{\\mathcal\{H\}\}denotes the RKHS norm ofg​\(x\)g\(x\)\. To derive the estimator, we first define:

hi:=𝔼​\[∫K​\(x,xi\)​f​\(x\|Z\)​𝑑x⋅𝔼​\[Y\|Z\]\];\\displaystyle h\_\{i\}:=\\mathbb\{E\}\\left\[\\int K\(x,x\_\{i\}\)f\(x\|Z\)dx\\cdot\\mathbb\{E\}\[Y\|Z\]\\right\]\\,;\(7\)K~i​j:=𝔼​\[∬K​\(x,xi\)​K​\(x′,xj\)​f​\(x\|Z\)​f​\(x′\|Z\)​𝑑x​𝑑x′\];\\displaystyle\\widetilde\{K\}\_\{ij\}:=\\mathbb\{E\}\\left\[\\iint K\(x,x\_\{i\}\)K\(x^\{\\prime\},x\_\{j\}\)f\(x\|Z\)f\(x^\{\\prime\}\|Z\)dxdx^\{\\prime\}\\right\]\\,;where the outer expectation is taken with respect toZZ\. Now we can rewrite the objective function as \(omitting the𝔼​\[Y\|Z\]2\\mathbb\{E\}\[Y\|Z\]^\{2\}term which is irrelevant from the RKHS coefficientαi\\alpha\_\{i\}\):

Sridge=∑i,jK~i​j​αi​αj−2​∑iαi​hi\+λ​∑i,jαi​αj​Ki​j,S\_\{\\text\{ridge\}\}=\\sum\_\{i,j\}\\widetilde\{K\}\_\{ij\}\\alpha\_\{i\}\\alpha\_\{j\}\-2\\sum\_\{i\}\\alpha\_\{i\}h\_\{i\}\+\\lambda\\sum\_\{i,j\}\\alpha\_\{i\}\\alpha\_\{j\}K\_\{ij\}\\,,\(8\)whereKi​jK\_\{ij\}is the standard Gram matrix elements, quite distinct fromK~i​j\\widetilde\{K\}\_\{ij\}which is the density smoothed version\. Taking the functional derivative with respect toαi\\alpha\_\{i\}and setting it to zero444For a detailed discussion of functional derivatives and how they act, we refer the readers to the appendix[D](https://arxiv.org/html/2606.00322#A4)\.:

δ​Sridgeδ​αi=0⇒∑jK~i​j​αj\+λ​Ki​j​αj=hi\.\\frac\{\\delta S\_\{\\text\{ridge\}\}\}\{\\delta\\alpha\_\{i\}\}=0\\Rightarrow\\sum\_\{j\}\\widetilde\{K\}\_\{ij\}\\alpha\_\{j\}\+\\lambda K\_\{ij\}\\alpha\_\{j\}=h\_\{i\}\\,\.\(9\)In matrix form\(K~\+λ​K\)​α=h\(\\widetilde\{K\}\+\\lambda K\)\\alpha=h\. The solution is the standard kernel ridge, denoted asα\(0\)\\alpha^\{\(0\)\}, the zeroth\-order solution in our perturbative approach\.

## 4Perturbative Approach

When condition numberκ\\kappais high or the effective rank is low, a well\-understood problem with kernel methods is that eigendirections with small eigenvalues become under\-represented in the solution space, which suggests a fast decay in the eigenvalues\. This would not be an issue if the target functionhhhappen to have features only in the eigendirections with large eigenvalues\. However, in the NPIV setting, in the presence of a convolution integral equation, for example in deep learning IV\(Wiltzeret al\.,[2024](https://arxiv.org/html/2606.00322#bib.bib45); Xuet al\.,[2023](https://arxiv.org/html/2606.00322#bib.bib36)\), the fast spectrum decay admits less and less effective eigenmodes\. High condition number can also occur when instrumental variables are weak where changes inZZcauses few to no changes ing​\(x\)g\(x\)\(Andrewset al\.,[2019](https://arxiv.org/html/2606.00322#bib.bib47); Stocket al\.,[2002b](https://arxiv.org/html/2606.00322#bib.bib46)\)\.

In order to tackle this, we need to induce mixing between different eigendirections by remixing their weights in the regression solution555We do the spectral analysis in appendix[A](https://arxiv.org/html/2606.00322#A1)\.To this end, we introduce a perturbative solution to kernel regression by incorporating additional higher\-order interactions666This is rather standard approach from physics, for a reference chapter\(de Faria and de Melo,[2010](https://arxiv.org/html/2606.00322#bib.bib44)\), and\(Peskin and Schroeder,[1995](https://arxiv.org/html/2606.00322#bib.bib37)\), where when perturbatively studying the scattering problem in quantum electrodynamics, a electron\-positron\-photon\-photon quartic interaction is added\. Also used in standard PDE solving and quantum mechanics\(Dyson,[1952](https://arxiv.org/html/2606.00322#bib.bib27); Volin,[2010](https://arxiv.org/html/2606.00322#bib.bib28); Magnen and Seneor,[1977](https://arxiv.org/html/2606.00322#bib.bib29); Lipatov,[1977](https://arxiv.org/html/2606.00322#bib.bib30)\), usually quoted as the WKB approximation method\.The simplest term one can add in the objective function is the triple moment, we define:

Ki​j​k=𝔼​\[∫K​\(x,xi\)​K​\(x,xj\)​K​\(x,xk\)​f​\(x\|Z\)​𝑑x\]K\_\{ijk\}=\\mathbb\{E\}\\left\[\\int K\(x,x\_\{i\}\)K\(x,x\_\{j\}\)K\(x,x\_\{k\}\)f\(x\|Z\)dx\\right\]\(10\)wherexxlabels an independent draw from the probability densityf​\(x\|Z\)f\(x\|Z\)\. Intuitively, this includes the potential influence of kernel functions from three different locations inℝn\\mathbb\{R\}^\{n\}on the regressed function when mapped to the Hilbert space\. We start with the objective function with a cubic term,

Spert=Sridge\+2​γ3\\displaystyle S\_\{\\text\{pert\}\}=S\_\{\\text\{ridge\}\}\+\\frac\{2\\gamma\}\{3\}∑i,j,kαi​αj​αk​Ki​j​k\.\\displaystyle\\sum\_\{i,j,k\}\\alpha\_\{i\}\\alpha\_\{j\}\\alpha\_\{k\}K\_\{ijk\}\\,\.\(11\)Taking the derivative with respect toαi\\alpha\_\{i\}, we obtain the minimization criteria for the objective function:

δ​Spertδ​αi=2𝔼\[λKi​jαj−\(𝔼\[Y\|Z\]−∫g\(x\)f\(x\|Z\)dx\)×\\displaystyle\\frac\{\\delta S\_\{\\text\{pert\}\}\}\{\\delta\\alpha\_\{i\}\}=2\\,\\mathbb\{E\}\\left\[\\lambda K\_\{ij\}\\alpha\_\{j\}\-\\left\(\\mathbb\{E\}\[Y\|Z\]\-\\int g\(x\)f\(x\|Z\)dx\\right\)\\times\\right\.∫K​\(x′,xi\)​f​\(x′\|Z\)​𝑑x′\+\\displaystyle\\qquad\\quad\\int K\(x^\{\\prime\},x\_\{i\}\)f\(x^\{\\prime\}\|Z\)dx^\{\\prime\}\+\(12\)γ∑j,kαjαk∫K\(x,xi\)K\(x,xj\)K\(x,xk\)f\(x\|Z\)dx\]\.\\displaystyle\\left\.\\gamma\\sum\_\{j,k\}\\alpha\_\{j\}\\alpha\_\{k\}\\int K\(x,x\_\{i\}\)K\(x,x\_\{j\}\)K\(x,x\_\{k\}\)f\(x\|Z\)dx\\right\]\\,\.We still expand the target functiong​\(x\)g\(x\)in the kernel Hilbert space, which traditionally in the ridge regression case relies on the optimal representer theorem\. We prove the generalized representer theorem in appendix[D](https://arxiv.org/html/2606.00322#A4)in the presence of the triple moment term\. Using the representer theorem like before, we simply set the variation \([4](https://arxiv.org/html/2606.00322#S4.Ex2)\) to zero and collecting terms, we arrive at the optimality condition on the kernel coefficientsαi\\alpha\_\{i\}that yields our perturbative expansion\.

∑j\(K~i​j\+λ​Ki​j\)​αj\+3​γ​∑j,kKi​j​k​αj​αk=hi,\\sum\_\{j\}\(\\widetilde\{K\}\_\{ij\}\+\\lambda K\_\{ij\}\)\\alpha\_\{j\}\+3\\gamma\\sum\_\{j,k\}K\_\{ijk\}\\alpha\_\{j\}\\alpha\_\{k\}=h\_\{i\}\\,,\(13\)where we have compactified the integrals,γ\\gammais a coupling parameter that controls the strength of the three\-point interaction \(three point moment\), we are essentially perturbing the regression problem around the zeroth order ”free” theory by adding this term\.

### 4\.1Perturbative Expansion

Due to the addition of triple moment term with scalarγ\\gamma, we expect the regressed resultαi\\alpha\_\{i\}to depend onγ\\gamma, which should reduce to the usual kernel ridge regression solution \([9](https://arxiv.org/html/2606.00322#S3.E9)\)\. A naive first ansatz one could use is simply the power series\. This is only a first guess so far, we are certainly not guaranteed to have an everywhere convergent series\. In fact the experience from physics is that they are almost never convergent, accordingly there are procedures to tame the divergent ones we shall introduce in the next two subsections\. Expand the coefficients as a power series inγ\\gamma

α=α\(0\)\+γ​α\(1\)\+γ2​α\(2\)\+…\.\\alpha=\\alpha^\{\(0\)\}\+\\gamma\\alpha^\{\(1\)\}\+\\gamma^\{2\}\\alpha^\{\(2\)\}\+\\ldots\\,\.\(14\)Substituting this into the optimality condition and collecting terms of the same power inγ\\gamma, we get:

∑j\(K~i​j\+λ​Ki​j\)​\(∑l=0∞γl​αj\(l\)\)=hi−γ​∑j,kKi​j​k​\(∑n=0∞γn​αj\(n\)\)​\(∑m=0∞γm​αk\(m\)\)\.\\sum\_\{j\}\(\\widetilde\{K\}\_\{ij\}\+\\lambda K\_\{ij\}\)\\left\(\\sum\_\{l=0\}^\{\\infty\}\\gamma^\{l\}\\alpha^\{\(l\)\}\_\{j\}\\right\)=\\\\ h\_\{i\}\-\\gamma\\sum\_\{j,k\}K\_\{ijk\}\\left\(\\sum\_\{n=0\}^\{\\infty\}\\gamma^\{n\}\\alpha^\{\(n\)\}\_\{j\}\\right\)\\left\(\\sum\_\{m=0\}^\{\\infty\}\\gamma^\{m\}\\alpha^\{\(m\)\}\_\{k\}\\right\)\\,\.\(15\)The benefit of expandingαi​\(γ\)\\alpha\_\{i\}\(\\gamma\)as a power series means that we can now use the simple trick of coefficient matching:

∂n∂γn​\(δ​Spertδ​αi​\(γ\)\)\|γ=0=0​forn∈ℕ\.\\left\.\\frac\{\\partial^\{n\}\}\{\\partial\\gamma^\{n\}\}\\left\(\\frac\{\\delta S\_\{\\text\{pert\}\}\}\{\\delta\\alpha\_\{i\}\}\(\\gamma\)\\right\)\\right\|\_\{\\gamma=0\}=0\\textbf\{ for $n\\in\\mathbb\{N\}$\}\\,\.\(16\)It is easy to see that we can solve for theα\(n\)\\alpha^\{\(n\)\}s order by order iteratively\. Looking at the coefficient ofγ0\\gamma^\{0\}term from both sides of the equation:

\(K~\+λ​K\)​α\(0\)=h\.\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(0\)\}=h\\,\.\(17\)This is the usual ridge regression term we talked about before\. Similarly for the coefficient ofγ1\\gamma^\{1\}, we have:

\(K~\+λ​K\)​α\(1\)=−3​∑j,kKi​j​k​αj\(0\)​αk\(0\),\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(1\)\}=\-3\\sum\_\{j,k\}K\_\{ijk\}\\alpha^\{\(0\)\}\_\{j\}\\alpha^\{\(0\)\}\_\{k\}\\,,\(18\)although notationally we omitted thexix\_\{i\}dependence, hopefully it is self\-explanatory777In full component form\(K~i​m\+λ​Ki​m\)​αm\(1\)=−3​∑j,kKi​j​k​αj\(0\)​αk\(0\)\(\\widetilde\{K\}\_\{im\}\+\\lambda K\_\{im\}\)\\alpha\_\{m\}^\{\(1\)\}=\-3\\sum\_\{j,k\}K\_\{ijk\}\\alpha^\{\(0\)\}\_\{j\}\\alpha^\{\(0\)\}\_\{k\}\.and subsequently the coefficient of genericγn\\gamma^\{n\}term:

\(K~\+λK\)α\(n\)=−3×\\displaystyle\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(n\)\}=\-3\\times\(19\)∑j,kKi​j​k​\(αj0​αkn−1\+αj1​αkn−2\+⋯\+αjn−1​αk0\)\.\\displaystyle\\sum\_\{j,k\}K\_\{ijk\}\\left\(\\alpha\_\{j\}^\{0\}\\alpha\_\{k\}^\{n\-1\}\+\\alpha\_\{j\}^\{1\}\\alpha\_\{k\}^\{n\-2\}\+\\dots\+\\alpha\_\{j\}^\{n\-1\}\\alpha\_\{k\}^\{0\}\\right\)\\,\.We see that this is an iterative algorithm, where the computation of each new perturbative order is a solving simple ridge regression problem on all the previous perturbative order coefficients888This is equivalent to a technique familiar in perturbative quantum field theory called the Feynman diagram, where the sum is computed by drawing all possible diagrams whose external legs are theαi\(n\)\\alpha\_\{i\}^\{\(n\)\}s\.\.

This creates non\-trivial nonlinear mixing among different frequency modes and rescue the representation power of the kernel in high dimensions when condition number is high by boosting the presence/occurrence of eigenmodes with small eigenvalues\. Explicitly, for cleaner presentation, the first two orders of contributions to the functiong​\(x\)g\(x\)can be written as a sum over different eigen modes of the kernel matrixϕm\\phi\_\{m\}and their corresponding eigenvalueλm\\lambda\_\{m\}\(in the case whenK~\\widetilde\{K\}is diagonalizable in the kernel eigenbasis\)999The case whereK~\\widetilde\{K\}is not diagonalizable in the kernel eigenbasis is slightly more involving, we present the easier version here for presentation purpose\.:

g\(0\)​\(x\)=∑i∑m=1Nλmηm\+λ​ϕmT​h​ϕm\\displaystyle g^\{\(0\)\}\(x\)=\\sum\_\{i\}\\sum\_\{m=1\}^\{N\}\\frac\{\\lambda\_\{m\}\}\{\\eta\_\{m\}\+\\lambda\}\\phi^\{T\}\_\{m\}h\\,\\phi\_\{m\}\(20\)g\(1\)​\(x\)=∑i∑m,n,p,q,rλm2​λn​λpηm\+λ​ϕqT​hηq\+λ​ϕrT​hηr\+λ×Φ~m​n​p​q​r,\\displaystyle g^\{\(1\)\}\(x\)=\\sum\_\{i\}\\sum\_\{m,n,p,q,r\}\\frac\{\\lambda^\{2\}\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}\}\{\\eta\_\{m\}\+\\lambda\}\\frac\{\\phi^\{T\}\_\{q\}h\}\{\\eta\_\{q\}\+\\lambda\}\\frac\{\\phi^\{T\}\_\{r\}h\}\{\\eta\_\{r\}\+\\lambda\}\\times\\tilde\{\\Phi\}\_\{mnpqr\}\\,,whereΦ~m​n​p​q​r\\tilde\{\\Phi\}\_\{mnpqr\}are various combinations of the eigenmode functions\.ηm\\eta\_\{m\}are square root of the eigenvalues of the kernel smoothedK~\\widetilde\{K\}in the kernel eigenbasis, which is slightly smaller thanλm\\lambda\_\{m\}due to the variance contraction of the conditional integral operator\. To simplify discussions, we have written the expression whenηm∼λm\\eta\_\{m\}\\sim\\lambda\_\{m\}, which essentially trivializes the conditional integral and brings us back to standard kernel ridge regression\. We expand further into the derivation of this in appendix[A](https://arxiv.org/html/2606.00322#A1)\. We note that generically,K~\\widetilde\{K\}would not be diagonalizable in the eigenbasis of the kernel matrix, which in some sense, induces spectral mixing between different eigenmodes of the kernel\. The spectral result we provide here is illustrating that on top of the linear mixing induced by the conditional expectation, the triple moment interaction term introduces additional non\-linear spectral mixing\.

In this simple setting, the order0contribution is the familiar spectral representation of standard kernel ridge regression\. We see that due to the presence of the regularization scaleλmλm\+λ\\frac\{\\lambda\_\{m\}\}\{\\lambda\_\{m\}\+\\lambda\}, the contributions tog​\(x\)g\(x\)only come from largerλm\\lambda\_\{m\}s while eigenmodes withλm≪λ\\lambda\_\{m\}\\ll\\lambdaare suppressed\. This is quantified by the condition numberκ=λmax/λmin\\kappa=\\lambda\_\{\\text\{max\}\}/\\lambda\_\{\\text\{min\}\}\. Highκ\\kappasignals the representation power of the kernel becomes concentrated to the first few large eigenvalues, effectively making other tail eigenmodes irrelevant \(due to the ridge regularization serving as a cut\-off threshold\)\. On the other hand, the first order contributiong\(1\)​\(x\)g^\{\(1\)\}\(x\)induces nonlinear spectral mixing, which even whenκ≫1\\kappa\\gg 1, reinstates the contribution from tail eigenmodes by pairing them with the larger eigenmodes101010We note that this is certainly not the only possible term on can add in the objective function, one can reverse engineer any term that does not violate the representer theorem starting with the desired spectral properties\.\. This remains true in the generic NPIV setting whereK~\\widetilde\{K\}is non\-diagonalizable in the eigenbasis of the kernel\.

### 4\.2Renormalization

As we mentioned earlier, the naive power series ansatz we assumed could be divergent, giving meaningless answers in the perturbative resultαi\\alpha\_\{i\}\. This can become dramatic especially when the inverse problem is more and more ill\-defined, the norms of the coefficient vectorsα\(k\)\\alpha^\{\(k\)\}can grow rapidly, leading to numerical instability\. Physics literature has explored multiple approaches to address this divergence\. One particular way, referred to as Wilsonian renormalization\(Wilson and Kogut,[1974](https://arxiv.org/html/2606.00322#bib.bib25)\), imposes a cut\-off/regularization parameter on physical systems \(the energy spectrum for example\) and rely on the belief that physical systems should not depend on such artificially imposed cut\-off\. This reveals the physically meaningful part of the divergent series which is stable under this process\.

Inspired by this, although the regularization parameter in the RKHS \(λ\\lambdaregularization strength in kernel ridge regression\) is not physical111111For readers familiar with the physics concept of renormalization might wonder if this is connected with the renormalization group equation, we discuss this further in appendix[C](https://arxiv.org/html/2606.00322#A3), asserting the differences with the usual renormalization and explaining why the physical assumption does not apply here\., we introduce renormalization parameters as a simple rescalingα=α\(0\)\+γ~1​α\(1\)\+γ~2​α\(2\)\+…\\alpha=\\alpha^\{\(0\)\}\+\\tilde\{\\gamma\}\_\{1\}\\alpha^\{\(1\)\}\+\\tilde\{\\gamma\}\_\{2\}\\alpha^\{\(2\)\}\+\\ldotswhereγ~k\\tilde\{\\gamma\}\_\{k\}are the renormalized coefficients\. Given by a simple rescaling:

γ~k=γk​min⁡\(‖α\(0\)‖‖α\(k\)‖,1\),\\tilde\{\\gamma\}\_\{k\}=\\gamma^\{k\}\\min\\left\(\\frac\{\\\|\\alpha^\{\(0\)\}\\\|\}\{\\\|\\alpha^\{\(k\)\}\\\|\},1\\right\)\\,,\(21\)as we shall further demonstrate in appendix[A\.4](https://arxiv.org/html/2606.00322#A1.SS4)using the spectral property of the kernel operators, the inevitable divergence of the naive power series, this simple rescaling tames the divergence\.

### 4\.3Resurgence

Yet another tool to tame potential divergence in the perturbative ansatz and extract useful information is resurgence\(Ecalle,[1981](https://arxiv.org/html/2606.00322#bib.bib26); Berry and Howls,[1991](https://arxiv.org/html/2606.00322#bib.bib41)\)\. In certain cases, the series could be an asymptotic series, which is factorially divergent in its coefficientsαi\(n\)\\alpha\_\{i\}^\{\(n\)\}\.

This happens simply because we have chosen the pointγ=0\\gamma=0to perturbatively expand \(which is assumed to be a saddle in the objective function landscapeSS\), this may very well be incorrect due to the presence of numerous other saddle points\. Luckily, resurgence comes to rescue, we give a much more detailed introduction on this subject practically in appendix[B](https://arxiv.org/html/2606.00322#A2)\. The belief of the resurgence trick is that divergent series are divergent for a reason, they diverge due to the existence of other nearby saddle points in the function landscape\. Naive perturbative expansion ”falls off the cliff” along the directions where the current saddle is a maxima\. Resurgence essentially extracts the nearby saddle point information from the divergent higher order terms which have factorially divergent coefficients using complex analysis, then rewriting the series to account for the existence of the nearby saddle points, we recover a trans\-series, which has finite radius of convergence\. Resurgence/Borel resummation finds ways to take factorially divergent series and resum them to an integral form\. Although this integral is not well\-defined around singularities in its domain, one can resolve this by adding additional series around those points suppressed by the ambiguity value\. This returns a convergent series with no ambiguity\. Curious readers are referred to appendix[B](https://arxiv.org/html/2606.00322#A2)for more detail\.

For now, on a practical level, we shall summerize things as an algorithm if we run into factorially divergent coefficient functions\. Check ifαi\(n\)∼n\!\\alpha\_\{i\}^\{\(n\)\}\\sim n\!, if false, then series not asymptotic, just stick to renormalized series\. If true, then series asymptotic, use Borel resummationg^​\(t​X\)=∑i∑n=0∞αi\(n\)n\!​γn​K​\(t​X,Xi\)\\hat\{g\}\(tX\)=\\sum\_\{i\}\\sum\_\{n=0\}^\{\\infty\}\\frac\{\\alpha\_\{i\}^\{\(n\)\}\}\{n\!\}\\gamma^\{n\}K\(tX,X\_\{i\}\)\. Then one examines singularity on complextt\-plane\. This gives a collection\{tk\}k=1m\\\{t\_\{k\}\\\}\_\{k=1\}^\{m\}singularities\. Computing the residue around those singularities returns a trans\-series with no ambiguity:

Φ​\(X\)=∑n=0∞α\(n\)​γn⏟original series\+∑k=1mCk​e−k​tk​\(X\)/X​gk​\(X\)⏟non\-perturbative near saddle,\\Phi\(X\)=\\underbrace\{\\sum\_\{n=0\}^\{\\infty\}\\alpha^\{\(n\)\}\\gamma\_\{n\}\}\_\{\\text\{original series\}\}\+\\underbrace\{\\sum\_\{k=1\}^\{m\}C\_\{k\}\\,e^\{\-k\\,t\_\{k\}\(X\)/X\}g\_\{k\}\(X\)\}\_\{\\text\{non\-perturbative near saddle\}\}\\,,\(22\)withCk=SC\+−SC−=Rest=tk​g^​\(t​X\)C\_\{k\}=S\_\{C\_\{\+\}\}\-S\_\{C\_\{\-\}\}=\\text\{Res\}\_\{t=t\_\{k\}\}\\hat\{g\}\(tX\)is the ambiguity of choosing complexified integration contour above or below the singularity named Stokes constant\.

The one line take away message of this part is that when the perturbative seriesαi​\(γ\)\\alpha\_\{i\}\(\\gamma\)is factorially divergent, one could use the method described above to obtain a series with finite \(nonzero\) radius of convergence and extract useful information from the originally wildly divergent series\.

## 5Implementation

To summerize the perturbative procedures we discussed so far, we practically have two stages\. First compute the kernel matrices needed, then iteratively compute the raw RKHS coefficientsαi​\(γ\)\\alpha\_\{i\}\(\\gamma\)as a power series inγ\\gammaup to orderNN\. Then we handle the case if the coefficientsαi\(n\)\\alpha\_\{i\}^\{\(n\)\}are divergent, we first use rescaling\. If they are factorially divergent, we further use resurgence\.

### 5\.1Computing the Kernel Matrices

For implementation, one can use the standard Nadaraya\-Watson estimator to estimate the two conditional densities𝔼​\[Y\|Z\]\\mathbb\{E\}\[Y\|Z\]andf​\(X\|Z\)f\(X\|Z\)after picking the kernel\. We do this for the triple moment kernel:

Ki​j​k≈1n​∑r=1nK​\(xr,xi\)​K​\(xr,xj\)​K​\(xr,xk\)⋅wr,K\_\{ijk\}\\approx\\frac\{1\}\{n\}\\sum\_\{r=1\}^\{n\}K\(x\_\{r\},x\_\{i\}\)K\(x\_\{r\},x\_\{j\}\)K\(x\_\{r\},x\_\{k\}\)\\cdot w\_\{r\}\\,,\(23\)wherewrw\_\{r\}is the weight that account for the conditional densityf​\(x\|Z\)f\(x\|Z\), which is estimated using standard Nadaraya\-Watson estimator\. And similarly for

𝔼​\[Y\|Z=z\]≈∑i=1nl​\(z,zi\)​yi∑i=1nl​\(z,zi\),\\mathbb\{E\}\[Y\|Z=z\]\\approx\\frac\{\\sum\_\{i=1\}^\{n\}l\(z,z\_\{i\}\)y\_\{i\}\}\{\\sum\_\{i=1\}^\{n\}l\(z,z\_\{i\}\)\}\\,,\(24\)wherel​\(z,z′\)l\(z,z^\{\\prime\}\)is a kernel function for the instrument space\. More traditionally in the usual kernel ridge regression, instead of directly estimating the densitybi​\(Z\)=𝔼​\[K​\(x,xi\)\|Z\]b\_\{i\}\(Z\)=\\mathbb\{E\}\[K\(x,x\_\{i\}\)\|Z\], people use the two stage method where the first stage is to compute the conditional mean embedding directly in RKHS via the covariance operator\. For eachz∈Zz\\in Z, the conditional feature vector isb​\(z\)=KX​\(L\+n​ρ​I\)−1​lzb\(z\)=K\_\{X\}\(L\+n\\rho I\)^\{\-1\}l\_\{z\}, where\(KX\)i​j=k​\(xi,xj\)\(K\_\{X\}\)\_\{ij\}=k\(x\_\{i\},x\_\{j\}\)andLk​l=l​\(zk,zl\)L\_\{kl\}=l\(z\_\{k\},z\_\{l\}\)are Gram matrices on X and Z,lz=\[l​\(z1,z\),…,l​\(zn,z\)\]Tl\_\{z\}=\[l\(z\_\{1\},z\),\\dots,l\(z\_\{n\},z\)\]^\{T\}denotes the row vector,nnis the dimension of the instrumental whileρ\\rhois some regularization parameter\. Then computing𝔼​\[Y\|Z\]=1n​lzT​\(L\+n​ρ​I\)−1​y\\mathbb\{E\}\[Y\|Z\]=\\frac\{1\}\{n\}\\,l\_\{z\}^\{T\}\(L\+n\\rho I\)^\{\-1\}y\. The desired quantities are computed usingbi​\(Z\)b\_\{i\}\(Z\)and𝔼​\[Y\|Z\]\\mathbb\{E\}\[Y\|Z\]:

hi=𝔼​\[bi​𝔼​\[Y\|Z\]\],Ki​j=𝔼​\[bi​bj\]=1n​∑l=1nbi​\(zl\)​bj​\(zl\)\.h\_\{i\}=\\mathbb\{E\}\[b\_\{i\}\\mathbb\{E\}\[Y\|Z\]\]\\,,K\_\{ij\}=\\mathbb\{E\}\[b\_\{i\}b\_\{j\}\]=\\frac\{1\}\{n\}\\sum\_\{l=1\}^\{n\}b\_\{i\}\(z\_\{l\}\)b\_\{j\}\(z\_\{l\}\)\\,\.

### 5\.2Algorithm

For prediction on new pointsxn​e​wx\_\{new\}, we compute:

g​\(xn​e​w\)=∑iαi\(0\)​K​\(xn​e​w,xi\)\+∑iγ1​αi\(1\)​K​\(xn​e​w,xi\)\+∑iγ2​αi\(2\)​K​\(xn​e​w,xi\)\+…,g\(x\_\{new\}\)=\\sum\_\{i\}\\alpha\_\{i\}^\{\(0\)\}K\(x\_\{new\},x\_\{i\}\)\+\\sum\_\{i\}\\gamma\_\{1\}\\alpha\_\{i\}^\{\(1\)\}K\(x\_\{new\},x\_\{i\}\)\\\\ \+\\sum\_\{i\}\\gamma\_\{2\}\\alpha\_\{i\}^\{\(2\)\}K\(x\_\{new\},x\_\{i\}\)\+\\ldots\\,,\(25\)whereαi\(n\)\\alpha\_\{i\}^\{\(n\)\}are the order by order RKHS coefficients induced by the triple moment andγn\\gamma\_\{n\}are the renormalized coupling powers\. The sum is truncated at a certain orderNN, where the marginal gain of computing a new order no longer beats the increase in computation time\.

The general procedure can be summerized in the following fashion\. We first compute the needed kernel matricesKi​jK\_\{ij\}, triple momentKi​j​kK\_\{ijk\}and conditional mean𝔼​\[Y\|Z\]\\mathbb\{E\}\[Y\|Z\]through conditional mean embedding and Nadaraya\-Watson estimator\. Then we compute the perturbative series in kernel coefficientsαi=∑n=0Nγn​αi\(n\)\\alpha\_\{i\}=\\sum\_\{n=0\}^\{N\}\\gamma\_\{n\}\\,\\alpha\_\{i\}^\{\(n\)\}order by order solving ridge regression problems up to a cut\-off orderNN\. If we detect divergence in the series, rescale \(renormalize\) the coupling strength adaptively, if the series is factorially divergent, use resurgence algorithm\. A detailed description for the first three steps is as in[1](https://arxiv.org/html/2606.00322#alg1), which has an time complexity of𝒪​\(n3×M\)\\mathcal\{O\}\(n^\{3\}\\times M\)withnndata points computed toMMorders, prior to this, computingKi​j​kK\_\{ijk\}requires a further𝒪​\(n4\)\\mathcal\{O\}\(n^\{4\}\)complexity computation\. We deferred the detail of last resurgence step to the appendix[2](https://arxiv.org/html/2606.00322#alg2)\.

## 6Experimental Results

We conducted extensive experiments to evaluate the performance of our perturbative renormalization approach across different dimensionality regimes\. The experiments were done withn=80n=80data points generated around the ground truth, which are the datasets described in appendix[F\.1](https://arxiv.org/html/2606.00322#A6.SS1), which create a series of challenging high\-dimensional NPIV problem with controlled properties including standard IV benchmarks\. Here we briefly summerize the overall performance of the algorithms across different kernel types: Gaussian RBF kernel and fractional Brownian kernel\. As a comparison, we also add the fitting RMSE of Deep IV method\(Hartfordet al\.,[2017](https://arxiv.org/html/2606.00322#bib.bib17)\)on the same dataset after500500epochs of training each regression stage with approximately 2k parameter neural nets\. We leave the complete tables of data across dimensions and renormalization strength in appendix[F\.3](https://arxiv.org/html/2606.00322#A6.SS3)\.

Table 1:NPIV best Performance Scaling with Dimension- •β\\betacontrols dimensiond=nβd=n^\{\\beta\}wheren=80n=80; values taken over best amongγ∈\{0\.4,0\.6,0\.8,0\.99\}\\gamma\\in\\\{0\.4,0\.6,0\.8,0\.99\\\}
- •”Order 0” is kernel ridge regression RMSE121212The results presented are all cross validated for different choices of the L2 regularization parameterλ\\lambda\.; ”Best Pert” refers to the best RMSE obtained using perturbative method; ”BestΔ%\\Delta\\%” represents the percentage gain in RMSE\.

![Refer to caption](https://arxiv.org/html/2606.00322v1/RBF_example.png)Figure 2:Example RBF kernel fitting![Refer to caption](https://arxiv.org/html/2606.00322v1/FB_example.png)Figure 3:Example fractional Brownian kernel fittingIn appendix[F\.4](https://arxiv.org/html/2606.00322#A6.SS4), we also include further standard NPIV datasets: Newey\-Powell\(Newey and Powell,[2003](https://arxiv.org/html/2606.00322#bib.bib21)\), weak/strong instrumental datasets, heteroscedastic dataset, nonlinear instrumental dataset and sparse signal dataset\. The performance result is summerized here, note that we only focus on the improvement against regularized kernel IV \(order0\) using fractional Brownian kernel, Deep IV method is there for completeness:

Table 2:RMSE Comparison on Alternative NPIV Datasets- •RMSE =MSE\\sqrt\{\\text\{MSE\}\}; “Order 0” is standard kernel ridge IV with fractional Brownian kernel; “Best Pert\.” is the best result acrossγ∈\{0\.4,0\.6,0\.8,0\.99\}\\gamma\\in\\\{0\.4,0\.6,0\.8,0\.99\\\}with fractional Brownian kernel\. Bold indicates best method per dataset\.

In addition to these, we further include empirical analysis on the performance of the algorithm on IV datasets with weak instrumentals in appendix[F\.5](https://arxiv.org/html/2606.00322#A6.SS5), a sensitivity analysis on the parameters used in the algorithm \(regularization parameterγ\\gamma, ridge regularization parameterλ\\lambdaand maximum order of perturbationNmaxN\_\{\\text\{max\}\}\) in appendix[F\.6](https://arxiv.org/html/2606.00322#A6.SS6)and experiments with larger sample size in appendix[F\.7](https://arxiv.org/html/2606.00322#A6.SS7)\.

Key Findings:The RBF kernel shows minimal improvement across different kernel widths\. If one plots the fitted curve on cross sections of the hypersurface as in figure[2](https://arxiv.org/html/2606.00322#S6.F2), where the cross section is taken with respect to dimension11, we see that RBF kernel almost zeros out due to its rotational invariance forces concentration in high dimensions, this almost trivial fitting gives a convenient baseline MSE for eachβ\\betacase\. We see in the table[1](https://arxiv.org/html/2606.00322#S6.T1)that perturbative method we have here also matches that baseline, it simply overlaps with the trivial order0kernel ridge regression fit in[2](https://arxiv.org/html/2606.00322#S6.F2)since the added interaction is also constrained by rotational invariance\. In appendix[A](https://arxiv.org/html/2606.00322#A1)we further discuss why the extra term we add works from a spectral perspective, while explaining why it respects rotational invariance hence contributing minimally to the fitted result\. We further propose another different kernel that breaks this invariance hence helping the performance of RBF kernels in lower dimensions in appendix[E](https://arxiv.org/html/2606.00322#A5)\. Although in high dimensions, the curse of dimensionality forces everything to fit trivially, making RBF kernels beyond rescue\.

The fractional Brownian kernel, on the other hand, shows superior performance over the trivial RBF kernel fit at order0when dimension is lowβ<0\.7\\beta<0\.7as in the first two rows of the table[1](https://arxiv.org/html/2606.00322#S6.T1)\. Atβ≥0\.7\\beta\\geq 0\.7, we see that it actually has worse RMSE than the trivial RBF fit, as plotted in[3](https://arxiv.org/html/2606.00322#S6.F3), this is due to the non\-trivial wiggle attempting to non\-trivially fit the curve\. The perturbative method effectively brings the fitted curve down to the true function, significantly improving the fitted curve\. The perturbatively revised fitting consistently outperforms the trivial RBF kernel baseline as well\. It is not obvious from the diagram whether resurgence was also used \(referred to as ”Borel” in the legend\), since it completely overlaps with the pure perturbative and renormalization \(”Renormalized”\) line, the mechanism only comes into effect when the series is factorially divergent, which is not the case in majority of the cases\.

We note that the two algorithms we wrote[1](https://arxiv.org/html/2606.00322#alg1),[2](https://arxiv.org/html/2606.00322#alg2)need to be comprehensively used to maximize their effectiveness\. In that[2](https://arxiv.org/html/2606.00322#alg2)is only useful when the perturbative series in kernel coefficientsα\\alphais asymptotic or factorially divergent\. The resummation hurts the performance when this is not the case\. It is not hard to see that this depends on the perturbative strength parameterγ\\gamma, it is easier to obtain an asymptotic series whenγ\\gammais closer to11\.γ→1\\gamma\\to 1suggests more significant contributions from higher order terms in[1](https://arxiv.org/html/2606.00322#alg1), which might hurt the its performance\. Hence there is a subtle empirical balance in the choice ofγ\\gammato balance them\.

Overall, solving the integral equation amplifies the ill\-conditioned\-ness of the kernel matrices in high dimensions \(highβ\\beta/dimension\-sample ratio\)\. Perturbative approach[1](https://arxiv.org/html/2606.00322#alg1)rescues this by signifying the contribution from eigendirections with small eigenvalues, hence enhancing the expressivity of the kernel\. Resurgence algorithm[2](https://arxiv.org/html/2606.00322#alg2)further corrects the behavior of the perturbative series when it is factorially divergent\. This is a novel method of extracting the most out of kernel ridge regression\. We note that the form of the additional term we added in the objective \([11](https://arxiv.org/html/2606.00322#S4.E11)\) can certainly be replaced with other ones, spectrally one could tailor the additional term needed in order for the induced perturbative solution to enhance a certain aspect of typical kernel ridge regression\. We demonstrate this in appendix[E](https://arxiv.org/html/2606.00322#A5), where a new term in the objective is introduced that explicitly breaks rotational invariance, tailored to tackle the extreme vanishing of rotationally invariant kernels in high dimensions\. We hope this paper–although as a crude first try–demonstrates the applicability of perturbative methods \(and the methods taming them\) in statistical/econometrics world131313Rigorous statistical proofs are not present, such as those finite sample bounds for standard kernel ridge regression\(Hang and Steinwart,[2018](https://arxiv.org/html/2606.00322#bib.bib4); Fischer and Steinwart,[2020](https://arxiv.org/html/2606.00322#bib.bib2); Caponnetto and de Vito,[2007](https://arxiv.org/html/2606.00322#bib.bib3)\), this is due to the proposed algorithm does not satisfy the standard assumptions under which most of the estimator efficiency proofs are constructed, we leave this for future work\.\.

## References

- \[1\]External Links:ISSN 00129682, 14680262,[Link](http://www.jstor.org/stable/2171753)Cited by:[§F\.5](https://arxiv.org/html/2606.00322#A6.SS5.p1.5)\.
- I\. Andrews, J\. H\. Stock, and L\. Sun \(2019\)Weak instruments in instrumental variables regression: theory and practice\.Annual Review of Economics11\(Volume 11, 2019\),pp\. 727–753\.External Links:[Document](https://dx.doi.org/https%3A//doi.org/10.1146/annurev-economics-080218-025643),[Link](https://www.annualreviews.org/content/journals/10.1146/annurev-economics-080218-025643),ISSN 1941\-1391Cited by:[§4](https://arxiv.org/html/2606.00322#S4.p1.4)\.
- F\. Bach \(2013\)Sharp analysis of low\-rank kernel matrix approximations\.External Links:1208\.2015,[Link](https://arxiv.org/abs/1208.2015)Cited by:[§1](https://arxiv.org/html/2606.00322#S1.p3.1)\.
- A\. Belloni, D\. Chen, V\. Chernozhukov, and C\. Hansen \(2015\)Sparse models and methods for optimal instruments with an application to eminent domain\.External Links:1010\.4345,[Link](https://arxiv.org/abs/1010.4345)Cited by:[§1\.1](https://arxiv.org/html/2606.00322#S1.SS1.p1.1)\.
- A\. Bennett, N\. Kallus, and T\. Schnabel \(2020\)Deep generalized method of moments for instrumental variable analysis\.External Links:1905\.12495,[Link](https://arxiv.org/abs/1905.12495)Cited by:[§1\.1](https://arxiv.org/html/2606.00322#S1.SS1.p1.1)\.
- M\. V\. Berry and C\. J\. Howls \(1991\)Hyperasymptotics for integrals with saddles\.Proceedings of the Royal Society of London\. Series A: Mathematical and Physical Sciences434\(1892\),pp\. 657–675\.External Links:[Document](https://dx.doi.org/10.1098/rspa.1991.0119),[Link](https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1991.0119),https://royalsocietypublishing\.org/doi/pdf/10\.1098/rspa\.1991\.0119Cited by:[Appendix B](https://arxiv.org/html/2606.00322#A2.p1.1),[§4\.3](https://arxiv.org/html/2606.00322#S4.SS3.p1.1)\.
- A\. Bhattacharya, J\. Cotler, A\. Dersy, and M\. D\. Schwartz \(2024\)Renormalons as Saddle Points\.External Links:2410\.07351Cited by:[§B\.1](https://arxiv.org/html/2606.00322#A2.SS1.p7.1)\.
- J\. Bound, D\. A\. Jaeger, and R\. M\. Baker \(1995\)Problems with instrumental variables estimation when the correlation between the instruments and the endogenous explanatory variable is weak\.Journal of the American Statistical Association90\(430\),pp\. 443–450\.External Links:[Document](https://dx.doi.org/10.1080/01621459.1995.10476536),[Link](https://doi.org/10.1080/01621459.1995.10476536),https://doi\.org/10\.1080/01621459\.1995\.10476536Cited by:[§F\.5](https://arxiv.org/html/2606.00322#A6.SS5.p1.5)\.
- A\. Caponnetto and E\. de Vito \(2007\)Optimal rates for the regularized least\-squares algorithm\.Foundations of Computational Mathematics7,pp\. 331–368\.External Links:[Link](https://api.semanticscholar.org/CorpusID:207063850)Cited by:[footnote 13](https://arxiv.org/html/2606.00322#footnote13)\.
- M\. Carrasco, J\. Florens, and E\. Renault \(2007\)Chapter 77 linear inverse problems in structural econometrics estimation based on spectral decomposition and regularization\.J\. J\. Heckman and E\. E\. Leamer \(Eds\.\),Handbook of Econometrics, Vol\.6,pp\. 5633–5751\.External Links:ISSN 1573\-4412,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/S1573-4412%2807%2906077-1),[Link](https://www.sciencedirect.com/science/article/pii/S1573441207060771)Cited by:[§1](https://arxiv.org/html/2606.00322#S1.p3.1)\.
- Z\. Chen, A\. Mustafi, P\. Glaser, A\. Korba, A\. Gretton, and B\. K\. Sriperumbudur \(2024\)\(De\)\-regularized maximum mean discrepancy gradient flow\.arXiv preprint arXiv:2409\.14980\.Cited by:[footnote 2](https://arxiv.org/html/2606.00322#footnote2)\.
- O\. Costin, G\. V\. Dunne, A\. Gruen, and S\. Gukov \(2023\)Going to the other side via the resurgent bridge\.External Links:2310\.12317,[Link](https://arxiv.org/abs/2310.12317)Cited by:[Appendix B](https://arxiv.org/html/2606.00322#A2.p1.1)\.
- S\. Darolles, Y\. Fan, J\. P\. Florens, and E\. Renault \(2011\)NONPARAMETRIC instrumental regression\.Econometrica79\(5\),pp\. 1541–1565\.External Links:ISSN 00129682, 14680262,[Link](http://www.jstor.org/stable/41237784)Cited by:[§1](https://arxiv.org/html/2606.00322#S1.p2.1)\.
- E\. de Faria and W\. de Melo \(2010\)Perturbative quantum field theory\.InMathematical Aspects of Quantum Field Theory,Cambridge Studies in Advanced Mathematics,pp\. 153–191\.Cited by:[§1](https://arxiv.org/html/2606.00322#S1.p5.3),[footnote 6](https://arxiv.org/html/2606.00322#footnote6)\.
- E\. Delabaere and F\. Pham \(1999\)Resurgent methods in semi\-classical asymptotics\.Annales de l’I\.H\.P\. Physique théorique71\(1\),pp\. 1–94\(en\)\.External Links:[Link](https://www.numdam.org/item/AIHPA_1999__71_1_1_0/),[MathReview Entry](https://www.ams.org/mathscinet-getitem?mr=1704654)Cited by:[Appendix B](https://arxiv.org/html/2606.00322#A2.p1.1)\.
- K\. Donhauser, M\. Wu, and F\. Yang \(2021\)How rotational invariance of common kernels prevents generalization in high dimensions\.InProceedings of the 38th International Conference on Machine Learning,M\. Meila and T\. Zhang \(Eds\.\),Proceedings of Machine Learning Research, Vol\.139,pp\. 2804–2814\.External Links:[Link](https://proceedings.mlr.press/v139/donhauser21a.html)Cited by:[§A\.3](https://arxiv.org/html/2606.00322#A1.SS3.p1.1),[§F\.1](https://arxiv.org/html/2606.00322#A6.SS1.p1.1),[§F\.2](https://arxiv.org/html/2606.00322#A6.SS2.p1.2),[§F\.4](https://arxiv.org/html/2606.00322#A6.SS4.p1.1),[§1\.1](https://arxiv.org/html/2606.00322#S1.SS1.p1.1),[§1](https://arxiv.org/html/2606.00322#S1.p2.1),[§1](https://arxiv.org/html/2606.00322#S1.p3.1),[§1](https://arxiv.org/html/2606.00322#S1.p4.1)\.
- D\. Dorigoni \(2019\)An introduction to resurgence, trans\-series and alien calculus\.Annals of Physics409,pp\. 167914\.Note:Accessible survey for physicists covering resurgent analysis and alien calculusExternal Links:[Document](https://dx.doi.org/10.1016/j.aop.2019.167914)Cited by:[Appendix B](https://arxiv.org/html/2606.00322#A2.p1.1)\.
- F\. J\. Dyson \(1952\)Divergence of perturbation theory in quantum electrodynamics\.Phys\. Rev\.85,pp\. 631–632\.External Links:[Document](https://dx.doi.org/10.1103/PhysRev.85.631)Cited by:[footnote 6](https://arxiv.org/html/2606.00322#footnote6)\.
- J\. Ecalle \(1981\)Les fonctions résurgentes: \(en trois parties\)\.Les fonctions résurgentes:,Université de Paris\-Sud, Département de Mathématique, Bât\. 425\.External Links:[Link](https://books.google.co.uk/books?id=oDfvAAAAMAAJ)Cited by:[Appendix B](https://arxiv.org/html/2606.00322#A2.p1.1),[§4\.3](https://arxiv.org/html/2606.00322#S4.SS3.p1.1)\.
- S\. Fischer and I\. Steinwart \(2020\)Sobolev norm learning rates for regularized least\-squares algorithm\.External Links:1702\.07254,[Link](https://arxiv.org/abs/1702.07254)Cited by:[footnote 13](https://arxiv.org/html/2606.00322#footnote13)\.
- O\. Hagrass, B\. K\. Sriperumbudur, and B\. Li \(2024\)Spectral regularized kernel two\-sample tests\.External Links:2212\.09201,[Link](https://arxiv.org/abs/2212.09201)Cited by:[footnote 2](https://arxiv.org/html/2606.00322#footnote2)\.
- P\. Hall and J\. L\. Horowitz \(2005\)Nonparametric methods for inference in the presence of instrumental variables\.The Annals of Statistics33\(6\)\.External Links:ISSN 0090\-5364,[Link](http://dx.doi.org/10.1214/009053605000000714),[Document](https://dx.doi.org/10.1214/009053605000000714)Cited by:[§1\.1](https://arxiv.org/html/2606.00322#S1.SS1.p1.1),[§1](https://arxiv.org/html/2606.00322#S1.p2.1)\.
- H\. Hang and I\. Steinwart \(2018\)Optimal learning with anisotropic gaussian svms\.External Links:1810\.02321,[Link](https://arxiv.org/abs/1810.02321)Cited by:[footnote 13](https://arxiv.org/html/2606.00322#footnote13)\.
- J\. Hartford, G\. Lewis, K\. Leyton\-Brown, and M\. Taddy \(2017\)Deep IV: a flexible approach for counterfactual prediction\.InProceedings of the 34th International Conference on Machine Learning,D\. Precup and Y\. W\. Teh \(Eds\.\),Proceedings of Machine Learning Research, Vol\.70,pp\. 1414–1423\.External Links:[Link](https://proceedings.mlr.press/v70/hartford17a.html)Cited by:[§1\.1](https://arxiv.org/html/2606.00322#S1.SS1.p1.1),[§6](https://arxiv.org/html/2606.00322#S6.p1.2)\.
- J\. Kim, D\. Meunier, A\. Gretton, T\. Suzuki, and Z\. Li \(2025\)Optimality and adaptivity of deep neural features for instrumental variable regression\.External Links:2501\.04898,[Link](https://arxiv.org/abs/2501.04898)Cited by:[§1\.1](https://arxiv.org/html/2606.00322#S1.SS1.p1.1),[§1](https://arxiv.org/html/2606.00322#S1.p4.1)\.
- L\. N\. Lipatov \(1977\)Divergence of the Perturbation Theory Series and the Quasiclassical Theory\.Sov\. Phys\. JETP45,pp\. 216–223\.Cited by:[footnote 6](https://arxiv.org/html/2606.00322#footnote6)\.
- J\. Magnen and R\. Seneor \(1977\)Phase Space Cell Expansion and Borel Summability for the Euclidean phi\*\*4 in Three\-Dimensions Theory\.Commun\. Math\. Phys\.56,pp\. 237\.External Links:[Document](https://dx.doi.org/10.1007/BF01614211)Cited by:[footnote 6](https://arxiv.org/html/2606.00322#footnote6)\.
- M\. Mariño \(2014\)Lectures on non‐perturbative effects in large n gauge theories, matrix models and strings\.Fortschritte der Physik62\(5–6\),pp\. 455–540\.External Links:ISSN 1521\-3978,[Link](http://dx.doi.org/10.1002/prop.201400005),[Document](https://dx.doi.org/10.1002/prop.201400005)Cited by:[Appendix B](https://arxiv.org/html/2606.00322#A2.p1.1)\.
- D\. Meunier, A\. Moulin, J\. Wornbard, V\. R\. Kostic, and A\. Gretton \(2025\)Demystifying spectral feature learning for instrumental variable regression\.External Links:2506\.10899,[Link](https://arxiv.org/abs/2506.10899)Cited by:[footnote 1](https://arxiv.org/html/2606.00322#footnote1)\.
- W\. Miao, Z\. Geng, and E\. T\. Tchetgen \(2018\)Identifying causal effects with proxy variables of an unmeasured confounder\.External Links:1609\.08816,[Link](https://arxiv.org/abs/1609.08816)Cited by:[footnote 3](https://arxiv.org/html/2606.00322#footnote3)\.
- K\. Muandet, A\. Mehrjou, S\. K\. Lee, and A\. Raj \(2020\)Dual instrumental variable regression\.External Links:1910\.12358,[Link](https://arxiv.org/abs/1910.12358)Cited by:[§1\.1](https://arxiv.org/html/2606.00322#S1.SS1.p1.1),[§1](https://arxiv.org/html/2606.00322#S1.p4.1)\.
- W\. K\. Newey and J\. L\. Powell \(2003\)Instrumental variable estimation of nonparametric models\.Econometrica71\(5\),pp\. 1565–1578\.External Links:ISSN 00129682, 14680262,[Link](http://www.jstor.org/stable/1555512)Cited by:[§F\.4](https://arxiv.org/html/2606.00322#A6.SS4.SSS0.Px1.p1.1),[§1\.1](https://arxiv.org/html/2606.00322#S1.SS1.p1.1),[§1](https://arxiv.org/html/2606.00322#S1.p2.1),[§6](https://arxiv.org/html/2606.00322#S6.p2.1)\.
- M\. E\. Peskin and D\. V\. Schroeder \(1995\)An Introduction to quantum field theory\.Addison\-Wesley,Reading, USA\.External Links:[Document](https://dx.doi.org/10.1201/9780429503559),ISBN 978\-0\-201\-50397\-5, 978\-0\-429\-50355\-9, 978\-0\-429\-49417\-8Cited by:[§1](https://arxiv.org/html/2606.00322#S1.p5.3),[footnote 6](https://arxiv.org/html/2606.00322#footnote6)\.
- M\. L\. Rizzo and G\. J\. Székely \(2016\)Energy distance\.WIREs Comput\. Stat\.8\(1\),pp\. 27–38\.External Links:ISSN 1939\-5108Cited by:[§F\.2](https://arxiv.org/html/2606.00322#A6.SS2.p2.6)\.
- B\. Schölkopf, R\. Herbrich, and A\. J\. Smola \(2001\)A generalized representer theorem\.InComputational Learning Theory,D\. Helmbold and B\. Williamson \(Eds\.\),Berlin, Heidelberg,pp\. 416–426\.External Links:ISBN 978\-3\-540\-44581\-4Cited by:[§3](https://arxiv.org/html/2606.00322#S3.p2.6)\.
- D\. Sejdinovic, B\. Sriperumbudur, A\. Gretton, and K\. Fukumizu \(2013\)Equivalence of distance\-based and rkhs\-based statistics in hypothesis testing\.The Annals of Statistics41\(5\)\.External Links:ISSN 0090\-5364,[Link](http://dx.doi.org/10.1214/13-AOS1140),[Document](https://dx.doi.org/10.1214/13-aos1140)Cited by:[§F\.2](https://arxiv.org/html/2606.00322#A6.SS2.p2.6),[§1](https://arxiv.org/html/2606.00322#S1.p6.1)\.
- R\. Singh, M\. Sahani, and A\. Gretton \(2020\)Kernel instrumental variable regression\.External Links:1906\.00232,[Link](https://arxiv.org/abs/1906.00232)Cited by:[§1\.1](https://arxiv.org/html/2606.00322#S1.SS1.p1.1),[§1](https://arxiv.org/html/2606.00322#S1.p2.1)\.
- I\. Steinwart and A\. Christmann \(2008\)Support vector machines\.Information science and statistics,Springer,New York, NY\(eng\)\.External Links:ISBN 978\-0\-387\-77241\-7 and 978\-1\-4899\-8963\-5 and 978\-6\-611\-92704\-2 and 978\-0\-387\-77242\-4Cited by:[§1](https://arxiv.org/html/2606.00322#S1.p3.1)\.
- J\. H\. Stock, J\. H\. Wright, and M\. Yogo \(2002a\)A survey of weak instruments and weak identification in generalized method of moments\.Journal of Business & Economic Statistics20\(4\),pp\. 518–529\.External Links:[Document](https://dx.doi.org/10.1198/073500102288618658),[Link](https://doi.org/10.1198/073500102288618658),https://doi\.org/10\.1198/073500102288618658Cited by:[§F\.5](https://arxiv.org/html/2606.00322#A6.SS5.p1.5)\.
- J\. H\. Stock, J\. H\. Wright, and M\. Yogo \(2002b\)A survey of weak instruments and weak identification in generalized method of moments\.Journal of Business & Economic Statistics20\(4\),pp\. 518–529\.External Links:[Document](https://dx.doi.org/10.1198/073500102288618658),[Link](https://doi.org/10.1198/073500102288618658),https://doi\.org/10\.1198/073500102288618658Cited by:[§4](https://arxiv.org/html/2606.00322#S4.p1.4)\.
- D\. Volin \(2010\)From the mass gap in O\(N\) to the non\-Borel\-summability in O\(3\) and O\(4\) sigma\-models\.Phys\. Rev\. D81,pp\. 105008\.External Links:0904\.2744,[Document](https://dx.doi.org/10.1103/PhysRevD.81.105008)Cited by:[footnote 6](https://arxiv.org/html/2606.00322#footnote6)\.
- K\. G\. Wilson and J\. B\. Kogut \(1974\)The Renormalization group and the epsilon expansion\.Phys\. Rept\.12,pp\. 75–199\.External Links:[Document](https://dx.doi.org/10.1016/0370-1573%2874%2990023-4)Cited by:[§4\.2](https://arxiv.org/html/2606.00322#S4.SS2.p1.2)\.
- H\. Wiltzer, J\. Farebrother, A\. Gretton, and M\. Rowland \(2024\)Foundations of multivariate distributional reinforcement learning\.External Links:2409\.00328,[Link](https://arxiv.org/abs/2409.00328)Cited by:[§4](https://arxiv.org/html/2606.00322#S4.p1.4)\.
- L\. Xu, Y\. Chen, S\. Srinivasan, N\. de Freitas, A\. Doucet, and A\. Gretton \(2023\)Learning deep features in instrumental variable regression\.External Links:2010\.07154,[Link](https://arxiv.org/abs/2010.07154)Cited by:[§1\.1](https://arxiv.org/html/2606.00322#S1.SS1.p1.1),[§1](https://arxiv.org/html/2606.00322#S1.p4.1),[§4](https://arxiv.org/html/2606.00322#S4.p1.4)\.

## Appendix ASpectral explanation of the triple moment term

### A\.1Mercer Decomposition and Covariance Operator

Every positive definite kernel admits the representation:

k​\(x,y\)=⟨φ​\(x\),φ​\(y\)⟩=∑m=0∞λm​gm​\(x\)​gm​\(y\),k\(x,y\)=\\langle\\varphi\(x\),\\varphi\(y\)\\rangle=\\sum\_\{m=0\}^\{\\infty\}\\lambda\_\{m\}g\_\{m\}\(x\)g\_\{m\}\(y\)\\,,\(26\)whereφ:𝒳→ℋ\\varphi:\\mathcal\{X\}\\to\\mathcal\{H\}is the feature map,\{gm:𝒳→ℝ\}\\\{g\_\{m\}:\\mathcal\{X\}\\to\\mathbb\{R\}\\\}are orthonormal eigenfunctions of the integral operatorB:L2​\(PX\)→L2​\(PX\)B:L^\{2\}\(P\_\{X\}\)\\to L^\{2\}\(P\_\{X\}\):

B​gm​\(x\)=∫𝒳k​\(x,y\)​gm​\(y\)​𝑑PX=λm​gm​\(x\)Bg\_\{m\}\(x\)=\\int\_\{\\mathcal\{X\}\}k\(x,y\)g\_\{m\}\(y\)dP\_\{X\}=\\lambda\_\{m\}g\_\{m\}\(x\)\(27\)wherePXP\_\{X\}denotes the marginal density ofXX\. We shall usePZP\_\{Z\}to denote the marginal density of the instrumentalZZto distinguish with the conditional integral operatorT:L2​\(PX\)→L2​\(PZ\)T:L^\{2\}\(P\_\{X\}\)\\to L^\{2\}\(P\_\{Z\}\)\.\{λm\}\\\{\\lambda\_\{m\}\\\}are eigenvalues of the covariance operator:

Σ=∫𝒳φ​\(x\)⊗φ​\(x\)​𝑑μ​\(x\)\.\\Sigma=\\int\_\{\\mathcal\{X\}\}\\varphi\(x\)\\otimes\\varphi\(x\)d\\mu\(x\)\\,\.\(28\)With finite data\{xi\}i=1n\\\{x\_\{i\}\\\}\_\{i=1\}^\{n\}, we approximate:

Ki​j\\displaystyle K\_\{ij\}=k​\(xi,xj\)=∑m=1nλm​ϕm​\[i\]​ϕm​\[j\]\\displaystyle=k\(x\_\{i\},x\_\{j\}\)=\\sum\_\{m=1\}^\{n\}\\lambda\_\{m\}\\phi\_\{m\}\[i\]\\phi\_\{m\}\[j\]\(29\)ϕm​\[i\]\\displaystyle\\phi\_\{m\}\[i\]≈gm​\(xi\),\\displaystyle\\approx g\_\{m\}\(x\_\{i\}\)\\,,\(30\)whereϕm∈ℝn\\phi\_\{m\}\\in\\mathbb\{R\}^\{n\}are eigenvectors of the kernel matrixKK\.

### A\.2Order\-by\-Order Construction

Recall the perturbative series we constructed

g​\(x\)=∑i∑n=0∞γn​α\(n\)​\(xi\)​K​\(x,xi\),g\(x\)=\\sum\_\{i\}\\sum\_\{n=0\}^\{\\infty\}\\gamma^\{n\}\\,\\alpha^\{\(n\)\}\(x\_\{i\}\)K\(x,x\_\{i\}\)\\,,\(31\)was solved iteratively:

Order 0:\(K~\+λ​K\)​α\(0\)=h\\displaystyle\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(0\)\}=hOrder 1:\(K~\+λ​K\)​α\(1\)=−3​∑j,kK3​\[i,j,k\]​αj\(0\)​αk\(0\)\\displaystyle\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(1\)\}=\-3\\sum\_\{j,k\}K\_\{3\}\[i,j,k\]\\alpha^\{\(0\)\}\_\{j\}\\alpha^\{\(0\)\}\_\{k\}Order n:\(K~\+λ​K\)​α\(n\)=RHSn​\(α\(0\),…,α\(n−1\)\),\\displaystyle\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(n\)\}=\\text\{RHS\}\_\{n\}\\left\(\\alpha^\{\(0\)\},\\ldots,\\alpha^\{\(n\-1\)\}\\right\)\\,,where thennth order right hand side is a function of various combinations of the previous coefficients\.

RHSn=−3​∑j,kK3​\[i,j,k\]×\(αj\(0\)​αk\(n−1\)\+αj\(1\)​αk\(n−2\)\+⋯\+αj\(n−1\)​αk\(0\)\)\\text\{RHS\}\_\{n\}=\-3\\sum\_\{j,k\}K\_\{3\}\[i,j,k\]\\times\\left\(\\alpha\_\{j\}^\{\(0\)\}\\alpha\_\{k\}^\{\(n\-1\)\}\+\\alpha\_\{j\}^\{\(1\)\}\\alpha\_\{k\}^\{\(n\-2\)\}\+\\dots\+\\alpha\_\{j\}^\{\(n\-1\)\}\\alpha\_\{k\}^\{\(0\)\}\\right\)\(32\)and the three\-point kernel tensor is:

K3​\(xi,xj,xk\)=∫K​\(x,xi\)​K​\(x,xj\)​K​\(x,xk\)​f​\(x\|Z\)​𝑑x,K\_\{3\}\(x\_\{i\},x\_\{j\},x\_\{k\}\)=\\int K\(x,x\_\{i\}\)K\(x,x\_\{j\}\)K\(x,x\_\{k\}\)f\(x\|Z\)dx\\,,\(33\)where we useK​\[i,j\]K\[i,j\]to shorthand entries of the Gram matrixK​\(xi,xj\)K\(x\_\{i\},x\_\{j\}\)\.

In order to see the result spectrally, we project bothK~\\widetilde\{K\}andKKonto the eigenbasis ofKK:

M:=ΦT​K~​Φ,K=Φ​Λ​ΦTM:=\\Phi^\{T\}\\widetilde\{K\}\\Phi\\,,\\quad K=\\Phi\\Lambda\\Phi^\{T\}\(34\)whereΛ\\Lambdais the diagonalized Gram matrixKKwith eigenvalues\{λm\}m=1n\\\{\\lambda\_\{m\}\\\}\_\{m=1\}^\{n\}and columns ofΦ\\Phiare the eigenvectorsϕm\\phi\_\{m\}\. The perturbative equations become:

\(M\+λ​Λ\)​\(ΦT​α\(n\)\)=ΦT​RHSn\(M\+\\lambda\\Lambda\)\(\\Phi^\{T\}\\alpha^\{\(n\)\}\)=\\Phi^\{T\}\\text\{RHS\}\_\{n\}\(35\)ifMMis diagonal, we would jump to conclusion and simply write the spectral version of this equation:

ΦT​α\(n\)=∑mϕmT​RHSnMm​m\+λ​λm\\Phi^\{T\}\\alpha^\{\(n\)\}=\\sum\_\{m\}\\frac\{\\phi\_\{m\}^\{T\}\\text\{RHS\}\_\{n\}\}\{M\_\{mm\}\+\\lambda\\lambda\_\{m\}\}\(36\)which in the standard IV kernel ridge regression settingMm​m=⟨T​\(λm​ϕm\),T​\(λm​ϕm\)⟩L2​\(PZ\)=ηm2<λm2M\_\{mm\}=\\langle T\(\\sqrt\{\\lambda\_\{m\}\}\\phi\_\{m\}\),T\(\\sqrt\{\\lambda\_\{m\}\}\\phi\_\{m\}\)\\rangle\_\{L^\{2\}\(P\_\{Z\}\)\}=\\eta\_\{m\}^\{2\}<\\lambda\_\{m\}^\{2\}\(due to the covariance reducing nature ofTT\) andRHS0=h\\text\{RHS\}\_\{0\}=h, which recovers the known expression for the0th order kernel ridge regression result:

g\(0\)​\(x\)=∑m=1Nλm2ηm2\+λ​λm​ϕmT​\[i\]​h​ϕm​\[X\]g^\{\(0\)\}\(x\)=\\sum\_\{m=1\}^\{N\}\\frac\{\\lambda\_\{m\}^\{2\}\}\{\\eta\_\{m\}^\{2\}\+\\lambda\\lambda\_\{m\}\}\\phi^\{T\}\_\{m\}\[i\]h\\,\\phi\_\{m\}\[X\]\(37\)where spectrally only eigenvectors with larger eigenvalues contribute\.

WhenMMis not diagonal in the IV setting, we whiten the equation usingΛ1/2\\Lambda^\{1/2\}as

\(Λ−1/2​M​Λ−1/2⏟A\+λ​I\)​\(Λ1/2​α\(n\)\)=Λ−1/2​ΦT​RHSn⏟RHSn~\(\\underbrace\{\\Lambda^\{\-1/2\}M\\Lambda^\{\-1/2\}\}\_\{A\}\+\\lambda I\)\(\\Lambda^\{1/2\}\\alpha^\{\(n\)\}\)=\\underbrace\{\\Lambda^\{\-1/2\}\\Phi^\{T\}\\text\{RHS\}\_\{n\}\}\_\{\\widetilde\{\\text\{RHS\}\_\{n\}\}\}\(38\)Diagonalizing the positive semi\-definiteA=U​Σ​UTA=U\\Sigma U^\{T\}gives

α\(n\)=Φ​Λ−1/2​U​\(Σ\+λ​I\)−1​UT​RHSn~\\alpha^\{\(n\)\}=\\Phi\\Lambda^\{\-1/2\}U\(\\Sigma\+\\lambda I\)^\{\-1\}U^\{T\}\\,\\widetilde\{\\text\{RHS\}\_\{n\}\}\(39\)which is shrinking along the eigenvector ofAAuniformly, not the original kernel eigenvector\. Note that the conditional expectation is also, in this sense, inducing spectral mixing by rotating different eigenmodes of the kernel matrix underUUdepending on how off\-diagonalMMis\. Although we shall see that in addition to this, the three\-point tensor introduces even further non\-linear spectral mixing\.

In addition, the three\-point tensor decomposes as:

K3​\[i,j,k\]=∑m,n,p=1Nλm​λn​λp​ϕm​\[i\]​ϕm​\[j\]​ϕn​\[j\]​ϕn​\[k\]​ϕp​\[k\]​ϕp​\[i\]\.K\_\{3\}\[i,j,k\]=\\sum\_\{m,n,p=1\}^\{N\}\\lambda\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}\\phi\_\{m\}\[i\]\\phi\_\{m\}\[j\]\\phi\_\{n\}\[j\]\\phi\_\{n\}\[k\]\\phi\_\{p\}\[k\]\\phi\_\{p\}\[i\]\\,\.\(40\)Define the three\-way coupling tensor:

Tm​n​p​\[i,j,k\]=ϕm​\[i\]​ϕm​\[j\]​ϕn​\[j\]​ϕn​\[k\]​ϕp​\[k\]​ϕp​\[i\]\.T\_\{mnp\}\[i,j,k\]=\\phi\_\{m\}\[i\]\\phi\_\{m\}\[j\]\\phi\_\{n\}\[j\]\\phi\_\{n\}\[k\]\\phi\_\{p\}\[k\]\\phi\_\{p\}\[i\]\\,\.\(41\)
This createscross\-eigenspace interactions:

K3​\[i,j,k\]=∑m,n,pλm​λn​λp​Tm​n​p​\[i,j,k\]\.K\_\{3\}\[i,j,k\]=\\sum\_\{m,n,p\}\\lambda\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}T\_\{mnp\}\[i,j,k\]\\,\.\(42\)The coupling strength between eigenspacesmm,nn,ppis weighted byλm​λn​λp\\lambda\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}\. For the first\-order correction:

RHS1​\[i\]\\displaystyle\\text\{RHS\}\_\{1\}\[i\]=−∑j,kK3​\[i,j,k\]​αj\(0\)​αk\(0\)\\displaystyle=\-\\sum\_\{j,k\}K\_\{3\}\[i,j,k\]\\alpha^\{\(0\)\}\_\{j\}\\alpha^\{\(0\)\}\_\{k\}\(43\)=−∑m,n,pλm​λn​λp​∑j,kTm​n​p​\[i,j,k\]​αj\(0\)​αk\(0\)\.\\displaystyle=\-\\sum\_\{m,n,p\}\\lambda\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}\\sum\_\{j,k\}T\_\{mnp\}\[i,j,k\]\\alpha^\{\(0\)\}\_\{j\}\\alpha^\{\(0\)\}\_\{k\}\\,\.
Sinceαj\(0\)=∑qαq\(0\)​ϕq​\[j\]\\alpha^\{\(0\)\}\_\{j\}=\\sum\_\{q\}\\alpha^\{\(0\)\}\_\{q\}\\phi\_\{q\}\[j\], this becomes:

RHS1​\[i\]=−∑m,n,p,q,rλm​λn​λp​αq\(0\)​αr\(0\)​∑j,kTm​n​p​\[i,j,k\]​ϕq​\[j\]​ϕr​\[k\]⏟Im​n​p​q​r\.\\text\{RHS\}\_\{1\}\[i\]=\-\\sum\_\{m,n,p,q,r\}\\lambda\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}\\alpha^\{\(0\)\}\_\{q\}\\alpha^\{\(0\)\}\_\{r\}\\underbrace\{\\sum\_\{j,k\}T\_\{mnp\}\[i,j,k\]\\phi\_\{q\}\[j\]\\phi\_\{r\}\[k\]\}\_\{I\_\{mnpqr\}\}\\,\.\(44\)The key coupling integral is:

Im​n​p​q​r\\displaystyle I\_\{mnpqr\}=∑j,kTm​n​p​\[i,j,k\]​ϕq​\[j\]​ϕr​\[k\]\\displaystyle=\\sum\_\{j,k\}T\_\{mnp\}\[i,j,k\]\\phi\_\{q\}\[j\]\\phi\_\{r\}\[k\]=∑j,kϕm​\[i\]​ϕm​\[j\]​ϕn​\[j\]​ϕn​\[k\]​ϕp​\[k\]​ϕp​\[i\]​ϕq​\[j\]​ϕr​\[k\]\\displaystyle=\\sum\_\{j,k\}\\phi\_\{m\}\[i\]\\phi\_\{m\}\[j\]\\phi\_\{n\}\[j\]\\phi\_\{n\}\[k\]\\phi\_\{p\}\[k\]\\phi\_\{p\}\[i\]\\phi\_\{q\}\[j\]\\phi\_\{r\}\[k\]=ϕm\[i\]ϕp\[i\]\(∑jϕm\[j\]ϕn\[j\]ϕq\[j\]\)×\\displaystyle=\\phi\_\{m\}\[i\]\\phi\_\{p\}\[i\]\\left\(\\sum\_\{j\}\\phi\_\{m\}\[j\]\\phi\_\{n\}\[j\]\\phi\_\{q\}\[j\]\\right\)\\times\(∑kϕn​\[k\]​ϕp​\[k\]​ϕr​\[k\]\)\.\\displaystyle\\left\(\\sum\_\{k\}\\phi\_\{n\}\[k\]\\phi\_\{p\}\[k\]\\phi\_\{r\}\[k\]\\right\)\\,\.We shall denote the triple products∑jϕa​\[j\]​ϕb​\[j\]​ϕc​\[j\]\\sum\_\{j\}\\phi\_\{a\}\[j\]\\phi\_\{b\}\[j\]\\phi\_\{c\}\[j\]evaluated at the sameXjX\_\{j\}asSa​b​cS\_\{abc\}\.Sa​b​cS\_\{abc\}depends on the specific kernel one uses, it symbolizes how much higher order interactions a kernel function can have\. Hence we have

RHS1​\[i\]=−∑m,n,p,q,rλm​λn​λp​αq\(0\)​αr\(0\)​ϕm​\[i\]​ϕp​\[i\]​Sm​n​q​Sn​p​r\.\\text\{RHS\}\_\{1\}\[i\]=\-\\sum\_\{m,n,p,q,r\}\\lambda\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}\\alpha^\{\(0\)\}\_\{q\}\\alpha^\{\(0\)\}\_\{r\}\\phi\_\{m\}\[i\]\\phi\_\{p\}\[i\]S\_\{mnq\}S\_\{npr\}\\,\.\(45\)Using the generic form of the spectral solution \([39](https://arxiv.org/html/2606.00322#A1.E39)\), we can substitute in theRHSn\\text\{RHS\}\_\{n\}to obtain the correspondingα\(n\)\\alpha^\{\(n\)\}\.

g\(0\)​\(x\)=∑iK​\(x,xi\)​Φ​Λ−1/2​U​\(Σ\+λ​I\)−1​UT​Λ−1/2​ΦT​h\\displaystyle g^\{\(0\)\}\(x\)=\\sum\_\{i\}K\(x,x\_\{i\}\)\\Phi\\Lambda^\{\-1/2\}U\(\\Sigma\+\\lambda I\)^\{\-1\}U^\{T\}\\Lambda^\{\-1/2\}\\Phi^\{T\}h\(46\)g\(1\)​\(x\)=∑iK​\(x,xi\)​Φ​Λ−1/2​U​\(Σ\+λ​I\)−1​UT​Λ−1/2​ΦT\\displaystyle g^\{\(1\)\}\(x\)=\\sum\_\{i\}K\(x,x\_\{i\}\)\\Phi\\Lambda^\{\-1/2\}U\(\\Sigma\+\\lambda I\)^\{\-1\}U^\{T\}\\Lambda^\{\-1/2\}\\Phi^\{T\}∑m,n,p,q,rλm​λn​λp​αq\(0\)​αr\(0\)​ϕm​\[i\]​ϕp​\[i\]​Sm​n​q​Sn​p​r\.\\displaystyle\\sum\_\{m,n,p,q,r\}\\lambda\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}\\alpha^\{\(0\)\}\_\{q\}\\alpha^\{\(0\)\}\_\{r\}\\phi\_\{m\}\[i\]\\phi\_\{p\}\[i\]S\_\{mnq\}S\_\{npr\}\\,\.Although we can already see that the first order contribution tog​\(x\)g\(x\)contains nonlinear mixing between the eigendirections throughλm​λn​λp\\lambda\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}, it is clearest in the simplified setting of kernel ridge regression\.

In order to see what is happening spectrally clearly, we work under the conditions thatMMis diagonal or almost diagonal in the Mercer basis, withMm​m=⟨T​\(λm​ϕm\),T​\(λm​ϕm\)⟩L2​\(PZ\)=ηm2<λm2M\_\{mm\}=\\langle T\(\\sqrt\{\\lambda\_\{m\}\}\\phi\_\{m\}\),T\(\\sqrt\{\\lambda\_\{m\}\}\\phi\_\{m\}\)\\rangle\_\{L^\{2\}\(P\_\{Z\}\)\}=\\eta\_\{m\}^\{2\}<\\lambda\_\{m\}^\{2\}, since the conditional operatorTTis contraction operator reducing the variance\. Summarizing the zeroth order and first order term:

g\(0\)​\(x\)=∑i∑m=1Nλm2ηm2\+λ​λm​ϕmT​\[i\]​h​ϕm​\[X\];\\displaystyle g^\{\(0\)\}\(x\)=\\sum\_\{i\}\\sum\_\{m=1\}^\{N\}\\frac\{\\lambda\_\{m\}^\{2\}\}\{\\eta\_\{m\}^\{2\}\+\\lambda\\lambda\_\{m\}\}\\phi^\{T\}\_\{m\}\[i\]h\\,\\phi\_\{m\}\[X\]\\,;\(47\)g\(1\)​\(x\)=∑i∑m,n,p,q,r\\displaystyle g^\{\(1\)\}\(x\)=\\sum\_\{i\}\\sum\_\{m,n,p,q,r\}λm3​λn​λpηm2\+λ​λm​αq\(0\)​αr\(0\)​ϕm​\[i\]​ϕp​\[i\]​ϕm​\[i\]​ϕm​\[X\]​Sm​n​q​Sn​p​r\\displaystyle\\frac\{\\lambda\_\{m\}^\{3\}\\lambda\_\{n\}\\lambda\_\{p\}\}\{\\eta\_\{m\}^\{2\}\+\\lambda\\lambda\_\{m\}\}\\alpha^\{\(0\)\}\_\{q\}\\alpha^\{\(0\)\}\_\{r\}\\phi\_\{m\}\[i\]\\phi\_\{p\}\[i\]\\phi\_\{m\}\[i\]\\phi\_\{m\}\[X\]S\_\{mnq\}S\_\{npr\}=\\displaystyle=∑i∑m,n,p,q,rλm3​λn​λpηm2\+λ​λmλq​ϕqT​hηq2\+λ​λqλr​ϕrT​hηr2\+λ​λr×\\displaystyle\\sum\_\{i\}\\sum\_\{m,n,p,q,r\}\\frac\{\\lambda^\{3\}\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}\}\{\\eta\_\{m\}^\{2\}\+\\lambda\\lambda\_\{m\}\}\\frac\{\\lambda\_\{q\}\\phi^\{T\}\_\{q\}h\}\{\\eta\_\{q\}^\{2\}\+\\lambda\\lambda\_\{q\}\}\\frac\{\\lambda\_\{r\}\\phi^\{T\}\_\{r\}h\}\{\\eta\_\{r\}^\{2\}\+\\lambda\\lambda\_\{r\}\}\\timesϕm​\[i\]​ϕp​\[i\]​ϕm​\[i\]​ϕm​\[X\]​Sm​n​q​Sn​p​r\.\\displaystyle\\phi\_\{m\}\[i\]\\phi\_\{p\}\[i\]\\phi\_\{m\}\[i\]\\phi\_\{m\}\[X\]S\_\{mnq\}S\_\{npr\}\\,\.
Here we introduce the quantity condition numberκ\\kappaas the ratio between the largest and smallest eigenvalues of the kernel matrix\. Since eigenvalues are labeled from11toNNin a descending order for positive definite matrices, we have the following definition:

κ=λ1/λN\.\\kappa=\\lambda\_\{1\}/\\lambda\_\{N\}\\,\.\(48\)
#### Scale separation\!

Comparing the zeroth and first order solutions, we see that condition numberκ\\kappaplays an important role in determining when the first order contribution matters\. For simplicity in the arguments, we assumeηm∼λm\\eta\_\{m\}\\sim\\lambda\_\{m\}\. Then

g\(0\)​\(x\)=∑i∑m=1Nλmλm\+λ​ϕmT​\[i\]​h​ϕm​\[X\];\\displaystyle g^\{\(0\)\}\(x\)=\\sum\_\{i\}\\sum\_\{m=1\}^\{N\}\\frac\{\\lambda\_\{m\}\}\{\\lambda\_\{m\}\+\\lambda\}\\phi^\{T\}\_\{m\}\[i\]h\\,\\phi\_\{m\}\[X\]\\,;\(49\)g\(1\)\(x\)=∑i∑m,n,p,q,rλm2​λn​λpλm\+λϕqT​hλq\+λϕrT​hλr\+λ×\\displaystyle g^\{\(1\)\}\(x\)=\\sum\_\{i\}\\sum\_\{m,n,p,q,r\}\\frac\{\\lambda^\{2\}\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}\}\{\\lambda\_\{m\}\+\\lambda\}\\frac\{\\phi^\{T\}\_\{q\}h\}\{\\lambda\_\{q\}\+\\lambda\}\\frac\{\\phi^\{T\}\_\{r\}h\}\{\\lambda\_\{r\}\+\\lambda\}\\timesϕm​\[i\]​ϕp​\[i\]​ϕm​\[i\]​ϕm​\[X\]​Sm​n​q​Sn​p​r\.\\displaystyle\\phi\_\{m\}\[i\]\\phi\_\{p\}\[i\]\\phi\_\{m\}\[i\]\\phi\_\{m\}\[X\]S\_\{mnq\}S\_\{npr\}\\,\.Whenκ\\kappais well\-conditioned,λm∼λp∼λq\\lambda\_\{m\}\\sim\\lambda\_\{p\}\\sim\\lambda\_\{q\}are of similar order of magnitude for allp,q,mp,q,m\. This simply putsα\(1\)\\alpha^\{\(1\)\}in the same footing asα\(0\)\\alpha^\{\(0\)\}, roughly counting the number of eigenvalues bigger than the regularization parameterλ\\lambda\. All eigenvectors contribute equally in this regime\. Whenκ≫1\\kappa\\gg 1,α\(0\)\\alpha^\{\(0\)\}only receives contributions from the first few eigenvectors, losing its representation power immensely\. However,α\(1\)\\alpha^\{\(1\)\}is still corrected by further eigenvalues since the coefficientλm​λp​λq\\lambda\_\{m\}\\lambda\_\{p\}\\lambda\_\{q\}can still remain large when we selectλm\\lambda\_\{m\}for largemmas long as we choosep,qp,qto be close to11\. This introduces additional representation power by mixing in the eigenvectors with small eigenvalues\. Ifκ\\kappais well\-conditioned, all eigenfunctions are represented similarly whileκ≫1\\kappa\\gg 1indicates only the first few eigenfunctions contribute\. Kernel functions used in the NPIV scenario has a incredibly large condition number due to the ill\-defined convolution integral inversion\. We observed that in regular kernel ridge regression, where condition number is small, the perturbative method \(addition of the third moment\) does not help the zeroth order solution\. It is only effective whenκ≫1\\kappa\\gg 1, we shall analyze this phenomenon from the spectral perspective\.

In the case when the condition number is low, all coupling weightsλm2​λn​λp≈λtypical4\\lambda^\{2\}\_\{m\}\\lambda\_\{n\}\\lambda\_\{p\}\\approx\\lambda\_\{\\text\{typical\}\}^\{4\}are similar\.α\(0\)\\alpha^\{\(0\)\}andα\(1\)\\alpha^\{\(1\)\}both receive ample equal contributions from all the eigenfunctions, soα\(1\)\\alpha^\{\(1\)\}only adds additional noise to the regressed function\.

However whenκ≫1\\kappa\\gg 1,α\(0\)​\(xi\)​K​\(x,xi\)\\alpha^\{\(0\)\}\(x\_\{i\}\)K\(x,x\_\{i\}\)only receives contributions from the eigenfunctions which have large eigenvalues compared toλ\\lambda\. This is where the triple moment helps, let us analyze case by case:

- •High\-to\-high coupling:m,n,pm,n,pall small⇒\\Rightarrowweightλm2​λn​λp\\lambda\_\{m\}^\{2\}\\lambda\_\{n\}\\lambda\_\{p\}is large
- •High\-to\-low coupling: Mixed indices create medium\-strength interactions
- •Low\-to\-low coupling:m,n,pm,n,pall large⇒\\Rightarrowweightλm2​λn​λp\\lambda\_\{m\}^\{2\}\\lambda\_\{n\}\\lambda\_\{p\}is small

We see thatα\(1\)\\alpha^\{\(1\)\}recovers contributions from eigenfunctions with larger eigenvalues as well\. This is the manifestation of improving the representation power of the kernel when condition number is high\.

### A\.3The failure of rotationally invariant kernels

Another issue we can address with the spectral perspective is why rotationally invariant kernels \(although having high condition numbers\) fit completely trivial lines in high dimensions\. The usual kernel ridge regression case is discussed extensively in\(Donhauseret al\.,[2021](https://arxiv.org/html/2606.00322#bib.bib15)\), we empirically observe that this is also the case even after adding the third moment term with RBF kernels\. The key is the triple eigenfunction integralSa​b​cS\_\{abc\}:

Sa​b​c=∑iϕa​\[i\]​ϕb​\[i\]​ϕc​\[i\]\.S\_\{abc\}=\\sum\_\{i\}\\phi\_\{a\}\[i\]\\phi\_\{b\}\[i\]\\phi\_\{c\}\[i\]\\,\.\(50\)For any rotationally invariant kernelK​\(x,y\)=k​\(‖x−y‖\)K\(x,y\)=k\(\\\|x\-y\\\|\): The eigenfunctions can be decomposed into the radial part and spherical harmonics:

ϕl,m​\(x\)=Rl​\(‖x‖\)​Yl,m​\(x/‖x‖\),\\phi\_\{l,m\}\(x\)=R\_\{l\}\(\\\|x\\\|\)Y\_\{l,m\}\(x/\\\|x\\\|\)\\,,\(51\)wherel∈ℕ0l\\in\\mathbb\{N\}\_\{0\}is the radial scaling homogeneity,m=1,…,dlm=1,\\ldots,d\_\{l\}wheredl=dim\(Harml​\(ℝn\)\)d\_\{l\}=\\dim\(\\text\{Harm\}\_\{l\}\(\\mathbb\{R\}^\{n\}\)\)andYl,mY\_\{l,m\}are orthonormal spherical harmonics of degreell\. The triple product has closed form

∑iϕa​\(xi\)​ϕb​\(xi\)​ϕc​\(xi\)=∑iRa​\(ri\)​Rb​\(ri\)​Rc​\(ri\)×Ca​bc​\(ω^i\),\\sum\_\{i\}\\phi\_\{a\}\(x\_\{i\}\)\\phi\_\{b\}\(x\_\{i\}\)\\phi\_\{c\}\(x\_\{i\}\)=\\sum\_\{i\}R\_\{a\}\(r\_\{i\}\)R\_\{b\}\(r\_\{i\}\)R\_\{c\}\(r\_\{i\}\)\\times C\_\{ab\}^\{c\}\(\\hat\{\\omega\}\_\{i\}\)\\,,\(52\)whereri=‖xi‖r\_\{i\}=\\\|x\_\{i\}\\\|,ω^i=xi/‖xi‖\\hat\{\\omega\}\_\{i\}=x\_\{i\}/\\\|x\_\{i\}\\\|andCa​bc​\(ω^\)=Ya​\(ω^\)​Yb​\(ω^\)​Yc​\(ω^\)C\_\{ab\}^\{c\}\(\\hat\{\\omega\}\)=Y\_\{a\}\(\\hat\{\\omega\}\)Y\_\{b\}\(\\hat\{\\omega\}\)Y\_\{c\}\(\\hat\{\\omega\}\)is the Clebsch\-Gordan coefficients\. The Clebsch\-Gordan coefficients satisfy strict selection rules:⟨l1​m1,l2​m2\|l3​m3⟩=0\\langle l\_\{1\}m\_\{1\},l\_\{2\}m\_\{2\}\|l\_\{3\}m\_\{3\}\\rangle=0unless\|l1−l2\|≤l3≤l1\+l2\|l\_\{1\}\-l\_\{2\}\|\\leq l\_\{3\}\\leq l\_\{1\}\+l\_\{2\}andl1\+l2\+l3​is evenl\_\{1\}\+l\_\{2\}\+l\_\{3\}\\text\{ is even\}andm1\+m2\+m3=0m\_\{1\}\+m\_\{2\}\+m\_\{3\}=0\.

These would set most of the combinations to0, essentially any combination that are not in the representation of the orthogonal groupO​\(d\)O\(d\)\. This explains the failure of RBF kernels even after adding the triple moment term\.

On the other hand for the polynomial kernelK​\(x,y\)=\(x⋅y\+c\)dK\(x,y\)=\(x\\cdot y\+c\)^\{d\}onℝn\\mathbb\{R\}^\{n\}, eigenfunctions are multinomial in the coordinates:

φa​\(x\)=λa​∏j=1nxjaj×\(normalization factor\),\\varphi\_\{a\}\(x\)=\\sqrt\{\\lambda\_\{a\}\}\\prod\_\{j=1\}^\{n\}x\_\{j\}^\{a\_\{j\}\}\\times\\text\{\(normalization factor\)\}\\,,\(53\)wherea=\(a1,…,an\)a=\(a\_\{1\},\\ldots,a\_\{n\}\)is a multi\-index with\|a\|=a1\+⋯\+an≤d\|a\|=a\_\{1\}\+\\cdots\+a\_\{n\}\\leq d\. The triple product is simply

∑iφa​\(xi\)​φb​\(xi\)​φc​\(xi\)\\displaystyle\\sum\_\{i\}\\varphi\_\{a\}\(x\_\{i\}\)\\varphi\_\{b\}\(x\_\{i\}\)\\varphi\_\{c\}\(x\_\{i\}\)\(54\)=\(normalization\)×∑i∏j=1nxi,jaj\+bj\+cj\\displaystyle=\\text\{\(normalization\)\}\\times\\sum\_\{i\}\\prod\_\{j=1\}^\{n\}x\_\{i,j\}^\{a\_\{j\}\+b\_\{j\}\+c\_\{j\}\}\(55\)=\(normalization\)×∑i∏j=1nxi,jdj,\\displaystyle=\\text\{\(normalization\)\}\\times\\sum\_\{i\}\\prod\_\{j=1\}^\{n\}x\_\{i,j\}^\{d\_\{j\}\}\\,,\(56\)wheredj=aj\+bj\+cjd\_\{j\}=a\_\{j\}\+b\_\{j\}\+c\_\{j\}\. This is just a power sum in each coordinate which does not impose any selection rule\.

### A\.4The need for renormalization

In this subsection, we perform some simple counting using the usual kernel ridge regression case to demonstrate the inevitable divergence of the naive perturbative series\. In the simplest setting whereK~\\widetilde\{K\}is diagonalizable in the kernel eigenbasis \([49](https://arxiv.org/html/2606.00322#A1.E49)\), from the recursive structure, we can derive that:

αm\(n\)∼λm3​n\(ηm\+λ\)2​n\+1×geometric factors\.\\alpha^\{\(n\)\}\_\{m\}\\sim\\frac\{\\lambda\_\{m\}^\{3n\}\}\{\\left\(\\eta\_\{m\}\+\\lambda\\right\)^\{2n\+1\}\}\\times\\text\{geometric factors\}\\,\.\(57\)
This leads to the crucial scaling:

Order 0:\|αm\(0\)\|\\displaystyle\\text\{Order 0:\}\\quad\|\\alpha^\{\(0\)\}\_\{m\}\|∼1ηm\+λ\\displaystyle\\sim\\frac\{1\}\{\\eta\_\{m\}\+\\lambda\}\(58\)Order 1:\|αm\(1\)\|\\displaystyle\\text\{Order 1:\}\\quad\|\\alpha^\{\(1\)\}\_\{m\}\|∼λm3\(ηm\+λ\)3\\displaystyle\\sim\\frac\{\\lambda\_\{m\}^\{3\}\}\{\(\\eta\_\{m\}\+\\lambda\)^\{3\}\}\(59\)Order 2:\|αm\(2\)\|\\displaystyle\\text\{Order 2:\}\\quad\|\\alpha^\{\(2\)\}\_\{m\}\|∼λm6\(ηm\+λ\)5\\displaystyle\\sim\\frac\{\\lambda\_\{m\}^\{6\}\}\{\(\\eta\_\{m\}\+\\lambda\)^\{5\}\}\(60\)Order n:\|αm\(n\)\|\\displaystyle\\text\{Order n:\}\\quad\|\\alpha^\{\(n\)\}\_\{m\}\|∼λm3​n\(ηm\+λ\)2​n\+1\.\\displaystyle\\sim\\frac\{\\lambda\_\{m\}^\{3n\}\}\{\(\\eta\_\{m\}\+\\lambda\)^\{2n\+1\}\}\\,\.\(61\)We shall assumeλm∼ηm\\lambda\_\{m\}\\sim\\eta\_\{m\}for simplicity\.

#### Well\-Conditioned Case \(κ\\kappanot divergent\)

When all eigenvalues are similar:λm≈λtypical\\lambda\_\{m\}\\approx\\lambda\_\{\\text\{typical\}\}for allmm\. Combining with the couplingγ\\gamma:

γn​\|α\(n\)\|γn−1​\|α\(n−1\)\|∼γ​λtypical3\(λtypical\+λ\)2​\{λtypical≫λ:γ<1λtypical∼λ:γ<4λλtypical≪λ:always convergent\\frac\{\\gamma^\{n\}\|\\alpha^\{\(n\)\}\|\}\{\\gamma^\{n\-1\}\|\\alpha^\{\(n\-1\)\}\|\}\\sim\\gamma\\,\\frac\{\\lambda\_\{\\text\{typical\}\}^\{3\}\}\{\(\\lambda\_\{\\text\{typical\}\}\+\\lambda\)^\{2\}\}\\begin\{cases\}\\lambda\_\{\\text\{typical\}\}\\gg\\lambda:\\gamma<1\\\\ \\lambda\_\{\\text\{typical\}\}\\sim\\lambda:\\gamma<\\frac\{4\}\{\\lambda\}\\\\ \\lambda\_\{\\text\{typical\}\}\\ll\\lambda:\\text\{always convergent\}\\end\{cases\}\(62\)We listed the convergence requirement forγ\\gamma, as long as we make appropriate choices ofλ\\lambdacompatible with the ridge regularization parameter, the series will be convergent\. Hence we would not need any renormalization anyways since all contributions vanish at high order\. The perturbative series terminates naturally, but provides no improvement over zeroth order\.

#### Ill\-Conditioned Case \(κ≫100\\kappa\\gg 100\)

Eigenvalues have extreme hierarchy: few largeλ1≫λ2≫⋯≫λ≫⋯≫λN\\lambda\_\{1\}\\gg\\lambda\_\{2\}\\gg\\dots\\gg\\lambda\\gg\\dots\\gg\\lambda\_\{N\}\. The convergence condition for the series becomes:

γn​\|α\(n\)\|γn−1​\|α\(n−1\)\|∼γ​λm​λp​λq\(λp\+λ\)2,\\frac\{\\gamma^\{n\}\|\\alpha^\{\(n\)\}\|\}\{\\gamma^\{n\-1\}\|\\alpha^\{\(n\-1\)\}\|\}\\sim\\gamma\\,\\frac\{\\lambda\_\{m\}\\lambda\_\{p\}\\lambda\_\{q\}\}\{\(\\lambda\_\{p\}\+\\lambda\)^\{2\}\}\\,,\(63\)which can go to infinity or0depending on the eigenvalues chosen\. The rescaling renormalization we used

γn=γn×min​\(\|α\(0\)\|\|α\(n\)\|,1\),\\gamma\_\{n\}=\\gamma^\{n\}\\times\\text\{min\}\\left\(\\frac\{\|\\alpha^\{\(0\)\}\|\}\{\|\\alpha^\{\(n\)\}\|\},1\\right\)\\,,\(64\)automatically implements spectral flow:

- •Stable eigenspaces: Contribute little to‖α\(n\)‖\\\|\\alpha^\{\(n\)\}\\\|⇒\\Rightarrowremain at full strength
- •Unstable eigenspaces: Dominate‖α\(n\)‖\\\|\\alpha^\{\(n\)\}\\\|⇒\\Rightarrowget suppressed by smallℛn\\mathcal\{R\}\_\{n\}
- •Critical eigenspaces: Drive the renormalization flow

## Appendix BBackground on resurgence

Using the mathematical tool of resurgence and Borel resummation\(Ecalle,[1981](https://arxiv.org/html/2606.00322#bib.bib26)\), we come up with a further enhancement of the renormalized perturbative series\. Note that the series \(even after renormalization\) is still an asymptotic series\. In standard QFT context, this is due to non\-perturbative effects such as field configurations that are trivial\(Berry and Howls,[1991](https://arxiv.org/html/2606.00322#bib.bib41); Mariño,[2014](https://arxiv.org/html/2606.00322#bib.bib39); Dorigoni,[2019](https://arxiv.org/html/2606.00322#bib.bib42); Delabaere and Pham,[1999](https://arxiv.org/html/2606.00322#bib.bib40); Costinet al\.,[2023](https://arxiv.org/html/2606.00322#bib.bib43)\), in our case, these are solutions that lie in the null space of the kernel matrix\. They are trivialized in the perturbative series yet contributes to the solution space of the ill\-defined inverse integral equation\.

The prevalence of saddle points in high dimensional optimization problems is the key issue here\. Several facts about them that makes the usual methods fail:

- •they are never correctly captured by finite truncations of any perturbative series\.
- •they are minimizer of the action/loss function that are ”trivial” in the sense that they are not solutions to the perturbative equation of motion truncated at any order\.
- •They live in the kernel of the perturbative operator, which is difficult to compute and usually the cause of asymptotic series in the perturbative solution\.

Luckily, Borel resummation is to rescue, it provides a systematic way of analyzing these saddles and incorporate them in the solution\.

You also run into these issues in neural networks, where they appear as saddles of the loss function, which grow in number as the dimension of the loss landscape increases\. This is usually tackled using optimizers and schedulers to create local perturbations that destabilizes the saddle configuration\.

### B\.1Introduction to resurgence

###### Definition 1\(Asymptotic Series\)\.

A formal power series∑n=0∞an​zn\\sum\_\{n=0\}^\{\\infty\}a\_\{n\}z^\{n\}is called an asymptotic expansion off​\(z\)f\(z\)asz→0z\\to 0in a sectorSSif for everyN≥0N\\geq 0:

f​\(z\)−∑n=0N−1an​zn=O​\(zN\)as​z→0​in​S\.f\(z\)\-\\sum\_\{n=0\}^\{N\-1\}a\_\{n\}z^\{n\}=O\(z^\{N\}\)\\quad\\text\{as \}z\\to 0\\text\{ in \}S\\,\.\(65\)We writef​\(z\)∼∑n=0∞an​znf\(z\)\\sim\\sum\_\{n=0\}^\{\\infty\}a\_\{n\}z^\{n\}\.

###### Example B\.1\(The Euler Integral\)\.

Consider the integral representation:

E​\(z\)=∫0∞e−t1\+z​t​𝑑t,\|arg⁡\(z\)\|<π2\.E\(z\)=\\int\_\{0\}^\{\\infty\}\\frac\{e^\{\-t\}\}\{1\+zt\}dt,\\quad\|\\arg\(z\)\|<\\frac\{\\pi\}\{2\}\\,\.\(66\)
Integration by parts yields the asymptotic expansion:

E​\(z\)∼∑n=0∞\(−1\)n​n\!​znas​z→0\.E\(z\)\\sim\\sum\_\{n=0\}^\{\\infty\}\(\-1\)^\{n\}n\!z^\{n\}\\quad\\text\{as \}z\\to 0\\,\.\(67\)
This series diverges for allz≠0z\\neq 0sincean=\(−1\)n​n\!a\_\{n\}=\(\-1\)^\{n\}n\!grows factorially, yet it provides the unique asymptotic expansion ofE​\(z\)E\(z\)\.

For factorially divergent series∑an​zn\\sum a\_\{n\}z^\{n\}with\|an\|∼n\!/rn\|a\_\{n\}\|\\sim n\!/r^\{n\}, the optimal truncation occurs at:

N∗≈r\|z\|N^\{\*\}\\approx\\frac\{r\}\{\|z\|\}\(68\)
The truncation error is exponentially small:

\|f​\(z\)−∑n=0N∗an​zn\|∼e−r/\|z\|\.\\left\|f\(z\)\-\\sum\_\{n=0\}^\{N^\{\*\}\}a\_\{n\}z^\{n\}\\right\|\\sim e^\{\-r/\|z\|\}\\,\.\(69\)
This exponential accuracy is the hallmark of asymptotic series, but it also reveals their fundamental limitation: exponentially small terms are completely invisible\.

###### Definition 2\(Borel Transform\)\.

For a formal power seriesϕ​\(z\)=∑n=0∞an​zn\\phi\(z\)=\\sum\_\{n=0\}^\{\\infty\}a\_\{n\}z^\{n\}, the Borel transform is:

ℬ​\[ϕ\]​\(ζ\)=∑n=0∞ann\!​ζn\.\\mathcal\{B\}\[\\phi\]\(\\zeta\)=\\sum\_\{n=0\}^\{\\infty\}\\frac\{a\_\{n\}\}\{n\!\}\\zeta^\{n\}\\,\.\(70\)

Whenℬ​\[ϕ\]​\(ζ\)\\mathcal\{B\}\[\\phi\]\(\\zeta\)can be analytically continued to the positive real axis, we can define:

###### Definition 3\(Borel Sum\)\.

The Borel sum ofϕ​\(z\)\\phi\(z\)is:

𝒮​\[ϕ\]​\(z\)=1z​∫0∞e−ζ/z​ℬ​\[ϕ\]​\(ζ\)​𝑑ζ,\\mathcal\{S\}\[\\phi\]\(z\)=\\frac\{1\}\{z\}\\int\_\{0\}^\{\\infty\}e^\{\-\\zeta/z\}\\mathcal\{B\}\[\\phi\]\(\\zeta\)d\\zeta\\,,\(71\)provided the integral converges\.

###### Example B\.2\(Borel Resummation of Euler Integral\)\.

ForE​\(z\)∼∑n=0∞\(−1\)n​n\!​znE\(z\)\\sim\\sum\_\{n=0\}^\{\\infty\}\(\-1\)^\{n\}n\!z^\{n\}, the Borel transform is:

ℬ​\[E\]​\(ζ\)=∑n=0∞\(−1\)n​n\!n\!​ζn=∑n=0∞\(−ζ\)n=11\+ζ\.\\mathcal\{B\}\[E\]\(\\zeta\)=\\sum\_\{n=0\}^\{\\infty\}\\frac\{\(\-1\)^\{n\}n\!\}\{n\!\}\\zeta^\{n\}=\\sum\_\{n=0\}^\{\\infty\}\(\-\\zeta\)^\{n\}=\\frac\{1\}\{1\+\\zeta\}\\,\.\(72\)
The Borel sum gives:

𝒮​\[E\]​\(z\)=1z​∫0∞e−ζ/z1\+ζ​𝑑ζ=E​\(z\)\.\\mathcal\{S\}\[E\]\(z\)=\\frac\{1\}\{z\}\\int\_\{0\}^\{\\infty\}\\frac\{e^\{\-\\zeta/z\}\}\{1\+\\zeta\}d\\zeta=E\(z\)\\,\.\(73\)
Thus, Borel resummation exactly recovers the original function\!

With a subtlety that we have secretly changed the position of an integral and a sum, which might not be convergent\. Forx<<1x<<1:

###### Example B\.3\.

∑n=0∞n\!​xn=∑n=0∞xn​∫0∞𝑑t​e−t​tn=∫0∞e−t​11−x​t\.\\sum\_\{n=0\}^\{\\infty\}n\!x^\{n\}=\\sum\_\{n=0\}^\{\\infty\}x^\{n\}\\int\_\{0\}^\{\\infty\}dt\\,e^\{\-t\}\\,t^\{n\}=\\int\_\{0\}^\{\\infty\}e^\{\-t\}\\,\\frac\{1\}\{1\-xt\}\\,\.\(74\)

This integral is singular on the positive real axis att=1xt=\\frac\{1\}\{x\}\. This forces us to analytically continue and change the contour\. There are naturally two choices to do that, which results in an ambiguity in the result\. They differ by the residue of the integrand at the singularity\. Surprisingly, this residue called Stokes constant also corresponds to the saddle point approximation of the integral representation of the series\. This is what we rely on to do the analysis on our series phrased as a minimization problem in high dimensions with multiple saddle points\.

###### Example B\.4\.

We can also show that if an integral representation exists for the equation we are trying to solve, the saddle point approximation of the integral is exactly what the residue at the poles on Borel plane of the Laplace transform is computing\. For the Airy equation:

d2​yd​x2−x​y=0,\\frac\{d^\{2\}y\}\{dx^\{2\}\}\-xy=0\\,,\(75\)perturbatively, it has the following solution:

y​\(x\)=∑n=0∞\(−34\)n​Γ​\(n\+5/6\)​Γ​\(n\+1/6\)n\!​\(x−3/2\)n\.y\(x\)=\\sum\_\{n=0\}^\{\\infty\}\\left\(\-\\frac\{3\}\{4\}\\right\)^\{n\}\\,\\frac\{\\Gamma\(n\+5/6\)\\Gamma\(n\+1/6\)\}\{n\!\}\\,\(x^\{\-3/2\}\)^\{n\}\\,\.\(76\)Its Borel resummation has a singularity att=−23​x3/2t=\-\\frac\{2\}\{3\}x^\{3/2\}\. The residue at the singularity is

Rest=−23​x3/2​𝑩​y​\[t\]=e2/3​x3/2\.\\text\{Res\}\_\{t=\-\\frac\{2\}\{3\}x^\{3/2\}\}\\boldsymbol\{B\}y\[t\]=e^\{2/3x^\{3/2\}\}\\,\.\(77\)Then we do a saddle point approximation for the Airy function’s integral representation:

y​\(x\)=∫0∞et3/3−x​t​𝑑t∼et3/3−x​t\|t=x​2​π​i1=e2/3​x3/2,y\(x\)=\\int\_\{0\}^\{\\infty\}e^\{t^\{3\}/3\-xt\}dt\\sim\\left\.e^\{t^\{3\}/3\-xt\}\\right\|\_\{t=\\sqrt\{x\}\}\\sqrt\{\\frac\{2\\pi i\}\{1\}\}=e^\{2/3x^\{3/2\}\}\\,,\(78\)which precisely matches the residue at the singularity on the Borel plane\. Usually it is difficult to solve for the saddle point of the minimizer in high dimensions, we can use the Borel analysis to bypass that\.

What was mentioned here is a generic pattern that we can formulate and prove as a theorem\(Bhattacharyaet al\.,[2024](https://arxiv.org/html/2606.00322#bib.bib38)\):

###### Theorem 1\.

The critical point contributions to the integral representation of an asymptotic series are in one\-to\-one correspondence with the poles on Borel plane\.

proof:A simple proof for this involves a trick in measure theory named the co\-area formula, which is used to reduce the dimension of an integral onto some lower dimensional domain level set\.

∫Ω∈ℝng​\(x\)​\|Jk​u​\(x\)\|​𝑑x=∫ℝk\(∫u−1​\(t\)g​\(x\)​𝑑Hn−k​\(x\)\)​𝑑t,\\int\_\{\\Omega\\in\\mathbb\{R\}^\{n\}\}g\(x\)\|J\_\{k\}u\(x\)\|dx=\\int\_\{\\mathbb\{R\}^\{k\}\}\\left\(\\int\_\{u^\{\-1\}\(t\)\}g\(x\)dH\_\{n\-k\}\(x\)\\right\)\\,dt\\,,\(79\)wherek<nk<n,u​\(x\)=tu\(x\)=tlabels then−kn\-kdimensional level set,\|Jk​u​\(x\)\|\|J\_\{k\}u\(x\)\|is itskk\-dim Jacobian andd​Hn−kdH\_\{n\-k\}represents the appropriate Hausdorff measure on the level set\. Multiplying both sides with a delta function to indicate the level setδ​\(t−u​\(x\)\)\\delta\(t\-u\(x\)\)\. Then we have

∫Ωg​\(x\)​δ​\(t−u​\(x\)\)​\|∇u​\(x\)\|​𝑑x=∫t=u​\(x\)g​\(x\)​𝑑H​\(x\)\.\\int\_\{\\Omega\}g\(x\)\\delta\(t\-u\(x\)\)\|\\nabla u\(x\)\|dx=\\int\_\{t=u\(x\)\}g\(x\)\\,dH\(x\)\\,\.\(80\)Now choosingg​\(x\)=1\|∇u​\(x\)\|g\(x\)=\\frac\{1\}\{\|\\nabla u\(x\)\|\}gives

∫Ωδ​\(t−u​\(x\)\)=∫t=u​\(x\)d​H​\(x\)\|∇u​\(x\)\|\.\\int\_\{\\Omega\}\\delta\(t\-u\(x\)\)=\\int\_\{t=u\(x\)\}\\frac\{dH\(x\)\}\{\|\\nabla u\(x\)\|\}\\,\.\(81\)We shall use this formula in the following proof\.

Writing the multivariant function as an asymptotic seriesf​\(g\)=∑n=0∞an​gnf\(g\)=\\sum\_\{n=0\}^\{\\infty\}a\_\{n\}g^\{n\}\. We denotes its Borel transform and the corresponding inverse transform:

B​\[f\]​\(t\)=∑n=0∞ann\!​\(t\)n\.B\[f\]\(t\)=\\sum\_\{n=0\}^\{\\infty\}\\frac\{a\_\{n\}\}\{n\!\}\(t\)^\{n\}\\,\.\(82\)And

ℬ​\(B​\[f\]\)​\(g\)=1g​∫0∞e−t/g​B​\(t\)\.\\mathcal\{B\}\(B\[f\]\)\(g\)=\\frac\{1\}\{g\}\\int\_\{0\}^\{\\infty\}e^\{\-t/g\}B\(t\)\\,\.\(83\)If an integral representation also exists:

f​\(g\)=∫e−S​\(x\)/g​𝑑x\.f\(g\)=\\int e^\{\-S\(x\)/g\}\\,dx\\,\.\(84\)Then we can write it on a level set then integrate over the choice of level:

f​\(g\)=1g​∫0∞𝑑t​e−t/g​\(g​∫𝑑x​δ​\(t−S​\(x\)\)\),f\(g\)=\\frac\{1\}\{g\}\\int\_\{0\}^\{\\infty\}dt\\,e^\{\-t/g\}\\left\(g\\int dx\\,\\delta\(t\-S\(x\)\)\\right\)\\,,\(85\)where we have already replaced the exponentiatedS​\(x\)S\(x\)withtt\. Using the co\-area trick \([81](https://arxiv.org/html/2606.00322#A2.E81)\), we have

f​\(g\)=1g​∫0∞𝑑t​e−t/g​\(g​∫t=S​\(x\)d​σ​\(x\)\|∇S​\(x\)\|\),f\(g\)=\\frac\{1\}\{g\}\\int\_\{0\}^\{\\infty\}dt\\,e^\{\-t/g\}\\,\\left\(g\\int\_\{t=S\(x\)\}\\frac\{d\\sigma\(x\)\}\{\|\\nabla S\(x\)\|\}\\right\)\\,,\(86\)whered​σ​\(x\)d\\sigma\(x\)is the appropriate measure\. Comparing this with the Laplace transform formula \([83](https://arxiv.org/html/2606.00322#A2.E83)\)\. We see that the Borel transform can be written as

B​\[1g​f\]​\(t\)=∫t=S​\(x\)d​σ​\(x\)\|∇S​\(x\)\|\.B\\left\[\\frac\{1\}\{g\}\\,f\\right\]\(t\)=\\int\_\{t=S\(x\)\}\\frac\{d\\sigma\(x\)\}\{\|\\nabla S\(x\)\|\}\\,\.\(87\)It is easy to see that singularities of the Borel transform occur precisely when∇S​\(x\)=0\\nabla S\(x\)=0, which are critical points of the integral representation\. This is a fact that can be extremely useful for more general usage of resurgence in neural networks141414To be elaborated in a separate section, off the the scope of this project\.\.

A further way to express the Borel transform derived from this is

B​\[f​\(g\)\]​\(t\)=∫dn​x​Θ​\(t−S​\(x\)\),B\[f\(g\)\]\(t\)=\\int d^\{n\}x\\,\\Theta\(t\-S\(x\)\)\\,,\(88\)where we see Borel transform offfis the volume of domain for whichS​\(x\)≤tS\(x\)\\leq t\. There are three possible sources of divergence or singularity on Borel plane\.

- •S\(x\) is critical atS​\(x\)=t∗S\(x\)=t\_\{\*\}
- •S\(x\) goes to some finitet∗t\_\{\*\}asxxgoes to infinity, the volume integral diverges
- •At somet→t∗t\\to t\_\{\*\},S​\(x\)S\(x\)has infinite volume\.

### B\.2The algorithm

We include the detailed algorithm for resurgence[2](https://arxiv.org/html/2606.00322#alg2)and full implementation of renormalization of perturbative series here, which asymptotically has a time complexity of𝒪​\(n2\)\\mathcal\{O\}\(n^\{2\}\), withnndata points\.

Algorithm 1Perturbative NPIV with RenormalizationInput:Data

\(X,Z,Y\)\(X,Z,Y\), kernel function

KK, ridge regularization coefficient

λ\\lambda, base perturbation parameter

γ\\gamma, max order

MM,

ϵ=10−4\\epsilon=10^\{\-4\}
Output:Coefficient vectors

α\(i\)\\alpha^\{\(i\)\}
Compute from dataset

Ki​jK\_\{ij\},

hih\_\{i\},

Ki​j​kK\_\{ijk\}
Solve

\(K~\+λ​K\)​α\(0\)=h\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(0\)\}=hfor zeroth\-order solution

‖α\(0\)‖←\\\|\\alpha^\{\(0\)\}\\\|\\leftarrownorm of zeroth\-order coefficients

for

m=1m=1to

MMdo

Compute right\-hand side:

r​h​si\(m\)=−3​∑j,kKi​j​k​\(αj\(0\)​αk\(m−1\)\+⋯\+αj\(m−1\)​αk\(0\)\)rhs\_\{i\}^\{\(m\)\}=\-3\\sum\_\{j,k\}K\_\{ijk\}\(\\alpha^\{\(0\)\}\_\{j\}\\alpha^\{\(m\-1\)\}\_\{k\}\+\\dots\+\\alpha^\{\(m\-1\)\}\_\{j\}\\alpha^\{\(0\)\}\_\{k\}\)
Solve

\(K~\+λ​K\)​α\(m\)=r​h​s\(m\)\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(m\)\}=rhs^\{\(m\)\}
‖α\(m\)‖←\\\|\\alpha^\{\(m\)\}\\\|\\leftarrownorm of

mm\-th order coefficients

ifrenormalization neededthen

Compute adaptive gamma:

γm=γm×min⁡\(1,‖α\(0\)‖/‖α\(m\)‖\)\\gamma\_\{m\}=\\gamma^\{m\}\\times\\min\(1,\\\|\\alpha^\{\(0\)\}\\\|/\\\|\\alpha^\{\(m\)\}\\\|\)
if

‖α\(0\)‖/‖α\(m\)‖<ϵ\\\|\\alpha^\{\(0\)\}\\\|/\\\|\\alpha^\{\(m\)\}\\\|<\\epsilonthen

γm=0\\gamma\_\{m\}=0\{Zero out correction if too unstable\}

endif

else

γm=γm\\gamma\_\{m\}=\\gamma^\{m\}\{Use standard perturbation parameter\}

endif

endfor

return

\{α\(n\)\}n=0M\\\{\\alpha^\{\(n\)\}\\\}\_\{n=0\}^\{M\}

Algorithm 2Resurgence algorithm for asymptotic seriesInput:the raw perturbative series

\{α\(i\)\}i=0M\\\{\\alpha^\{\(i\)\}\\\}\_\{i=0\}^\{M\},

\{γi\}i=1M\\\{\\gamma\_\{i\}\\\}\_\{i=1\}^\{M\}
ifNot asymptotic series:then

returnrenormalised series:

g​\(xn​e​w\)=∑n=0M∑iγn​αi\(n\)​K​\(xn​e​w,xi\)g\(x\_\{new\}\)=\\sum\_\{n=0\}^\{M\}\\sum\_\{i\}\\gamma\_\{n\}\\alpha\_\{i\}^\{\(n\)\}K\(x\_\{new\},x\_\{i\}\)
else

for

i=1i=1to

\|α\|\|\\alpha\|do

Compute Borel transform:

α^i​\(s\)=∑n=0Mαi\(n\)n\!​sn\\hat\{\\alpha\}\_\{i\}\(s\)=\\sum\_\{n=0\}^\{M\}\\frac\{\\alpha\_\{i\}^\{\(n\)\}\}\{n\!\}s^\{n\}
endfor

Construct Borel transformed function:

g^​\(t​X\)=∑iα^i​\(γ​t\)​K​\(t​X,Xi\)\\hat\{g\}\(tX\)=\\sum\_\{i\}\\hat\{\\alpha\}\_\{i\}\(\\gamma t\)K\(tX,X\_\{i\}\)
Initialize singularity set:

𝒮=∅\\mathcal\{S\}=\\emptyset
for

tton complex contour around origindo

if

\|g^​\(t​X\)\|\>1δ\|\\hat\{g\}\(tX\)\|\>\\frac\{1\}\{\\delta\}or

g^​\(t​X\)\\hat\{g\}\(tX\)undefinedthen

𝒮←𝒮∪\{t\}\\mathcal\{S\}\\leftarrow\\mathcal\{S\}\\cup\\\{t\\\}

endif

endfor

Let

\{tk\}k=1m\\\{t\_\{k\}\\\}\_\{k=1\}^\{m\}be the dominant singularities from

𝒮\\mathcal\{S\}
for

k=1k=1to

mmdo

Compute Stokes constant via residue theorem:

Ck←Rest=tk​g^​\(t​X\)C\_\{k\}\\leftarrow\\text\{Res\}\_\{t=t\_\{k\}\}\\hat\{g\}\(tX\)\{Residue at singularity

tkt\_\{k\}\}

for

i=1i=1to

\|α\|\|\\alpha\|do

gk,i​\(X\)←g\_\{k,i\}\(X\)\\leftarrowanalytic continuation of

α^i\\hat\{\\alpha\}\_\{i\}around

tkt\_\{k\}
endfor

endfor

forprediction point

xn​e​wx\_\{new\}do

Base perturbative series:

gp​e​r​t​\(xn​e​w\)←∑n=0M∑iγn​αi\(n\)​K​\(xn​e​w,xi\)g\_\{pert\}\(x\_\{new\}\)\\leftarrow\\sum\_\{n=0\}^\{M\}\\sum\_\{i\}\\gamma\_\{n\}\\alpha\_\{i\}^\{\(n\)\}K\(x\_\{new\},x\_\{i\}\)
Trans\-series corrections:

gt​r​a​n​s​\(xn​e​w\)←∑k=1mCk​exp⁡\(−k​tk​\(xn​e​w\)xn​e​w\)​∑igk,i​\(xn​e​w\)​K​\(xn​e​w,xi\)g\_\{trans\}\(x\_\{new\}\)\\leftarrow\\sum\_\{k=1\}^\{m\}C\_\{k\}\\exp\\left\(\-k\\frac\{t\_\{k\}\(x\_\{new\}\)\}\{x\_\{new\}\}\\right\)\\sum\_\{i\}g\_\{k,i\}\(x\_\{new\}\)K\(x\_\{new\},x\_\{i\}\)

Combined prediction:

Φ​\(xn​e​w\)←gp​e​r​t​\(xn​e​w\)\+gt​r​a​n​s​\(xn​e​w\)\\Phi\(x\_\{new\}\)\\leftarrow g\_\{pert\}\(x\_\{new\}\)\+g\_\{trans\}\(x\_\{new\}\)
endfor

return

Φ​\(X\)\\Phi\(X\)
endif

return

\{α\(i\)\}i=0M\\\{\\alpha^\{\(i\)\}\\\}\_\{i=0\}^\{M\},

\{γi\}i=1M\\\{\\gamma\_\{i\}\\\}\_\{i=1\}^\{M\},

Φ​\(x\)\\Phi\(x\)

## Appendix CThe irrelevance of physical renormalization

In physics, the Wilsonian paradigm of renormalization is based on the following philosophy\. If there is a scaleΛ\\Lambdaone introduces in a physical system beyond which the degrees of freedom is not considered\. The physical observables of the system should not depend on the magnitude ofΛ\\Lambda\. This requirement creates a spectral flow automatically called the renormalization group flow, adjusting the coupling constant accordingly\. In our case, there is a beautiful analogy\. Recall the series ansatz for the regressed solution:

g​\(x\)=∑n=0∞∑iγn​α\(n\)​\(Xi\)​K​\(Xi,X\)\.g\(x\)=\\sum\_\{n=0\}^\{\\infty\}\\sum\_\{i\}\\gamma^\{n\}\\alpha^\{\(n\)\}\(X\_\{i\}\)K\(X\_\{i\},X\)\\,\.\(89\)We shall argue that the regularization scaleλ\\lambdaautomatically serves as the cut\-off scaleΛ\\Lambda\. To see this, we define the following quantity:

ρ=λmλm\+λ,\\rho=\\frac\{\\lambda\_\{m\}\}\{\\lambda\_\{m\}\+\\lambda\}\\,,\(90\)which whenλ≫λm\\lambda\\gg\\lambda\_\{m\},ρ→0\\rho\\to 0and whenλ≪λm\\lambda\\ll\\lambda\_\{m\},ρ→1\\rho\\to 1\. Hence it roughly is measuring the amount of eigenvalues larger thanλ\\lambda, the ones belowλ\\lambdado not really contribute\. Taking the target functiong​\(X\)g\(X\), a physical renormalization process is triggered by recognizing thatg​\(X\)g\(X\)should not depend on the choice ofλ\\lambda151515This is not quite the case in regression, sinceλ\\lambdacontrols the bias/variance trade off, its scale physically impacts the regressed result\.\. Although it is intriguing to make the exact identification with physics, it is not quite correct due the aforementioned reason\. We think it is nevertheless important to include this in the script\.

Nevertheless, if we follow this logic, we see that scale independence ofg​\(X\)g\(X\)implies the Euler vector action:

λ​dd​λ​g​\(X\)=0\.\\lambda\\frac\{d\}\{d\\lambda\}\\,g\(X\)=0\\,\.\(91\)This can not be achieved ifγ\\gammais independent ofλ\\lambda, hence this physical requirement and the existence of imbalanced condition number creates a run on the couplingγ\\gamma, demanding it to depend onλ\\lambda\. More specifically, after some computations, one finds:

β​\(λ\)=λ​dd​λ=−∑nγn​∑mλλm\+λ​\(∂λ\(ϕmT​RHSn\)−αm\(n\)​\(λ\)\)​ϕm∑nγn−1​∑mαm\(n\)​\(λ\)​ϕm,\\beta\(\\lambda\)=\\lambda\\frac\{d\}\{d\\lambda\}=\-\\frac\{\\sum\_\{n\}\\gamma\_\{n\}\\sum\_\{m\}\\frac\{\\lambda\}\{\\lambda\_\{m\}\+\\lambda\}\\left\(\\partial\_\{\\lambda\}\(\\phi\_\{m\}^\{T\}\\,\\text\{RHS\}\_\{n\}\)\-\\alpha^\{\(n\)\}\_\{m\}\(\\lambda\)\\right\)\\phi\_\{m\}\}\{\\sum\_\{n\}\\gamma\_\{n\-1\}\\sum\_\{m\}\\alpha^\{\(n\)\}\_\{m\}\(\\lambda\)\\phi\_\{m\}\}\\,,\(92\)where all regimes ofλm\\lambda\_\{m\}contributes\. Our method of rescaling the couplingγ\\gammabased on the norm ratio of theα\(n\)\\alpha^\{\(n\)\}s is inspired by the physical method of renormalization\.

## Appendix DGeneralized Representer Theorem for Perturbative NPIV

We extend the classical representer theorem to accommodate higher\-order interaction terms\. This section provides a derivation of why our series expansion method maintains a finite\-dimensional representation at each order, even with cubic interaction terms\.

#### Problem Formulation and Notation

Let us begin by precisely defining the mathematical framework:

###### Definition 4\(Reproducing Kernel Hilbert Space\)\.

A Reproducing Kernel Hilbert Space \(RKHS\)ℋ\\mathcal\{H\}with kernelK:𝒳×𝒳→ℝK:\\mathcal\{X\}\\times\\mathcal\{X\}\\rightarrow\\mathbb\{R\}is a Hilbert space of functionsf:𝒳→ℝf:\\mathcal\{X\}\\rightarrow\\mathbb\{R\}such that:

1. 1\.For allx∈𝒳x\\in\\mathcal\{X\}, the functionK​\(⋅,x\)K\(\\cdot,x\)belongs toℋ\\mathcal\{H\}\.
2. 2\.For allx∈𝒳x\\in\\mathcal\{X\}and allf∈ℋf\\in\\mathcal\{H\}, the reproducing property holds:f​\(x\)=⟨f,K​\(⋅,x\)⟩ℋf\(x\)=\\langle f,K\(\\cdot,x\)\\rangle\_\{\\mathcal\{H\}\}\.

###### Definition 5\(Conditional Expectation Operators\)\.

LetXXandZZbe random variables with joint distributionPX,ZP\_\{X,Z\}\. We define:

1. 1\.The conditional expectation operatorT:ℋ→L2​\(PZ\)T:\\mathcal\{H\}\\rightarrow L^\{2\}\(P\_\{Z\}\)as\(T​f\)​\(z\)=𝔼​\[f​\(X\)\|Z=z\]\(Tf\)\(z\)=\\mathbb\{E\}\[f\(X\)\|Z=z\]\.
2. 2\.The adjoint operatorT∗:L2​\(PZ\)→ℋT^\{\*\}:L^\{2\}\(P\_\{Z\}\)\\rightarrow\\mathcal\{H\}satisfies⟨T​f,g⟩L2​\(PZ\)=⟨f,T∗​g⟩ℋ\\langle Tf,g\\rangle\_\{L^\{2\}\(P\_\{Z\}\)\}=\\langle f,T^\{\*\}g\\rangle\_\{\\mathcal\{H\}\}for allf∈ℋf\\in\\mathcal\{H\}andg∈L2​\(PZ\)g\\in L^\{2\}\(P\_\{Z\}\)\.

We consider the nonparametric instrumental variable \(NPIV\) problem with a cubic interaction term:

S​\[f\]\\displaystyle S\[f\]=S0​\[f\]\+γ​S1​\[f\]\\displaystyle=S\_\{0\}\[f\]\+\\gamma S\_\{1\}\[f\]\(93\)=𝔼Z​\[\(𝔼​\[Y\|Z\]−𝔼​\[f​\(X\)\|Z\]\)2\]\+λ​‖f‖ℋ2\+γ​S1​\[f\],\\displaystyle=\\mathbb\{E\}\_\{Z\}\\left\[\\left\(\\mathbb\{E\}\[Y\|Z\]\-\\mathbb\{E\}\[f\(X\)\|Z\]\\right\)^\{2\}\\right\]\+\\lambda\\\|f\\\|^\{2\}\_\{\\mathcal\{H\}\}\+\\gamma S\_\{1\}\[f\]\\,,\(94\)whereS1​\[f\]S\_\{1\}\[f\]represents the non\-independent three\-point interaction:

S1​\[f\]=23​𝔼Z​\[𝔼​\[f​\(X\)3\|Z\]\]\.\\displaystyle S\_\{1\}\[f\]=\\frac\{2\}\{3\}\\mathbb\{E\}\_\{Z\}\\left\[\\mathbb\{E\}\\left\[f\(X\)^\{3\}\\,\\bigg\|\\,Z\\right\]\\right\]\\,\.\(95\)The expectation overZZ\(notated as𝔼Z\\mathbb\{E\}\_\{Z\}\) is needed since the moment condition𝔼​\[Y−g​\(X\)\|Z\]=0\\mathbb\{E\}\[Y\-g\(X\)\|Z\]=0requires considering all possible values ofZZ, it is also what we compute in the acutal implementation\.

Our goal is to find the functionf∗∈ℋf^\{\*\}\\in\\mathcal\{H\}that minimizesS​\[f\]S\[f\]\. We approach this using a perturbative expansion:

f​\(x\)=f0​\(x\)\+γ​f1​\(x\)\+γ2​f2​\(x\)\+…\.\\displaystyle f\(x\)=f\_\{0\}\(x\)\+\\gamma f\_\{1\}\(x\)\+\\gamma^\{2\}f\_\{2\}\(x\)\+\\ldots\\,\.\(96\)

#### Functional Derivatives and Variational Calculus

To apply variational methods, we must understand functional derivatives:

###### Definition 6\(Functional Derivative\)\.

The functional derivative of a functionalS​\[f\]S\[f\]with respect to the functionffat pointxx, denotedδ​S​\[f\]δ​f​\(x\)\\frac\{\\delta S\[f\]\}\{\\delta f\(x\)\}, represents the rate of change ofS​\[f\]S\[f\]whenffis perturbed infinitesimally at the pointxx:

δ​S​\[f\]δ​f​\(x\)=limϵ→0S​\[f\+ϵ​δx\]−S​\[f\]ϵ,\\displaystyle\\frac\{\\delta S\[f\]\}\{\\delta f\(x\)\}=\\lim\_\{\\epsilon\\to 0\}\\frac\{S\[f\+\\epsilon\\delta\_\{x\}\]\-S\[f\]\}\{\\epsilon\}\\,,\(97\)whereδx\\delta\_\{x\}is the Dirac delta function centered atxx\.

The adjoint operatorT∗T^\{\*\}plays a crucial role in our derivations:

###### Proposition 2\(Adjoint ofTTin an RKHS\)\.

For\(T​f\)​\(z\)=𝔼​\[f​\(X\)∣Z=z\]\(Tf\)\(z\)=\\mathbb\{E\}\[f\(X\)\\mid Z=z\]with RKHS kernelKK, the adjointT∗:L2​\(PZ\)→ℋT^\{\*\}:L^\{2\}\(P\_\{Z\}\)\\to\\mathcal\{H\}is

T∗​g=𝔼​\[g​\(Z\)​K​\(⋅,X\)\]\.T^\{\*\}g\\;=\\;\\mathbb\{E\}\\\!\\big\[g\(Z\)\\,K\(\\cdot,X\)\\big\]\.

###### Proof\.

For anyf∈ℋf\\in\\mathcal\{H\}andg∈L2​\(PZ\)g\\in L^\{2\}\(P\_\{Z\}\),

⟨T​f,g⟩L2​\(PZ\)=𝔼​\[𝔼​\[f​\(X\)∣Z\]​g​\(Z\)\]=𝔼​\[⟨f,K​\(⋅,X\)⟩ℋ​g​\(Z\)\]=⟨f,𝔼​\[g​\(Z\)​K​\(⋅,X\)\]⟩ℋ\.\\langle Tf,g\\rangle\_\{L^\{2\}\(P\_\{Z\}\)\}=\\mathbb\{E\}\\\!\\big\[\\mathbb\{E\}\[f\(X\)\\mid Z\]\\;g\(Z\)\\big\]=\\mathbb\{E\}\\\!\\big\[\\langle f,K\(\\cdot,X\)\\rangle\_\{\\mathcal\{H\}\}\\;g\(Z\)\\big\]=\\langle f,\\,\\mathbb\{E\}\[g\(Z\)\\,K\(\\cdot,X\)\]\\rangle\_\{\\mathcal\{H\}\}\.\(98\)∎

###### Proposition 3\(Explicit Form of the Adjoint Operator\)\.

For the conditional expectation operatorTTdefined above, the adjoint operatorT∗:L2​\(PZ\)→ℋT^\{\*\}:L^\{2\}\(P\_\{Z\}\)\\rightarrow\\mathcal\{H\}can be explicitly written as:

\(T∗​g\)​\(x\)=𝔼​\[g​\(Z\)\|X=x\]\.\\displaystyle\(T^\{\*\}g\)\(x\)=\\mathbb\{E\}\[g\(Z\)\|X=x\]\\,\.\(99\)

###### Proof\.

For anyf∈ℋf\\in\\mathcal\{H\}andg∈L2​\(PZ\)g\\in L^\{2\}\(P\_\{Z\}\), we have:

⟨T​f,g⟩L2​\(PZ\)\\displaystyle\\langle Tf,g\\rangle\_\{L^\{2\}\(P\_\{Z\}\)\}=∫\(T​f\)​\(z\)​g​\(z\)​𝑑PZ​\(z\)\\displaystyle=\\int\(Tf\)\(z\)g\(z\)dP\_\{Z\}\(z\)\(100\)=∫𝔼​\[f​\(X\)\|Z=z\]​g​\(z\)​𝑑PZ​\(z\)\\displaystyle=\\int\\mathbb\{E\}\[f\(X\)\|Z=z\]g\(z\)dP\_\{Z\}\(z\)\(101\)=∬f​\(x\)​𝑑PX\|Z​\(x\|z\)​g​\(z\)​𝑑PZ​\(z\)\\displaystyle=\\iint f\(x\)dP\_\{X\|Z\}\(x\|z\)g\(z\)dP\_\{Z\}\(z\)\(102\)=∬f​\(x\)​g​\(z\)​𝑑PX,Z​\(x,z\)\\displaystyle=\\iint f\(x\)g\(z\)dP\_\{X,Z\}\(x,z\)\(103\)=∫f​\(x\)​\(∫g​\(z\)​𝑑PZ\|X​\(z\|x\)\)​𝑑PX​\(x\)\\displaystyle=\\int f\(x\)\\left\(\\int g\(z\)dP\_\{Z\|X\}\(z\|x\)\\right\)dP\_\{X\}\(x\)\(104\)=∫f​\(x\)​𝔼​\[g​\(Z\)\|X=x\]​𝑑PX​\(x\)\\displaystyle=\\int f\(x\)\\mathbb\{E\}\[g\(Z\)\|X=x\]dP\_\{X\}\(x\)\(105\)=⟨f,T∗​g⟩ℋ\\displaystyle=\\langle f,T^\{\*\}g\\rangle\_\{\\mathcal\{H\}\}\(106\)Thus,\(T∗​g\)​\(x\)=𝔼​\[g​\(Z\)\|X=x\]\(T^\{\*\}g\)\(x\)=\\mathbb\{E\}\[g\(Z\)\|X=x\]\. ∎

#### Zeroth\-Order Solution

Forγ=0\\gamma=0, we consider the standard regularized NPIV problem:

minf∈ℋ⁡S0​\[f\]=minf∈ℋ⁡\{𝔼Z​\[\(𝔼​\[Y\|Z\]−𝔼​\[f​\(X\)\|Z\]\)2\]\+λ​‖f‖ℋ2\}\.\\min\_\{f\\in\\mathcal\{H\}\}S\_\{0\}\[f\]=\\min\_\{f\\in\\mathcal\{H\}\}\\left\\\{\\mathbb\{E\}\_\{Z\}\\left\[\\left\(\\mathbb\{E\}\[Y\|Z\]\-\\mathbb\{E\}\[f\(X\)\|Z\]\\right\)^\{2\}\\right\]\+\\lambda\\\|f\\\|^\{2\}\_\{\\mathcal\{H\}\}\\right\\\}\\,\.\(107\)
Using operator notation, we can write this as:

S0​\[f\]=‖h−T​f‖L2​\(PZ\)2\+λ​‖f‖ℋ2,\\displaystyle S\_\{0\}\[f\]=\\\|h\-Tf\\\|\_\{L^\{2\}\(P\_\{Z\}\)\}^\{2\}\+\\lambda\\\|f\\\|^\{2\}\_\{\\mathcal\{H\}\}\\,,\(108\)whereh​\(z\)=𝔼​\[Y\|Z=z\]h\(z\)=\\mathbb\{E\}\[Y\|Z=z\]\.

###### Theorem 4\(Zeroth\-Order Representer Theorem\)\.

The minimizerf0f\_\{0\}ofS0​\[f\]S\_\{0\}\[f\]has the form:

f0​\(x\)=∑i=1nαi\(0\)​K​\(x,xi\),\\displaystyle f\_\{0\}\(x\)=\\sum\_\{i=1\}^\{n\}\\alpha\_\{i\}^\{\(0\)\}K\(x,x\_\{i\}\)\\,,\(109\)where\{x1,x2,…,xn\}\\\{x\_\{1\},x\_\{2\},\\ldots,x\_\{n\}\\\}are the observed data points\. The coefficientsα\(0\)\\alpha^\{\(0\)\}are the solution to:

\(K~\+λ​K\)​α\(0\)=h,\\displaystyle\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(0\)\}=h\\,,\(110\)whereK~i​j=𝔼​\[∫∫K​\(x,xi\)​K​\(x′,xj\)​f​\(x\|Z\)​f​\(x′\|Z\)​𝑑x​𝑑x′\]\\widetilde\{K\}\_\{ij\}=\\mathbb\{E\}\[\\int\\int K\(x,x\_\{i\}\)K\(x^\{\\prime\},x\_\{j\}\)f\(x\|Z\)f\(x^\{\\prime\}\|Z\)dxdx^\{\\prime\}\],Ki​j=K​\(xi,xj\)K\_\{ij\}=K\(x\_\{i\},x\_\{j\}\)andhi=\(T∗​h\)​\(xi\)h\_\{i\}=\(T^\{\*\}h\)\(x\_\{i\}\)\.

###### Proof\.

Taking the functional derivative ofS0​\[f\]S\_\{0\}\[f\]and setting it to zero:

δ​S0​\[f\]δ​f​\(x\)\\displaystyle\\frac\{\\delta S\_\{0\}\[f\]\}\{\\delta f\(x\)\}=−⟨h−T​f,T​δx⟩L2​\(PZ\)\+λ​⟨f,δx⟩ℋ\\displaystyle=\-\\langle h\-Tf,T\\delta\_\{x\}\\rangle\_\{L^\{2\}\(P\_\{Z\}\)\}\+\\lambda\\langle f,\\delta\_\{x\}\\rangle\_\{\\mathcal\{H\}\}\(111\)=−\(T∗​\(h−T​f\)\)​\(x\)\+λ​f​\(x\),\\displaystyle=\-\(T^\{\*\}\(h\-Tf\)\)\(x\)\+\\lambda f\(x\)\\,,\(112\)where we’ve used the definition of the adjoint operatorT∗T^\{\*\}\. Setting this to zero:

\(T∗​\(h−T​f\)\)​\(x\)=λ​f​\(x\)\\displaystyle\(T^\{\*\}\(h\-Tf\)\)\(x\)=\\lambda f\(x\)\(113\)
By the standard representer theorem argument decomposing the function in basis in the RKHS and orthogonal, minimizing such objective easily gives us the form of the solution, the solution has the formf0​\(x\)=∑i=1nαi\(0\)​K​\(x,xi\)f\_\{0\}\(x\)=\\sum\_\{i=1\}^\{n\}\\alpha\_\{i\}^\{\(0\)\}K\(x,x\_\{i\}\)\. Substituting this and evaluating at a finite collection of data points\{xi\}i=1n\\\{x\_\{i\}\\\}\_\{i=1\}^\{n\}gives the linear system:

\(K~\+λ​K\)​α\(0\)=h,\\displaystyle\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(0\)\}=h\\,,\(114\)wherehi=\(T∗​h\)​\(xi\)=𝔼​\[Y​K​\(X,xi\)\|Z\]h\_\{i\}=\(T^\{\*\}h\)\(x\_\{i\}\)=\\mathbb\{E\}\[Y\\,K\(X,x\_\{i\}\)\|Z\]\. ∎

#### First\-Order Correction

For the first\-order correction, we need to findf1f\_\{1\}by expandingS​\[f\]S\[f\]aroundf0f\_\{0\}:

S​\[f0\+γ​f1\]=S0​\[f0\]\+γ​\(δ​S0​\[f\]δ​f\|f=f0⋅f1\+S1​\[f\]\|f=f0\)\+O​\(γ2\)\.S\[f\_\{0\}\+\\gamma f\_\{1\}\]=S\_\{0\}\[f\_\{0\}\]\+\\gamma\\left\(\\frac\{\\delta S\_\{0\}\[f\]\}\{\\delta f\}\\bigg\|\_\{f=f\_\{0\}\}\\cdot f\_\{1\}\+S\_\{1\}\[f\]\\bigg\|\_\{f=f\_\{0\}\}\\right\)\+O\(\\gamma^\{2\}\)\\,\.\(115\)Setting the variational derivative with respect tof1f\_\{1\}to zero, using

δ​S0​\[f\]δ​f\|f=f0=T∗​\(h−T​f0\)\+λ​f0\.\\displaystyle\\frac\{\\delta S\_\{0\}\[f\]\}\{\\delta f\}\\bigg\|\_\{f=f\_\{0\}\}=T^\{\*\}\(h\-Tf\_\{0\}\)\+\\lambda f\_\{0\}\\,\.\(116\)Using operator notation, the equation forf1f\_\{1\}becomes:

\(T∗​T\+λ​I\)​f1=−δ​S1​\[f\]δ​f\|f=f0\.\\displaystyle\(T^\{\*\}T\+\\lambda I\)f\_\{1\}=\-\\frac\{\\delta S\_\{1\}\[f\]\}\{\\delta f\}\\bigg\|\_\{f=f\_\{0\}\}\\,\.\(117\)Then for the second termS1​\[f\]=23​𝔼Z​\[𝔼​\[f​\(X\)3\|Z\]\]S\_\{1\}\[f\]=\\frac\{2\}\{3\}\\mathbb\{E\}\_\{Z\}\[\\mathbb\{E\}\[f\(X\)^\{3\}\|Z\]\]:

δ​S1​\[f\]δ​f​\(x\)\\displaystyle\\frac\{\\delta S\_\{1\}\[f\]\}\{\\delta f\(x\)\}=δδ​f​\(x\)​\(23​𝔼Z​\[𝔼​\[f​\(X\)3\|Z\]\]\)\\displaystyle=\\frac\{\\delta\}\{\\delta f\(x\)\}\\left\(\\frac\{2\}\{3\}\\mathbb\{E\}\_\{Z\}\[\\mathbb\{E\}\[f\(X\)^\{3\}\|Z\]\]\\right\)\(118\)=23​𝔼Z​\[δδ​f​\(x\)​𝔼​\[f​\(X\)3\|Z\]\]\.\\displaystyle=\\frac\{2\}\{3\}\\mathbb\{E\}\_\{Z\}\\left\[\\frac\{\\delta\}\{\\delta f\(x\)\}\\mathbb\{E\}\[f\(X\)^\{3\}\|Z\]\\right\]\\,\.\(119\)For any functiong​\(X\)g\(X\), the functional derivative with respect tof​\(x\)f\(x\)is:

δ​g​\(X\)δ​f​\(x\)=∂g​\(X\)∂f​\(X\)​δ​f​\(X\)δ​f​\(x\)=∂g​\(X\)∂f​\(X\)​δ​\(X−x\)\.\\displaystyle\\frac\{\\delta g\(X\)\}\{\\delta f\(x\)\}=\\frac\{\\partial g\(X\)\}\{\\partial f\(X\)\}\\frac\{\\delta f\(X\)\}\{\\delta f\(x\)\}=\\frac\{\\partial g\(X\)\}\{\\partial f\(X\)\}\\delta\(X\-x\)\\,\.\(120\)Therefore:

δδ​f​\(x\)​𝔼​\[f​\(X\)3\|Z\]\\displaystyle\\frac\{\\delta\}\{\\delta f\(x\)\}\\mathbb\{E\}\[f\(X\)^\{3\}\|Z\]=𝔼​\[δ​f​\(X\)3δ​f​\(x\)\|Z\]\\displaystyle=\\mathbb\{E\}\\left\[\\frac\{\\delta f\(X\)^\{3\}\}\{\\delta f\(x\)\}\\bigg\|Z\\right\]\(121\)=𝔼​\[3​f​\(X\)2​δ​\(X−x\)\|Z\]\.\\displaystyle=\\mathbb\{E\}\\left\[3f\(X\)^\{2\}\\delta\(X\-x\)\\bigg\|Z\\right\]\\,\.\(122\)The presence of the delta functionδ​\(X−x\)\\delta\(X\-x\)inside the expectation requires careful treatment\. For a continuous random variableXXwith densityp​\(X\|Z\)p\(X\|Z\):

𝔼​\[δ​\(X−x\)​g​\(X\)\|Z\]\\displaystyle\\mathbb\{E\}\[\\delta\(X\-x\)g\(X\)\|Z\]=∫δ​\(y−x\)​g​\(y\)​p​\(y\|Z\)​𝑑y\\displaystyle=\\int\\delta\(y\-x\)g\(y\)p\(y\|Z\)dy\(123\)=g​\(x\)​p​\(x\|Z\)\.\\displaystyle=g\(x\)p\(x\|Z\)\\,\.\(124\)Thus:

𝔼​\[3​f​\(X\)2​δ​\(X−x\)\|Z\]=3​f​\(x\)2​p​\(x\|Z\)\.\\displaystyle\\mathbb\{E\}\[3f\(X\)^\{2\}\\delta\(X\-x\)\|Z\]=3f\(x\)^\{2\}p\(x\|Z\)\\,\.\(125\)Substituting back:

δ​S1​\[f\]δ​f​\(x\)\\displaystyle\\frac\{\\delta S\_\{1\}\[f\]\}\{\\delta f\(x\)\}=23​𝔼Z​\[3​f​\(x\)2​p​\(x\|Z\)\]\\displaystyle=\\frac\{2\}\{3\}\\mathbb\{E\}\_\{Z\}\[3f\(x\)^\{2\}p\(x\|Z\)\]\(126\)=2​f​\(x\)2​𝔼Z​\[p​\(x\|Z\)\]\\displaystyle=2f\(x\)^\{2\}\\mathbb\{E\}\_\{Z\}\[p\(x\|Z\)\]\(127\)=2​f​\(x\)2​p​\(x\),\\displaystyle=2f\(x\)^\{2\}p\(x\)\\,,\(128\)wherep​\(x\)p\(x\)is the marginal density ofXX\.

However, in our RKHS setting, we work with functions rather than explicitly with densities\. In particular, when evaluating at our data points, we can use the connection to the adjoint operator and kernel functions\.

Evaluating the equation at a data pointxix\_\{i\}:

\(T∗​T\+λ​I\)​f1​\(xi\)\\displaystyle\(T^\{\*\}T\+\\lambda I\)f\_\{1\}\(x\_\{i\}\)=−f0​\(xi\)2​p​\(xi\)\.\\displaystyle=\-f\_\{0\}\(x\_\{i\}\)^\{2\}p\(x\_\{i\}\)\\,\.\(129\)Substituting the representer form off0f\_\{0\}:

f0​\(xi\)2=\(∑j=1nαj\(0\)​K​\(xi,xj\)\)2=∑j=1n∑k=1nαj\(0\)​αk\(0\)​K​\(xi,xj\)​K​\(xi,xk\)\.f\_\{0\}\(x\_\{i\}\)^\{2\}=\\left\(\\sum\_\{j=1\}^\{n\}\\alpha\_\{j\}^\{\(0\)\}K\(x\_\{i\},x\_\{j\}\)\\right\)^\{2\}=\\sum\_\{j=1\}^\{n\}\\sum\_\{k=1\}^\{n\}\\alpha\_\{j\}^\{\(0\)\}\\alpha\_\{k\}^\{\(0\)\}K\(x\_\{i\},x\_\{j\}\)K\(x\_\{i\},x\_\{k\}\)\\,\.\(130\)When we apply this in the RKHS context, the termp​\(xi\)p\(x\_\{i\}\)is effectively handled by kernel methods\. In particular, the right\-hand side of our equation becomes:

r1,i=−T∗​\[𝔼​\[∑j=1n∑k=1nαj\(0\)​αk\(0\)​K​\(X,xj\)​K​\(X,xk\)\|Z\]\]​\(xi\),r\_\{1,i\}=\-T^\{\*\}\\left\[\\mathbb\{E\}\\left\[\\sum\_\{j=1\}^\{n\}\\sum\_\{k=1\}^\{n\}\\alpha\_\{j\}^\{\(0\)\}\\alpha\_\{k\}^\{\(0\)\}K\(X,x\_\{j\}\)K\(X,x\_\{k\}\)\\,\\bigg\|\\,Z\\right\]\\right\]\(x\_\{i\}\)\\,,\(131\)whereZ=⋅Z=\\cdotis selected to be the appropriate values ofZZfor whichT∗T^\{\*\}reproduces the correct Hilbert space function\. By the properties of the adjoint operatorT∗T^\{\*\}, this is equivalent to averaging over the function conditioned on all possible instrumentalZZs, using properties of expectations, we can rearrange:

r1,i=−∑j=1n∑k=1nαj\(0\)αk\(0\)×\\displaystyle r\_\{1,i\}=\-\\sum\_\{j=1\}^\{n\}\\sum\_\{k=1\}^\{n\}\\alpha\_\{j\}^\{\(0\)\}\\alpha\_\{k\}^\{\(0\)\}\\times\(132\)𝔼Z​\[𝔼​\[K​\(X,xi\)​K​\(X,xj\)​K​\(X,xk\)\|Z\]\]\.\\displaystyle\\mathbb\{E\}\_\{Z\}\\left\[\\mathbb\{E\}\\left\[K\(X,x\_\{i\}\)K\(X,x\_\{j\}\)K\(X,x\_\{k\}\)\\,\\bigg\|\\,Z\\right\]\\right\]\\,\.\(133\)Due to the symmetries in the three\-point interaction and the specific form of our cubic term, an additional combinatorial factor of 3 appears, giving us:

r1,i=−3×\\displaystyle r\_\{1,i\}=\-3\\times\(134\)∑j=1n∑k=1nαj\(0\)​αk\(0\)​𝔼Z​\[𝔼​\[K​\(X,xi\)​K​\(X,xj\)​K​\(X,xk\)\|Z\]\]\.\\displaystyle\\sum\_\{j=1\}^\{n\}\\sum\_\{k=1\}^\{n\}\\alpha\_\{j\}^\{\(0\)\}\\alpha\_\{k\}^\{\(0\)\}\\mathbb\{E\}\_\{Z\}\\left\[\\mathbb\{E\}\\left\[K\(X,x\_\{i\}\)K\(X,x\_\{j\}\)K\(X,x\_\{k\}\)\\,\\bigg\|\\,Z\\right\]\\right\]\\,\.Finally, defining the non\-independent three\-point kernel tensor:

Ki​j​k=𝔼Z​\[𝔼​\[K​\(X,xi\)​K​\(X,xj\)​K​\(X,xk\)\|Z\]\]\.\\displaystyle K\_\{ijk\}=\\mathbb\{E\}\_\{Z\}\\left\[\\mathbb\{E\}\\left\[K\(X,x\_\{i\}\)K\(X,x\_\{j\}\)K\(X,x\_\{k\}\)\\,\\bigg\|\\,Z\\right\]\\right\]\\,\.\(135\)We arrive at:

r1,i\\displaystyle r\_\{1,i\}=−3​∑j=1n∑k=1nαj\(0\)​αk\(0\)​Ki​j​k\.\\displaystyle=\-3\\sum\_\{j=1\}^\{n\}\\sum\_\{k=1\}^\{n\}\\alpha\_\{j\}^\{\(0\)\}\\alpha\_\{k\}^\{\(0\)\}K\_\{ijk\}\\,\.\(136\)This matches exactly the computation in our algorithm\.

###### Theorem 5\(First\-Order Representer Theorem\)\.

The first\-order correctionf1f\_\{1\}has the representer form:

f1​\(x\)=∑i=1nαi\(1\)​K​\(x,xi\),\\displaystyle f\_\{1\}\(x\)=\\sum\_\{i=1\}^\{n\}\\alpha\_\{i\}^\{\(1\)\}K\(x,x\_\{i\}\)\\,,\(137\)whereα\(1\)\\alpha^\{\(1\)\}is the solution to the linear system:

\(K~\+λ​K\)​α\(1\)=r1,\\displaystyle\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(1\)\}=r\_\{1\}\\,,\(138\)withr1,i=−3​∑j=1n∑k=1nαj\(0\)​αk\(0\)​Ki​j​kr\_\{1,i\}=\-3\\sum\_\{j=1\}^\{n\}\\sum\_\{k=1\}^\{n\}\\alpha\_\{j\}^\{\(0\)\}\\alpha\_\{k\}^\{\(0\)\}K\_\{ijk\}\.

###### Proof\.

Since the operator\(T∗​T\+λ​I\)\(T^\{\*\}T\+\\lambda I\)maps the RKHS to itself while preserving the finite\-dimensional subspace spanned by the data points, and the right\-hand sider1r\_\{1\}lies in this subspace, the solutionf1f\_\{1\}must also have the representer form\. ∎

Essentially, the representer theorem does not incur any difficulty here either since the addition of anything orthogonal to the Hilbert space basis still projects to0, hence not improving the minization of the entire functional\.

#### Second\-Order and Higher Corrections

Following the same procedure, we can derive the second\-order correction:

S​\[f0\+γ​f1\+γ2​f2\]=S0​\[f0\]\+𝒪​\(γ\)\+γ2​\(δ​S0​\[f\]δ​f\|f=f0⋅f2\+12​δ2​S0​\[f\]δ​f2\|f=f0​\(f1,f1\)\+δ​S1​\[f\]δ​f\|f=f0⋅f1\)\+O​\(γ3\)\.S\[f\_\{0\}\+\\gamma f\_\{1\}\+\\gamma^\{2\}f\_\{2\}\]=S\_\{0\}\[f\_\{0\}\]\+\\mathcal\{O\}\(\\gamma\)\+\\\\ \\gamma^\{2\}\\left\(\\frac\{\\delta S\_\{0\}\[f\]\}\{\\delta f\}\\bigg\|\_\{f=f\_\{0\}\}\\cdot f\_\{2\}\+\\frac\{1\}\{2\}\\frac\{\\delta^\{2\}S\_\{0\}\[f\]\}\{\\delta f^\{2\}\}\\bigg\|\_\{f=f\_\{0\}\}\(f\_\{1\},f\_\{1\}\)\+\\frac\{\\delta S\_\{1\}\[f\]\}\{\\delta f\}\\bigg\|\_\{f=f\_\{0\}\}\\cdot f\_\{1\}\\right\)\+O\(\\gamma^\{3\}\)\\,\.\(139\)Setting the variational derivative with respect tof2f\_\{2\}to zero, we have the constraint\.

###### Theorem 6\(Second\-Order Representer Theorem\)\.

The second\-order correctionf2f\_\{2\}also has the representer form:

f2​\(x\)=∑i=1nαi\(2\)​K​\(x,xi\),\\displaystyle f\_\{2\}\(x\)=\\sum\_\{i=1\}^\{n\}\\alpha\_\{i\}^\{\(2\)\}K\(x,x\_\{i\}\)\\,,\(140\)whereα\(2\)\\alpha^\{\(2\)\}is the solution to the linear system:

\(K~\+λ​K\)​α\(2\)=r2,\\displaystyle\(\\widetilde\{K\}\+\\lambda K\)\\alpha^\{\(2\)\}=r\_\{2\}\\,,\(141\)withr2,i=−3​∑j=1n∑k=1n\(αj\(0\)​αk\(1\)\+αj\(1\)​αk\(0\)\)​Ki​j​kr\_\{2,i\}=\-3\\sum\_\{j=1\}^\{n\}\\sum\_\{k=1\}^\{n\}\(\\alpha\_\{j\}^\{\(0\)\}\\alpha\_\{k\}^\{\(1\)\}\+\\alpha\_\{j\}^\{\(1\)\}\\alpha\_\{k\}^\{\(0\)\}\)K\_\{ijk\}\.

#### General Result for All Orders

We can generalize these results to all orders in the perturbative expansion:

###### Theorem 7\(Series Expansion Representer Theorem\)\.

For the NPIV problem with a cubic interaction term, the solution can be expressed as a perturbative series:

f∗​\(x\)=∑k=0∞γk​fk​\(x\)=∑k=0∞γk​∑i=1nαi\(k\)​K​\(x,xi\),\\displaystyle f^\{\*\}\(x\)=\\sum\_\{k=0\}^\{\\infty\}\\gamma^\{k\}f\_\{k\}\(x\)=\\sum\_\{k=0\}^\{\\infty\}\\gamma^\{k\}\\sum\_\{i=1\}^\{n\}\\alpha\_\{i\}^\{\(k\)\}K\(x,x\_\{i\}\)\\,,\(142\)where eachα\(k\)\\alpha^\{\(k\)\}is determined by solving a linear system that depends on\{α\(j\)\}j=0k−1\\\{\\alpha^\{\(j\)\}\\\}\_\{j=0\}^\{k\-1\}\.

This is exactly what we have had in the algorithm\.

## Appendix EInteraction for rotationally invariant kernels

The original triple moment interaction we added was:

Ki​j​k=∫K​\(x,xi\)​K​\(x,xj\)​K​\(x,xk\)​f​\(x\|Z\)​𝑑x\.K\_\{ijk\}=\\int K\(x,x\_\{i\}\)K\(x,x\_\{j\}\)K\(x,x\_\{k\}\)f\(x\|Z\)dx\\,\.\(143\)
For rotationally invariant kernelsK​\(x,y\)=k​\(‖x−y‖\)K\(x,y\)=k\(\\\|x\-y\\\|\)with Mercer decomposition:

K​\(x,y\)=∑αλα​ϕα​\(x\)​ϕα​\(y\),K\(x,y\)=\\sum\_\{\\alpha\}\\lambda\_\{\\alpha\}\\phi\_\{\\alpha\}\(x\)\\phi\_\{\\alpha\}\(y\)\\,,\(144\)whereϕα​\(x\)=ϕl,m​\(x\)=Rl​\(‖x‖\)​Ylm​\(x/‖x‖\)\\phi\_\{\\alpha\}\(x\)=\\phi\_\{l,m\}\(x\)=R\_\{l\}\(\\\|x\\\|\)Y\_\{l\}^\{m\}\(x/\\\|x\\\|\)withYlm​\(x/‖x‖\)Y\_\{l\}^\{m\}\(x/\\\|x\\\|\)the usual spherical harmonics\.

As we have seen, these triple moments contribute to the final integral as the Clebsch\-Gordan or Wigner 3\-j symbols which has a strict selection rule for the labels, namely it is only non\-zero if

m1\+m2\+m3=0,l1,l2,l3​even and form triangle\.m\_\{1\}\+m\_\{2\}\+m\_\{3\}=0,\\quad l\_\{1\},l\_\{2\},l\_\{3\}\\text\{ even and form triangle\}\\,\.\(145\)This is rather sparse and zeros out most of the contributions\. The most restrictive conditions are actually the delta function inmmquantum number and the fact that thells satisfy the triangle inequality\. We would like to explicitly overcome this, we simply introduce a factor that breaks rotational invariance:

Ki​j​kaniso=∫K​\(x,xi\)​K​\(x,xj\)​K​\(x,xk\)​\[f​\(x\|Z\)\]3​WA​\(x\)​𝑑x,K\_\{ijk\}^\{\\text\{aniso\}\}=\\int K\(x,x\_\{i\}\)K\(x,x\_\{j\}\)K\(x,x\_\{k\}\)\[f\(x\|Z\)\]^\{3\}W\_\{A\}\(x\)dx\\,,\(146\)
where the anisotropic weight is:

WA​\(x\)=exp⁡\(−12​xT​A−1​x\),W\_\{A\}\(x\)=\\exp\\left\(\-\\frac\{1\}\{2\}x^\{T\}A^\{\-1\}x\\right\)\\,,\(147\)
withA=diag​\(a1,a2,…,an\)A=\\text\{diag\}\(a\_\{1\},a\_\{2\},\\ldots,a\_\{n\}\)andai\>0a\_\{i\}\>0all distinct\. With the anisotropic weightWA​\(x\)=exp⁡\(−12​xT​A−1​x\)W\_\{A\}\(x\)=\\exp\(\-\\frac\{1\}\{2\}x^\{T\}A^\{\-1\}x\):

𝒯α​β​γaniso=∫ϕα​\(x\)​ϕβ​\(x\)​ϕγ​\(x\)​\[f​\(x\|Z\)\]3​exp⁡\(−12​xT​A−1​x\)​𝑑x\.\\mathcal\{T\}\_\{\\alpha\\beta\\gamma\}^\{\\text\{aniso\}\}=\\int\\phi\_\{\\alpha\}\(x\)\\phi\_\{\\beta\}\(x\)\\phi\_\{\\gamma\}\(x\)\[f\(x\|Z\)\]^\{3\}\\exp\\left\(\-\\frac\{1\}\{2\}x^\{T\}A^\{\-1\}x\\right\)dx\\,\.\(148\)The anisotropic factor \*\*breaks rotational symmetry\*\*, making the integral generically non\-zero even when the standard Clebsch\-Gordan coefficient vanishes\.

Forϕα​\(x\)=Rl​\(‖x‖\)​Ylm​\(x/‖x‖\)\\phi\_\{\\alpha\}\(x\)=R\_\{l\}\(\\\|x\\\|\)Y\_\{l\}^\{m\}\(x/\\\|x\\\|\):

𝒯α​β​γaniso=∫ℝnRl1​\(‖x‖\)​Yl1m1​\(x/‖x‖\)​Rl2​\(‖x‖\)​Yl2m2​\(x/‖x‖\)​Rl3​\(‖x‖\)​Yl3m3​\(x/‖x‖\)​\[f​\(‖x‖\|Z\)\]3​exp⁡\(−12​xT​A−1​x\)​𝑑x\.\\mathcal\{T\}\_\{\\alpha\\beta\\gamma\}^\{\\text\{aniso\}\}=\\int\_\{\\mathbb\{R\}^\{n\}\}R\_\{l\_\{1\}\}\(\\\|x\\\|\)Y\_\{l\_\{1\}\}^\{m\_\{1\}\}\(x/\\\|x\\\|\)R\_\{l\_\{2\}\}\(\\\|x\\\|\)Y\_\{l\_\{2\}\}^\{m\_\{2\}\}\(x/\\\|x\\\|\)R\_\{l\_\{3\}\}\(\\\|x\\\|\)Y\_\{l\_\{3\}\}^\{m\_\{3\}\}\(x/\\\|x\\\|\)\[f\(\\\|x\\\|\|Z\)\]^\{3\}\\exp\\left\(\-\\frac\{1\}\{2\}x^\{T\}A^\{\-1\}x\\right\)dx\\,\.\(149\)
In spherical coordinatesx=r​ω^x=r\\hat\{\\omega\}:

𝒯α​β​γaniso=∫0∞∫𝕊n−1Rl1​\(r\)​Rl2​\(r\)​Rl3​\(r\)​\[f​\(r\|Z\)\]3​Yl1m1​\(ω^\)​Yl2m2​\(ω^\)​Yl3m3​\(ω^\)​exp⁡\(−12​r2​ω^T​A−1​ω^\)​rn−1​𝑑r​𝑑ω^\.\\mathcal\{T\}\_\{\\alpha\\beta\\gamma\}^\{\\text\{aniso\}\}=\\int\_\{0\}^\{\\infty\}\\int\_\{\\mathbb\{S\}^\{n\-1\}\}R\_\{l\_\{1\}\}\(r\)R\_\{l\_\{2\}\}\(r\)R\_\{l\_\{3\}\}\(r\)\[f\(r\|Z\)\]^\{3\}Y\_\{l\_\{1\}\}^\{m\_\{1\}\}\(\\hat\{\\omega\}\)Y\_\{l\_\{2\}\}^\{m\_\{2\}\}\(\\hat\{\\omega\}\)Y\_\{l\_\{3\}\}^\{m\_\{3\}\}\(\\hat\{\\omega\}\)\\exp\\left\(\-\\frac\{1\}\{2\}r^\{2\}\\hat\{\\omega\}^\{T\}A^\{\-1\}\\hat\{\\omega\}\\right\)r^\{n\-1\}drd\\hat\{\\omega\}\\,\.\(150\)The key angular integral is:

𝒜l1​m1,l2​m2,l3​m3=∫𝕊n−1Yl1m1​\(ω^\)​Yl2m2​\(ω^\)​Yl3m3​\(ω^\)​exp⁡\(−12​r2​ω^T​A−1​ω^\)​𝑑ω^\.\\mathcal\{A\}\_\{l\_\{1\}m\_\{1\},l\_\{2\}m\_\{2\},l\_\{3\}m\_\{3\}\}=\\int\_\{\\mathbb\{S\}^\{n\-1\}\}Y\_\{l\_\{1\}\}^\{m\_\{1\}\}\(\\hat\{\\omega\}\)Y\_\{l\_\{2\}\}^\{m\_\{2\}\}\(\\hat\{\\omega\}\)Y\_\{l\_\{3\}\}^\{m\_\{3\}\}\(\\hat\{\\omega\}\)\\exp\\left\(\-\\frac\{1\}\{2\}r^\{2\}\\hat\{\\omega\}^\{T\}A^\{\-1\}\\hat\{\\omega\}\\right\)d\\hat\{\\omega\}\\,\.\(151\)IfA−1=α​IA^\{\-1\}=\\alpha I\(isotropic\), then:

exp⁡\(−12​r2​ω^T​A−1​ω^\)=exp⁡\(−12​α​r2\)=const\.\.\\exp\\left\(\-\\frac\{1\}\{2\}r^\{2\}\\hat\{\\omega\}^\{T\}A^\{\-1\}\\hat\{\\omega\}\\right\)=\\exp\\left\(\-\\frac\{1\}\{2\}\\alpha r^\{2\}\\right\)=\\text\{const\.\}\\,\.\(152\)
This gives:

𝒜l1​m1,l2​m2,l3​m3=exp⁡\(−12​α​r2\)​∫𝕊n−1Yl1m1​\(ω^\)​Yl2m2​\(ω^\)​Yl3m3​\(ω^\)​𝑑ω^\.\\mathcal\{A\}\_\{l\_\{1\}m\_\{1\},l\_\{2\}m\_\{2\},l\_\{3\}m\_\{3\}\}=\\exp\\left\(\-\\frac\{1\}\{2\}\\alpha r^\{2\}\\right\)\\int\_\{\\mathbb\{S\}^\{n\-1\}\}Y\_\{l\_\{1\}\}^\{m\_\{1\}\}\(\\hat\{\\omega\}\)Y\_\{l\_\{2\}\}^\{m\_\{2\}\}\(\\hat\{\\omega\}\)Y\_\{l\_\{3\}\}^\{m\_\{3\}\}\(\\hat\{\\omega\}\)d\\hat\{\\omega\}\\,\.\(153\)
The integral is just the Clebsch\-Gordan coefficient, subject to selection rules\. WhenA−1=diag​\(a1−1,a2−1,…,an−1\)A^\{\-1\}=\\text\{diag\}\(a\_\{1\}^\{\-1\},a\_\{2\}^\{\-1\},\\ldots,a\_\{n\}^\{\-1\}\)with distinct values:

exp⁡\(−12​r2​ω^T​A−1​ω^\)=exp⁡\(−12​r2​∑j=1naj−1​ωj2\),\\exp\\left\(\-\\frac\{1\}\{2\}r^\{2\}\\hat\{\\omega\}^\{T\}A^\{\-1\}\\hat\{\\omega\}\\right\)=\\exp\\left\(\-\\frac\{1\}\{2\}r^\{2\}\\sum\_\{j=1\}^\{n\}a\_\{j\}^\{\-1\}\\omega\_\{j\}^\{2\}\\right\)\\,,\(154\)whereωj\\omega\_\{j\}are the components ofω^\\hat\{\\omega\}in Cartesian coordinates\. We consider the key angular integral \([151](https://arxiv.org/html/2606.00322#A5.E151)\)\.

The integral is subject to a strict parity selection rule\. The spherical harmonic has the propertyYlm​\(−ω^\)=\(−1\)l​Ylm​\(ω^\)Y\_\{l\}^\{m\}\(\-\\hat\{\\omega\}\)=\(\-1\)^\{l\}Y\_\{l\}^\{m\}\(\\hat\{\\omega\}\)\. The exponential term is an even function with respect to the transformationω^→−ω^\\hat\{\\omega\}\\to\-\\hat\{\\omega\}\. The parity of the integrand is therefore determined solely by the product of the three spherical harmonics:

Integrand​\(−ω^\)=\(−1\)l1\+l2\+l3×Integrand​\(ω^\)\\text\{Integrand\}\(\-\\hat\{\\omega\}\)=\(\-1\)^\{l\_\{1\}\+l\_\{2\}\+l\_\{3\}\}\\times\\text\{Integrand\}\(\\hat\{\\omega\}\)As the integral is over a symmetric domain, the integral of any odd function is zero\.Conclusion:The integralIIis zero unless the suml1\+l2\+l3l\_\{1\}\+l\_\{2\}\+l\_\{3\}is an even integer\.

Assumingl1\+l2\+l3l\_\{1\}\+l\_\{2\}\+l\_\{3\}is even, the integral can be computed by expanding the exponential term in a power series:

exp⁡\(−12​ρ2​ω^T​A​ω^\)=∑k=0∞1k\!​\(−ρ22\)k​\(ω^T​A​ω^\)k\.\\exp\\left\(\-\\frac\{1\}\{2\}\\rho^\{2\}\\hat\{\\omega\}^\{T\}A\\hat\{\\omega\}\\right\)=\\sum\_\{k=0\}^\{\\infty\}\\frac\{1\}\{k\!\}\\left\(\-\\frac\{\\rho^\{2\}\}\{2\}\\right\)^\{k\}\(\\hat\{\\omega\}^\{T\}A\\hat\{\\omega\}\)^\{k\}\\,\.\(155\)Substituting this into the integral and interchanging summation and integration yields a series representation forII:

I=∑k=0∞\(−ρ2/2\)kk\!​ℳk,I=\\sum\_\{k=0\}^\{\\infty\}\\frac\{\(\-\\rho^\{2\}/2\)^\{k\}\}\{k\!\}\\mathcal\{M\}\_\{k\}\\,,\(156\)whereℳk\\mathcal\{M\}\_\{k\}are the moment integrals:

ℳk=∫𝕊n−1Yl1m1​\(ω^\)​Yl2m2​\(ω^\)​Yl3m3​\(ω^\)​\(∑i=1nλi​ωi2\)k​𝑑ω^\.\\mathcal\{M\}\_\{k\}=\\int\_\{\\mathbb\{S\}^\{n\-1\}\}Y\_\{l\_\{1\}\}^\{m\_\{1\}\}\(\\hat\{\\omega\}\)Y\_\{l\_\{2\}\}^\{m\_\{2\}\}\(\\hat\{\\omega\}\)Y\_\{l\_\{3\}\}^\{m\_\{3\}\}\(\\hat\{\\omega\}\)\\left\(\\sum\_\{i=1\}^\{n\}\\lambda\_\{i\}\\omega\_\{i\}^\{2\}\\right\)^\{k\}d\\hat\{\\omega\}\\,\.\(157\)These moments can be computed by expressing the integrand as a homogeneous polynomial in the Cartesian components ofω^\\hat\{\\omega\}and integrating term\-by\-term\.

#### Selection Rules in Three Dimensions

The presence of the anisotropic matrixAAbreaks the continuous rotational symmetry of the sphere, which relaxes the strict selection rules associated with the conservation of angular momentum\. To derive the new rules, we expand the anisotropic termg​\(ω^\)=exp⁡\(−12​ρ2​ω^T​A​ω^\)g\(\\hat\{\\omega\}\)=\\exp\(\-\\frac\{1\}\{2\}\\rho^\{2\}\\hat\{\\omega\}^\{T\}A\\hat\{\\omega\}\)in spherical harmonics:

g​\(ω^\)=∑L=0∞∑M=−LLCL,M​YLM​\(ω^\)\.g\(\\hat\{\\omega\}\)=\\sum\_\{L=0\}^\{\\infty\}\\sum\_\{M=\-L\}^\{L\}C\_\{L,M\}Y\_\{L\}^\{M\}\(\\hat\{\\omega\}\)\\,\.\(158\)The integralIIcan then be written as a sum over integrals of four harmonics:

I=∑L,MCL,M​∫𝕊2Yl1m1​\(ω^\)​Yl2m2​\(ω^\)​Yl3m3​\(ω^\)​YLM​\(ω^\)​𝑑Ω\.I=\\sum\_\{L,M\}C\_\{L,M\}\\int\_\{\\mathbb\{S\}^\{2\}\}Y\_\{l\_\{1\}\}^\{m\_\{1\}\}\(\\hat\{\\omega\}\)Y\_\{l\_\{2\}\}^\{m\_\{2\}\}\(\\hat\{\\omega\}\)Y\_\{l\_\{3\}\}^\{m\_\{3\}\}\(\\hat\{\\omega\}\)Y\_\{L\}^\{M\}\(\\hat\{\\omega\}\)d\\Omega\\,\.\(159\)The selection rules forIIare determined by the properties of the coefficientsCL,MC\_\{L,M\}and the four\-harmonic integrals\.

The selection rules for the magnetic quantum numbersmim\_\{i\}are determined by the azimuthal symmetry ofg​\(ω^\)g\(\\hat\{\\omega\}\)\. In 3D spherical coordinates, the argument of the exponential is:

ω^T​A​ω^\\displaystyle\\hat\{\\omega\}^\{T\}A\\hat\{\\omega\}=λ1​sin2⁡θ​cos2⁡ϕ\+λ2​sin2⁡θ​sin2⁡ϕ\+λ3​cos2⁡θ\\displaystyle=\\lambda\_\{1\}\\sin^\{2\}\\theta\\cos^\{2\}\\phi\+\\lambda\_\{2\}\\sin^\{2\}\\theta\\sin^\{2\}\\phi\+\\lambda\_\{3\}\\cos^\{2\}\\theta=\(λ1\+λ22​sin2⁡θ\+λ3​cos2⁡θ\)\\displaystyle=\\left\(\\frac\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\{2\}\\sin^\{2\}\\theta\+\\lambda\_\{3\}\\cos^\{2\}\\theta\\right\)\+\\displaystyle\+\(λ1−λ22​sin2⁡θ\)​cos⁡\(2​ϕ\)\.\\displaystyle\\left\(\\frac\{\\lambda\_\{1\}\-\\lambda\_\{2\}\}\{2\}\\sin^\{2\}\\theta\\right\)\\cos\(2\\phi\)\\,\.The dependence on the angleϕ\\phiis entirely contained within acos⁡\(2​ϕ\)\\cos\(2\\phi\)term\. Any function ofcos⁡\(2​ϕ\)\\cos\(2\\phi\)has a Fourier series containing only frequencies that are multiples of 2 \(i\.e\., terms likecos⁡\(2​k​ϕ\)\\cos\(2k\\phi\)\)\. The expansion coefficientsCL,M=∫g​\(ω^\)​YLM⁣∗​\(ω^\)​𝑑ΩC\_\{L,M\}=\\int g\(\\hat\{\\omega\}\)Y\_\{L\}^\{M\*\}\(\\hat\{\\omega\}\)d\\Omegainvolve an integral overϕ\\phiof the form:

∫02​πF​\(cos⁡\(2​ϕ\)\)​e−i​M​ϕ​𝑑ϕ\.\\int\_\{0\}^\{2\\pi\}F\(\\cos\(2\\phi\)\)e^\{\-iM\\phi\}d\\phi\\,\.This integral is non\-zero only ifMMmatches a frequency present in the function ofcos⁡\(2​ϕ\)\\cos\(2\\phi\), which meansMMmust be an even integer\. The integral of four harmonics is non\-zero only if the sum of the magnetic quantum numbers is zero:m1\+m2\+m3\+M=0m\_\{1\}\+m\_\{2\}\+m\_\{3\}\+M=0\. Since only terms with evenMMcontribute to the sum forII, we arrive at the relaxed selection rule:

m1\+m2\+m3=−M⟹m1\+m2\+m3=an even integer\.m\_\{1\}\+m\_\{2\}\+m\_\{3\}=\-M\\implies m\_\{1\}\+m\_\{2\}\+m\_\{3\}=\\text\{an even integer\.\}\(160\)
The triangle inequality for\(l1,l2,l3\)\(l\_\{1\},l\_\{2\},l\_\{3\}\)is lifted because the anisotropic field provides access to multiple angular momentum channelsLL\. The integral of four harmonics is non\-zero if\(l1,l2,l3,L\)\(l\_\{1\},l\_\{2\},l\_\{3\},L\)can be coupled to a total angular momentum of zero\. Sinceg​\(ω^\)g\(\\hat\{\\omega\}\)is an even function, its expansion only contains even values ofLL\.

#### Counterexample:

Consider the triplet\(l1,l2,l3\)=\(1,5,1\)\(l\_\{1\},l\_\{2\},l\_\{3\}\)=\(1,5,1\)\. This violates the standard triangle inequality\|l1−l2\|≤l3≤l1\+l2\|l\_\{1\}\-l\_\{2\}\|\\leq l\_\{3\}\\leq l\_\{1\}\+l\_\{2\}\(since1∉\[4,6\]1\\notin\[4,6\]\), so its Gaunt integral is zero\. For the anisotropic integral, we check if a coupling pathway exists through an evenLL\.

1. 1\.Couple\(l1,l2\)=\(1,5\)⟹\(l\_\{1\},l\_\{2\}\)=\(1,5\)\\impliesintermediate angular momentuml12∈\{4,5,6\}l\_\{12\}\\in\\\{4,5,6\\\}\.
2. 2\.We now seek an evenLLthat allows coupling\(l3,L\)=\(1,L\)\(l\_\{3\},L\)=\(1,L\)to one of thesel12l\_\{12\}values\. Let’s targetl12=4l\_\{12\}=4\.
3. 3\.The coupling rule requires\|l3−L\|≤l12≤l3\+L\|l\_\{3\}\-L\|\\leq l\_\{12\}\\leq l\_\{3\}\+L, which for our values becomes\|1−L\|≤4≤1\+L\|1\-L\|\\leq 4\\leq 1\+L\.
4. 4\.The inequality4≤1\+L4\\leq 1\+LimpliesL≥3L\\geq 3\. The inequality\|1−L\|≤4\|1\-L\|\\leq 4implies−3≤L≤5\-3\\leq L\\leq 5\.

The conditions requireLLto be in the range\[3,5\]\[3,5\]\. We can choose the even valueL=4L=4\. The expansion ofg​\(ω^\)g\(\\hat\{\\omega\}\)contains a non\-zeroC4,MC\_\{4,M\}term\. Therefore, a coupling pathway exists via theL=4L=4channel, and the integral can be non\-zero\. This demonstrates that the triangle inequality on\(l1,l2,l3\)\(l\_\{1\},l\_\{2\},l\_\{3\}\)is lifted\.

#### Summary:

Hence we see, we have managed to lift two of the most restrictive selection rule on the Clebsch\-Gordan symbol, now with the new interaction term, it has non\-zero contribution as long as the sum of bothmmandllare even, which is not such a restrictive condition\. By no means is this choice unique, we could in principle choose any interaction term that integrates to object with little selection rule, hence improving the performance of the rotationally invariant kernels\.

### E\.1results

Here we list some of the results we obtained in various dimensions, how this anisotropic term helps rescuing the performance of the Gaussian RBF kernel\.

![Refer to caption](https://arxiv.org/html/2606.00322v1/anisotropic.png)Figure 4:MSE comparison across dimension and RBF kernel bandwidth\.As shown in figure[4](https://arxiv.org/html/2606.00322#A5.F4), it does help in lower dimension when kernels themselves do not vanish, even this depends on the length scale of the RBF kernel\. For some choices the additional term even hurts the performance\. Although there is a lift in coefficients, due to the vanishing of the kernels themselves in high dimensions, this still fails inevitably\.

## Appendix FExperiments

We introduce several datasets and summerize the performance of the perturbative method against kernel IV and Deep IV161616All results on Deep IV was done using a two layered approximately 2k parameters neural net trained for 500 epochs for each regression stage\.\. We note that all tests are done repeatedly 1000 times and the average RMSE was recorded here as the result for comparison\.

### F\.1Primary experiment dataset

The primary dataset is constructed to match the experimental setup in\(Donhauseret al\.,[2021](https://arxiv.org/html/2606.00322#bib.bib15)\), with an instrumental variable twist\.

For a given sample sizennand dimensionality parameterβ\\beta, we setd=⌊nβ⌋d=\\lfloor n^\{\\beta\}\\rfloorand generate the following data\. First the instrumental variableZi∼𝒩​\(0,Id\)Z\_\{i\}\\sim\\mathcal\{N\}\(0,I\_\{d\}\), then the confounding term with correlation across dimensions:UX,i∼𝒩​\(0,Id\)U\_\{X,i\}\\sim\\mathcal\{N\}\(0,I\_\{d\}\)\. Endogenous variablesXiX\_\{i\}with nonlinear dependence onZZ:

Xi=0\.7​Zi\+0\.3​Zi2\+0\.2​UX,i,X\_\{i\}=0\.7Z\_\{i\}\+0\.3Z\_\{i\}^\{2\}\+0\.2U\_\{X,i\}\\,,\(161\)where the operations are element\-wise\. Ground truth function combines nonlinear effects from multiple dimensions:

g​\(Xi\)=0\.5​sin⁡\(∑j=110wj​Xi,j\)\+0\.5​e−0\.1​∑j=1120wj​Xi,j2,g\(X\_\{i\}\)=0\.5\\sin\\left\(\\sum\_\{j=1\}^\{10\}w\_\{j\}X\_\{i,j\}\\right\)\+0\.5e^\{\-0\.1\\sum\_\{j=11\}^\{20\}w\_\{j\}X\_\{i,j\}^\{2\}\}\\,,\(162\)wherewj=1/jw\_\{j\}=1/\\sqrt\{j\}are dimension\-specific weights\. Effectively the useful feature dimensions are the first2020dimensions\. Dimensions beyond2020are essentially noise dimensions, which creates a challenging scenario with high condition number for kernel ridge regression in NPIV, essentially kernel function needs to learn to ignore these noise\. Finally we add error term with endogeneity:

ϵY,i=0\.8⋅U¯X,i\+0\.2⋅νi,\\epsilon\_\{Y,i\}=0\.8\\cdot\\bar\{U\}\_\{X,i\}\+0\.2\\cdot\\nu\_\{i\}\\,,\(163\)whereU¯X,i=1d​∑j=1dUX,i,j\\bar\{U\}\_\{X,i\}=\\frac\{1\}\{d\}\\sum\_\{j=1\}^\{d\}U\_\{X,i,j\}is the average error across dimensions andνi∼𝒩​\(0,1\)\\nu\_\{i\}\\sim\\mathcal\{N\}\(0,1\)with finally outcome:

Yi=g​\(Xi\)\+ϵY,i\.Y\_\{i\}=g\(X\_\{i\}\)\+\\epsilon\_\{Y,i\}\\,\.\(164\)This construction creates a challenging high\-dimensional regression problem where the true function depends nonlinearly on multiple dimensions, endogeneity is present through correlated errors and the dimensionality parameterβ\\betacontrols how quickly dimension grows with sample size\.

We tested across various dimensionality regimes by setting dimensiond≈nβd\\approx n^\{\\beta\}forβ∈\{0\.3,0\.5,0\.7,0\.9,1\.1\}\\beta\\in\\\{0\.3,0\.5,0\.7,0\.9,1\.1\\\}\. This allowed us to observe how performance varies as the dimension increases relative to sample size\. We tested the following parameters to identify optimal values for base perturbation parameterγ∈\{0\.4,0\.6,0\.8,1\.0\}\\gamma\\in\\\{0\.4,0\.6,0\.8,1\.0\\\}\.

### F\.2Kernels

Gaussian RBF kernels are simply:

KRBF​\(x,x′\)=exp⁡\(−‖x−x′‖22​σ2\),K\_\{\\text\{RBF\}\}\(x,x^\{\\prime\}\)=\\exp\\left\(\-\\frac\{\\\|x\-x^\{\\prime\}\\\|^\{2\}\}\{2\\sigma^\{2\}\}\\right\)\\,,\(165\)whereσ\>0\\sigma\>0is the bandwidth parameter which we test in logarithmic scale0\.01,0\.1,1,100\.01,0\.1,1,10\. It is rotationally invariant and has been proved to suffer from curse of dimensionality in\(Donhauseret al\.,[2021](https://arxiv.org/html/2606.00322#bib.bib15)\)\.

The Fractional Brownian kernel is based on fractional Brownian motion covariance functions\. For pointsx,x′∈ℝdx,x^\{\\prime\}\\in\\mathbb\{R\}^\{d\}, the kernel reads:

Kq​\(x,x′\)=12​\(‖x‖q\+‖x′‖q−‖x−x′‖q\),K\_\{q\}\(x,x^\{\\prime\}\)=\\frac\{1\}\{2\}\\left\(\\\|x\\\|^\{q\}\+\\\|x^\{\\prime\}\\\|^\{q\}\-\\\|x\-x^\{\\prime\}\\\|^\{q\}\\right\)\\,,\(166\)where0<q≤20<q\\leq 2is the fractional parameter\. It explicitly breaks rotational invariance by choosing a reference point \(here we chose the origin\)\. This kernel belongs to the class of fractional Brownian kernels and corresponds to the covariance structure of fractional Brownian motion with Hurst parameterH=q/2H=q/2\. The kernel captures long\-range dependence whenq\>1q\>1and exhibits different smoothness properties compared to standard RBF kernels\(Sejdinovicet al\.,[2013](https://arxiv.org/html/2606.00322#bib.bib31); Rizzo and Székely,[2016](https://arxiv.org/html/2606.00322#bib.bib32)\)\. In our implementation, we useq=1q=1\(corresponding toH=0\.5H=0\.5, standard Brownian motion\)\. We shall refer to this kernel as the fractional Brownian kernel in the experiment section\.

### F\.3high dimensional results

Here we include the performance of our algorithm across different kernels, dimensionalities and renormalization strength in table[3](https://arxiv.org/html/2606.00322#A6.T3)and[4](https://arxiv.org/html/2606.00322#A6.T4)\.

Table 3:Order 0 RMSE \(Baseline\) Across All Kernels and Deep IVTable 4:Renormalized RMSE Across All Kernels
### F\.4Further datasets results

We evaluate our methods on several alternative NPIV datasets beyond the one from\(Donhauseret al\.,[2021](https://arxiv.org/html/2606.00322#bib.bib15)\)described in appendix[F\.1](https://arxiv.org/html/2606.00322#A6.SS1)\. Each dataset is generated withn=80n=80samples with results summerized in table[5](https://arxiv.org/html/2606.00322#A6.T5)\.

#### Newey\-Powell Dataset\.

A classic NPIV benchmark following\(Newey and Powell,[2003](https://arxiv.org/html/2606.00322#bib.bib21)\):

Zi\\displaystyle Z\_\{i\}∼Uniform​\(−3,3\),Vi∼𝒩​\(0,1\)\\displaystyle\\sim\\text\{Uniform\}\(\-3,3\),\\quad V\_\{i\}\\sim\\mathcal\{N\}\(0,1\)\(167\)Xi\\displaystyle X\_\{i\}=0\.6​Zi\+0\.4​Vi\+0\.2​νi,νi∼𝒩​\(0,1\)\\displaystyle=0\.6Z\_\{i\}\+0\.4V\_\{i\}\+0\.2\\nu\_\{i\},\\quad\\nu\_\{i\}\\sim\\mathcal\{N\}\(0,1\)\(168\)g​\(X\)\\displaystyle g\(X\)=sin⁡\(X\)\+0\.5​X\\displaystyle=\\sin\(X\)\+0\.5X\(169\)Yi\\displaystyle Y\_\{i\}=g​\(Xi\)\+0\.5​Vi\+0\.5​ϵi,ϵi∼𝒩​\(0,1\)\\displaystyle=g\(X\_\{i\}\)\+0\.5V\_\{i\}\+0\.5\\epsilon\_\{i\},\\quad\\epsilon\_\{i\}\\sim\\mathcal\{N\}\(0,1\)\(170\)

#### Weak/Strong Instrument Dataset\.

Tests robustness to instrument strengthρ∈\{0\.1,0\.3,0\.7\}\\rho\\in\\\{0\.1,0\.3,0\.7\\\}withd=5d=5:

Zi,Ui\\displaystyle Z\_\{i\},U\_\{i\}∼𝒩​\(0,Id\)\\displaystyle\\sim\\mathcal\{N\}\(0,I\_\{d\}\)\(171\)Xi\\displaystyle X\_\{i\}=ρ​Zi\+\(1−ρ\)​Ui\+0\.1​νi,νi∼𝒩​\(0,Id\)\\displaystyle=\\rho Z\_\{i\}\+\(1\-\\rho\)U\_\{i\}\+0\.1\\nu\_\{i\},\\quad\\nu\_\{i\}\\sim\\mathcal\{N\}\(0,I\_\{d\}\)\(172\)g​\(X\)\\displaystyle g\(X\)=0\.5​sin⁡\(∑j=1dwj​Xj\)\+0\.3​cos⁡\(∑j=12wj​Xj\)\\displaystyle=0\.5\\sin\\left\(\\sum\_\{j=1\}^\{d\}w\_\{j\}X\_\{j\}\\right\)\+0\.3\\cos\\left\(\\sum\_\{j=1\}^\{2\}w\_\{j\}X\_\{j\}\\right\)\(173\)Yi\\displaystyle Y\_\{i\}=g​\(Xi\)\+0\.6​U¯i\+0\.4​ϵi\\displaystyle=g\(X\_\{i\}\)\+0\.6\\bar\{U\}\_\{i\}\+0\.4\\epsilon\_\{i\}\(174\)wherewj=1/jw\_\{j\}=1/\\sqrt\{j\}andU¯i=d−1​∑j=1dUi​j\\bar\{U\}\_\{i\}=d^\{\-1\}\\sum\_\{j=1\}^\{d\}U\_\{ij\}\.

#### Heteroscedastic Dataset\.

Features variance that depends onXXwithd=10d=10:

Zi,Ui\\displaystyle Z\_\{i\},U\_\{i\}∼𝒩​\(0,Id\),Xi=0\.7​Zi\+0\.3​Ui\\displaystyle\\sim\\mathcal\{N\}\(0,I\_\{d\}\),\\quad X\_\{i\}=0\.7Z\_\{i\}\+0\.3U\_\{i\}\(175\)g​\(X\)\\displaystyle g\(X\)=sin⁡\(∑j=15wj​Xj\)\+0\.3​exp⁡\(−0\.1​∑j=6dXj2\)\\displaystyle=\\sin\\left\(\\sum\_\{j=1\}^\{5\}w\_\{j\}X\_\{j\}\\right\)\+0\.3\\exp\\left\(\-0\.1\\sum\_\{j=6\}^\{d\}X\_\{j\}^\{2\}\\right\)\(176\)σ2​\(X\)\\displaystyle\\sigma^\{2\}\(X\)=0\.5\+0\.5​\|X¯\|,Yi=g​\(Xi\)\+0\.5​U¯i\+σ​\(Xi\)​ϵi\\displaystyle=0\.5\+0\.5\|\\bar\{X\}\|,\\quad Y\_\{i\}=g\(X\_\{i\}\)\+0\.5\\bar\{U\}\_\{i\}\+\\sigma\(X\_\{i\}\)\\epsilon\_\{i\}\(177\)

#### Nonlinear Instrument Dataset\.

Tests nonlinear first\-stage relationships withd=10d=10:

Zi,Ui\\displaystyle Z\_\{i\},U\_\{i\}∼𝒩​\(0,Id\)\\displaystyle\\sim\\mathcal\{N\}\(0,I\_\{d\}\)\(178\)Xi\\displaystyle X\_\{i\}=0\.5​Zi\+0\.3​Zi2\+0\.2​sin⁡\(Zi\)\+0\.2​Ui\\displaystyle=0\.5Z\_\{i\}\+0\.3Z\_\{i\}^\{2\}\+0\.2\\sin\(Z\_\{i\}\)\+0\.2U\_\{i\}\(179\)g​\(X\)\\displaystyle g\(X\)=0\.4​tanh⁡\(∑j=15wj​Xj\)\+0\.3​cos⁡\(∑j=610wj​Xj\)\\displaystyle=0\.4\\tanh\\left\(\\sum\_\{j=1\}^\{5\}w\_\{j\}X\_\{j\}\\right\)\+0\.3\\cos\\left\(\\sum\_\{j=6\}^\{10\}w\_\{j\}X\_\{j\}\\right\)\(180\)Yi\\displaystyle Y\_\{i\}=g​\(Xi\)\+0\.6​U¯i\+0\.4​ϵi\\displaystyle=g\(X\_\{i\}\)\+0\.6\\bar\{U\}\_\{i\}\+0\.4\\epsilon\_\{i\}\(181\)

#### Sparse Signal Dataset\.

High\-dimensional setting \(d=50d=50\) with onlys=5s=5active dimensions:

Zi,Ui\\displaystyle Z\_\{i\},U\_\{i\}∼𝒩​\(0,Id\),Xi=0\.7​Zi\+0\.3​Ui\\displaystyle\\sim\\mathcal\{N\}\(0,I\_\{d\}\),\\quad X\_\{i\}=0\.7Z\_\{i\}\+0\.3U\_\{i\}\(182\)g​\(X\)\\displaystyle g\(X\)=sin⁡\(1s​∑j=1sXj\)\+0\.3​cos⁡\(2​X1\)\\displaystyle=\\sin\\left\(\\frac\{1\}\{\\sqrt\{s\}\}\\sum\_\{j=1\}^\{s\}X\_\{j\}\\right\)\+0\.3\\cos\(2X\_\{1\}\)\(183\)Yi\\displaystyle Y\_\{i\}=g​\(Xi\)\+0\.5​U¯1:s,i\+0\.5​ϵi\\displaystyle=g\(X\_\{i\}\)\+0\.5\\bar\{U\}\_\{1:s,i\}\+0\.5\\epsilon\_\{i\}\(184\)whereU¯1:s,i=s−1​∑j=1sUi​j\\bar\{U\}\_\{1:s,i\}=s^\{\-1\}\\sum\_\{j=1\}^\{s\}U\_\{ij\}averages only the active dimensions\.

Table 5:RMSE Comparison on Alternative NPIV Datasets- •RMSE =MSE\\sqrt\{\\text\{MSE\}\}; “Order 0” is standard kernel ridge IV; “Best Pert\.” is the best result acrossγ∈\{0\.4,0\.6,0\.8,0\.99\}\\gamma\\in\\\{0\.4,0\.6,0\.8,0\.99\\\}and kernels∈\\in\{FB, RBF\}\. Bold indicates best method per dataset\.

In table[5](https://arxiv.org/html/2606.00322#A6.T5), we note the improvement of the perturbative method on order0base kernel IV across all these datasets\. Although DeepIV performs well on most datasets, the perturbative method shows competitive or superior performance in two notable cases: \(1\) the weak instrument setting \(ρ=0\.1\\rho=0\.1\), where both neural network training and kernel IV struggle but the perturbative correction provides meaningful improvement; and \(2\) the sparse high\-dimensional setting \(d=50d=50\), where the perturbative approach’s ability to enhance kernel expressivity in specific eigendirections proves advantageous\. The Order 0 baseline \(standard kernel ridge regression\) performs poorly across all datasets, highlighting the importance of addressing the ill\-posedness inherent in NPIV problems\.

### F\.5Weak Instrument Analysis

A fundamental challenge in instrumental variable estimation is theweak instrument problem, which occurs when instrumentsZZhave low correlation with the endogenous variableXX\. Weak instruments lead to biased estimates, inflated standard errors, and unreliable inference\([J\. Bound, D\. A\. Jaeger, and R\. M\. Baker \(1995\)](https://arxiv.org/html/2606.00322#bib.bib6);[1](https://arxiv.org/html/2606.00322#bib.bib7)\)\. The first\-stageR2R^\{2\}—measuring the proportion of variance inXXexplained byZZ—serves as a diagnostic: values below 0\.1 are conventionally considered indicative of weak instruments\(Stocket al\.,[2002a](https://arxiv.org/html/2606.00322#bib.bib5)\)\.

We systematically vary the instrument strength parameterρ∈\[0\.05,0\.9\]\\rho\\in\[0\.05,0\.9\]in our weak instrument dataset \(Section[F\.4](https://arxiv.org/html/2606.00322#A6.SS4)\) to examine how different NPIV methods degrade as instruments weaken\. The first\-stage relationship isX=ρ​Z\+\(1−ρ\)​U\+noiseX=\\rho Z\+\(1\-\\rho\)U\+\\text\{noise\}, soρ\\rhodirectly controls instrument relevance\. We summerize the results in table[6](https://arxiv.org/html/2606.00322#A6.T6)\.

Table 6:RMSE Across Instrument Strength Levels- •RMSE =MSE\\sqrt\{\\text\{MSE\}\}\. First\-stageR2<0\.1R^\{2\}<0\.1indicates weak instruments \(shaded regionρ≤0\.2\\rho\\leq 0\.2\)\. Bold indicates best method\.

As shown in table[6](https://arxiv.org/html/2606.00322#A6.T6), perhaps unsurprisingly, DeepIV dominates across all instrument strengths\. The neural network’s flexible function approximation proves advantageous even with weak instruments, though its RMSE increases from 0\.312 \(ρ=0\.9\\rho=0\.9\) to 0\.882 \(ρ=0\.1\\rho=0\.1\)—a 2\.8×\\timesdegradation\. Compared to order0kernel IV with regularization, perturbative corrections helps to reduce Order 0 RMSE by 18–37% across settings, demonstrating the value of higher\-order corrections even in challenging regimes\. However, it cannot match DeepIV’s performance\.

### F\.6Sensitivity Analysis

The perturbative renormalization framework introduces several hyperparameters\. We conduct a systematic sensitivity analysis on the Donhauser dataset withβ=1\.0\\beta=1\.0\(d=60d=60,n=60n=60\) to understand their effects and identify robust default values\.

We study parameters across the range, base perturbation coupling constantγ∈\{0\.1,0\.2,…,0\.9,0\.95,0\.99\}\\gamma\\in\\\{0\.1,0\.2,\\ldots,0\.9,0\.95,0\.99\\\}, Tikhonov regularization strengthλ∈\{0\.01,0\.05,0\.1,0\.3,0\.6,1\.0,2\.0\}\\lambda\\in\\\{0\.01,0\.05,0\.1,0\.3,0\.6,1\.0,2\.0\\\}, maximum perturbative orderNmax∈\{2,3,4,5,6,7\}N\_\{\\max\}\\in\\\{2,3,4,5,6,7\\\}, summerized in table[7](https://arxiv.org/html/2606.00322#A6.T7)\.

Table 7:Sensitivity Analysis: Optimal Parameters and RMSE Ranges- •Sensitivity rated as: High \(\>\>50% RMSE variation\), Medium \(20–50%\), Low \(<<20%\)\.

As shown in table[7](https://arxiv.org/html/2606.00322#A6.T7), the coupling constantγ\\gammais the most sensitive parameter, as is confirmed in the result in table[4](https://arxiv.org/html/2606.00322#A6.T4)\. The optimal valueγ∗≈0\.6\\gamma^\{\*\}\\approx 0\.6reflects a bias\-variance tradeoff: smallγ\\gammaunderweights higher\-order corrections \(high bias\), whileγ→1\\gamma\\to 1amplifies potentially divergent terms \(high variance\)\.

Moderate sensitivity with Regularizationλ\\lambdawith optimalλ∗=0\.6\\lambda^\{\*\}=0\.6\. Too little regularization \(λ<0\.1\\lambda<0\.1\) allows ill\-conditioned kernel inversions; too much \(λ\>1\\lambda\>1\) over\-smooths the solution\.

Performance is stable forNmax≥4N\_\{\\max\}\\geq 4, indicating that most information is captured in the first few orders\. This aligns with asymptotic theory: optimal truncation typically occurs atN∗=O​\(γ−1\)N^\{\*\}=O\(\\gamma^\{\-1\}\)\.

#### Recommended Defaults\.

Based on this analysis, we recommend:γ=0\.6\\gamma=0\.6,λ=0\.6\\lambda=0\.6,Nmax=5N\_\{\\max\}=5\. These values perform well across the datasets in Sections[6](https://arxiv.org/html/2606.00322#S6)–[F\.4](https://arxiv.org/html/2606.00322#A6.SS4)without dataset\-specific tuning\.

### F\.7Larger sample experiments

Here we include a further set of experiments with larger sample size\(n=120\)\(n=120\)insteadn=80n=80\. We summerize the results here in table[8](https://arxiv.org/html/2606.00322#A6.T8)\.

Table 8:RMSE Comparison on Alternative NPIV Datasets- •RMSE =MSE\\sqrt\{\\text\{MSE\}\}; “Order 0” is standard kernel ridge IV; “Best Pert\.” is the best result acrossγ∈\{0\.4,0\.6,0\.8,0\.99\}\\gamma\\in\\\{0\.4,0\.6,0\.8,0\.99\\\}and kernels∈\\in\{FB, RBF\}\. Bold indicates best method per dataset\.