Stabilized Higher-Order Influence Functions: Statistical Theory of a Class of Bilinear Forms

bibunit[plainnat] \affil[1]{School of Mathematical Sciences, Shanghai Jiao Tong University} \affil[2]{Department of Statistics, University of Virginia} \affil[3]{Department of Statistics and Data Science, Tsinghua University} \affil[4]{Institute of Natural Sciences, MOE--LSC, CMA--Shanghai, SJTU--Yale Joint Center for Biostatistics and Data Science, Shanghai Jiao Tong University} \begin{abstract} Higher-order influence functions, introduced in a series of articles robins2008higher, robins2009quadratic, van2014higher, robins2016technical, robins2023minimax, liu2017semiparametric, are a unified framework for constructing rate-optimal point estimates of a class of statistical functionals under various complexity-reducing assumptions on the posited statistical model that generates the observed data. Although higher-order (influence functions) estimators are theoretically appealing, they have very limited practical uptake compared to their first-order counterparts. The original higher-order estimators proposed in robins2008higher and robins2017minimax involve nonparametric density estimation of multi-dimensional covariates, a highly nontrivial statistical and computational problem on its own. The density estimator is, in turn, used in the evaluation of the inverse population Gram matrix $\Omega$ of a set of $k$-dimensional basis transformations of covariates. There, $k$ is allowed to be as large as $o (n^{2})$. To partially address this potential shortcoming, liu2017semiparametric restrict $k$ to $o (n)$ and instead estimates $\Omega$ directly using the inverse sample Gram matrix estimator, but computed from an independent sample often obtained by sample-splitting. liu2017semiparametric refer to this alternative estimator as the empirical higher-order estimator. Although the empirical higher-order estimator bypasses density estimation, it suffers from numerical instability due to potentially inverting a large-dimensional sample Gram matrix. In this article, for a class of bilinear forms/functionals that often appear in substantive fields such as economics, epidemiology, and clinical medicine, we propose a new stabilized higher-order estimator without sample splitting, which exhibits more stable finite-sample performance compared to the empirical higher-order estimator. More importantly, we prove that this new class of higher-order estimators enjoys similar statistical guarantees to those of liu2017semiparametric. \end{abstract} { Keywords: Causal Inference, Functional Estimation, Higher-Order Influence Functions, M\"{o}bius Inversion, Enumerative Combinatorics} \onehalfspacing \allowdisplaybreaks \section{Introduction} One of the unique features of modern statistics, which distinguishes itself from other related areas such as machine learning or AI, is the enormous interest in learning about smooth (statistical) functionals of the possibly infinite-dimensional probabilistic model that generates the observed data, instead of the model itself bickel1988estimating, ritov1990achieving, van1991differentiable, bickel1998efficient, robins2008higher. In this article, a functional is a mapping $\psi: {\mathcal{P}} \to {\mathbb{R}}$, from the underlying statistical model, denoted by ${\mathcal{P}}$, to the reals ${\mathbb{R}}$. A statistical model ${\mathcal{P}}$ contains all possible observed-data-generating probability distributions, posited by a statistician. A functional $\psi$ is said to be smooth in the sense of van1991differentiable, that is, the pathwise derivative of $\psi ({\mathbb{P}})$, along any parametric submodel $\{{\mathbb{P}}_{t}: {\mathbb{P}}_{0} = {\mathbb{P}}\} \subseteq {\mathcal{P}}$, allows the following representation: \begin{align*} \left. \frac{{\mathrm d}}{{\mathrm d} t} \right|_{t = 0} \psi ({\mathbb{P}}_{t}) = {\mathbb{E}} \{\mathsf{IF}_{\psi} \cdot g (O)\}, \end{align*} where $g$ is the score function associated with the parametric submodel ${\mathbb{P}}_{t}$, and $\mathsf{IF}_{\psi} \equiv \mathsf{IF}_{\psi, {\mathbb{P}}}$ is the (first-order) efficient influence function (IF) (or canonical gradient) of $\psi$ locally at ${\mathbb{P}} \in {\mathcal{P}}$ fisher2021visually, hines2022demystifying. It is also required that $\mathsf{IF}_{\psi}$ has mean zero at ${\mathbb{P}}$. Examples of smooth functionals abound: in causal inference, common target parameters of interest, such as the average treatment effect, the average treatment effect on the treated, and the quantile treatment effect, are all smooth functionals under standard causal identification conditions (consistency, positivity, and ignorability) robins1994estimation, hahn1998role, hahn2004functional, van2006targeted, abadie2018econometric; in (conditional) independence testing, dependence measures such as the generalized covariance measure shah2020hardness, niu2024reconciling and $f$-divergence kandasamy2015nonparametric, are also smooth functionals. This article specifically tackles the problem of constructing “good” estimators for smooth functionals, which we abbreviate as the problem of functional estimation. A natural attempt to estimate $\psi$ is to start with the “plug-in” estimator $\widehat{\psi}_{0} = \psi (\widehat{{\mathbb{P}}})$, where $\widehat{{\mathbb{P}}}$ is some estimator of ${\mathbb{P}}$. However, a common theme in the functional estimation literature tells us that the plug-in estimator $\widehat{\psi}_{0}$ has a sub-optimal convergence rate in many settings robins2009semiparametric, balakrishnan2026fundamental. The sub-optimality of the plug-in estimator is often resulting from its large bias. A popular (and almost dominating) paradigm in the current statistics literature is to use the IF of $\psi$, $\mathsf{IF}_{\psi}$, to de-bias the plug-in estimator $\widehat{\psi}_{0}$ scharfstein1999adjusting, van2006targeted, chernozhukov2018double, ray2020semiparametric, breunig2025double. We refer to these debiased estimators based solely on $\mathsf{IF}_{\psi}$ as first-order estimators liu2026asymptotic, which include popular methods in applications such as double machine learning/Neyman orthogonal scores chernozhukov2018double and targeted maximum likelihood estimation (TMLE) van2006targeted. In many settings, however, first-order estimators are still sub-optimal in terms of convergence rates liu2024assumption, bonvini2024doubly, liu2023root. To resolve the potential sub-optimality of $\widehat{\psi}_{1}$, building upon von Mises functional expansions and higher-order scores mises1947asymptotic, pfanzagl1983asymptotic, pfanzagl1990estimation, pfanzagl2011parametric, small1989projection, waterman1996projected, bobkov2024fisher, villani2025fisher, robins2008higher, robins2009quadratic, robins2016technical develop a general framework called higher-order influence functions (HOIFs) that generalize the concept of IF from first-order to higher-orders, for constructing (nearly) rate-optimal estimators in various settings. We also refer to bonhomme2026higher for related development in higher-order Neyman orthogonal scores and to diaz2016second, van2021higher for related development in higher-order TMLE (HOTMLE). TMLE-related methodologies generally enjoy favorable finite sample performance. The HOIF framework has also been used to construct estimators in related infinite-dimensional problems kennedy2024minimax, bonvini2022fast and to understand the statistical properties of irregular estimators of causal parameters bonvini2024doubly. One key insight of robins2008higher, robins2016technical, robins2023minimax is to find an approximation of the target functional $\psi$ by a particular bilinear form $\widetilde{\psi}_{k} = \mu^{\top} \Sigma^{-1} \eta$, where $\Sigma = {\mathbb{E}} (X X^{\top})$ is the $k \times k$ population Gram matrix of some random vector $X$, and $\mu$ and $\eta$ are two $k$-dimensional vectors that can be written respectively as $\mu = {\mathbb{E}} (X A)$ and $\eta = {\mathbb{E}} (X Y)$ for some random variables $A$ and $Y$ (see Section (ref) for details). Once this step is accomplished, HOIFs offer a unified scheme of constructing rate-optimal estimators of the bilinear form $\widetilde{\psi}_{k}$, and the resulting estimators are higher-order $U$-statistics. Fortunately, many of the aforementioned examples of smooth functionals indeed admit such a bilinear form approximation; again, see Section (ref) for concrete examples (Examples (ref)--(ref)). As will be clear in Section (ref), in this article, we will directly take the bilinear form $\widetilde{\psi}_{k}$ as the target parameter $\psi$ without worrying about the bias due to this bilinear approximation. The HOIF estimators proposed in robins2008higher, robins2016technical, robins2023minimax allow the dimension $k$ to be as large as of order $o (n^{2})$, but require a nonparametric density estimation step when estimating $\Sigma$ from data. Given the difficulty of nonparametric density estimation even in moderate dimensions, the original HOIF estimators have not been routinely deployed in practice. When the dimension $k$ is of order $o (n)$ so $\Sigma^{-1}$ can be consistently estimated by the inverse of the sample Gram matrix $\widehat{\Sigma}^{-1}$, liu2017semiparametric proposed the so-called empirical HOIF estimators, simply estimating $\Sigma^{-1}$ by $\widehat{\Sigma}^{-1}$ from a separate sample independent of the main sample used to estimate $\psi$. To our knowledge, the empirical HOIF estimator remains the only $\sqrt{n}$-consistent and asymptotic normal ($\sqrt{n}$-CAN) estimator of $\psi$ when $k = o (n)$, without imposing any assumption on the covariate density. zhang2026higher extend both versions of HOIF estimators to parameters defined implicitly via $Z$/$M$-estimation problems, such as quantile treatment effects and expected shortfalls. More recently, newey2018cross initiated the research program on constructing estimators motivated by but much simpler than HOIFs, with follow-up work in various directions kennedy2023towards, mcgrath2026nuisance, mcclean2026double. Finally, we also mention in passing that similar bias correction ideas have also been independently developed in the econometric and general mathematical statistics literature newey2004twicing, cattaneo2018kernel, cattaneo2018inference, cattaneo2019two, breunig2024adaptive, cavaliere2024bootstrap, koltchinskii2022bootstrap, koltchinskii2025estimation. Although empirical HOIF estimators neither estimate nor impose any complexity-reducing assumptions on the density of $X$, inverting the sample Gram matrix $\widehat{\Sigma}$ may easily lead to numerical instability when $k$ is relatively large compared to $n$. This potential instability has been documented in the simulation studies conducted in liu2020nearly, liu2017semiparametric, liu2024assumption, zhang2026higher, being a primary reason for the limited practical uptake of empirical HOIF estimators. However, it is less well known that liu2020nearly also proposed alternative empirical HOIF estimators (at orders $2$ and $3$, in retrospect) that still estimate the population Gram matrix $\Sigma$ by its sample analog $\widehat{\Sigma}$ but from the same sample used to compute the final $U$-statistic estimator. Since sample splitting is not used, liu2020nearly did not prove that this new alternative HOIF estimator works in theory; interestingly, for the same reason, these alternative HOIF estimators exhibit much improved finite-sample performance compared to the original ones proposed in liu2017semiparametric, in particular in terms of their numerical stability, even allowing practitioners to choose $k$ very close to $n$ (see Remark (ref) for further explanations). For the sake of completeness, this is demonstrated in Figure (ref) in Section (ref), which display the numerical results of a simple simulation study, the setup of which is described in Appendix (ref). \subsection{Our contributions} The main contribution of this article is to offer theoretical guarantees for the aforementioned alternative HOIF estimators, which we refer to as numerically stable HOIF estimators. The main technical difficulty arises from the dependence of the $U$-statistic kernel on the entire sample through $\widehat{\Sigma}^{-1}$ when sample splitting is not employed. To overcome this challenge, we have to deviate from the analysis strategy for the original empirical HOIF estimators taken in liu2017semiparametric and instead perform a more meticulous analysis that involves various complex expansions and nontrivial counting stanley2011enumerative. We obtain results similar to those for the empirical HOIF estimators of liu2017semiparametric, in the sense that the new HOIF estimators are also $\sqrt{n}$-CAN for the bilinear forms $\psi$, as long as $k = o (n)$ without any further complexity-reducing assumptions on the density of $X$. Specifically, we bring in tools from enumerative combinatorics and graph theory lauritzen1996graphical, chen2010mobius, stanley2011enumerative, shpitser2011efficient, richardson2023nested to prove the bias and variance bounds for this new class of HOIF estimators. These tools were recently exploited in chen2025computing to design efficient algorithms for the exact computation of higher-order $U$-statistics. In addition, schafer2026mobius also uses these tools to give a new combinatorial interpretation of the iterative bootstrap procedure. However, to our knowledge, these tools have not been used to establish statistical properties for estimators that involve higher-order $U$-statistics. The second article of this series will further delineate the connection between our new stabilized HOIF estimators and various other higher-order bias correction schemes in mathematical statistics at large, together with a more comprehensive set of simulation studies to benchmark the finite-sample performance of different higher-order bias correction methods. \subsection{Notation} Throughout the article, ${\mathbb{U}}_{n,j}$ denotes the $j$-th order $U$-statistic operator: for any measurable $h: {\mathcal{O}}_{1} \times \cdots \times {\mathcal{O}}_{j} \to {\mathbb{R}}$, \begin{align*} {\mathbb{U}}_{n, j} \{h (O_{1}, \cdots, O_{j})\} \coloneqq \frac{(n - j)!}{n!} \sum_{1 \leq i_{1} \neq \cdots \neq i_{j} \leq n} h (O_{i_{1}}, \cdots, O_{i_{j}}). \end{align*} We reserve $\Sigma$ and $\widehat{\Sigma}$ for the population and sample Gram matrices of $X$, and write $\Omega \coloneqq \Sigma^{-1}$ and $\widehat{\Omega} \coloneqq \widehat{\Sigma}^{-1}$ for their inverses whenever these exist ($\widehat{\Sigma}$ being invertible almost surely under our assumptions). The identity matrix is denoted by $I$. For a random variable $W$ and $p \geq 1$, $\Vert W \Vert_{p} \coloneqq \{{\mathbb{E}} (|W|^{p})\}^{1/p}$ denotes the $L^{p}({\mathbb{P}})$-norm of $W$. To lighten notation, for any sample-index subset $S \subseteq [n]$, we write $O_{S} \coloneqq \{O_{i} : i \in S\}$, and given any positive integer $\ell$, we let $[\ell] \coloneqq \{1, \cdots, \ell\}$. We write $\{i_1, \cdots, i_k\}$ as a set including elements $i_1, \cdots, i_k$ and write $(i_1, \cdots, i_k)$ as an ordered tuple, in which all elements are distinct and are assigned a particular ordering (mostly a canonical ordering). \subsection{Organizations} The remainder of this article is structured as follows. Section (ref) sets the stage by describing the problem setting, regularity assumptions, and providing a brief review of the empirical HOIF estimator of liu2017semiparametric. In Section (ref), we present the main result of this article, in which we first introduce the new numerically stable HOIF estimators and then characterize their bias, variance, and asymptotic distribution. The theoretical results are all encapsulated in Theorem (ref), the main theorem in our article. Section (ref) provides a proof sketch of Theorem (ref), with technical details deferred to the Appendix. Section (ref) concludes the article with a discussion of future topics. \section{Problem Setting and A Brief Review of Existing HOIF Estimators} Let $O \coloneqq (X, A, Y)$ denote a triple of the observed random vector, where $X \in {\mathcal{X}} \subseteq {\mathbb{R}}^{k}$ is a $k$-dimensional vector, $A \in {\mathcal{A}} \subset {\mathbb{R}}$ and $Y \in {\mathcal{Y}} \subseteq {\mathbb{R}}$ denote some outcomes of interest. We assume access to $n$ i.i.d. observations ${\mathcal{D}} \coloneqq (O_{1}, \cdots, O_{n})$, drawn from a common data-generating distribution ${\mathbb{P}} \in {\mathcal{P}}$, where ${\mathcal{P}}$ denotes the statistical model restricted by the following regularity conditions. \begin{assumption} The distribution of $X$ satisfies the following: \begin{align} & {\mathbb{E}} (X^{\top} X) = O (k), \\ & \Vert X^{\top} X \Vert_{\infty} = O (k), \end{align} and the eigenvalues of $\Sigma$ are strictly bounded away from $0$ and $\infty$. \end{assumption} In addition, in this article, we restrict to the case $k = o (n)$. But we will state the more precise condition on $k$ in the statement of related theoretical claims. We also need to impose the following $L_{\infty}$-stability assumption on the projection on the span of $X$, as commonly done in previous work on HOIFs robins2008higher, robins2016technical, robins2017minimax, robins2023minimax, liu2017semiparametric, liu2024assumption. \begin{assumption} For every bounded measurable function $h: {\mathcal{X}} \to {\mathbb{R}}$, define the following integral operator: \begin{equation*} (\Pi h) (x) \coloneqq x^{\top} \Omega {\mathbb{E}} \{X h(X)\}, \ \forall \ x \in {\mathcal{X}}. \end{equation*} We assume that $\Pi$ is uniformly bounded as an operator on $L_\infty ({\mathcal{X}})$: there exists a strictly bounded constant $C_{\Pi} < \infty$, independent of $k$ and $n$, such that \begin{equation} \|\Pi h\|_{\infty} \leq C_\Pi \|h\|_{\infty}. \end{equation} \end{assumption} Finally, for convenience, we further impose the following condition on $A$ and $Y$. \begin{assumption} Both $A$ and $Y$ are bounded almost surely. \end{assumption} \begin{remark} The above assumptions are made for technical convenience. For example, if we relax Assumption (ref) from boundedness to light-tailed assumptions, we need to further develop exponential and moment inequalities for higher-order $U$-statistics with unbounded kernels, which is an important research topic in applied probability on its own chakrabortty2025tail. \end{remark} For ease of exposition, throughout the article we consider the following functional of ${\mathbb{P}}$ as the target parameter: \begin{equation} \psi \equiv \psi ({\mathbb{P}}) \coloneqq \mu^{\top} \Omega \eta, \quad where \mu \coloneqq {\mathbb{E}} (X A), \eta \coloneqq {\mathbb{E}} (X Y), \Sigma \coloneqq {\mathbb{E}} (X X^{\top}) and \Omega \coloneqq \Sigma^{-1}. \end{equation} Although $\psi$ takes a very simple bilinear form, it encapsulates many substantively important smooth functionals that appear in the literature. We use several examples to demonstrate the ubiquity of $\psi$. \begin{example}[Quadratic functional of a density] Suppose that $Y \sim {\mathbb{P}}$ with $p$ being the probability density function of ${\mathbb{P}}$, the target functional is $\psi = \int_{{\mathcal{Y}}} p (y)^{2} {\mathrm d} y$, and $p$ can be represented as a linear combination of $\bar{\phi}$, assumed to be orthonormal with respect to the Lebesgue measure over ${\mathcal{Y}}$. Thus, there exists $\eta \in {\mathbb{R}}^{k}$ such that $p (\cdot) = \eta^{\top} \bar{\phi} (\cdot)$. We identify $X \coloneqq \bar{\phi} (Y)$ and $A \equiv Y$. Then given $O = (X, A, Y)$, $\psi = \eta^{\top} \eta$ with $\Sigma = I$. This quadratic functional of a density is one of the most well-studied smooth functionals in the statistics literature bickel1988estimating. \end{example} \begin{example}[Signal-to-noise ratio] Suppose that $(X, Y) \sim {\mathbb{P}}$, and the target functional is $\psi = {\mathbb{E}} \{b (X)^{2}\}$ where $b (\cdot) \coloneqq {\mathbb{E}} (Y \mid X = \cdot)$. We further assume that $b (\cdot) = \beta^{\top} (\cdot)$ for some $\beta \in {\mathbb{R}}^{k}$. Then $\psi = \beta^{\top} \Sigma \beta = \eta^{\top} \Omega \eta$, with $\beta = \Omega \eta$. Similar parameters have been extensively studied in the past decade in the context of high-dimensional (generalized) linear models verzelen2018adaptive, chen2024method. \end{example} \begin{example}[Treatment-specific counterfactual mean] Suppose that $(Z, A, Y) \sim {\mathbb{P}}$ constitutes the observed data of an unconfounded observational study, in which $A$ is the binary treatment variable, $Y$ is an outcome of interest, and $Z$ is the baseline covariates that contain all confounders between $A$ and $Y$. The target parameter is the treatment-specific counterfactual mean $\psi = {\mathbb{E}} Y (1) = {\mathbb{E}} \{A a (Z) Y\} = {\mathbb{E}} \{b (Z)\}$, where $a (\cdot) \coloneqq {\mathbb{E}}^{-1} (A \mid Z = \cdot)$ and $b (\cdot) \coloneqq {\mathbb{E}} (Y \mid Z = \cdot, A = 1)$. Let $X = A \bar{\phi} (Z)$. As shown in robins2007comment, liu2017semiparametric, bruns2026augmented, if we posit that $a (\cdot) \equiv \alpha^{\top} \bar{\phi} (\cdot)$ and $b (\cdot) \equiv \beta^{\top} \bar{\phi} (\cdot)$, where $\alpha, \beta \in {\mathbb{R}}^{k}$, then $\psi = \alpha^{\top} \Sigma \beta = \mu^{\top} \Omega \eta$, where $\alpha = \Omega \mu$ and $\beta = \Omega \eta$ with $\mu = {\mathbb{E}} (X a (X))$ and $\eta = {\mathbb{E}} (X A Y)$. For implicitly defined parameters such as the quantile treatment effect and the $\alpha$-expected shortfall, zhang2026higher also showed how to represent the estimating equation of the parameter of interest in this bilinear form. \end{example} \begin{example}[Generalized covariance measure] When testing the conditional independence between $A$ and $Y$ given $Z$, shah2020hardness proposed to construct test statistics based on the generalized covariance measure $\tau = {\mathbb{E}} \{(A - a (Z)) (Y - b (Z))\}$, where $a (\cdot) \coloneqq {\mathbb{E}} (A \mid Z = \cdot)$ and $b (\cdot) \coloneqq {\mathbb{E}} (Y \mid Z = \cdot)$. To estimate $\tau$, the most difficult component is $\psi = {\mathbb{E}} \{a (Z) b (Z)\}$. In liu2020nearly, it was shown that if both $a (\cdot) = \alpha^{\top} \bar{\phi} (\cdot)$ and $b (\cdot) = \beta^{\top} \bar{\phi} (\cdot)$ are linear combinations of $\bar{\phi}$, then by identifying $X = \bar{\phi} (Z)$, $\psi = \alpha^{\top} \Sigma \beta = \mu^{\top} \Omega \eta$, once we set $\alpha = \Omega \mu$ and $\beta = \Omega \eta$. \end{example} More related examples can also be found in robins2008higher, rotnitzky2021characterization, chernozhukov2022locally, rotnitzky2026note. For all of the above examples, when $\Omega$ is known (referred to as the oracle case in liu2020nearly), $\psi$ can be unbiasedly estimated by its oracle second-order influence function, which is the following second-order $U$-statistic: \begin{equation} \widehat{\psi}_{2, k} (\Omega) \coloneqq \widehat{\mathbb{IF}}_{2, 2, k} (\Omega) = {\mathbb{U}}_{n, 2} \{\widehat{\mathsf{IF}}_{2, 2, k} (\Omega)\}. \end{equation} In contrast to the settings of robins2008higher and liu2017semiparametric, we consider a slightly more simplified setting in which the first-order estimator $\widehat{\psi}_{1} = 0$; otherwise $\widehat{\psi}_{2, k} (\Omega) = \widehat{\psi}_{1} + \widehat{\mathbb{IF}}_{2, 2, k} (\Omega)$. When $\Omega$ is unknown, one can construct the so-called empirical HOIF estimators taking the following form liu2017semiparametric: \begin{align} & \widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}}) \coloneqq \sum_{j = 2}^{m} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}_{\mathrm{nuis}}), \nonumber \\ & \text{where } \widehat{\mathbb{IF}}_{2, 2, k} (\widetilde{\Omega}) \coloneqq {\mathbb{U}}_{n, 2} \{\widehat{\mathsf{IF}}_{2, 2, k} (\widetilde{\Omega})\} \text{ with } \widehat{\mathsf{IF}}_{2, 2, k} (\widetilde{\Omega}) \coloneqq A_{1} X_{1}^{\top} \widetilde{\Omega} X_{2} Y_{2}, \text{ and for $j = 3, 4, \cdots$} \\ & \widehat{\mathbb{IF}}_{j, j, k} (\widetilde{\Omega}) \coloneqq (-1)^{j} {\mathbb{U}}_{n, j} \{\widehat{\mathsf{IF}}_{j, j, k} (\widetilde{\Omega})\} \text{ with } \widehat{\mathsf{IF}}_{j, j, k} (\widetilde{\Omega}) \coloneqq A_{1} X_{1}^{\top} \widetilde{\Omega} \Big\{ \prod_{s = 3}^{j} (X_{s} X_{s}^{\top} - \widetilde{\Sigma}) \widetilde{\Omega} \Big\} X_{2} Y_{2}. \nonumber \end{align} Here, $\widetilde{\Sigma}$ and $\widetilde{\Omega}$ denote, respectively, some generic estimators of $\Sigma$ and $\Omega$. Furthermore, $\widehat{\Omega}_{\mathrm{nuis}} = \widehat{\Sigma}_{\mathrm{nuis}}^{-1}$, with $\widehat{\Sigma}_{\mathrm{nuis}}$ the sample Gram matrix estimator computed from a separate sample ${\mathcal{D}}_{\mathrm{nuis}}$ independent of our main sample ${\mathcal{D}}$. \begin{remark} We choose the above notation convention to strictly follow earlier works on HOIFs robins2008higher, robins2016technical, robins2023minimax, liu2017semiparametric, liu2024assumption. For example, robins2008higher reserves the notation $\widehat{\mathbb{IF}}_{j, k} (\widetilde{\Omega})$ for $\widehat{\mathbb{IF}}_{j, k} (\widetilde{\Omega}) \coloneqq \sum_{l = 2}^{j} \widehat{\mathbb{IF}}_{l, l, k} (\widetilde{\Omega})$. We also choose to use $\widehat{\mathsf{IF}}$ and $\widehat{\mathbb{IF}}$ instead of $\mathsf{IF}$ and $\mathbb{IF}$ throughout to keep the notation more aligned with the scenario in which all $A, Y, X$ may in fact depend on some first-step nuisance estimates. \end{remark} In particular, liu2017semiparametric established the following results on $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$. Here, we only provide the simplified version of their results and liu2017semiparametric in fact provide more comprehensive characterizations of both the bias and variance bounds of $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$. \begin{proposition} Under Assumptions (ref)--(ref), the following results hold. \begin{enumerate}[label = (\arabic*)] • The bias of $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$ can be bounded as follows: \begin{align*} |{\mathbb{E}} \{\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}}) - \psi\}| \lesssim \Vert A \Vert_{2} \cdot \Vert Y \Vert_{2} \cdot \Big( \frac{k}{n} \Big)^{m / 2}. \end{align*} • The variance of $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$ can be bounded as follows if $k \lesssim \frac{n}{\log^{3} n}$ and $m \asymp \log n$: \begin{align*} \mathrm{var} \{\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})\} \lesssim \frac{1}{n} + \frac{k}{n^{2}}. \end{align*} • Under the same additional conditions in (2), let $\nu_{\mathrm{nuis}}^{2} \coloneqq \lim_{n \rightarrow \infty} n \mathrm{var} \{\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})\}$. $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$ is $\sqrt{n}$-CAN, that is: \begin{align*} \sqrt{n} \{\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}}) - \psi\} \rightsquigarrow_{{\mathbb{P}}} {\mathcal{N}} (0, \nu^{2}_{\mathrm{nuis}}). \end{align*} \end{enumerate} \end{proposition} \section{The New HOIF Estimators, Statistical Guarantees, and M\"obius Inversion} \subsection{The new HOIF estimators and statistical guarantees} As alluded to in the Introduction, although the empirical HOIF estimator $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$ dispenses with the need of a (nonparametric) density estimator $\widehat{g}$ of $g$, it can be numerically unstable when the dimension $k$ is large compared to the sample size $n$. As demonstrated in simulation studies shown in recent work liu2017semiparametric, zhang2026higher, the finite-sample performance of $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$ indeed degrades as the condition number $\rho = \rho (n) \coloneqq k / n$ increases with $k$. To resolve the numerical instability of $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$, we instead construct the following HOIF estimator: \begin{equation} \widehat{\psi}_{m, k} (\widehat{\Omega}) \coloneqq \sum_{j = 2}^{m} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}). \end{equation} As mentioned, the 2nd- and 3rd-order versions of $\widehat{\psi}_{m, k} (\widehat{\Omega})$ have appeared in the previous work of the last author of this article liu2020nearly, but there was no theoretical proof. The sole difference between our new HOIF estimator $\widehat{\psi}_{m, k} (\widehat{\Omega})$ and the empirical HOIF estimator $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$ is that we now estimate $\Omega = \Sigma^{-1}$ by the inverse sample Gram matrix estimator $\widehat{\Omega}$ not from another independent sample ${\mathcal{D}}_{\mathrm{nuis}}$, but from the same sample ${\mathcal{D}}$ used to construct the HOIF estimator. Due to the correlation induced by $\widehat{\Omega}$, it is more challenging to analyze the statistical properties of $\widehat{\psi}_{m, k} (\widehat{\Omega})$, compared to $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$ in liu2017semiparametric. Overcoming this technical challenge to obtain theoretical guarantees parallel to those in Proposition (ref) is the main contribution of this article. \begin{remark} We explain why $\widehat{\psi}_{m, k} (\widehat{\Omega})$ has improved stability compared to $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$. Intuitively, since $\widehat{\Sigma}$ contains the same sample ${\mathcal{D}}$ and enters $\widehat{\psi}_{m, k} (\widehat{\Omega})$ as a “denominator”, it exhibits a self-normalization phenomenon not shared by $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$, as $\widehat{\Omega}_{\mathrm{nuis}} = \widehat{\Sigma}_{\mathrm{nuis}}^{-1}$ is computed from a different sample. We refer readers to Section S4.3 of liu2020nearly for further explanations. \end{remark} \begin{remark} chen2025computing develop an algorithm for the exact computation of $\widehat{\psi}_{m, k} (\widehat{\Omega})$. In particular, they showed that the exact time complexity arora2009computational of computing $\widehat{\psi}_{m, k} (\widehat{\Omega})$ is $O (n^{\kappa})$, where $\kappa$ is the treewidth of an undirected graph associated with the $U$-statistic kernel of $\widehat{\psi}_{m, k} (\widehat{\Omega})$. If one is willing to sacrifice some efficiency, it is entirely possible to compute each $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$ as an incomplete higher-order $U$-statistic with almost the same complexity as $j$ matrix multiplications kong2018estimating. \end{remark} Next, we present Theorem (ref), the main and most advanced result of this article. \begin{theorem} Under Assumptions (ref)--(ref), the following results hold. \begin{enumerate}[label = (\arabic*)] • The bias of $\widehat{\psi}_{m, k} (\widehat{\Omega})$ can be characterized as follows: \begin{align*} |{\mathbb{E}} \{\widehat{\psi}_{m, k} (\widehat{\Omega}) - \psi\}| \lesssim (\Vert A \Vert_{2} \cdot \Vert Y \Vert_{2} + \Vert A \Vert_{\infty} \cdot \Vert Y \Vert_{2} + \Vert A \Vert_{2} \cdot \Vert Y \Vert_{\infty}) \Big( \frac{k m}{n} \Big)^{ \lceil \frac{m - 1}{4} \rceil \vee 1}. \end{align*} • The variance of $\widehat{\psi}_{m, k} (\widehat{\Omega})$ can be bounded as follows if $k \lesssim \frac{n}{\log^{3} n}$ and $m \asymp \log n$: \begin{align*} \mathrm{var} \{\widehat{\psi}_{m, k} (\widehat{\Omega})\} \lesssim \frac{1}{n} + \frac{k}{n^{2}}. \end{align*} • If $m \lesssim \log n$ and $k \lesssim \frac{n}{\log^{3} n}$, let $\nu^{2} \coloneqq \lim_{n \rightarrow \infty} n \mathrm{var} \{\widehat{\psi}_{m, k} (\widehat{\Omega})\}$. $\widehat{\psi}_{m, k} (\widehat{\Omega})$ is $\sqrt{n}$-CAN, that is: \begin{align*} \sqrt{n} \{\widehat{\psi}_{m, k} (\widehat{\Omega}) - \psi\} \rightsquigarrow_{{\mathbb{P}}} {\mathcal{N}} (0, \nu^{2}). \end{align*} \end{enumerate} \end{theorem} In Section (ref) below, we will provide a proof sketch of the above theorem, to illustrate the main steps. The details of the proof are delegated to the Appendix. \begin{remark} In fact, once the bias of $\widehat{\psi}_{m, k} (\widehat{\Omega})$ can be shown to be $o (n^{-1 / 2})$, it is straightforward to establish the $\sqrt{n}$-CAN of $\widehat{\psi}_{m, k} (\widehat{\Omega})$ because $\widehat{\psi}_{2, k} (\Omega)$ is an unbiased and $\sqrt{n}$-CAN estimator of $\psi$, following bhattacharya1992class; see liu2020nearly for a proof and bobkov2019higher, gotze1984expansions, dobler2022functional, chakrabortty2025tail for some recent related progress on the probability theory side. \end{remark} To demonstrate the better finite-sample performance of $\widehat{\psi}_{m, k} (\widehat{\Omega})$ compared to $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$, a simple simulation study is conducted, with the setup described in Appendix (ref). Specifically, Figure (ref) compares the performance between $\widehat{\psi}_{m, k} (\widehat{\Omega})$ and $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$ when $m = 3$, varying $\rho = k / n$. All summary statistics are computed based on 250 Monte Carlo runs. It is evident that the performance of $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$ starts to break down as $\rho$ increases, whereas $\widehat{\psi}_{m, k} (\widehat{\Omega})$ maintains a very stable performance even when $\rho$ is near $1$. In particular, based on Figure (ref)(a), the RMSEs of $\widehat{\psi}_{m, k} (\widehat{\Omega})$ track those of $\widehat{\psi}_{m, k} (I)$ quite well even when $\rho$ is as large as $0.7$. In a follow-up paper, we will report numerical results from a set of more comprehensive simulation studies. \begin{figure} \caption{Finite-sample comparison between the sample-split empirical HOIF estimator $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$ and the same-sample stabilized HOIF estimator $\widehat{\psi}_{m, k} (\widehat{\Omega})$ at order $m = 3$. Panel (a) reports the RMSE on a logarithmic scale as $\rho = k / n$ varies. Panel (b) decomposes the error into absolute bias and standard deviation. The sample-split estimator becomes unstable as $\rho$ increases, whereas the stabilized estimator remains numerically stable.} \end{figure} \subsection{The M\"{o}bius inversion decomposition} Before proving our main theorem, we record an (interesting) observation regarding $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$. \begin{lemma} Write $H_{i} \coloneqq (X_{i} X_{i}^{\top} - \widehat{\Sigma}) \widehat{\Omega} = X_{i} X_{i}^{\top} \widehat{\Omega} - I$ for $i \in [n]$ (note that $H_{i}$'s appear repeatedly in the $U$-statistic kernel of $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$). Then the following elementary identity holds. \begin{equation} \sum_{i = 1}^{n} H_{i} \equiv 0. \end{equation} \end{lemma} With Lemma (ref), by exploiting a classical tool in enumerative combinatorics, \emph{M\"obius inversion} on partition lattices lauritzen1996graphical, stanley2011enumerative, mccullagh2018tensor, we can then decompose $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$ into a finite sum of lower-order $U$-statistics, which will not only be useful in the proof of Theorem (ref) to be presented in Section (ref), but also shed some light on more detailed bias reduction mechanisms of each $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$ for $j \in [m]$. Before presenting this \emph{M\"{o}bius inversion decomposition}, we introduce some additional notation. Fix any $j \geq 3$. Let $\iota = j - 2$ and ${\mathcal{R}}_{i_{1} i_{2}} \coloneqq H_{i_{1}} + H_{i_{2}}$. Let ${\mathbb{B}}_{\iota}$ consist of all finite collections ${\mathcal{B}} = \{B_{1}, \cdots, B_{r}\}$ of pairwise disjoint subsets of $[\iota]$ such that $|B_{\nu}| \ge 2$ for every $\nu \in [r]$. The collection ${\mathcal{B}}$ is allowed to be empty and is not required to cover $[\iota]$. For any ${\mathcal{B}} \in {\mathbb{B}}_{\iota}$, order its elements according to their smallest elements and define \begin{align} & K_{{\mathcal{B}}} (i_{1}, i_{2}; a_{1}, \cdots, a_{r}) \coloneqq A_{i_{1}} X_{i_{1}}^{\top} \widehat{\Omega} \Big\{ \prod_{l = 1}^{\iota} G_{{\mathcal{B}}, l}^{i_{1} i_{2}} (a_{1}, \cdots, a_{r}) \Big\} X_{i_{2}} Y_{i_{2}}, \\ & \text{ where } G_{{\mathcal{B}}, l}^{i_{1} i_{2}} (a_{1}, \cdots, a_{r}) \coloneqq \begin{cases} H_{a_{\nu}}, & l \in \bigcup_{\nu = 1}^{r} B_{\nu},\\ {\mathcal{R}}_{i_{1} i_{2}}, & l \notin \bigcup_{\nu = 1}^{r} B_{\nu}. \nonumber \end{cases} \end{align} When ${\mathcal{B}} = \emptyset$, we let $K_{\emptyset} (i_{1}, i_{2}) \coloneqq A_{i_1} X_{i_1}^\top \widehat{\Omega} {\mathcal{R}}_{i_{1} i_{2}}^{\iota} X_{i_{2}} Y_{i_{2}}$. We are now ready to present the following lemma, a proof of which is deferred to Appendix (ref). \begin{lemma} $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$ can be decomposed as follows: \begin{equation} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) = \sum_{{\mathcal{B}} \in {\mathbb{B}}_{\iota}} c_{{\mathcal{B}}, n} {\mathbb{U}}_{n, 2 + |{\mathcal{B}}|} (K_{{\mathcal{B}}}), \ \text{where} \ c_{{\mathcal{B}}, n} \coloneqq (-1)^{|{\mathcal{B}}|} \Big\{ \prod_{B \in {\mathcal{B}}}(|B| - 1) \Big\} \frac{(n - j) !}{(n - 2 - |{\mathcal{B}}|) !}. \end{equation} The coefficients $c_{{\mathcal{B}},n}$ are the so-called Möbius coefficients. In particular, every term in the expansion is a $U$-statistic of order at most $2 + \lfloor \frac{\iota}{2} \rfloor = 2 + \lfloor \frac{j - 2}{2} \rfloor$. \end{lemma} \begin{remark} We illustrate Lemma (ref) with the cases $j = 3$ and $j = 4$. \begin{itemize} • When $j = 3$ , we have $\iota = 1$. Since no non-singleton element can be formed from the singleton set $\{1\}$, the only element family is ${\mathcal{B}} = \emptyset$. Hence, \begin{align*} \widehat{\mathbb{IF}}_{3, 3, k} (\widehat{\Omega}) & = \frac{1}{n - 2} {\mathbb{U}}_{n, 2} (A_{1} X_{1}^{\top} \widehat{\Omega} {\mathcal{R}}_{1 2} X_{2} Y_{2}) = \frac{1}{n - 2} {\mathbb{U}}_{n,2} (A_{1} X_{1}^{\top} \widehat{\Omega} (H_{1} + H_{2}) X_{2} Y_{2}) \\ & = \frac{1}{n - 2} {\mathbb{U}}_{n, 2} (A_{1} X_{1}^{\top} \widehat{\Omega} X_{1} X_{1}^{\top} \widehat{\Omega} X_{2} Y_{2}) + \frac{1}{n - 2} {\mathbb{U}}_{n, 2} (A_{1} X_{1}^{\top} \widehat{\Omega} X_{2} X_{2}^{\top} \widehat{\Omega} X_{2} Y_{2}) \\ & \quad - \frac{2}{n - 2} \underbrace{{\mathbb{U}}_{n, 2} (A_{1} X_{1}^{\top} \widehat{\Omega} X_{2} Y_{2})}_{\equiv \, \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})}. \end{align*} In particular, it is not difficult to see that the dominating terms in $\widehat{\mathbb{IF}}_{3, 3, k} (\widehat{\Omega})$, corresponding to the first two terms in the last equality of the above display, match the dominating bias terms of $\widehat{\psi}_{2, k} (\widehat{\Omega}) = \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$, except that $\Omega$ is replaced by $\widehat{\Omega}$. It is also worth noting that the monomials of the leverage scores (terms of the form $X_{i}^{\top} \widehat{\Omega} X_{j}$) up to degree $2$ appear in $\widehat{\mathbb{IF}}_{3, 3, k} (\widehat{\Omega})$. • When $j = 4$, we have $\iota = 2$. There are two possible element families: ${\mathcal{B}} = \emptyset$ and ${\mathcal{B}} = \{\{1, 2\}\}$. Hence, \begin{align*} \widehat{\mathbb{IF}}_{4, 4, k}(\widehat{\Omega}) = \frac{1}{(n - 2) (n - 3)} {\mathbb{U}}_{n, 2} (A_{1} X_{1}^{\top} \widehat{\Omega} {\mathcal{R}}_{1 2}^{2} X_{2} Y_{2}) - \frac{1}{n - 3} {\mathbb{U}}_{n, 3} (A_{1} X_{1}^{\top} \widehat{\Omega} H_{3}^{2} X_{2} Y_{2}). \end{align*} By elementary algebra, we have the following: \begin{align*} {\mathcal{R}}_{1 2}^{2} & = (X_{1} X_{1}^{\top} \widehat{\Omega} + X_{2} X_{2}^{\top} \widehat{\Omega} - 2 I)^2 \\ & = X_{1} X_{1}^{\top} \widehat{\Omega} X_{1} X_{1}^{\top} \widehat{\Omega} + X_{1} X_{1}^{\top} \widehat{\Omega} X_{2} X_{2}^{\top} \widehat{\Omega} + X_{2} X_{2}^{\top} \widehat{\Omega} X_{1} X_{1}^{\top} \widehat{\Omega} \\ & \quad + X_{2} X_{2}^{\top} \widehat{\Omega} X_{2} X_{2}^{\top} \widehat{\Omega} - 4 X_{1} X_{1}^{\top} \widehat{\Omega} - 4 X_{2} X_{2}^{\top} \widehat{\Omega} + 4 I, \end{align*} and \begin{align*} H_{3}^{2} = (X_{3} X_{3}^{\top} \widehat{\Omega} - I)^{2} = X_{3} X_{3}^{\top} \widehat{\Omega} X_{3} X_{3}^{\top} \widehat{\Omega} - 2 X_{3} X_{3}^{\top} \widehat{\Omega} + I. \end{align*} Therefore, $\widehat{\mathbb{IF}}_{4, 4, k} (\widehat{\Omega})$ reads as follows: \begin{align*} & \ \widehat{\mathbb{IF}}_{4, 4, k} (\widehat{\Omega}) \\ = & \ \frac{1}{(n - 2)(n - 3)} \Big\{ {\mathbb{U}}_{n,2} (A_{1} (X_{1}^{\top} \widehat{\Omega} X_{1})^{2} X_{1}^{\top} \widehat{\Omega} X_{2} Y_{2}) + {\mathbb{U}}_{n, 2} (A_{1} X_{1}^{\top} \widehat{\Omega} X_{2} (X_{2}^{\top} \widehat{\Omega} X_{2})^{2} Y_{2}) \Big\} \\ & + \frac{1}{(n - 2) (n - 3)} \Big\{ {\mathbb{U}}_{n,2} (A_{1} (X_{1}^{\top} \widehat{\Omega} X_{2})^{3} Y_{2}) + {\mathbb{U}}_{n,2} (A_{1} X_{1}^{\top} \widehat{\Omega} X_{1} X_{1}^{\top} \widehat{\Omega} X_{2} X_{2}^{\top} \widehat{\Omega} X_{2} Y_{2}) \Big\} \\ & - \frac{6}{(n - 2)(n - 3)} \Big\{ {\mathbb{U}}_{n,2} (A_{1} X_{1}^{\top} \widehat{\Omega} X_{1} X_{1}^{\top} \widehat{\Omega} X_{2} Y_{2}) + {\mathbb{U}}_{n,2} (A_{1} X_{1}^{\top} \widehat{\Omega} X_{2} X_{2}^{\top} \widehat{\Omega} X_{2} Y_{2}) \Big\} \\ & + \frac{n + 6}{(n - 2) (n - 3)} \underbrace{{\mathbb{U}}_{n,2} (A_{1} X_{1}^{\top} \widehat{\Omega} X_{2} Y_{2})}_{\equiv \, \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})} - \frac{1}{n - 3} {\mathbb{U}}_{n, 3} (A_{1} X_{1}^{\top} \widehat{\Omega} X_{3} X_{3}^{\top} \widehat{\Omega} X_{3} X_{3}^{\top} \widehat{\Omega} X_{2} Y_{2}). \end{align*} Similarly, $\widehat{\mathbb{IF}}_{4, 4, k} (\widehat{\Omega})$ matches the dominating bias terms of $\widehat{\psi}_{3, k} (\widehat{\Omega}) = \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) + \widehat{\mathbb{IF}}_{3, 3, k} (\widehat{\Omega})$, except that $\Omega$ is replaced by $\widehat{\Omega}$. It is also straightforward to see that the monomials of the leverage scores up to degree $3$ appear in $\widehat{\mathbb{IF}}_{4, 4, k} (\widehat{\Omega})$. \end{itemize} \end{remark} \section{Proof Sketch of Theorem (ref)} In this section, we sketch the proof of Theorem (ref). We focus only on the first two statements of Theorem (ref), as we have argued in Remark (ref) how to prove that $\widehat{\psi}_{m, k} (\widehat{\Omega})$ is $\sqrt{n}$-CAN. Specifically, Section (ref) below provides a sketch of the bias analysis establishing part (1) of Theorem (ref), whereas Section (ref) sketches the proof of variance bound in Theorem (ref). Before embarking on the proof sketch, in Section (ref), we first introduce a useful proof device, which we refer to as the \emph{graph-counting lemma} (Lemma (ref)). Lemma (ref) turns the problem of controlling moment bounds of certain $U$-statistic kernels into an enumerative combinatorics problem on graphs, drastically simplifying the proof. Throughout the bias and variance analyses, we impose $\Sigma = \Omega = I$ without loss of generality by Assumption (ref). \subsection{A graph-counting lemma} The following graph-counting lemma gives the required bound in terms of the first Betti number (or equivalently, the circuit rank) of $G$ stanley2011enumerative. \begin{lemma} Let $G = (V, E)$ be a fixed undirected graph, where $V$ is a collection of observation labels and each edge $e = (u,v) \in E$ represents a \emph{bilinear structure} $X_{u}^{\top} B_{e} X_{v}$, in the sense that two vertices $u$ and $v$ are contracted by an edge induced by this bilinear structure. Self-loops are admissible and each self-loop contributes two half-edges at the same vertex. Let \begin{equation*} v(G) \coloneqq|V|, \ e(G) \coloneqq|E|, \ \kappa(G) \coloneqq\text{the number of connected components in }G, \end{equation*} and $\mathfrak{r} (G) \coloneqq e (G) + \kappa (G) - v (G)$ is the first Betti number of $G$. Assume that the matrices $\{B_{e}: e\in E\}$ are independent of the vectors $\{X_{v} : v\in V\}$ and satisfy \begin{equation*} \max_{e \in E} \|B_{e}\|_{\mathrm{op}}\le C, \end{equation*} almost surely. Suppose that Assumptions (ref) and (ref) hold, we have \begin{equation} \Big| {\mathbb{E}} \Big( \prod_{e = (u, v) \in E} X_{u}^{\top} B_{e} X_{v} \Big) \Big| \lesssim k^{\mathfrak{r} (G)}. \end{equation} The implicit constant depends only on the fixed graph $G$, moments of the observed data $O$, and the uniform operator-norm bound, but not on $n$ or $k$. \end{lemma} A proof of this result can be found in Appendix (ref). Lemma (ref) associates $U$-statistic kernels only involving products in the form of $\prod_{e = (u,v) \in E} X_{u}^{\top} B_{e} X_{v}$ (the integrand in (ref)), which we refer to as \emph{multiplicative-kernels}, with an (undirected) graph $G$, with which controlling moment bounds in the form of (ref) can be conveniently translated into counting the first Betti number $\mathfrak{r} (G)$ of the graph $G$. \subsection{Bias analysis} Since $\widehat{\psi}_{2, k} (\Omega)$, as defined in (ref), is unbiased for $\psi$, we can represent the bias of $\widehat{\psi}_{m, k} (\widehat{\Omega})$ as: \begin{equation} {\mathcal{B}}_{m, k} \coloneqq {\mathbb{E}} \{\widehat{\psi}_{m, k} (\widehat{\Omega}) - \psi\} = {\mathbb{E}} \{\widehat{\psi}_{m, k} (\widehat{\Omega}) - \widehat{\psi}_{2, k} (\Omega)\} = {\mathbb{E}} \{\widehat{\psi}_{m, k} (\widehat{\Omega}) - \widehat{\psi}_{2, k} (I)\}. \end{equation} We divide the bias analysis into the following steps. The detailed proofs can be found in Appendix (ref). \begin{enumerate}[label = \roman*.] • The first step rewrites ${\mathcal{B}}_{m, k}$ by applying Lemma (ref) and Lemma (ref) presented later in this subsection in a sequence, up to the point that ${\mathcal{B}}_{m, k}$ can be decomposed into a remainder ${\mathcal{R}}_{m, k, J}$ of the form in (ref) and a summation of terms ${\mathcal{M}}_{c}^{(J)}$ defined in (ref). The essential idea is to “linearize” $\widehat{\Omega}$ by the Neumann series expansion (Lemma (ref) in Appendix (ref)). • In the second step, we further refine the representation of ${\mathcal{B}}_{m, k}$ obtained in \textbf{Step i}. Specifically, Lemma (ref), to be presented later in this subsection, demonstrates that many ${\mathcal{M}}_{c}^{(J)}$'s obtained in \textbf{Step i} are zero when $c$ is sufficiently small in the decomposition. This critical observation results from a couple of intermediate results (Lemma (ref) and Lemma (ref)), which we detail in the proof of Lemma (ref) in Appendix (ref). As will be clear in the proof, these intermediate results are used to show that the terms in ${\mathcal{M}}_{c}^{(J)}$ cancel each other meticulously when $c$ is below a certain threshold (denoted by ${\mathsf{c}}_{m} \coloneqq \lceil (m - 1) / 2 \rceil$). • We next bound all relevant terms from \textbf{Step ii} by applying the graph-counting lemma (Lemma (ref)) introduced in Section (ref), culminating in Lemma (ref). Finally, the remainder term ${\mathcal{R}}_{m, k, J}$ is controlled by Lemma (ref), which completes the analysis of the bias bound. \end{enumerate} \paragraph{Step i.} We first represent ${\mathcal{B}}_{m, k}$ in a particular form as stated in the following lemma; see its proof at the beginning of Appendix (ref). \begin{lemma} ${\mathcal{B}}_{m, k}$ admits the following alternative representations: \begin{align} {\mathcal{B}}_{m, k} & = \sum_{j = 1}^{m - 1} (-1)^{j + 1} \binom{m - 1}{j} {\mathbb{E}} \Big\{ A_{m - 1} X_{m - 1}^{\top} \Big( \prod_{s = 0}^{j - 1} X_{s} X_{s}^{\top} \widehat{\Omega} - I \Big) X_{m} Y_{m} \Big\} \notag \\ & = \sum_{j = 1}^{m - 1} (-1)^{j + 1} \binom{m - 1}{j} \sum_{\emptyset \neq S \subseteq [j - 1] \cup \{0\}} {\mathbb{E}} \Big\{ A_{m - 1} X_{m - 1}^{\top} \prod_{s = 0}^{j - 1} \Big( X_{s} X_{s}^{\top} (\widehat{\Omega} - I)^{\mathbbm{1} \{s \in S\}} \Big) X_{m} Y_{m} \Big\}. \end{align} Here, we use the convention that $s = 0$ corresponds to the identity matrix $I$. \end{lemma} By Lemma (ref), ${\mathcal{B}}_{m, k}$ can be expressed as a binomially weighted sum of ordered product expectations indexed by nonempty subsets of positions at which the factor $\widehat{\Omega} - I$ is inserted. The Neumann series expansion (Lemma (ref) in Appendix (ref)) gives: \begin{equation} \widehat{\Omega} - I = \sum_{j = 1}^{J} \Delta_{n}^{j} + {\mathsf{R}}_{J}, \quad \Delta_{n} \coloneqq I - \widehat{\Sigma}, \quad {\mathsf{D}}_{J} \coloneqq \sum_{j = 1}^{J} \Delta_{n}^{j}, \quad {\mathsf{R}}_{J} \coloneqq \Delta_{n}^{J + 1} \widehat{\Omega}. \end{equation} \begin{remark} The identity (ref) is exact for every $J$. However, in the proof, to avoid the last term ${\mathsf{R}}_{J}$ as it involves the nonlinear $\widehat{\Omega}$, we take $J = J (n) = \lceil C_{0} \log n \rceil$ for some sufficiently large constant $C_{0}$. This choice of $J$ makes ${\mathsf{R}}_J$ negligible: on the event $\|\Delta_{n}\|_{\mathrm{op}} \le r_{n}$ and $\|\widehat{\Omega}\|_{\mathrm{op}} \le C$ for some large enough constant $C > 0$, $\|{\mathsf{R}}_{J}\|_{\mathrm{op}} \le \|\Delta_{n}\|_{\mathrm{op}}^{J + 1} \|\widehat{\Omega}\|_{\mathrm{op}} \lesssim r_{n}^{J + 1}$. At the same time, under the regime $m \asymp \log n$ and $k \lesssim n / \log^{3} n$, this choice satisfies \begin{equation} \frac{m J k}{n} \lesssim \frac{\log^{2} n}{\log^{3} n} = o(1). \end{equation} The condition (ref) is needed in various places in the proof details; e.g., Lemma (ref) in Appendix (ref). \end{remark} We next state a lemma that further decomposes ${\mathcal{B}}_{m, k}$ into components that share the same multiplicity of $\Delta_{n} = I - \widehat{\Sigma}$, after the Neumann series expansion of $\widehat{\Omega} - I$. The proof is delegated to Appendix (ref). \begin{lemma} For an integer $J \geq 1$, \begin{equation} {\mathcal{B}}_{m, k} = \sum_{c = 1}^{(m - 1) J} {\mathcal{M}}_{c}^{(J)} + {\mathcal{R}}_{m, k, J}, \end{equation} where ${\mathcal{R}}_{m, k, J}$ is the collection of all terms containing at least one occurrence of ${\mathsf{R}}_{J}$, namely \begin{align} & {\mathcal{R}}_{m, k, J} \coloneqq \\ & \sum_{j = 1}^{m - 1} (-1)^{j + 1} \binom{m - 1}{j} \sum_{\emptyset \neq S \subseteq [j - 1] \cup \{0\}} \sum_{\emptyset \neq T \subseteq S} {\mathbb{E}} \Big\{ A_{m - 1} X_{m - 1}^{\top} \prod_{s = 0}^{j - 1} (X_{s} X_{s}^{\top} {\mathsf{D}}_{J}^{\mathbbm{1}\{s \in S \setminus T\}} {\mathsf{R}}_{J}^{\mathbbm{1} \{s \in T\}}) X_{m} Y_{m} \Big\}, \nonumber \end{align} and ${\mathcal{M}}_{c}^{(J)}$ is defined as: \begin{align} & {\mathcal{M}}_{c}^{(J)} \coloneqq \\ & \sum_{j = 1}^{m - 1}(-1)^{j + 1} \binom{m - 1}{j} \sum_{r = 1}^{c \wedge j} \sum_{\substack{S \subseteq [j - 1] \cup \{0\} \\ |S| = r}} \sum_{\substack{(\ell_{s'})_{s' \in S} \in [J]^{r} \\ \sum_{s' \in S} \ell_{s'} = c}} {\mathbb{E}} \Big\{ A_{m - 1} X_{m - 1}^{\top} \prod_{s = 0}^{j - 1} (X_{s} X_{s}^{\top} \Delta_{n}^{\ell_{s} \mathbbm{1} \{s \in S\}}) X_{m} Y_{m} \Big\}. \nonumber \end{align} Here, we use the convention that $\ell_{s} = 0$ for $s \notin S$. \end{lemma} Equivalently, ${\mathcal{M}}_{c}^{(J)}$ sums up all terms for which the multiplicity $\Delta_{n} = I - \widehat{\Sigma}$ equals $c$. When the truncation level $J$ is fixed, we write ${\mathcal{M}}_{c}$ for ${\mathcal{M}}_{c}^{(J)}$ to simplify the notation. This step reduces the analysis to each term ${\mathcal{M}}_{c}^{(J)}$ for $c \in [(m - 1) J]$ and the remainder ${\mathcal{R}}_{m, k, J}$. \paragraph{Step ii.} Recall that, by Lemma (ref), we have $\sum_{s = 0}^{j - 1} \ell_{s} = c$. We further refine ${\mathcal{M}}_{c}^{(J)}$ for $c \in [(m - 1) J]$ by showing that ${\mathcal{M}}_{c}^{(J)} = 0$ when $c$ is sufficiently small. More concretely, we establish Lemma (ref) below. \begin{lemma} Under the notation of Lemma (ref), let $ {\mathsf{c}}_{m} \coloneqq \left \lceil \frac{m - 1}{2} \right \rceil$. Then, for every integer $J \ge 1$, \begin{equation} {\mathcal{M}}_{c}^{(J)} = 0, \quad 1 \leq c < {\mathsf{c}}_{m}. \end{equation} Equivalently, when $J$ is fixed and we write ${\mathcal{M}}_{c}$ for ${\mathcal{M}}_{c}^{(J)}$, one has ${\mathcal{M}}_{c} = 0$ for all $1 \le c < {\mathsf{c}}_{m}$. \end{lemma} The proof of Lemma (ref) is deferred to Appendix (ref). As mentioned, showing that ${\mathcal{M}}_{c}^{(J)}$ is exactly zero demands a careful calculation to demonstrate that all terms involved in ${\mathcal{M}}_{c}^{(J)}$ cancel each other out. To achieve this, in the proof, we first establish Lemma (ref) and Lemma (ref), based on which Lemma (ref) can be proved. \paragraph{Step iii.} We now bound the remainder term ${\mathcal{R}}_{m, k, J}$ and the non-zero ${\mathcal{M}}_{c}^{(J)}$'s after \textbf{Step ii}. Define \begin{equation*} s_{c} \coloneqq \Big\lceil \frac{c}{2} \Big\rceil \vee 1, \ \rho_{j} \coloneqq j \rho = \frac{j k}{n},\ \zeta_{A, Y}\coloneqq \|A\|_2 \|Y\|_2 + \|A\|_\infty \|Y\|_2 + \|A\|_2 \|Y\|_\infty. \end{equation*} First, Lemma (ref) below exhibits the order of ${\mathcal{M}}_{c}^{(J)}$ when it is not identically zero. \begin{lemma} Let $ {\mathsf{c}}_{m} \coloneqq \left \lceil \frac{m - 1}{2} \right \rceil, \ s_{c} \coloneqq \left \lceil \frac{c}{2} \right \rceil \vee 1 $. Suppose that $ \frac{C m J k}{n} \leq \eta < 1 $. Then, for every ${\mathsf{c}}_{m} \leq c \leq (m - 1)J$, \begin{equation} \left| {\mathcal{M}}_{c}^{(J)} \right| \lesssim_{\eta} \zeta_{A, Y} \Big( \frac{C (m \vee c) k}{n} \Big)^{s_{c}}. \end{equation} \end{lemma} Then, Lemma (ref) below controls the order of the remainder term ${\mathcal{R}}_{m, k, J}$. \begin{lemma} Let $J = \lceil C_{0} \log n \rceil$ for some sufficiently large universal constant $C_{0}$. Suppose that $m \asymp \log n$ and $\frac{C m J k}{n} \leq \eta < 1 $. Then \begin{equation*} |{\mathcal{R}}_{m, k, J}| \lesssim \zeta_{A, Y} \Big( \frac{C m k}{n} \Big)^{s_{{\mathsf{c}}_{m}}} . \end{equation*} \end{lemma} Again, we defer the proofs of the above two lemmas to Appendix (ref). In particular, the proofs of both results rely on the graph-counting Lemma (ref) by associating $U$-statistic kernels emerged from rewriting ${\mathcal{B}}_{m, k}$ with undirected graphs. Specifically, bounding the mean of these $U$-statistic kernels will be reduced to counting the first Betti number of the associated undirected graph. By Lemma (ref), \begin{equation*} |{\mathcal{M}}_{c}^{(J)}| \lesssim \zeta_{A, Y} \Big(\frac{C (m \vee c) k}{n} \Big)^{s_{c}},\ {\mathsf{c}}_{m} \leq c \leq (m - 1) J. \end{equation*} We then divide our analysis into two scenarios. \begin{itemize} • For ${\mathsf{c}}_{m} \leq c \leq m$, we have $m \vee c = m,\ \frac{C (m \vee c) k}{n} = C \rho_{m}$. Therefore, \begin{equation*} \sum_{c = {\mathsf{c}}_{m}}^{m} |{\mathcal{M}}_{c}^{(J)}| \lesssim \zeta_{A, Y} \sum_{c = {\mathsf{c}}_{m}}^{m} (C \rho_{m})^{s_{c}}. \end{equation*} Since $s_{c} = \lceil c / 2 \rceil \vee 1$, pairing adjacent values of $c$ shows that each exponent $s_{c}$ occurs at most twice. Thus, under $C \rho_{m} < 1$, \begin{equation*} \sum_{c = {\mathsf{c}}_{m}}^{m} (C \rho_{m})^{s_{c}} \leq 2 \sum_{\ell = s_{{\mathsf{c}}_{m}}}^{s_{m}} (C \rho_{m})^{\ell} \lesssim (C \rho_{m})^{s_{{\mathsf{c}}_{m}}}. \end{equation*} Consequently, \begin{equation*} \sum_{c = {\mathsf{c}}_{m}}^{m} |{\mathcal{M}}_{c}^{(J)}| \lesssim \zeta_{A, Y} (C \rho_{m})^{s_{{\mathsf{c}}_{m}}}. \end{equation*} • For $c > m$, set $g_{c} \coloneqq (C \rho_{c})^{c / 2}$. Since $C \rho_{c} < 1$ and $s_{c} \geq c / 2$, \begin{equation*} |{\mathcal{M}}_{c}^{(J)}| \lesssim \zeta_{A, Y} g_{c}. \end{equation*} There exists $C > 0$ such that $\sqrt{C e \rho m J} \leq q < 1$. Then \begin{equation*} \frac{g_{c + 1}}{g_{c}} = \sqrt{C \rho_{c + 1}} \Big(1 + \frac{1}{c} \Big)^{c/2} \leq \sqrt{C e\rho_{c + 1}} \leq q, \end{equation*} so the sequence $\{g_{c}\}_{c > m}$ decreases to zero at a geometric rate. Therefore, \begin{equation*} \sum_{c = m + 1}^{(m - 1) J} |{\mathcal{M}}_{c}^{(J)}| \lesssim \zeta_{A, Y} g_{m} = \zeta_{A, Y} (C \rho_{m})^{m / 2} \leq \zeta_{A, Y} (C \rho_{m})^{s_{{\mathsf{c}}_{m}}}. \end{equation*} \end{itemize} Integrating the above two scenarios has the following consequence: \begin{equation*} \sum_{c = {\mathsf{c}}_{m}}^{(m - 1) J} |{\mathcal{M}}_{c}^{(J)}| \lesssim \zeta_{A, Y} \Big( \frac{C m k}{n} \Big)^{s_{{\mathsf{c}}_{m}}}. \end{equation*} Combining (ref), (ref) and Lemma (ref) yields the following: \begin{equation*} |{\mathcal{B}}_{m, k}| \lesssim \zeta_{A, Y} \Big( \frac{C m k}{n} \Big)^{s_{{\mathsf{c}}_{m}}}, \ \text{where} \ s_{{\mathsf{c}}_{m}} = \Big\lceil \frac{{\mathsf{c}}_{m}}{2} \Big\rceil = \Big\lceil \frac{m - 1}{4} \Big\rceil. \end{equation*} This completes of the proof of the bias bound. \subsection{Variance analysis} The variance analysis is much more complicated than that of $\widehat{\psi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}})$ in liu2017semiparametric, because we can no longer use Hoeffding decomposition. We divide the variance analysis into the following steps: \begin{enumerate}[label = \roman*.] • We first apply Minkowski's inequality to reduce the task of bounding the variance of $ \widehat{\psi}_{m, k} (\widehat{\Omega})$ to the task of bounding the variance of each $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$ for $j = 2, \cdots, m$; see Lemma (ref). We then invoke Lemma (ref) (through M\"obius inversion) to rewrite each $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$ as a finite sum of lower-order $U$-statistics. • Starting from the lower-order $U$-statistics obtained in \textbf{Step i}, we further expand each $U$-statistic kernel into kernels involving only products of bilinear forms $X_{i}^{\top} M X_{j}$ for $i, j \in [n]$ (abbreviated as \emph{multiplicative-kernels}), with $M$ being some square matrix of size $k$. We then associate each multiplicative-kernel with an undirected graph, whose vertices correspond to all sample indices $i, j$ involved in the aforementioned bilinear forms $X_{i}^{\top} M X_{j}$ and whose edges describe whether a pair of indices $i, j$ are present in any of these bilinear forms. We then prove a generic variance bound for these $U$-statistics by combining several technical ingredients: \begin{enumerate}[label = (\arabic*)] • a standard decomposition of the variance of a $U$-statistic into a sum of terms organized by the size of overlapped indices; • a counting argument based on the first Betti number of the graph associated with the kernel, as stated previously in Lemma (ref); • the Neumann series expansion of $\widehat{\Omega}$ and leave-*-out analysis; and finally • the Efron--Stein inequality efron1981jackknife, rajendran2023concentration. \end{enumerate} • Finally, we combine the expansion based on M\"obius inversion in Lemma (ref) in \textbf{Step i} and the results in \textbf{Step ii} to obtain the desired variance bound for $\widehat{\mathbb{IF}}_{j, j, k}(\widehat{\Omega})$ for each $j = 2, 3, \cdots, m$. \end{enumerate} The logical flow of the argument is summarized in Figure (ref). \begin{figure} \begin{tikzpicture}[ >={Latex[length=2mm]}, roadbox/.style={ draw=black, rounded corners=2pt, align=center, text width=0.435\textwidth, inner xsep=6pt, inner ysep=5pt, font=\scriptsize }, arr/.style={ ->, line width=0.7pt }, downmark/.style={ font= } ] \node[roadbox] (s1) { \textbf{1. Minkowski inequality}\\left[1mm] $\displaystyle \mathrm{var}^{1 / 2} \Big\{ \sum_{j = 2}^{m}\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) \Big\} \le \sum_{j = 2}^{m} \mathrm{var}^{1 / 2} \Big\{ \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) \Big\} $\\left[1mm] Lemma (ref) }; \node[roadbox, right=0.55cm of s1] (s2) { \textbf{2. M\"obius-inversion expansion}\\left[1mm] $\displaystyle \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) = \sum_{{\mathcal{B}} \in {\mathbb{B}}_{\iota}} c_{{\mathcal{B}}, n}\, {\mathbb{U}}_{n, 2 + |{\mathcal{B}}|} (K_{{\mathcal{B}}}), \ \iota = j-2 $\\left[1mm] Lemma (ref) }; \node[roadbox, below=0.55cm of s1] (s3) { \textbf{3. Reduction to effective chains}\\left[1mm] $\displaystyle \begin{aligned} &\quad\qquad{\mathbb{U}}_{n, 2 + r} (K_{{\mathcal{B}}}) = \sum_{\varepsilon \in {\mathcal{E}} ({\mathcal{B}})} d_{{\mathcal{B}}, \varepsilon} T_{a_{\varepsilon}, \ell_{\varepsilon}, \gamma_{\varepsilon}},\\left[1mm] &T_{a, \ell, \gamma} = {\mathbb{U}}_{n, a} \Big\{ A_{i_{1}} X_{i_{1}}^{\top} \widehat{\Omega} \Big( \prod_{s = 1}^{\ell} X_{i_{\gamma (s)}} X_{i_{\gamma (s)}}^{\top} \widehat{\Omega} \Big) X_{i_2} Y_{i_2} \Big\}. \end{aligned} $\\left[1mm] Lemma (ref) }; \node[roadbox, right=0.55cm of s3] (s4) { \textbf{4. Covariance decomposition}\\left[1mm] $\displaystyle \mathrm{var} (T_{a, \ell, \gamma}) \le \sum_{\alpha = 1}^{a} V_{\alpha} + V_{0}^{\mathrm{cross}} + V_{0}^{\mathrm{loc}} $\\left[1mm] $\displaystyle \begin{aligned} & V_{\alpha} :\ \alpha\text{-overlap covariance terms},\quad 1 \le \alpha \le a,\\ & V_{0}^{\mathrm{cross}},\,V_{0}^{\mathrm{loc}} :\ \text{two zero-overlap covariance terms}. \end{aligned} $\\left[1mm] Lemma (ref) in Appendix (ref) }; \node[roadbox, below=0.55cm of s3] (s5) { \textbf{5. Leave-*-out expansion and “graph lemma”}\\left[1mm] $\begin{aligned} \widehat{\Omega} = B_{S}& - M_{S}, \ B_{S} = \widehat{\Omega}_{-S},\\ k_{a, \ell, \gamma} ({\mathbf{i}}) k_{a, \ell, \gamma} ({\mathbf{i}}') & = n^{-q} \prod_{e = (u, v) \in E(G)} X_{u}^{\top} B_{e,S} X_{v} . \end{aligned}$\\left[1mm] $\begin{aligned} | {\mathbb{E}} \prod_{e = (u, v) \in E(G)} X_{u}^{\top} B_{e} X_{v} | &\lesssim k^{\mathfrak{r} (G)}. \end{aligned}$\\left[1mm] Lemma (ref) (graph-counting lemma) and Lemma (ref) in Appendix (ref) }; \node[roadbox, right=0.55cm of s5] (s6) { \textbf{6. Counting and summation}\\left[1mm] $\displaystyle \begin{aligned} \mathrm{var} (T_{a, \ell, \gamma}) & \lesssim \frac{k^{2 \ell - 2 a + 4}}{n} \Gamma_{a, \ell, n},\\left[1mm] \mathrm{var} \{\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})\} & \lesssim \frac{j^{2}}{n} \frac{ \exp (C_{\eta} j^{2} k/n + C j^{2}/n)} {(1 - C j k/n)^2} \\ & \quad \times \Big(C j \frac{k}{n}\Big)^{2\lfloor (j-1) / 2 \rfloor}. \end{aligned} $\\left[1mm] Lemmas (ref)--(ref) }; \draw[arr] (s1.east) -- (s2.west); \draw[arr] (s3.east) -- (s4.west); \draw[arr] (s5.east) -- (s6.west); \draw[arr] (s2.south) -- ++(0,-0.28) -| (s3.north); \draw[arr] (s4.south) -- ++(0,-0.28) -| (s5.north); \end{tikzpicture} \caption{Schematic overview of the variance analysis. The diagram displays the main algebraic reductions, the covariance decomposition for the generic multiplicative-kernel $U$-statistic, the leave-*-out expansion followed by graph-counting bounds, and the final summation over collection levels and correction orders.} \end{figure} \paragraph{Step i.} We have the following result, which is a direct consequence of Minkowski's inequality. \begin{lemma} The following inequality holds. \begin{align} \mathrm{var}^{1 / 2} \Big\{ \sum_{j = 2}^{m} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) \Big\} \le \sum_{j = 2}^{m} \mathrm{var}^{1 / 2} \{\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})\}. \end{align} \end{lemma} Thus, by Lemma (ref), the variance analysis of $ \widehat{\psi}_{m, k}(\widehat{\Omega})$ reduces to bounding each fixed-order term $\widehat{\mathbb{IF}}_{j, j, k}(\widehat{\Omega})$, while keeping track of the dependence on $j$ (and eventually on $m$), $k$, and $n$. \paragraph{Step ii.} In this part, we recall all the notations defined in Section (ref). We bound the variance of $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$ by using the M\"{o}bius inversion decomposition (ref) presented in Lemma (ref): \begin{align*} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) = \sum_{{\mathcal{B}} \in {\mathbb{B}}_{\iota}} c_{{\mathcal{B}}, n} {\mathbb{U}}_{n, 2 + |{\mathcal{B}}|} (K_{{\mathcal{B}}}), \end{align*} where the form of $K_{{\mathcal{B}}}$ is recorded in (ref). For each $l \in [\iota]$ ($\iota = j - 2$), define \begin{equation*} {\mathcal{E}}_{l} ({\mathcal{B}}) \coloneqq \begin{cases} \{0, \nu\}, & l \in B_{\nu} \text{ for some } \nu \in \{1, \cdots, r\},\\ \{0, 1, 2\}, & l \notin \bigcup_{\nu = 1}^{r} B_{\nu}. \end{cases} \end{equation*} Here, the value $0$ corresponds to $- I$ from an $H = X X^{\top} \widehat{\Omega} - I$ or $- 2 I$ from a ${\mathcal{R}}_{1 2} = X_{1} X_{1}^{\top} \widehat{\Omega} + X_{2} X_{2}^{\top} \widehat{\Omega} - 2 I$. The nonzero values correspond to terms that involve $X_{u} X_{u}^{\top} \widehat{\Omega}$. Let ${\mathcal{E}} ({\mathcal{B}}) \coloneqq \prod_{l = 1}^{\iota} {\mathcal{E}}_{l} ({\mathcal{B}})$. For $\varepsilon = (\varepsilon_{1}, \cdots, \varepsilon_{\iota}) \in {\mathcal{E}} ({\mathcal{B}})$, define $d_{{\mathcal{B}}, \varepsilon} \coloneqq (-1)^{N_{H} (\varepsilon)} (-2)^{N_{R} (\varepsilon)}$, where \begin{equation*} N_{H} (\varepsilon) \coloneqq \Big| \Big\{ l \in \bigcup_{\nu = 1}^{r} B_{\nu}: \varepsilon_{l} = 0 \Big\} \Big|, \quad N_{R} (\varepsilon) \coloneqq \Big| \Big\{ l \notin \bigcup_{\nu = 1}^{r} B_{\nu}: \varepsilon_{l} = 0 \Big\} \Big|. \end{equation*} Now let \begin{equation*} {\mathcal{A}}_{\varepsilon} \coloneqq \left\{ \nu \in \{1, \cdots, r\}: \varepsilon_{l} = \nu \text{ for at least one } l \in B_{\nu}\right\}. \end{equation*} Thus ${\mathcal{A}}_{\varepsilon}$ records those $\nu$ for which the corresponding sample index $a_{\nu}$ appears through a term $X_{a_{\nu}} X_{a_{\nu}}^{\top} \widehat{\Omega}$. For $l \in [\iota] \setminus \bigcup_{\nu = 1}^{r} B_{\nu}$, any term of the form $X_{u} X_{u}^{\top} \widehat{\Omega}$ involves only $u = i_{1}$ or $u = i_{2}$. These two endpoint indices remain in the resulting kernel and are not included in ${\mathcal{A}}_{\varepsilon}$. Set $b_{\varepsilon} \coloneqq |{\mathcal{A}}_{\varepsilon}|$. Write $ {\mathcal{A}}_{\varepsilon} = \{\nu_{1}, \cdots, \nu_{b_{\varepsilon}}\}, \ \nu_{1} < \cdots < \nu_{b_{\varepsilon}}$. For every $\nu \notin {\mathcal{A}}_{\varepsilon}$, the index $a_{\nu}$ does not appear in the displayed kernel and can therefore be summed out exactly. After this summation, the original $U$-statistic of order $2 + r$ reduces to a $U$-statistic of order $ a_{\varepsilon} = 2 + b_{\varepsilon}$, with remaining displayed indices ordered as ${\mathbf{i}}_{a_{\varepsilon}} = (i_{1}, i_{2}, a_{\nu_{1}}, \cdots, a_{\nu_{b_{\varepsilon}}})$. Let $\ell_{\varepsilon} \coloneqq |\{l \in [\iota] : \varepsilon_{l} \neq 0\}|$. Writing the positions $l \in [\iota]$ with $\varepsilon_{l} \neq 0$ in increasing order defines an index assignment $\gamma_{\varepsilon} : \{1, \cdots, \ell_{\varepsilon} \} \to [a_{\varepsilon}]$: for the $h$-th non-identity position $l_{h}$, if $l_{h} \in B_{\nu_{t}}$ for some $t \in \{1, \cdots, b_{\varepsilon}\}$, then $\gamma_{\varepsilon} (h) = 2 + t$; otherwise, $\varepsilon_{l_{h}} \in \{1, 2\}$ and $\gamma_{\varepsilon} (h) = \varepsilon_{l_{h}}$. With the above preparation, we are ready to present the following lemma, which further decomposes ${\mathbb{U}}_{n, 2 + r} (K_{{\mathcal{B}}})$ into $U$-statistics with multiplicative-kernels. A proof can be found in Appendix (ref). \begin{lemma} Each summand in the M\"{o}bius inversion decomposition of $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$ in (ref) admits the following decomposition: \begin{equation*} {\mathbb{U}}_{n, 2 + r} (K_{{\mathcal{B}}}) = \sum_{\varepsilon \in {\mathcal{E}}({\mathcal{B}})} d_{{\mathcal{B}}, \varepsilon} T_{a_{\varepsilon}, \ell_{\varepsilon}, \gamma_{\varepsilon}}, \end{equation*} where, for $a \ge 2$, $\ell \ge 0$, and $\gamma : \{1, \cdots, \ell\} \to [a]$, \begin{equation*} T_{a, \ell, \gamma} = {\mathbb{U}}_{n, a} \Big\{ A_{1} X_{1}^{\top} \widehat{\Omega} \Big( \prod_{s = 1}^{\ell} X_{\gamma (s)} X_{\gamma (s)}^{\top} \widehat{\Omega} \Big) X_{2} Y_{2} \Big\}. \end{equation*} In particular, in $T_{a, \ell, \gamma}$, the following constraint holds: $\{1, 2, \gamma (1), \cdots, \gamma (\ell)\} = [a]$. Moreover, for every $\varepsilon \in {\mathcal{E}} ({\mathcal{B}})$, $\ell_{\varepsilon} - b_{\varepsilon} \le \iota - r$. For $T_{2 + b, q, \gamma}$, we also have: \begin{equation*} q - b \le \iota - r. \end{equation*} \end{lemma} By Lemma (ref), each summand ${\mathbb{U}}_{n, 2 + |{\mathcal{B}}|} (K_{{\mathcal{B}}})$ is a finite linear combination of $U$-statistics $T_{a, \ell, \gamma}$. It remains to control the variance of $T_{a, \ell, \gamma}$ uniformly in $(a, \ell, \gamma)$. Throughout the variance analysis, we use \begin{equation*} {\mathcal{I}}_{n, a} \coloneqq \{(i_{1}, \cdots, i_{a}) \in [n]^{a}: i_{s} \ne i_{t} \text{ for } s \ne t\} \end{equation*} to denote the set of ordered tuples of pairwise distinct sample indices. For ${\mathbf{i}} = (i_{1}, \cdots, i_{a}) \in {\mathcal{I}}_{n, a}$, write the corresponding kernel as \begin{equation*} k_{a, \ell, \gamma} (O_{{\mathbf{i}}}) \coloneqq A_{i_{1}} X_{i_{1}}^{\top} \widehat{\Omega} \Big( \prod_{s = 1}^{\ell} X_{i_{\gamma (s)}} X_{i_{\gamma (s)}}^{\top} \widehat{\Omega} \Big) X_{i_{2}} Y_{i_{2}}. \end{equation*} Here $\gamma : \{1, \cdots, \ell\} \to [a]$, and the condition $\{1, 2, \gamma (1), \cdots, \gamma (\ell)\} = [a]$ means that every entry of ${\mathbf{i}} = (i_{1}, \cdots, i_{a})$ appears in the kernel (when spelling out the $U$-statistic operator), either as one of the endpoint indices $i_{1}, i_{2}$ or through some $i_{\gamma (s)}$. For set operations, we write $\operatorname{ind} ({\mathbf{i}}) \coloneqq \{i_{1}, \cdots, i_{a}\}$ for the unordered set of sample indices appearing in the tuple ${\mathbf{i}}$. To control the variance of $T_{a, \ell, \gamma}$, we analyze the covariance between the kernels of $T_{a, \ell, \gamma}$ indexed by the ordered tuples ${\mathbf{i}} = (i_{1}, \cdots, i_{a})$ and ${\mathbf{i}}' = (i'_{1}, \cdots, i'_{a})$. We group the covariances by the number of the shared sample indices: \begin{equation*} \alpha ({\mathbf{i}}, {\mathbf{i}}') \coloneqq |\operatorname{ind} ({\mathbf{i}}) \cap \operatorname{ind} ({\mathbf{i}}')|. \end{equation*} Thus, $\mathrm{var} (T_{a, \ell, \gamma})$ decomposes into a summation of covariances indexed by $\alpha = 0, 1, \cdots, a$. More precisely, Lemma (ref) in Appendix (ref) bounds $\mathrm{var} (T_{a, \ell, \gamma})$ as follows: \begin{equation} \mathrm{var} (T_{a, \ell, \gamma}) \le \sum_{\alpha = 1}^{a} V_{\alpha} + V_{0}^{\mathrm{cross}} + V_{0}^{\mathrm{loc}}. \end{equation} The term $V_{\alpha}$ collects all covariances between kernels that share exactly $\alpha$ indices with $\alpha \geq 1$. The terms $V_{0}^{\mathrm{cross}}$ and $V_{0}^{\mathrm{loc}}$ collect the cases with $\alpha = 0$, and the two different terms arise from leave-*-out expansion of $V_{0}$, which we describe next. We next bound these terms by the graph-counting Lemma (ref). To this end, we first record the following result, which is proved in Appendix (ref). \begin{lemma} Given any $S \subseteq [n]$, define \begin{equation*} \widehat{\Sigma}_{-S} \coloneqq \widehat{\Sigma} - \frac{1}{n} \sum_{r \in S} X_{r} X_{r}^{\top}, \ B_{S} \coloneqq \widehat{\Omega}_{-S} \coloneqq \widehat{\Sigma}_{-S}^{-1}. \end{equation*} Then \begin{equation*} \widehat{\Omega} = B_{S} - M_{S},\ \text{where} \ M_{S} \coloneqq \sum_{q = 1}^{\infty} \frac{(-1)^{q - 1}}{n^{q}} \sum_{r_{1}, \cdots, r_{q}\in S} B_{S} X_{r_{1}} X_{r_{1}}^{\top} B_{S} X_{r_{2}} X_{r_{2}}^{\top} B_{S} \cdots X_{r_{q}} X_{r_{q}}^{\top} B_{S}. \end{equation*} \end{lemma} Fix a covariance pair indexed by ordered tuples ${\mathbf{i}}$ and ${\mathbf{i}}'$, and set $S \coloneqq \operatorname{ind} ({\mathbf{i}}) \cup \operatorname{ind} ({\mathbf{i}}')$. By Lemma (ref), expanding each occurrence of $\widehat{\Omega}$ around the leave-*-out inverse $B_{S} = \widehat{\Omega}_{-S}$ rewrites every expanded covariance term as \begin{equation*} n^{-q} \prod_{e = (u, v)\in E(G)} X_{u}^{\top} B_{e,S} X_{v}, \end{equation*} up to endpoint factors ($A$ and $Y$), where $q$ is the number of “inserted” $X X^{\top}$. Each insertion contributes a factor $n^{-1}$ and adds an edge to the associated undirected graph $G$. We now apply Lemma (ref), together with Lemma (ref), to the three types of covariances in (ref). Figure (ref) provides a graphical illustration of the three types of terms in (ref). The bounds for these three types of terms are proved in Lemma (ref) in Appendix (ref), but we provide some heuristic explanations below. \begin{figure} \tikzset{ vtx/.style={ circle, draw=blue, line width=0.65pt, minimum size=6.8mm, inner sep=0pt, font= }, sharedvtx/.style={ vtx, line width=0.85pt, fill=gray!10 }, gedge/.style={ draw=black!75, line width=0.65pt }, dashededge/.style={ draw=black, dashed, line width=0.75pt, dash pattern=on 3pt off 2pt }, elab/.style={ font=\scriptsize, fill=white, inner sep=1pt } } \captionsetup[subfigure]{justification=centering} \begin{subfigure}[t]{0.31\textwidth} \makebox[\linewidth][c]{ \resizebox{0.78\linewidth}{!}{ \begin{tikzpicture} \path[use as bounding box] (-0.45,-1.55) rectangle (2.75,1.55); \node[vtx] (a1) at (0.0,0.0) {$1$}; \node[sharedvtx] (a3) at (1.2,0.0) {$3$}; \node[vtx] (a2) at (1.2,-1.2) {$2$}; \node[vtx] (a4) at (1.2,1.2) {$4$}; \node[vtx] (a5) at (2.4,1.2) {$5$}; \draw[gedge] (a1) -- node[above, elab] {$a$} (a3); \draw[gedge] (a3) -- node[right, elab] {$b$} (a2); \draw[gedge] (a3) -- node[left, elab] {$c$} (a4); \draw[gedge] (a4) -- node[above, elab] {$d$} (a5); \end{tikzpicture} }} \caption{$\alpha \geq 1$} \resizebox{\linewidth}{!}{ \begin{tikzpicture} \node[align=left, font=\scriptsize, inner sep=1pt] {$ \begin{aligned} k_{a, \ell, \gamma} (O_{{\mathbf{i}}}) & = A_{1} X_{1}^{\top} \widehat{\Omega} (X_{3} X_{3}^{\top} \widehat{\Omega})X_{2} Y_{2}, \\ k_{a,\ell,\gamma'} (O_{{\mathbf{i}}'}) & = A_{3} X_{3}^{\top} \widehat{\Omega} (X_{4} X_{4}^{\top} \widehat{\Omega}) X_{5} Y_{5} . \end{aligned} $}; \end{tikzpicture} } \end{subfigure} \begin{subfigure}[t]{0.31\textwidth} \makebox[\linewidth][c]{ \resizebox{0.78\linewidth}{!}{ \begin{tikzpicture} \path[use as bounding box] (-0.45,-1.55) rectangle (2.75,1.55); \node[vtx] (b1) at (0.0,1.0) {$1$}; \node[vtx] (b2) at (1.2,1.0) {$2$}; \node[vtx] (b3) at (2.4,1.0) {$3$}; \node[vtx] (b4) at (0.0,-1.0) {$4$}; \node[vtx] (b5) at (1.2,-1.0) {$5$}; \node[vtx] (b6) at (2.4,-1.0) {$6$}; \draw[gedge] (b1) -- node[above, elab] {$a$} (b2); \draw[gedge] (b2) -- node[above, elab] {$b$} (b3); \draw[gedge] (b4) -- node[below, elab] {$c$} (b5); \draw[gedge] (b5) -- node[below, elab] {$d$} (b6); \draw[dashededge] (b2) -- node[right, elab] {$e$} (b5); \end{tikzpicture} }} \caption{$V_{0}^{\mathrm{cross}}$} \resizebox{\linewidth}{!}{ \begin{tikzpicture} \node[align=left, font=\scriptsize, inner sep=1pt] {$ \begin{aligned} k_{a, \ell, \gamma} (O_{{\mathbf{i}}}) & = A_{1} X_{1}^{\top} \widehat{\Omega} (X_{2} X_{2}^{\top} \widehat{\Omega}) X_{3} Y_{3}, \\ k_{a,\ell,\gamma'} (O_{{\mathbf{i}}'}) & = A_{4} X_{4}^{\top} \widehat{\Omega} (X_{5} X_{5}^{\top} \widehat{\Omega}) X_{6} Y_{6} . \end{aligned} $}; \end{tikzpicture} } \end{subfigure} \begin{subfigure}[t]{0.31\textwidth} \makebox[\linewidth][c]{ \resizebox{0.78\linewidth}{!}{ \begin{tikzpicture} \path[use as bounding box] (-0.45,-1.55) rectangle (2.75,1.55); \node[vtx] (c1) at (0.0,1.0) {$1$}; \node[vtx] (c2) at (1.2,1.0) {$2$}; \node[vtx] (c3) at (2.4,1.0) {$3$}; \node[vtx] (c4) at (0.0,-1.0) {$4$}; \node[vtx] (c5) at (1.2,-1.0) {$5$}; \node[vtx] (c6) at (2.4,-1.0) {$6$}; \node[vtx] (cr) at (1.2,0.0) {$r$}; \draw[gedge] (c1) -- node[above, elab] {$a$} (c2); \draw[gedge] (c2) -- node[above, elab] {$b$} (c3); \draw[gedge] (c4) -- node[below, elab] {$c$} (c5); \draw[gedge] (c5) -- node[below, elab] {$d$} (c6); \draw[dashededge] (c2) -- (cr); \draw[dashededge] (cr) -- (c5); \end{tikzpicture} }} \caption{$V_{0}^{\mathrm{loc}}$} \resizebox{\linewidth}{!}{ \begin{tikzpicture} \node[align=left, font=\scriptsize, inner sep=1pt] {$ \begin{aligned} k_{a, \ell, \gamma} (O_{{\mathbf{i}}}) & = A_{1} X_{1}^{\top} \widehat{\Omega} (X_{2} X_{2}^{\top} \widehat{\Omega}) X_{3} Y_{3}, \\ k_{a,\ell,\gamma'} (O_{{\mathbf{i}}'}) & = A_{4} X_{4}^{\top} \widehat{\Omega} (X_{5} X_{5}^{\top} \widehat{\Omega}) X_{6} Y_{6} . \end{aligned} $}; \end{tikzpicture} } \end{subfigure} \caption{Three graph structures in the covariance decomposition. In each panel, the digits on vertices denote the sample indices appearing in the two kernels $k_{a, \ell, \gamma} (O_{{\mathbf{i}}})$ and $k_{a, \ell, \gamma} (O_{{\mathbf{i}}'})$, and solid edges denote bilinear forms already present before leave-*-out expansion in these kernels. Panels (a)--(c) correspond respectively to $V_{\alpha}$ for $\alpha \geq 1$, $V_{0}^{\mathrm{cross}}$, and $V_{0}^{\mathrm{loc}}$. Dashed edges denote the additional edge introduced either by the leave-*-out expansion (for (b)) or by the introduction of an independent copy when applying the Efron--Stein inequality (for (c)). Below each panel, we exhibit the kernel-pair formulae corresponding to the undirected graphs.} \end{figure} \begin{enumerate}[label = (\roman*)] • \textit{$V_{\alpha}$ for $\alpha \ge 1$:} After replacing $\widehat{\Omega}$ by its leave-*-out expansion as in Lemma (ref), the shared indices ensure that the associated undirected graph is connected, as illustrated in Figure (ref). The pure leave-*-out term, in which every inverse is replaced by $B_{S}$, gives a connected graph. For this leading graph, \begin{equation*} e = 2 (\ell + 1), \ v = 2 a - \alpha,\ \kappa = 1. \end{equation*} Hence, Lemma (ref) gives the factor \begin{equation*} k^{\mathfrak{r}} = k^{2 (\ell + 1) - (2 a - \alpha) + 1}. \end{equation*} The remaining terms in the leave-*-out expansion insert additional $X X^{\top}$'s. Each such insertion adds an edge to the graph and contributes one factor $1 / n$ from the expansion; and hence, it leads to an additional factor of order $k / n$ after graph counting. Summing all insertion patterns only changes the bound by the factor $\Gamma^{\mathrm{ov}}_{a, \ell, n}$ of $O (1)$ depending on $a, \ell, n$ (see Lemma (ref) in Appendix (ref) for its explicit form). Therefore, \begin{equation*} \sum_{\alpha = 1}^{a} V_{\alpha} \lesssim \frac{k^{2 \ell- 2 a + 4}}{n} \Gamma^{\mathrm{ov}}_{a, \ell, n}. \end{equation*} • \textit{$V_{0}^{\mathrm{cross}}$:} When $\alpha = 0$, as in the leave-*-out expansion described in Lemma (ref), some terms contain explicit $X X^{\top}$-insertions that connect the two undirected graphs associated with the two kernels in the covariance, as illustrated in Figure (ref). The graph simply adds a new edge between existing vertices, and Lemma (ref) applies in the same way as in the case with $\alpha \geq 1$ just discussed. Summing over all such insertion patterns gives \begin{equation*} V_{0}^{\mathrm{cross}} \lesssim \frac{k^{2 \ell - 2 a + 4}}{n} \Gamma^{\mathrm{cross}}_{a, \ell, n}. \end{equation*} Here $\Gamma^{\mathrm{cross}}_{a, \ell, n}$ collects the connected insertion patterns and the geometric summation over their insertion orders; its explicit form is given in Lemma (ref) in Appendix (ref). • \textit{$V_{0}^{\mathrm{loc}}$:} $V_{0}^{\mathrm{loc}}$ collects the remaining terms in the case $\alpha = 0$ with the two graphs corresponding to the kernel pair not connected even after leave-*-out expansion. Conditional on $B_{S}$, the kernels are independent, so their covariance is reduced to the covariance of their conditional means. This term is controlled by the Efron--Stein inequality (Lemma (ref) in Appendix (ref)). When applying the Efron--Stein inequality, observations not in $S$ will be replaced by an independent copy, introducing a shared vertex that connects the originally disconnected graphs corresponding to the two kernels. We then apply the graph-counting Lemma (ref) to the newly connected graph (see Figure (ref) for an illustration). This yields \begin{equation*} V_{0}^{\mathrm{loc}} \lesssim \frac{k^{2 \ell - 2 a + 4}}{n} \Gamma^{\mathrm{loc}}_{a, \ell, n}, \end{equation*} where $\Gamma^{\mathrm{loc}}_{a, \ell, n}$ depends on $a, \ell$ and $\rho = k / n$; its explicit form is given in Lemma (ref). \end{enumerate} Combining the three contributions gives the generic multiplicative-kernel variance bound \begin{equation*} \mathrm{var} (T_{a, \ell, \gamma}) \lesssim \frac{k^{2 \ell - 2 a + 4}}{n} \Gamma_{a, \ell, n}, \quad \Gamma_{a, \ell, n} = \Gamma^{\mathrm{ov}}_{a, \ell, n} + \Gamma^{\mathrm{cross}}_{a, \ell, n} + \Gamma^{\mathrm{loc}}_{a, \ell, n}. \end{equation*} The exact form of $\Gamma_{a, \ell, n}$ is given in Lemma (ref) in Appendix (ref). \paragraph{Step iii.} We now combine the M\"{o}bius-inversion expansion in Lemma (ref) with the generic multiplicative-kernel bound obtained in \textbf{Step ii}. Let $r_{\iota}^{*} \coloneqq \lfloor \frac{\iota}{2} \rfloor$. For $0 \le r \le r_{\iota}^{*}$, let $ {\mathbb{B}}_{\iota,r} \coloneqq \{{\mathcal{B}} \in {\mathbb{B}}_{\iota}: |{\mathcal{B}}| = r\}$. Equivalently, $r$ is the number of sets in the collection ${\mathcal{B}}$. The M\"obius-inversion expansion of $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$ can then be rewritten as \begin{equation*} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) = \sum_{r = 0}^{r_{\iota}^{*}} Z_{\iota, r}, \quad Z_{\iota, r} \coloneqq \sum_{{\mathcal{B}} \in {\mathbb{B}}_{\iota, r}} c_{{\mathcal{B}}, n} {\mathbb{U}}_{n, 2 + r} (K_{{\mathcal{B}}}) \end{equation*} The maximal possible value of $r$ is $r_{\iota}^{*}$ because every set $B_{\nu}$ in ${\mathcal{B}}$ has cardinality of at least two. The next lemma first controls the contribution from a given $r$. \begin{lemma} Let $j \ge 3$, and for $0 \le r \le r_{\iota}^{*}$, let $\Gamma_{\iota, r, n} \coloneqq \max_{\substack{0 \le b \le r\\ 0 \le q \le \iota}} \Gamma_{2 + b, q, n}$. When $n \ge 2 (\iota + 2)$, \begin{equation*} \mathrm{var} (Z_{\iota, r}) \lesssim 2^{2 (\iota - r)} w_{\iota, r}^{2} \Gamma_{\iota, r, n} \frac{k^{2 (\iota - r)}}{n^{2 (\iota - r) + 1}}. \end{equation*} $w_{\iota, r}$ is defined as follows. For ${\mathcal{B}} \in {\mathbb{B}}_{\iota, r}$, define $D ({\mathcal{B}}) \coloneqq \sum_{\nu = 1}^{r} |B_{\nu}|$, with the convention $D (\emptyset) = 0$. Then: \begin{equation*} w_{\iota, r} \coloneqq \sum_{{\mathcal{B}} \in {\mathbb{B}}_{\iota, r}} \Big\{ \prod_{B \in {\mathcal{B}}}(|B| - 1) \Big\} 2^{D ({\mathcal{B}})} 4^{\iota - D ({\mathcal{B}})}. \end{equation*} \end{lemma} In Lemma (ref), $w_{\iota, r}$ bounds the absolute sum of the coefficients in the expansion of all level-$r$ terms $c_{{\mathcal{B}} ,n} {\mathbb{U}}_{n, 2 + r} (K_{{\mathcal{B}}})$, ${\mathcal{B}} \in {\mathbb{B}}_{\iota, r}$, into multiplicative-kernel $U$-statistics, up to the common factor $(n - 2 - \iota) ! / (n - 2 - r) !$; see Lemma (ref) in Appendix (ref). The M\"obius-inversion expansion of $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$ can also be represented by $Z_{\iota, r}$'s: $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) = \sum_{r = 0}^{r_{\iota}^{*}} Z_{\iota, r}$. The next lemma gives the variance bound of $\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$ after summing over $\mathrm{var} (Z_{\iota, r})$ for $r = 0, \cdots, r_{\iota}^{*}$. \begin{lemma} Let $j \ge 3$. Suppose that $n \ge 2 j$ and $C j k/n \le \eta < 1$. Then \begin{equation*} \mathrm{var} \{\widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})\} \lesssim \frac{j^{2}}{n} \frac{ \exp (C_{\eta} j^{2} k/n + C j^{2}/n)} {(1 - C j k/n)^2} \Big(C j \frac{k}{n}\Big)^{2 \lfloor (j - 1) / 2\rfloor}. \end{equation*} \end{lemma} The proofs of Lemmas (ref) and (ref) are deferred to Appendix (ref). The term $j = 2$ is controlled by Lemma (ref) in Appendix (ref), which gives \begin{equation*} \mathrm{var}^{1 / 2}\{\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})\} \lesssim \frac{1}{\sqrt n} \Big(1 + \frac{k}{n}\Big)^{1 / 2}. \end{equation*} Therefore, \begin{align*} \mathrm{var}^{1 / 2} \Big\{ \sum_{j = 2}^{m} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) \Big\} \lesssim \frac{1}{\sqrt n} \Big( 1 + \frac{k}{n} \Big)^{1 / 2} + \frac{1}{\sqrt n} \sum_{j = 3}^{m} j \frac{\exp (C_{\eta} j^{2} k/n + C j^{2}/n)}{1 - C j k/n} \Big( C j\frac{k}{n} \Big)^{\lfloor (j - 1) / 2 \rfloor}. \end{align*} Recall that $\rho = k / n$. Under $C m \rho < 1$, pairing adjacent orders $j = 2 \ell + 1$ and $j = 2 \ell + 2$ yields \begin{align*} \sum_{j = 3}^{m} j \frac{ \exp (C_{\eta} j^{2} \rho + C j^{2} / n)}{1 - C j \rho} (C j \rho)^{\lfloor (j - 1) / 2 \rfloor} & \lesssim \frac{\exp (C_{\eta} m^{2} \rho + C m^{2} / n)}{1 - C m \rho} \sum_{\ell = 1}^{\lfloor (m - 1) / 2 \rfloor} \ell (C \ell \rho)^{\ell} \\ & \lesssim \frac{\rho \exp (C_{\eta} m^{2} \rho + C m^{2} / n)}{(1 - C m \rho)^{4}}. \end{align*} Consequently, \begin{equation*} \mathrm{var}^{1 / 2} \Big\{ \sum_{j = 2}^{m} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) \Big\} \lesssim \frac{1}{\sqrt{n}} \Big\{ (1 + \rho)^{1 / 2} + \rho \frac{\exp (C_{\eta} m^{2} \rho + C m^{2} / n)}{(1 - C m \rho)^{4}} \Big\}, \end{equation*} which immediately implies that \begin{equation*} \mathrm{var} \Big\{ \sum_{j = 2}^{m} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) \Big\} \lesssim \frac{1}{n} \Big\{ 1 + \frac{k}{n} \frac{\exp \{(C_{\eta} m^{2} k + C m^{2}) / n\}}{(1 - C m \rho)^{4}} \Big\}^2. \end{equation*} The proof of the variance bound is now complete. \section{Concluding Remarks} We conclude our article by mentioning several future research directions. \begin{enumerate}[label = (\arabic*)] • It will be interesting to study if one can extend the idea developed in this article to the case $k \gtrsim n$ by, for instance, estimating $\Omega$ via shrinkage or regularized methods. As conjectured in robins2016technical, the optimal convergence rate of the functionals studied in this article may depend on the regularity of the density of $X$. It is then reasonable to conjecture that the shrinkage or regularization also depends on the density of $X$. Simulation studies in liu2020nearly suggest the nonlinear shrinkage covariance matrix estimators ledoit2012nonlinear, ledoit2020analytical could be a viable option. It will also be interesting to investigate the statistical theoretical guarantees when $\Omega$ is estimated by the inverse of the ridge penalized estimator cheng2024dimension in the proportional asymptotic regime ($k \asymp n$) chen2024method. • We expect to see the analysis strategy developed here to be further generalized to more complex problems, such as assumption-lean estimands vansteelandt2022assumption, vansteelandt2025towards, functionals of NPIV models breunig2024adaptive, functionals beyond bilinear forms lin2024worthwhile, zhang2026higher, moment-condition models \citetext{bonhomme2026higher; robins2016technical}, multi-index models damian2025generative, joshi2026learning, and other related problems wein2019kikuchi, lasserre2024moment, liu2025quantum. \end{enumerate} \section*{Acknowledgments} Lin Liu thanks the Isaac Newton Institute (INI) of Mathematical Sciences at the University of Cambridge, the School of Mathematics and Statistics at the University College Dublin, and the Center of Data Science at Zhejiang University for hospitality during the completion of this work. The authors thank Rohit Bhattacharya, Kwun Chuen Gary Chan, Fengnan Gao, Zhenyu Liao, Rajarshi Mukherjee, Jamie Robins, Andrea Rotnitzky, Eric Tchetgen Tchetgen, Aad van der Vaart, Cheng Wang, and participants in the \href{https://www.newton.ac.uk/event/cifw04/}{Causality and Machine Learning Workshop} held at INI for helpful discussions. This research is supported by the National Key R&D Program of China Project Number 2025YFA1016700, NSFC Grant No.12471274, and Science and Technology Talent and Platform Program of Yunnan Province Grant No.202605AF35007. \putbib[Master.bib]

bibunit[plainnat]