The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
100,490 characters
New $n$-consistent, numerically stable higher-order influence function estimators
\title{New $\sqrt{n}$-consistent, numerically stable higher-order influence function estimators}
\author{Lin Liu\thanks{Correspondence: \href{[email removed]}{[email removed]}. The authors would like to thank \href{https://gaofn.xyz/}{Fengnan Gao}, \href{https://zhenyu-liao.github.io/}{Zhenyu Liao}, \href{https://scholar.harvard.edu/rajarshi/home}{Rajarshi Mukherjee}, \href{https://www.hsph.harvard.edu/james-robins/}{Jamie Robins}, and \href{https://sites.google.com/view/zheng-zhang}{Zheng Zhang} for invaluable discussions on this paper. The authors gratefully acknowledges funding support by NSFC Grant No.12101397 and No.12090024, Shanghai Municipal Science and Technology Grant No.2021SHZDZX0102, Shanghai Science and Technology Commission Grant No.21JC1402900, Shanghai Natural Science Foundation Grant No.21ZR1431000.}}\affil[2]{Institute of Natural Sciences, MOE-LSC, School of Mathematical Sciences, CMA-Shanghai, SJTU-Yale Joint Center for Biostatistics and Data Science, Shanghai Jiao Tong University and Shanghai Artificial Intelligence Laboratory, Shanghai, China}
\author{Chang Li}\affil[1]{Department of Statistics, University of Virginia, Charlottesville, VA, USA}
\date{\today}
\maketitle
\begin{abstract}
Higher-Order Influence Functions (HOIFs) provide a unified theory for constructing rate-optimal estimators for a large class of low-dimensional (smooth) statistical functionals/parameters (and sometimes even infinite-dimensional functions) that arise in substantive fields including epidemiology, economics, and the social sciences. Since the introduction of HOIFs by \citet{robins2008higher} or \citet{robins2016technical}\footnote{\citet{robins2016technical} is the complete version of \citet{robins2008higher}, including more results and proofs. We therefore only refer to \citet{robins2016technical} in the sequel.}, they have been viewed mostly as a theoretical benchmark rather than a useful tool for statistical practice. Works aimed to flip the script are scant, but a few recent papers \citet{liu2017semiparametric, liu2021assumption} make some partial progress. In this paper, we take a fresh attempt at achieving this goal by constructing new, numerically stable HOIF estimators (or sHOIF estimators for short with ``s'' standing for ``stable'') with provable statistical and computational guarantees. This new class of sHOIF estimators (up to the 2nd order) was foreshadowed in synthetic experiments conducted by \citet{liu2020nearly}.
\end{abstract}
{\footnotesize \textbf{Keywords:} Causal Inference, Functional Estimation, Higher-Order Influence Functions, Semiparametric Theory, Combinatorics}
\newpage
\section{Introduction}
\label{sec:intro}
{\it Higher-Order Influence Functions} (HOIFs) \citep{robins2016technical} are higher-order generalizations of the first-order influence functions (IFs), a staple in semiparametric statistical theory \citep{newey1990semiparametric, bickel1998efficient, van2002part}. HOIFs are a powerful and unified approach to constructing minimax rate-optimal estimators for a class of statistical functionals/parameters (and sometimes even functions; see \citet{kennedy2022minimax}) that arise in (bio)statistics, epidemiology, economics, and the social sciences. HOIF estimators originally proposed in \citet{robins2016technical, robins2017minimax}\footnote{See \citet{robins2022corrigenda} for corrections of the proofs in \citet{robins2017minimax}.} remain the only known minimax rate-optimal estimators for statistical functionals/parameters with substantive interests in the above disciplines, including the Average Treatment Effect (ATE) under the strong ignorability assumption\footnote{In \citet{liu2021assumption}, we derived the HOIFs for the ATE functional even when the strong ignorability assumption fails to hold, provided that we have access to valid proxies for both the treatment and outcome, following a series of works on proximal causal learning \citep{tchetgen2020introduction}.} and the expected conditional covariance of two random variables $A$ and $Y$ given a third random variable $X$, even after highly active research by the statistics and econometrics communities in recent years \citep{newey2018cross, kennedy2020optimal, hirshberg2021augmented, yu2020treatment}. More recent works \citep{kennedy2022minimax, bonvini2022fast} also initiated the application of the HOIF machinery to the minimax optimal estimation of Conditional Average Treatment Effect (CATE) function or dose response curves. Their results lay important theoretical foundation for individualized decision making problems, e.g. personalized medicine. This is the first instance when HOIF estimators are also shown to be effective, at least in theory, for function estimation problems, or more precisely, ``hybrid function and functional estimation problems''. Similar idea has also been applied to dose-response curve estimation \citep{bonvini2022fast}. For an introductory level review of HOIFs, we refer the interested readers to \citet{van2014higher} and Section 1 of \citet{liu2020rejoinder}. A relatively more technical review of HOIFs is delegated to Section \ref{sec:review}.
Over the past decade, Robins and colleagues initiated the research program of establishing theoretical foundations for HOIFs and estimators based on HOIFs \citep{robins2004optimal, van2014higher, robins2016technical, robins2017minimax, liu2017semiparametric} for a class of statistical functionals/parameters recently characterized in \citet{rotnitzky2021characterization}, which are heretofore termed as {\it Doubly Robust Functionals} (DRF) in this paper. We adopt this terminology to reflect the fact that their nonparametric first-order IFs give rise to doubly robust estimators \citep{scharfstein1999adjusting, robins2001comments, chernozhukov2018double}. This class of DRFs subsumes the class of functionals studied in \cite{robins2016technical} and \cite{chernozhukov2018riesz}. Under the standard \text{H\"{o}lder}{}-regularity assumptions on the nuisance parameters (abbreviated as \text{H\"{o}lder}{} nuisance models), \citet{robins2017minimax} constructed minimax optimal but non-adaptive HOIF estimators for a sub-class of DRFs. \citet{liu2021adaptive} constructed adaptive second-order IF estimators for DRFs using the celebrated Lepski\v{i}'s adaptation scheme \citep{lepskii1991problem}, within a strict submodel of the \text{H\"{o}lder}{} nuisance models. But both estimators require estimating the density of the potentially high-dimensional covariates $X$, even in $\sqrt{n}$-estimable regimes. When the dimension $d$ of the covariates is only moderately large (e.g. $d = 10$), nonparametric density estimation is already a daunting computational and statistical task.
To overcome the above issues, \citet{liu2017semiparametric} introduced {\it empirical} HOIF (eHOIF for short) estimators that obviate multi-dimensional density estimation by inverting the sample/empirical Gram matrix of vector-valued basis transformation of the covariates computed using a separate sample independent of the sample used to construct the estimator of the DRF. This sample-splitting strategy is adopted mainly for simplifying the mathematical analysis, leading to rather straightforward analysis of the statistical properties of the eHOIF estimators. In particular, the eHOIF estimators are still the only class of estimators that achieves $\sqrt{n}$-consistency and semiparametric efficiency for DRFs under the minimal \text{H\"{o}lder}{}-regularity assumptions \citep{robins2009semiparametric}. These nice statistical properties of the eHOIF estimators also motivate the development of a class of assumption-lean hypothesis tests statistic that is designed to falsify if the standard $(1 - \alpha) \times 100\%$ Wald confidence interval of the DRF has the claimed coverage probability \citep{liu2020nearly, liu2021assumption}. At this point, astute readers must wonder why we need a new class of empirical HOIF estimators at all, which is what this article is all about.
\subsection{Motivation and main contributions}
\label{sec:novelty}
Despite the effort in \citet{liu2017semiparametric}, from our past experience of using eHOIF estimators in practice \citep{liu2017semiparametric, liu2020nearly, liu2021assumption, wanis2023machine}, several singular issues of their finite-sample performance were unveiled by large-scale simulation experiments\footnote{For interested readers, these simulation experiments have also been used to expose the gap between the (nonparametric) statistical theory deep neural networks (DNNs) and their practice in \citet{xu2022deepmed}. One can access computer codes of generating such simulations \href{https://github.com/siqixu/DeepMed}{here}.}:
\begin{enumerate}[label = (\roman*)]
\item \underline{Numerical instability}: In \citet{liu2017semiparametric}, although eHOIF estimators exhibit better finite-sample performance than the original HOIF estimators in \citet{robins2017minimax}, the simulations were restricted to very low condition number $k / n$: e.g. $k \approx 1,000$ and $n \approx 10,000$. In the simulation studies of \citet{liu2020nearly}, when $k$ gets near $n$, eHOIF estimators blow up numerically (see Section S3.1 of \citet{liu2020nearly}) already at order two. What is more striking is that the eHOIF estimators at higher orders, though supposed to be correcting the bias, can only exacerbate the numeric blow-up.
\item \underline{Non-monotone bias reduction}: Theoretical results in \citet{liu2017semiparametric} hint that increasing the orders of the estimator should {\it in principle} reduce the bias. However, we found that this is not usually the case for eHOIF estimators in practice (e.g. see Section 5 of \citep{liu2021assumption}). Interestingly, sHOIF estimators do not seem to suffer from this problem in simulations, elevating the theoretical results from {\it mere principles} closer to {\it empirical facts}; see \citet{liu2020nearly} or \citet{wanis2023machine} for simulations at orders 2 or 3.
\end{enumerate}
Our contributions are three-fold.
\begin{itemize}
\item \underline{Methodology and practical relevance}: This article proposes a new class of numerically stable sHOIF estimators for DRFs, that overcomes the above two major limitations of eHOIF estimators. The stable Second-Order IF (SOIF) estimators first appeared in the simulation studies of \citet{liu2020nearly}, but their statistical properties remain elusive.
\item \underline{Theory and the proof strategy}: Obtaining a deeper theoretical underpinning of this phenomenon mandates meticulous calculations rather than crude upper bounds. This is the critical technical innovation vis-\`{a}-vis other HOIF-related works. In particular, we intensively use the following proof techniques: leave-out analysis, matrix-valued Taylor expansion, and combinatorial calculations (i.e. corollaries of the binomial identity). The proof strategy developed in this paper may be of independent interest.
\item \underline{Extensions of sHOIFs beyond ATE settings}: We also generalize sHOIF estimators to all the DRFs, allowing us to handle more structural parameters in the current causal inference (or econometrics) literature.
\end{itemize}
\subsection{Notation}
\label{sec:notation}
Before proceeding, we gather some frequently used notation throughout the paper. We denote the observed data random vector as $O \in {\mathcal{O}}$, where ${\mathcal{O}}$ is its corresponding sample space. Let $\bar{{\mathsf{z}}}_{k} \coloneqq (z_{1}, \cdots, z_{k})^{\top}$ denote a collection of $k$ different functions, each of which has input domain ${\mathcal{X}}$. Fix some $\theta' \in \Theta$. ${\mathbb{E}}_{\theta'}$, $\mathsf{var}_{\theta'}$, and $\mathsf{cov}_{\theta'}$ are, respectively, the expectation, variance, and covariance operators under the probability law ${\mathbb{P}}_{\theta'}$. For any measurable function $h: {\mathcal{X}} \rightarrow {\mathbb{R}}$, let $\Vert h \Vert_{\infty} \coloneqq \mathrm{ess } \sup_{x \in {\mathcal{X}}} h (x)$ and $\Vert h \Vert_{\theta', p} \coloneqq \left\{ {\mathbb{E}}_{\theta'} \left[ h (X)^{p} \right] \right\}^{1 / p}$ for any $p \geq 1$. We adopt standard (stochastic) asymptotic notation $\lesssim$, $\gtrsim$, $\asymp$, $\gg$, $\ll$, $o (\cdot)$, $\omega (\cdot)$, $O (\cdot)$, $\Omega (\cdot)$, $o_{{\mathbb{P}}_{\theta'}} (\cdot)$, $O_{{\mathbb{P}}_{\theta'}} (\cdot)$. For any real-valued vector $\bar{v}$ and any $q \in {\mathbb{R}}$, let $v^{q}$ be the element-wise $q$-th power of $v$.
Furthermore, define $\Sigma_{\theta'} \coloneqq {\mathbb{E}}_{\theta'} [Q]$, where $Q \coloneqq A \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}$, as the ($A$-weighted) population Gram matrix of $\bar{{\mathsf{z}}}_{k} (X)$, till Section \ref{sec:drf}, in which we generalize all our results from $\psi (\theta) \coloneqq {\mathbb{E}}_{\theta} [Y (a = 1)]$\footnote{We use the potential outcome notation without introducing it, which will not affect the understanding of the main theme of this work.}, the mean of outcome $Y$ in the treated group under strong ignorability, to all members of the DRFs. Similarly, define $\widehat{\Sigma} \coloneqq {\mathbb{P}}_{n} [Q_{k}] \equiv n^{-1} \sum_{i = 1}^{n} Q_{i}$ as the ($A$-weighted) sample Gram matrix of $\bar{{\mathsf{z}}}_{k} (X)$, again till Section \ref{sec:drf}. Here ${\mathbb{P}}_{n} [\cdot]$ denotes the sample mean operator.
To further lighten the notation, we let $Q_{I} \coloneqq \sum_{i \in I} Q_{i}$ for any multi-index set $I \subseteq [n]$. For convenience, we also denote multi-index set $\{i_{1}, i_{2}, \cdots, i_{j}\} \subseteq [n]$ as $\bar{i}_{j}$ for $j \leq n$. $\Omega_{\theta'} \equiv \Sigma_{\theta'}^{-1}$ and $\widehat{\Omega} \equiv \widehat{\Sigma}^{-1}$, when they exist, are respectively the inverse of the population and sample Gram matrices. The kernels constructed from $\bar{{\mathsf{z}}}_{k}$ are denoted as $K_{\theta', k} (x, x') \coloneqq \bar{{\mathsf{z}}}_{k} (x)^{\top} \Omega_{\theta'} \bar{{\mathsf{z}}}_{k} (x')$ and $\widehat{K}_{k} (x, x') = \bar{{\mathsf{z}}}_{k} (x)^{\top} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (x')$. Given a set of functions $\bar{{\mathsf{v}}}_{k}: {\mathcal{X}} \rightarrow {\mathbb{R}}^{k}$ and any $L_{2}$ function $h: {\mathcal{O}} \rightarrow {\mathbb{R}}$, $\Pi [h | \bar{{\mathsf{v}}}_{k}]$ denotes the linear projection operator of projecting $h$ onto the linear span of $\bar{{\mathsf{v}}}_{k}$: formally,
$$
\Pi [h | \bar{{\mathsf{v}}}_{k}] (x) \coloneqq \bar{{\mathsf{v}}}_{k} (x)^{\top} {\mathbb{E}} [\bar{{\mathsf{v}}}_{k} (X) \bar{{\mathsf{v}}}_{k} (X)^{\top}]^{-1} {\mathbb{E}} [\bar{{\mathsf{v}}}_{k} (X) h (O)].
$$
We use ${\mathbb{P}}_{\theta}$ to denote the true data generating law, unless stated otherwise. When the reference measure is the true law ${\mathbb{P}}_{\theta}$, we often drop the dependence on $\theta$: for example, we write $\Vert h \Vert_{p} \equiv \Vert h \Vert_{\theta, p}$, $\Sigma \equiv \Sigma_{\theta}$, $\Omega \equiv \Omega_{\theta}$ and ${\mathbb{P}}$, ${\mathbb{E}}$, $\mathsf{var}$, $\mathsf{cov}$ correspond to ${\mathbb{P}}_{\theta}$, ${\mathbb{E}}_{\theta}$, $\mathsf{var}_{\theta}$, $\mathsf{cov}_{\theta}$. Note that $\Omega$ should not be confused with the asymptotic notation $\Omega (\cdot)$ and this will be clear from the context. A statistic is said to be ``oracle'' whenever it depends on some part(s) of the unknown {\it true} data generating law ${\mathbb{P}}_{\theta}$ (such as $\Omega$); otherwise it is said to be ``feasible''. We also introduce $\mathsf{Diag}$ as the operator of extracting the diagonal elements of a matrix.
Finally, let ${\mathbb{U}}_{n, m} [\cdot]$ denote the $m$-th order $U$-statistic operator: for any function $h: {\mathbb{R}}^{m} \rightarrow {\mathbb{R}}$
\begin{align*}
{\mathbb{U}}_{n, m} [h (O_{1}, \cdots, O_{m})] \coloneqq \frac{(n - m)!}{n!} \sum_{1 \leq i_{1} \neq \cdots \neq i_{m} \leq n} h (O_{i_{1}}, \cdots, O_{i_{m}}).
\end{align*}
When $m = 1$, ${\mathbb{U}}_{n, m} [\cdot]$ reduces to the sample mean operator ${\mathbb{P}}_{n} [\cdot]$. Similarly, let ${\mathbb{V}}_{n, m}$ be the corresponding $V$-statistic operator\footnote{Here we use the scaling $\frac{(n - m)!}{n!}$ instead of the more conventional $\frac{1}{n^{m}}$ for notational convenience.}:
\begin{align*}
{\mathbb{V}}_{n, m} [h (O_{1}, \cdots, O_{m})] \coloneqq \frac{(n - m)!}{n!} \sum_{i_{1} = 1}^{n} \cdots \sum_{i_{m} = 1}^{n} h (O_{i_{1}}, \cdots, O_{i_{m}}).
\end{align*}
Later in the paper, for $m \geq 2$, we will define ``oracle'' $m$-th order influence function estimators constructed using the dictionary $\bar{{\mathsf{z}}}_{k}$, denoted as $\widehat{\mathbb{IF}}_{m, m, k} (\Omega) \equiv \widehat{\mathbb{IF}}_{m, m, k} = {\mathbb{U}}_{n, m} [\widehat{\mathsf{IF}}_{m, m, k, \bar{i}_{m}}]$ with $U$-statistic kernel $\mathsf{IF}_{m, m, k, \bar{i}_{m}} \equiv \widehat{\mathsf{IF}}_{m, m, k, \bar{i}_{m}} (\Omega)$. Its stable feasible version is denoted by $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ with the corresponding kernel $\widehat{\mathsf{IF}}_{m, m, k, \bar{i}_{m}} (\widehat{\Omega})$.
\subsection{The setup and a review of the theory of HOIFs}
\label{sec:review}
With the notation just introduced, we are poised to state the problem setup and briefly review the theory of HOIFs relevant for this paper, in particular the theory of eHOIFs.
Suppose that we are given $N$ i.i.d. observations $\{O_{i}\}_{i = 1}^{N} \sim {\mathbb{P}}_{\theta}$, where $\theta \in \Theta$ is the so-called nuisance parameter and $\Theta$ is its underlying parameter space. Let ${\mathcal{P}} \coloneqq \left\{ {\mathbb{P}}_{\theta}: \theta \in \Theta \right\}$ be the space of data generating probability measures. Our primary interest is to estimate and draw statistical inference on a smooth statistical functional $\psi (\theta): \rightarrow {\mathbb{R}}$, in the sense of \citet{van1991differentiable}. We restrict $\psi (\theta)$ to be the DRFs defined in \citet{rotnitzky2021characterization}. As mentioned, our running example is $\psi (\theta) = {\mathbb{E}}_{\theta} [Y (a = 1)]$ the mean of an outcome $Y$ in the treated group $A = 1$. Here the observed data specializes to $O = (X, A, Y)$: respectively the $d$-dimensional covariates belonging to a compact subset ${\mathcal{X}} \equiv [-B, B]^{d}$ of ${\mathbb{R}}^{d}$, the binary treatment assignment, and the bounded outcome variable. Under unconfoundedness assumption (that can be relaxed by using the HOIFs of $\psi (\theta)$ under the proximal causal inference setting \citep{liu2021assumption}), $\psi (\theta)$ can be identified by either of the two statistical functionals of the observed data distribution:
\begin{equation}
\label{fnl}
\psi (\theta) \equiv {\mathbb{E}} \left[ A a (X) Y \right] \equiv {\mathbb{E}} [b (X)]
\end{equation}
where $a (x) \coloneqq \{{\mathbb{E}} [A | X = x]\}^{-1}$ and $b (x) \coloneqq {\mathbb{E}} [Y | X = x, A = 1]$ except Section \ref{sec:drf}. For this functional $\psi (\theta)$, the nuisance parameter is $\theta \equiv (a, b, g)$ where $g (x)$ is the probability density/mass function of the covariates $X$ conditional on $A = 1$. Hence the nuisance parameter space $\Theta = {\mathcal{A}} \times {\mathcal{B}} \times {\mathcal{G}}$, where ${\mathcal{A}}, {\mathcal{B}}, {\mathcal{G}}$ are, respectively, the space where $a, b, g$ lie. We further divide the whole $N$ data points into two parts: one with sample size $n$, called the estimation sample, and the other with sample size $N - n$, called the nuisance sample used to estimate the nuisance parameter $\theta$. Throughout this paper, we condition on the nuisance sample data by treating it or any quantity computed from it as fixed.
For a smooth statistical functional $\psi (\theta): \Theta \rightarrow {\mathbb{R}}$ in the sense of \citet{van1991differentiable}, its first-order influence function $\mathbb{IF}_{1} (\theta)$ is a mean-zero first-order $U$-statistic satisfying the following functional equation
\begin{align*}
\left. \frac{{\mathrm d} \psi (\theta_{t})}{{\mathrm d} t} \right\vert_{t = 0} = {\mathbb{E}} \left[ \mathbb{IF}_{1} (\theta) \cdot {\mathbb{S}}_{1} \right]
\end{align*}
where ${\mathbb{P}}_{\theta_{t}}$ is any parametric submodels in $\{{\mathbb{P}}_{\theta}, \theta \in \Theta\}$, such that when $t = 0$, ${\mathbb{P}}_{\theta_{t}} \equiv {\mathbb{P}}_{\theta}$, the true data generating law, and ${\mathbb{S}}_{1}$ is its first-order score vector, as defined in \citet{waterman1996projected}; also see \citet{robins2016technical}. Here $\mathbb{IF}_{1} (\theta)$ has the following form \citep{robins1994estimation}:
\begin{equation}
\label{if1}
\mathbb{IF}_{1} (\theta) \equiv \frac{1}{n} \sum_{i = 1}^{n} \mathsf{IF}_{1, i} (\theta), \text{ where } \mathsf{IF}_{1} (\theta) = A a (X) (Y - b (X)) + b (X) - \psi (\theta).
\end{equation}
Typically, classical semiparametric theory \citep{newey1990semiparametric, bickel1998efficient} constructs semiparametric efficient first-order estimators $\widehat{\psi}_{1}$ of $\psi (\theta)$ based on its first-order influence function follows:
$$
\widehat{\psi}_{1} = \frac{1}{n} \sum_{i = 1}^{n} A_{i} \widehat{a} (X_{i}) (Y_{i} - \widehat{b} (X_{i})) + \widehat{b} (X_{i})
$$
where $\widehat{a}, \widehat{b}$ are nuisance parameter estimates computed from the nuisance sample. In particular, $\widehat{\psi}_{1}$ has bias
\begin{equation}
\label{first-order bias}
\begin{split}
\mathsf{bias} (\widehat{\psi}_{1}) & = {\mathbb{E}} [\widehat{\psi}_{1} - \psi (\theta)] = {\mathbb{E}} [\mathsf{IF}_{1} (\widehat{\theta}) - \mathsf{IF}_{1} (\theta)] \\
& = {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X)}{a (X)} - 1 \right) (b (X) - \widehat{b} (X)) \right].
\end{split}
\end{equation}
Formally, $\mathsf{bias} (\widehat{\psi}_{1})$ is a {\it product of two nuisance estimation errors}\footnote{\citet{rotnitzky2021characterization} actually define the general class of statistical functionals that permit doubly-robust estimators based on this second-order bias property; see Section \ref{sec:drf}.}, and hence {\it doubly-robust} \citep{scharfstein1999adjusting}.
Despite being doubly-robust, the veracity of inference based on first-order estimators like $\widehat{\psi}_{1}$ may nonetheless be questionable when the nuisance parameter $\theta$ is of high complexity: e.g. functions with low smoothness or without sparsity. For example, when $a, b$ belong to \text{H\"{o}lder}{} functions with smoothness $s_{a}, s_{b}$ and $g$ arbitrarily complex, by far no first-order estimators are known to be $\sqrt{n}$-consistency for estimating $\psi (\theta)$ throughout the entire range
\begin{equation}
\label{minimal}
\{(s_{a}, s_{b}): (s_{a} + s_{b}) / 2 \geq d / 4\}
\end{equation}
but the eHOIF estimators of \citet{liu2017semiparametric} or the original HOIF estimators of \citet{robins2016technical} if additionally assuming $g$ to be \text{H\"{o}lder}{} with smoothness $s_{g} > 0$. In fact, \citet{robins2009semiparametric} also showed that \eqref{minimal} is the minimal condition for the existence of $\sqrt{n}$-consistent estimators of $\psi (\theta)$ under the \text{H\"{o}lder}{} nuisance modeling assumption. Outside \eqref{minimal}, $\psi (\theta)$ is non $\sqrt{n}$-estimable and the only known estimator with the optimal rate of convergence in minimax sense is again the HOIF estimator \citep{robins2016technical, robins2017minimax, robins2022corrigenda}. When restricting to highly smooth $g$, \citet{liu2021adaptive} construct minimax optimal and adaptive estimator of $\psi (\theta)$ by combining the HOIF estimators with the celebrated Lepskii's adaptation scheme \citep{lepskii1991problem}.
This article is about the $\sqrt{n}$-estimable regime \eqref{minimal}, so we will focus our attention on the eHOIF estimators. First, we choose a set of $k$-dimensional functions $\bar{{\mathsf{z}}}_{k} \equiv (z_{1}, \cdots, z_{k})^{\top}: {\mathcal{X}} \rightarrow {\mathbb{R}}^{k}$ satisfying certain regularity conditions to be given later in Section \ref{sec:soif}. The Second-Order Influence Function (SOIF) estimator of $\psi (\theta)$ is the following second-order $U$-statistic:
\begin{equation}
\label{soif}
\begin{split}
& \widehat{\psi}_{2, k} (\Omega) \coloneqq \widehat{\psi}_{1} + \widehat{\mathbb{IF}}_{2, 2, k} \text{ where } \widehat{\mathbb{IF}}_{2, 2, k} \equiv \widehat{\mathbb{IF}}_{2, 2, k} (\Omega) \coloneqq {\mathbb{U}}_{n, 2} \left[ \mathsf{IF}_{2, 2, k; 1, 2} \right]
\end{split}
\end{equation}
and
\begin{align*}
\mathsf{IF}_{2, 2, k; 1, 2} & \equiv \mathsf{IF}_{2, 2, k; 1, 2} (\Omega) \coloneqq \left( A_{1} \widehat{a} (X_{1}) - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \Omega \bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} (Y_{2} - \widehat{b} (X_{2})) \\
& \equiv \left( A_{1} \widehat{a} (X_{1}) - 1 \right) K_{k} (X_{1}, X_{2}) A_{2} (Y_{2} - \widehat{b} (X_{2})).
\end{align*}
Based on the definition of HOIFs \citep{robins2016technical}, $- \widehat{\mathbb{IF}}_{2, 2, k}$ is in fact the SOIF of $\mathsf{bias} (\widehat{\psi}_{1})$\footnote{The difference in the signs in $\widehat{\mathbb{IF}}_{2, 2, k}$ between here and \citet{robins2016technical} is non-essential.}. A more intuitively appealing explanation goes as follows: $- \widehat{\mathbb{IF}}_{2, 2, k}$ is an unbiased estimator of the following quantity:
\begin{equation}
\label{bias_k}
\begin{split}
\mathsf{bias}_{k} (\widehat{\psi}_{1}) & = {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X)}{a (X)} - 1 \right) \bar{{\mathsf{z}}}_{k} (X)^{\top} \right] \Omega {\mathbb{E}} \left[ A \bar{{\mathsf{z}}}_{k} (X) (b (X) - \widehat{b} (X)) \right] \\
& = {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X_{1})}{a (X_{1})} - 1 \right) K_{k} (X_{1}, X_{2}) A_{2} (b (X_{2}) - \widehat{b} (X_{2})) \right]
\end{split}
\end{equation}
which is simply replacing the estimation errors $\widehat{a} / a - 1$ and $b - \widehat{b}$ in \eqref{first-order bias} by
$$
\Pi \left[ \left. \frac{\widehat{a}}{a} - 1 \right\vert \bar{{\mathsf{z}}}_{k} \right] \text{ and } \Pi \left[ \left. b - \widehat{b} \right\vert A \bar{{\mathsf{z}}}_{k} \right].
$$
Hence $\widehat{\mathbb{IF}}_{2, 2, k}$ can be interpreted as a bias correction term that {\it partially} debiases $\mathsf{bias} (\widehat{\psi}_{1})$.
However, evaluating $\Omega$ in practice relies on the knowledge of $g$, which is generally unknown to the analyst. The initial attempt by \citet{robins2016technical} and \citet{robins2017minimax} was to estimate $g$ from the nuisance sample by $\widehat{g}$, leading to statistical properties affected by $g - \widehat{g}$ and thus complexity-reducing assumptions on ${\mathcal{G}} \ni g$. To completely resolve this reliance, \citet{liu2017semiparametric} choose to estimate $\Omega$ by its empirical analogue using the {\it nuisance sample}, denoted as $\widehat{\Omega}_{\mathrm{nuis}} = \widehat{\Sigma}_{\mathrm{nuis}}^{-1}$. The resulting estimated kernel is denoted as $\widehat{K}_{k}^{\mathrm{nuis}} (x, x')$, similar to $\widehat{K}_{k}$ defined in Section \ref{sec:notation}. Then the empirical SOIF (eSOIF) estimator $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}})$ of $\mathsf{bias}_{k} (\widehat{\psi}_{1})$ is
\begin{align*}
\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}}) \equiv {\mathbb{U}}_{n, 2} \left[ \mathsf{IF}_{2, 2, k; 1, 2} (\widehat{\Omega}_{\mathrm{nuis}}) \right],
\end{align*}
which, unlike $\widehat{\mathbb{IF}}_{2, 2, k}$, incurs a kernel estimation bias
\begin{align*}
{\mathbb{E}} [\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}}) - \widehat{\mathbb{IF}}_{2, 2, k}] & = {\mathbb{E}} [(A_{1} \widehat{a} (X_{1}) - 1) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top}] (\widehat{\Omega}_{\mathrm{nuis}} - {\mathbb{I}}) {\mathbb{E}} [\bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} (Y_{2} - \widehat{b} (X_{2}))] \\
& = {\mathbb{E}} [(A_{1} \widehat{a} (X_{1}) - 1) (\widehat{K}_{k}^{\mathrm{nuis}} (X_{1}, X_{2}) - K_{k} (X_{1}, X_{2})) A_{2} (Y_{2} - \widehat{b} (X_{2}))],
\end{align*}
shown to be of order at most $\sqrt{k \log k / n}$ in \citet{liu2017semiparametric}. To further reduce the kernel estimation bias, one can consider the following $m$-th order eHOIF estimator, which is an $m$-th order $U$-statistic:
\begin{align*}
& \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}_{\mathrm{nuis}}) \coloneqq \sum_{j = 2}^{m} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}_{\mathrm{nuis}}) \\
\text{where } & \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}_{\mathrm{nuis}}) \coloneqq {\mathbb{U}}_{n, j} \left[ \widehat{\mathsf{IF}}_{j, j, k; 1, \cdots, j} (\widehat{\Omega}_{\mathrm{nuis}}) \right]
\end{align*}
and
\begin{align*}
\widehat{\mathsf{IF}}_{j, j, k; 1, \cdots, j} (\widehat{\Omega}_{\mathrm{nuis}}) = (-1)^{j} (A_{1} \widehat{a} (X_{1}) - 1) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega}_{\mathrm{nuis}} \prod_{s = 3}^{j} \left\{ (Q_{s} - \widehat{\Sigma}_{\mathrm{nuis}}) \widehat{\Omega}_{\mathrm{nuis}} \right\} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} (Y_{2} - \widehat{b} (X_{2})).
\end{align*}
\citet{liu2017semiparametric} showed that the kernel estimation bias of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}_{\mathrm{nuis}})$ is of order at most $(k \log k / n)^{m / 2}$ and variance of order at most $1 / n \vee k / n^{2}$. Hence by taking $m \asymp \sqrt{\log n}$ and $k \asymp n / \log^{c} n$ for some absolute constant $c > 0$, we could estimate $\mathsf{bias}_{k} (\widehat{\psi}_{1})$ with essentially no bias without inflating the order of the variance of $\widehat{\psi}_{1}$. Furthermore, under \text{H\"{o}lder}{} nuisance models on ${\mathcal{A}} \times {\mathcal{B}}$, \citet{liu2017semiparametric} demonstrate that the sHOIF estimator $\widehat{\psi}_{1} + \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}_{\mathrm{nuis}})$, with said choices of $m$ and $k$, is $\sqrt{n}$-consistent in \eqref{minimal} and semiparametric efficient in the interior of \eqref{minimal} under some additional mild assumptions. In this paper, the sHOIF estimators to be introduced in Section \ref{sec:shoif} simply replace $\widehat{\Sigma}_{\mathrm{nuis}}$ and $\widehat{\Omega}_{\mathrm{nuis}}$ in the eHOIF estimators by $\widehat{\Sigma}$ and $\widehat{\Omega}$, the empirical analogues of $\Sigma$ and $\Omega$ computed from the estimation sample. One can easily see that, due to the correlation between $\widehat{\Omega}$ and the estimation sample, the analysis of the statistical properties of sHOIF estimators becomes significantly more challenging.
\subsection{Plan} \label{sec:outline}
The rest of the paper is organized as follows. Section \ref{sec:soif} defines the stable Second-Order IF (sSOIF) estimators and studies their statistical and numerical properties as a warm-up. Section \ref{sec:shoif} presents the full version of sHOIF estimators, together with their statistical, numerical, and computational properties. We then apply sHOIF estimators and their statistical properties to two concrete problems Section \ref{sec:applications}: one is to show that sHOIF estimators for $\psi (\theta)$ achieve semiparametric efficiency under the minimal conditions within the classical \text{H\"{o}lder}{} nuisance models; the other is to use sHOIF estimators to test if the nominal $(1 - \alpha) \times 100\%$ Wald confidence interval centered at the first-order DML estimator has the claimed coverage, a novel assumption-lean statistical procedure recently proposed in \citet{liu2020nearly}, and further developed in \citet{liu2021assumption}. To demonstrate the generality of sHOIF estimators, Section \ref{sec:extensions} extends results heretofore in several directions.
Finally, Section \ref{sec:discussion} concludes the paper and discusses several open problems and possible future directions. Appendix contains technical details that provide insights on the proof strategy. The remaining technical details are deferred to Supplementary Materials \citep{shoif_supp}.
\section{Assumptions and warm-up: Stable second-order influence function estimators}
\label{sec:soif}
In this section, we disclose the main assumptions, accompanied with an illustration of the main results using the stable second-order influence function (sSOIF) estimator $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ as a warm-up of what follows.
The assumptions below are imposed throughout the paper unless stated otherwise.
\begin{assumption}[Conditions on initial first-step nuisance parameter estimates.]
\label{cond:nuis}
Nuisance parameter estimators $\widehat{a}$ and $\widehat{b}$ are attained from a separate independent frozen nuisance sample. For simplicity, we assume this sample to also have size $n$. $\widehat{a}$ and $\widehat{b}$ further satisfy the following properties until otherwise noticed:
\begin{enumerate}[label = (\roman*)]
\item $\Vert \widehat{a} - a \Vert_{2} = o (1)$ and $\Vert \widehat{b} - b \Vert_{2} = o (1)$, i.e. both nuisance parameter estimators are $L_{2}$-consistent;
\item $\Vert a \Vert_{\infty}$, $\Vert b \Vert_{\infty}$, $\Vert \widehat{a} \Vert_{\infty}$ and $\Vert \widehat{b} \Vert_{\infty}$ are bounded by some absolute constant $B > 0$.
\item In the case of $\psi (\theta) = {\mathbb{E}} [Y (1)]$ under strong ignorability, we additionally need $1 / a$ and $1 / \widehat{a}$ to be bounded between $(c, 1 - c)$ for some absolute constant $0 < c < 0.5$.
\end{enumerate}
\end{assumption}
\begin{assumption}[Conditions on $\bar{{\mathsf{z}}}_{k}$ related quantities.]
\label{cond:b}
The following are assumed on the basis functions $\bar{{\mathsf{z}}}_{k}$ and the corresponding (inverse) Gram matrices $\Sigma, \widehat{\Sigma}, \Omega, \widehat{\Omega}$ and projection kernels $K_{k}$ and $\widehat{K}_{k}$:
\begin{enumerate}[label = (\roman*)]
\item There exists an absolute constant $B > 0$ such that $\sup_{x \in {\mathcal{X}}} K_{k} (x, x) \leq B k$ and $\sup_{x \in {\mathcal{X}}} \widehat{K}_{k} (x, x) \leq B k$;
\item Both $\Sigma$ and $\widehat{\Sigma}$ have bounded spectra;
\item The projection kernel satisfies the following $L_{\infty}$-stability condition: for any measurable function $h: {\mathcal{X}} \rightarrow {\mathbb{R}}$,
\begin{equation}
\label{l_inf_stability}
\left\Vert \Pi \left[ h | \bar{{\mathsf{z}}}_{k} \right] (\cdot) \right\Vert_{\infty} \lesssim \Vert h \Vert_{\infty}.
\end{equation}
\end{enumerate}
\end{assumption}
\begin{remark}[Comments on Assumptions \ref{cond:nuis} and \ref{cond:b}] \leavevmode
\begin{enumerate}[label = (\roman*)]
\item Given Assumption \ref{cond:b}(ii), there is no loss of generality by assuming $\Sigma \equiv \Omega \equiv {\mathbb{I}}$, the identity matrix of the same size as $\Sigma$ or $\Omega$. We make such a simplification throughout the paper, unless stated otherwise.
\item The assumptions on the nuisance parameters and their estimators in Assumption \ref{cond:nuis} are quite mild. In particular, we do not assume $\widehat{a}$, $\widehat{b}$ converge to $a$, $b$ at any algebraic rate in $L_{2}$-norm. In fact, if content with $\sqrt{n}$-consistency instead of semiparametric efficiency, $\Vert \widehat{a} - a \Vert_{2} = o (1)$ and $\Vert \widehat{b} - b \Vert_{2} = o (1)$ can be even relaxed to $\Vert \widehat{a} - a \Vert_{2} = O (1)$ and $\Vert \widehat{b} - b \Vert_{2} = O (1)$; see \citet{liu2017semiparametric}.
\item Assumption \ref{cond:b} on the dictionary $\bar{{\mathsf{z}}}_{k}$ also appeared in \citet{robins2017minimax, liu2017semiparametric, liu2020nearly, liu2021assumption}; also see comments in \citet{liu2020rejoinder}. The $L_{\infty}$-stability condition (iii) have been established for Cohen-Daubechies-Vial wavelets, B-splines, and local polynomial partition series \citep{belloni2015some}. It is possible to relax such a condition to a high-probability version, which we decide not to further pursue in this paper.
\end{enumerate}
\end{remark}
The following result on the sSOIF estimator $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ is a special case of Theorem \ref{thm:properties} to be revealed in Section \ref{sec:shoif}.
\begin{proposition}[Bias and variance bounds of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$.]
\label{prop:soif}
Under Assumptions \ref{cond:nuis} -- \ref{cond:b}, with $k = o (n)$, one has the following:
\begin{enumerate}[label = {\normalfont(\roman*)}]
\item The kernel estimation bias of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ satisfies
\begin{equation} \label{EB2}
\begin{split}
& \mathsf{kern\mbox{-}bias}_{2, k} (\widehat{\psi}_{1}) \coloneqq {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) \right] - \mathsf{bias}_{\theta, k} (\widehat{\psi}_{1}) \equiv {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \\
& \lesssim \frac{k}{n} \left\{ \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2} + \left\Vert \frac{\widehat{a} - a}{a} \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2} + \left( \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{\infty} \left\Vert \widehat{b} - b \right\Vert_{2} \right) \right\}.
\end{split}
\end{equation}
\item The variance of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ satisfies
\begin{equation} \label{var2}
\begin{split}
\mathsf{var} \left[ \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) \right] \lesssim \frac{1}{n} \left\{ \frac{k}{n} + \left( \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{\infty} \left\Vert \widehat{b} - b \right\Vert_{2} \right) \right\}.
\end{split}
\end{equation}
\end{enumerate}
\end{proposition}
\begin{remark}
The dependence on the condition number $k / n$ in the kernel estimation bias upper bound of the eSOIF estimator $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}})$ in \citet{liu2017semiparametric} ($\sqrt{k \log k / n}$) is worse than that of the sSOIF estimator $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ reported here ($k / n$).
\end{remark}
\allowdisplaybreaks
\subsection{Proof sketch of Proposition \ref{prop:soif}}
\label{sec:proof_idea_soif}
\subsubsection{Kernel estimation bias bound}
\label{sec:proof_idea_soif_eb}
$\mathsf{kern\mbox{-}bias}_{2, k} (\widehat{\psi}_{1})$ can be controlled by repeatedly using the matrix identity $(A - B)^{-1} - A^{-1} = - A^{-1} B (A - B)^{-1}$ with $A \equiv \widehat{\Sigma}^{\dag} \coloneqq n^{-1} \sum_{i = 3}^{n} Q_{i}$ and $B = n^{-1} Q_{1, 2}$:
\begin{align*}
& \ \mathsf{kern\mbox{-}bias}_{2, k} (\widehat{\psi}_{1}) \\
= & \ {\mathbb{E}} \left[ (A_{1} \widehat{a} (X_{1}) - 1) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} (\widehat{\Omega}^{-1} - {\mathbb{I}}) \bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} (Y_{2} - \widehat{b} (X_{2})) \right] \\
= & \ \sum_{j = 1}^{J - 1} {\mathbb{E}} \left[ (A_{1} \widehat{a} (X_{1}) - 1) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \left( {\mathbb{I}} - \widehat{\Sigma}^{\dag} - \frac{Q_{1, 2}}{n} \right)^{j} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} (Y_{2} - \widehat{b} (X_{2})) \right] \\
& + {\mathbb{E}} \left[ (A_{1} \widehat{a} (X_{1}) - 1) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \left( {\mathbb{I}} - \widehat{\Sigma}^{\dag} - \frac{Q_{1, 2}}{n} \right)^{J} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} (Y_{2} - \widehat{b} (X_{2})) \right].
\end{align*}
By choosing $J \asymp \log n$, the second term of the above display can be shown to be $o (n^{- 1 / 2})$.
For the first term, we only look at $j = 1, 2$ in the main text and the remaining analysis is a special case of the proof of Theorem \ref{thm:properties} in Appendix \ref{app:main}.
For $j = 1$, we have
\allowdisplaybreaks
\begin{align}
& \ {\mathbb{E}} \left[ (A_{1} \widehat{a} (X_{1}) - 1) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \left( {\mathbb{I}} - \widehat{\Sigma}^{\dag} - \frac{Q_{1, 2}}{n} \right) \bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} (Y_{2} - \widehat{b} (X_{2})) \right] \notag \\
= & \ \frac{2}{n} {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X_{1})}{a (X_{1})} - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) (b (X_{2}) - \widehat{b} (X_{2})) \right] \notag \\
& - \frac{1}{n} {\mathbb{E}} \left[ \left( A_{1} \widehat{a} (X_{1}) - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} Q_{1, 2} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} \left( Y_{2} - \widehat{b} (X_{2}) \right) \right] \label{bias_j1} \\
= & \ O \left( \frac{1}{n} \right) - \frac{1}{n} {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X_{1}) - 1}{a (X_{1})} \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) (b (X_{2}) - \widehat{b} (X_{2})) \right] \notag \\
& - \frac{1}{n} {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X_{1})}{a (X_{1})} - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) \bar{{\mathsf{z}}}_{k} (X_{2})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) (b (X_{2}) - \widehat{b} (X_{2})) \right] \notag \\
\lesssim & \ \frac{1}{n} + \frac{k}{n} \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \Vert b - \widehat{b} \Vert_{2} + \frac{k}{n} \left\Vert \frac{\widehat{a}}{a} - 1 \right\Vert_{2} \Vert b - \widehat{b} \Vert_{2} \notag
\end{align}
where the last line follows from triangle inequality, Cauchy-Schwarz inequality and Assumptions \ref{cond:nuis}, \ref{cond:b}(i) and \ref{cond:b}(ii).
For $j = 2$, we have
\begin{align}
& \ {\mathbb{E}} \left[ (A_{1} \widehat{a} (X_{1}) - 1) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \left( {\mathbb{I}} - \widehat{\Sigma}^{\dag} - \frac{Q_{1, 2}}{n} \right)^{2} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} (Y_{2} - \widehat{b} (X_{2})) \right] \notag \\
= & \ {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X_{1})}{a (X_{1})} - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} \left[ \left( {\mathbb{I}} - \widehat{\Sigma}^{\dag} \right)^{2} \right] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) (b (X_{2}) - \widehat{b} (X_{2})) \right] \notag \\
& - \frac{2}{n^{2}} {\mathbb{E}} \left[ \left( A_{1} \widehat{a} (X_{1}) - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} Q_{1, 2} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} \left( Y_{2} - \widehat{b} (X_{2}) \right) \right] \notag \\
& + \frac{1}{n^{2}} {\mathbb{E}} \left[ \left( A_{1} \widehat{a} (X_{1}) - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} Q_{1, 2}^{2} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} \left( Y_{2} - \widehat{b} (X_{2}) \right) \right] \label{bias_j2_1} \\
\eqqcolon & \ (\mathrm{I}) + (\mathrm{II}) + (\mathrm{III}). \notag
\end{align}
Since $(\mathrm{II})$ is dominated by the term for $j = 1$, we only need to further analyze $(\mathrm{I})$ and $(\mathrm{III})$. $\mathrm{(III)}$ can be bounded by
\begin{equation}
\label{first-non-commutative}
(\mathrm{III}) \lesssim \left( \frac{k}{n} \right)^{2} \left\{ \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \Vert b - \widehat{b} \Vert_{2} + \left\Vert \frac{\widehat{a}}{a} - 1 \right\Vert_{2} \Vert b - \widehat{b} \Vert_{2} + \left( \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{\infty} \Vert b - \widehat{b} \Vert_{2} \right) \wedge \left( \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \Vert b - \widehat{b} \Vert_{\infty} \right) \right\}
\end{equation}
where the first two terms are due to the first three terms in the (non-commutative) expansion of
\begin{equation}
\label{q12_expand}
Q_{1, 2}^{2} = Q_{1}^{2} + Q_{2}^{2} + Q_{1} Q_{2} + Q_{2} Q_{1},
\end{equation}
and the third term comes from the last term in the above expansion. The appearance of the estimation error in $L_{\infty}$-norm is due to the opposite order of sample points indexed by $1$ and $2$ between the ``meat'' $Q_{2} Q_{1}$ and the ``bread slices'' of the ``sandwich'' structure $\bar{{\mathsf{z}}}_{k} (X_{1})^{\top} [\cdots] \bar{{\mathsf{z}}}_{k} (X_{2})$.
For $(\mathrm{I})$, we need to expand $({\mathbb{I}} - \widehat{\Sigma}^{\dag})^{2}$.
\begin{align}
(\mathrm{I}) = & \ {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X_{1})}{a (X_{1})} - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} \left[ {\mathbb{I}} - 2 \widehat{\Sigma}^{\dag} + \widehat{\Sigma}^{\dag}{}^{2} \right] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) (b (X_{2}) - \widehat{b} (X_{2})) \right] \notag \\
= & \ \left( \frac{4}{n} - 1 \right) {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X_{1})}{a (X_{1})} - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) (b (X_{2}) - \widehat{b} (X_{2})) \right] \notag \\
& + \frac{(n - 2) (n - 3)}{n^{2}} {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X_{1})}{a (X_{1})} - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) (b (X_{2}) - \widehat{b} (X_{2})) \right] \notag \\
& + \frac{n - 2}{n^{2}} {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X_{1})}{a (X_{1})} - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} \left[ Q_{3}^{2} \right] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) (b (X_{2}) - \widehat{b} (X_{2})) \right] \notag \\
= & \ \left( \frac{6}{n^{2}} - \frac{1}{n} \right) {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X_{1})}{a (X_{1})} - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) (b (X_{2}) - \widehat{b} (X_{2})) \right] \notag \\
& + \left( \frac{1}{n} - \frac{2}{n^{2}} \right) {\mathbb{E}} \left[ \left( \frac{\widehat{a} (X_{1})}{a (X_{1})} - 1 \right) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} \left[ Q_{3}^{2} \right] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) (b (X_{2}) - \widehat{b} (X_{2})) \right]. \label{bias_j2_2}
\end{align}
It is straightforward to see the first term in the last equality of the above display is dominated by the term for $j = 1$, whereas the second term can be shown to be bounded by
\begin{align*}
\frac{k}{n} \left\Vert \frac{\widehat{a}}{a} - 1 \right\Vert_{2} \Vert b - \widehat{b} \Vert_{2}.
\end{align*}
Taken together, the terms for $j = 1$ and $j = 2$ give the desired bound for $\mathsf{kern\mbox{-}bias}_{2, k} (\widehat{\psi}_{1})$ in \eqref{EB2}. It remains to prove the terms for $j \geq 3$ are of smaller order, which is deferred to Appendix \ref{app:main_bias} for the general case. For $j \geq 3$, the corresponding term is of order
\begin{align*}
\left( \frac{k}{n} \right)^{j - 1} \left\Vert \frac{\widehat{a}}{a} - 1 \right\Vert_{2} \Vert b - \widehat{b} \Vert_{2}.
\end{align*}
\begin{remark}
\label{rem:compare}
Now is a perfect time to compare how the analysis of the kernel estimation bias of sSOIF differs from that of eSOIF of \citet{liu2017semiparametric}. The only difference between the eSOIF and sSOIF estimators are the samples used to estimate $\Omega$. Using the nuisance sample instead, the (conditional) kernel estimation bias of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}})$ conditioning on the nuisance sample data is
\begin{align*}
& \ {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \\
= & \ {\mathbb{E}} \left[ (A_{1} \widehat{a} (X_{1}) - 1) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] \left( \widehat{\Omega}_{\mathrm{nuis}} - {\mathbb{I}} \right) {\mathbb{E}} \left[ A_{2} \bar{{\mathsf{z}}}_{k} (X_{2}) (Y_{2} - \widehat{b} (X_{2})) \right].
\end{align*}
From this, we can conclude
\begin{align*}
\left\vert {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \right\vert \lesssim \left( \frac{k \log k}{n} \right)^{1 / 2} \left\Vert \frac{\widehat{a}}{a} - 1 \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2}
\end{align*}
by using matrix Bernstein or Khintchine inequality \citep{rudelson1999random, bandeira2021matrix}; also see \citet{couillet2022random}. \citet{liu2017semiparametric} further show that
\begin{align*}
\left\vert {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}_{\mathrm{nuis}}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \right\vert \lesssim \left( \frac{k \log k}{n} \right)^{m / 2} \left\Vert \frac{\widehat{a}}{a} - 1 \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2}.
\end{align*}
However, as pointed out in Section \ref{sec:novelty}, the finite-sample performance of eHOIF estimators is not well-reflected by these upper bounds, prompting the need of developing sHOIF estimators.
\end{remark}
\subsubsection{Variance bound}
\label{sec:proof_idea_soif_var}
The variance bound is technically involved. The missing steps can be found in Appendix \ref{app:soif_var}. The key step is to show
\begin{equation}
\label{var-bound-key}
{\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{3} \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{4}) \right] - \left( {\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{2}) \right] \right)^{2} = O \left( \frac{1}{n} \right).
\end{equation}
To prove \eqref{var-bound-key}, it is sufficient to exhibit
\begin{equation}
\label{var-bound-key-1}
\begin{split}
& \ {\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{3} \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} \left( \widehat{\Omega} - {\mathbb{I}} \right) \bar{{\mathsf{z}}}_{k} (X_{4}) \right] \\
& - {\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) \right] {\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \left( \widehat{\Omega} - {\mathbb{I}} \right) \bar{{\mathsf{z}}}_{k} (X_{2}) \right] = O \left( \frac{1}{n} \right)
\end{split}
\end{equation}
and
\begin{equation}
\label{var-bound-key-2}
\begin{split}
& \ {\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \left( \widehat{\Omega} - {\mathbb{I}} \right) \bar{{\mathsf{z}}}_{k} (X_{2}) A_{3} \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} \left( \widehat{\Omega} - {\mathbb{I}} \right) \bar{{\mathsf{z}}}_{k} (X_{4}) \right] \\
& - {\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \left( \widehat{\Omega} - {\mathbb{I}} \right) \bar{{\mathsf{z}}}_{k} (X_{2}) \right] {\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \left( \widehat{\Omega} - {\mathbb{I}} \right) \bar{{\mathsf{z}}}_{k} (X_{2}) \right] = O \left( \frac{1}{n} \right).
\end{split}
\end{equation}
Recall $\widehat{\Omega} = \widehat{\Sigma}^{-1}$ and $\widehat{\Sigma} = n^{-1} \sum_{i = 1}^{n} Q_{i}$. We introduce independent ``ghost copies'' $Q_{1}', Q_{2}'$ of $Q_{1}, Q_{2}$ and denote $\widehat{\Omega}' = (\widehat{\Sigma}')^{-1}$ and $\widehat{\Sigma}' = n^{-1} \sum_{i = 3}^{n} Q_{i} + n^{-1} (Q_{1}' + Q_{2}')$. Then \eqref{var-bound-key-1} is equivalent to
\begin{equation}
\label{var-bound-key-1-1}
{\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{3} \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} \left( \widehat{\Omega} - \widehat{\Omega}' \right) \bar{{\mathsf{z}}}_{k} (X_{4}) \right] = O \left( \frac{1}{n} \right).
\end{equation}
Let $\bar{\Omega} \equiv \bar{\Sigma}$ and $\bar{\Sigma} \equiv n^{-1} \sum_{i = 5}^{n} Q_{i}$. Repeating the matrix identity $(A + B)^{-1} - A^{-1} = - A^{-1} B (A + B)^{-1}$ on \eqref{var-bound-key-1-1} by setting $A = \bar{\Sigma}$ and $B = n^{-1} (Q_{1, 2} + Q_{3, 4})$ or $B = n^{-1} (Q_{1, 2}' + Q_{3, 4})$, we have
\begin{align*}
& \ {\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{3} \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} \left( \widehat{\Omega} - \widehat{\Omega}' \right) \bar{{\mathsf{z}}}_{k} (X_{4}) \right] \\
= & \ \sum_{j = 1}^{J - 1} (-1)^{j} {\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{3} \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} \left\{ \left( \bar{\Omega} \frac{Q_{1, 2} + Q_{3, 4}}{n} \right)^{j} - \left( \bar{\Omega} \frac{Q_{1, 2}' + Q_{3, 4}}{n} \right)^{j} \right\} \bar{\Omega} \bar{{\mathsf{z}}}_{k} (X_{4}) \right] \\
& + (-1)^{J} {\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{3} \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} \left\{ \left( \bar{\Omega} \frac{Q_{1, 2} + Q_{3, 4}}{n} \right)^{J} - \left( \bar{\Omega} \frac{Q_{1, 2}' + Q_{3, 4}}{n} \right)^{J} \right\} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{4}) \right].
\end{align*}
Let $J \asymp \log n$. The second term of the above display can be shown to be $o (1 / n)$. Proceeding to the first term, it is easy to see from Assumption \ref{cond:b}(iii) that the term corresponding to $j = 1$:
\begin{align*}
& \ \frac{1}{n} {\mathbb{E}} \left[ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{3} \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} \bar{\Omega} (Q_{1, 2} - Q_{1, 2}') \bar{\Omega} \bar{{\mathsf{z}}}_{k} (X_{4}) \right] \\
= & \ \frac{1}{n} \left( \begin{array}{c}
{\mathbb{E}} \left\{ A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} {\mathbb{E}} [\bar{{\mathsf{z}}}_{k} (X_{2})] {\mathbb{E}} [A_{3} \bar{{\mathsf{z}}}_{k} (X_{3})]^{\top} \bar{\Omega} \bar{{\mathsf{z}}}_{k} (X_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{\Omega} {\mathbb{E}} [\bar{{\mathsf{z}}}_{k} (X_{4})] \right\} \\
+ \ {\mathbb{E}} \left\{ {\mathbb{E}} [A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})]^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathbb{E}} [A_{3} \bar{{\mathsf{z}}}_{k} (X_{3})]^{\top} \bar{\Omega} A_{2} \bar{{\mathsf{z}}}_{k} (X_{2}) \bar{{\mathsf{z}}}_{k} (X_{2})^{\top} \bar{\Omega} {\mathbb{E}} [\bar{{\mathsf{z}}}_{k} (X_{4})] \right\} \\
- \ {\mathbb{E}} \left\{ {\mathbb{E}} [A_{1} \bar{{\mathsf{z}}}_{k} (X_{1})]^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathbb{E}} [A_{3} \bar{{\mathsf{z}}}_{k} (X_{3})]^{\top} \bar{\Omega} {\mathbb{E}} [Q_{1, 2}'] \bar{\Omega} {\mathbb{E}} [\bar{{\mathsf{z}}}_{k} (X_{4})] \right\}
\end{array} \right) = O \left( \frac{1}{n} \right).
\end{align*}
Similarly, the terms corresponding to $j \geq 2$ can be shown to be $O \left( \frac{1}{n} \left( \frac{2 k}{n} \right)^{j - 1} \right)$, a consequence of Lemma \ref{lem:nc} below. Note that the extra factor $2$ appears because there are $O (2^{j})$ terms in total by expanding out $\left( \frac{Q_{1, 2} + Q_{3, 4}}{n} \right)^{j}$ and $\left( \frac{Q_{1, 2}' + Q_{3, 4}}{n} \right)^{j}$.
\begin{lemma}
\label{lem:nc}
Given a positive integer $j$. Given any pair of integers $j_{1} \geq 0, j_{2} > 0$ such that $j_{1} + j_{2} = j$, for any subset $\mathfrak{j} \subseteq \{0, 1\}^{j}$ of the $j$-dimensional Boolean hypercube with $\Vert \mathfrak{j} \Vert_{1} = j_{1}$, we have
\begin{equation}
\label{nc-good}
{\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} \left( \prod_{\ell = 1}^{j} Q_{1, 2}^{\mathfrak{j}_{\ell}} Q_{3, 4}^{(1 - \mathfrak{j}_{\ell})} \right) \bar{{\mathsf{z}}}_{k} (X_{4}) \right] \lesssim k^{j - 1}.
\end{equation}
However, if $j_{1} = 0$ and $j_{2} = j$, we have
\begin{equation}
\label{nc-bad}
{\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} Q_{3, 4}^{j} \bar{{\mathsf{z}}}_{k} (X_{4}) \right] \lesssim k^{j}.
\end{equation}
\end{lemma}
The proof of Lemma \ref{lem:nc} can be found in Appendix \ref{app:nc}. Finally, we defer the proof of \eqref{var-bound-key-2} to the online supplements, which can be proved in a similar fashion.
\subsection{Numerical stability and time complexity of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$}
\label{sec:numerical_soif}
As discussed in Section \ref{sec:intro}, the key motivation for proposing sHOIF estimators is the numerical instability observed for eHOIF estimators. As a warm-up, we rigorously prove the numerical stability and calculate the time complexity of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ in this section. The reason why $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ can be numerically unstable is that when $k$ is near $n$, it is highly likely $\lambda_{\min} (\widehat{\Sigma}) \approx 0$ and hence $\lambda_{\max} (\widehat{\Omega})$ is close to infinity. But:
\begin{proposition} \label{prop:stability}
$\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ does not depend on the eigenvalues of $\widehat{\Omega}$.
\end{proposition}
For ease of exposition, in what follows we let
\begin{itemize}
\item $\bar{{\mathsf{Z}}}_{n, k} \coloneqq \left( \bar{{\mathsf{z}}}_{k} (X_{1}), \cdots, \bar{{\mathsf{z}}}_{k} (X_{n}) \right)^{\top}$ as the $n \times k$-matrix of the dictionary vectors for all $n$ samples;
\item $\bar{{\mathsf{Z}}}_{n, k}^{A} \coloneqq \left( A_{1} \bar{{\mathsf{z}}}_{k} (X_{1}), \cdots, A_{n} \bar{{\mathsf{z}}}_{k} (X_{n}) \right)^{\top}$ as the $n \times k$-matrix of the $A$-weighted dictionary vectors for all $n$ samples;
\item $\bm{{\mathcal{E}}}_{n, a} (\widehat{a}) \coloneqq \left( {\mathcal{E}}_{a} (\widehat{a}, O_{1}), \cdots, {\mathcal{E}}_{a} (\widehat{a}, O_{n}) \right)^{\top}$ and $\bm{{\mathcal{E}}}_{n, b} (\widehat{b}) \coloneqq \left( {\mathcal{E}}_{b} (\widehat{b}, O_{1}), \cdots, {\mathcal{E}}_{b} (\widehat{b}, O_{n}) \right)^{\top}$.
\end{itemize}
\begin{proof}
We can rewrite $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ as
\begin{equation}
\label{soif_v_stats}
\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) = \frac{1}{n - 1} \bm{{\mathcal{E}}}_{n, a} (\widehat{a})^{\top} \left[ {\mathbb{I}} - \mathsf{Diag} \right] \left\{ \bar{{\mathsf{Z}}}_{n, k} \left( \bar{{\mathsf{Z}}}_{n, k}^{\top} \bar{{\mathsf{Z}}}_{n, k}^{A} \right)^{-1} \bar{{\mathsf{Z}}}_{n, k}^{A \top} \right\} \bm{{\mathcal{E}}}_{n, b} (\widehat{b}).
\end{equation}
Now apply Singular Value Decomposition (SVD) on the matrices $\bar{{\mathsf{Z}}}_{n, k}^{A}$:
\begin{align*}
\bar{{\mathsf{Z}}}_{n, k}^{A} = U^{A} \mathsf{Diag} (D^{A}) V^{A \top}.
\end{align*}
Then
\begin{align*}
\bar{{\mathsf{Z}}}_{n, k} \left( \bar{{\mathsf{Z}}}_{n, k}^{\top} \bar{{\mathsf{Z}}}_{n, k}^{A} \right)^{-1} \bar{{\mathsf{Z}}}_{n, k}^{A \top} = \bar{{\mathsf{Z}}}_{n, k}^{A} \left( \bar{{\mathsf{Z}}}_{n, k}^{A \top} \bar{{\mathsf{Z}}}_{n, k}^{A} \right)^{-1} \bar{{\mathsf{Z}}}_{n, k}^{A \top} = U^{A} U^{A \top}.
\end{align*}
So
\begin{align*}
\eqref{soif_v_stats} = \frac{1}{n - 1} \bm{{\mathcal{E}}}_{n, a} (\widehat{a})^{\top} \left[ {\mathbb{I}} - \mathsf{Diag} \right] \left\{ U^{A} U^{A \top} \right\} \bm{{\mathcal{E}}}_{n, b} (\widehat{b}),
\end{align*}
which is completely independent of the eigenvalues of $\widehat{\Omega}$ ($(D^{A})^{2}$ up to constant).
\end{proof}
Hence it is not surprising that $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ is numerically stable even when $k \rightarrow n$.
\begin{remark}
\label{rem:self normalization}
In a sense, $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ can be viewed as a {\it self-normalized} version of $\widehat{\mathbb{IF}}_{2, 2, k}$. It is generally expected that self-normalized statistics could have better statistical properties than the non-self-normalized ones \citep{pena2008self}. However, whether the perspective of self-normalization is useful for establishing statistical properties of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ is still unclear to us and is worth pursuing as a research problem.
\end{remark}
Furthermore, not only does the alternative formula \eqref{soif_v_stats} of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ directly imply its numerical stability, but also it hints at the complexity of computing $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$. Barring the time complexity of SVD ($O (n k^{2})$), the time complexity of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ scales with $k n$ at a linear instead of a quadratic rate. This can be seen from \eqref{soif_v_stats}, in which only two vector-matrix products are involved, each taking $O (n k)$ operations. Thus we have
\begin{proposition}
The time complexity of computing $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ is $O (n k^{2})$, dominated by that of SVD.
\end{proposition}
\begin{remark}
Another alternative way of arriving at the above conclusion is to observe that the $U$-statistic kernel of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$, denoted as $\widehat{\mathsf{IF}}_{2, 2, k, \bar{i}_{2}} (\widehat{\Omega})$, is separable, in the following sense: there exists a pair (but not necessarily a unique pair) of functions $h_{1}, h_{2}$ such that
\begin{align*}
\widehat{\mathsf{IF}}_{2, 2, k, \bar{i}_{2}} (\widehat{\Omega}) \equiv h_{1} (O_{i_{1}}, \widehat{\Omega}) \cdot h_{2} (O_{i_{2}}, \widehat{\Omega}).
\end{align*}
\end{remark}
\section{The hierarchy of sHOIF estimators}
\label{sec:shoif}
As indicated in Section \ref{sec:review}, the $m$-th order sHOIF estimator $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ takes the same form as the $m$-th order eHOIF estimator $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega}_{\mathrm{nuis}})$, with the sole difference that $\widehat{\Omega}_{\mathrm{nuis}}$ is replaced by $\widehat{\Omega}$. Formally, the $m$-th order sHOIF and the corresponding $m$-th order estimator of $\psi (\theta)$ read as follows:
\begin{equation}
\label{shoif}
\begin{split}
\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega}) & \coloneqq (-1)^{m} {\mathbb{U}}_{n, m} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} \prod_{s = 3}^{m} \left\{ \left( Q_{s} - \widehat{\Sigma} \right) \widehat{\Omega} \right\} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right], \\
\widehat{\psi}_{m, k} (\widehat{\Omega}) & \coloneqq \widehat{\psi}_{1} + \sum_{j = 2}^{m} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega}) \equiv \widehat{\psi}_{1} + \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}).
\end{split}
\end{equation}
\begin{remark} \label{rem:form}
The above sHOIF statistics are the same as the eHOIF statistics except that $\widehat{\Omega}$ is constructed from the estimation sample instead of the nuisance sample.
\end{remark}
In this section, we first explain heuristically why $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ is enough to correct for the kernel estimation bias (see Section \ref{sec:explain}), after which the statistical, numerical and computational properties of $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ are stated formally.
\subsection{Heuristic explanation}
\label{sec:explain}
In what follows we explain heuristically why the kernel estimation bias $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ can be further corrected by adding:
\begin{align*}
\widehat{\mathbb{IF}}_{3, 3, k} (\widehat{\Omega}) \coloneqq \frac{n - 2}{n} {\mathbb{U}}_{n, 3} \left[ - \ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} \left( Q_{3} - \widehat{\Sigma} \right) \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right]
\end{align*}
and
\begin{align*}
\widehat{\mathbb{IF}}_{4, 4, k} (\widehat{\Omega}) \coloneqq \frac{n - 3}{n} {\mathbb{U}}_{n, 4} \left[ (A_{1} \widehat{a} (X_{1}) - 1) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} \left( Q_{3} - \widehat{\Sigma} \right) \widehat{\Omega} \left( Q_{4} - \widehat{\Sigma} \right) \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{2}) A_{2} (Y_{2} - \widehat{b} (X_{2})) \right].
\end{align*}
For short, we define $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (3, 3), k} (\widehat{\Omega}) \coloneqq \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) + \widehat{\mathbb{IF}}_{3, 3, k} (\widehat{\Omega})$ and $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (4, 4), k} (\widehat{\Omega}) \coloneqq \sum_{j = 2}^{4} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$. Also note that the majority of this section is written for mathematical rigor.
Simple algebra gives
\begin{align*}
& \ \widehat{\mathbb{IF}}_{3, 3, k} (\widehat{\Omega}) \\
\equiv & - \frac{n - 2}{n} \frac{(n - 3)!}{n!} \sum_{1 \leq i_{1} \neq i_{2} \neq i_{3} \leq n} {\mathcal{E}}_{a} (\widehat{a}; O_{i_{1}}) \bar{{\mathsf{z}}}_{k} (X_{i_{1}})^{\top} \widehat{\Omega} \left( Q_{i_{3}} - \widehat{\Sigma} \right) \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{i_{2}}) {\mathcal{E}}_{b} (\widehat{b}; O_{i_{2}}) \\
= & \ \frac{1}{n} {\mathbb{U}}_{n, 2} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} Q_{1, 2} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right] - \frac{2}{n} \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) \\
\approx & \ \frac{1}{n} {\mathbb{U}}_{n, 2} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} Q_{1, 2} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right] \eqqcolon \widetilde{\widehat{\mathbb{IF}}}_{3, 3, k} (\widehat{\Omega})
\end{align*}
and
\allowdisplaybreaks
\begin{align*}
& \ \widehat{\mathbb{IF}}_{4, 4, k} (\widehat{\Omega}) \\
\equiv & \ \frac{n - 3}{n} \frac{(n - 4)!}{n!} \sum_{1 \leq i_{1} \neq i_{2} \neq i_{3} \neq i_{4} \leq n} {\mathcal{E}}_{a} (\widehat{a}; O_{i_{1}}) \bar{{\mathsf{z}}}_{k} (X_{i_{1}})^{\top} \widehat{\Omega} \prod_{s = 3}^{4} \left[ \left( Q_{i_{s}} - \widehat{\Sigma} \right) \widehat{\Omega} \right] \bar{{\mathsf{z}}}_{k} (X_{i_{2}}) {\mathcal{E}}_{b} (\widehat{b}; O_{i_{2}}) \\
= & \ \frac{1}{n (n - 2)} {\mathbb{U}}_{n, 2} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} Q_{1, 2} \widehat{\Omega} Q_{1, 2} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right] \\
& - \frac{1}{n} {\mathbb{U}}_{n, 3} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} Q_{3} \widehat{\Omega} Q_{3} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right] \\
& - \frac{1}{n - 2} \left( 6 \widehat{\mathbb{IF}}_{3, 3, k} (\widehat{\Omega}) - \left( 1 - \frac{6}{n} \right) \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) \right) \\
\approx & \ \frac{1}{n^{2}} {\mathbb{U}}_{n, 2} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} Q_{1, 2} \widehat{\Omega} Q_{1, 2} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right] \\
& - \frac{1}{n} {\mathbb{U}}_{n, 3} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} Q_{3} \widehat{\Omega} Q_{3} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right] \\
\eqqcolon & \ \widetilde{\widehat{\mathbb{IF}}}_{4, 4, k} (\widehat{\Omega}).
\end{align*}
First, observe that the expectation of the oracle version of $\widetilde{\widehat{\mathbb{IF}}}_{3, 3, k} (\widehat{\Omega})$
\begin{align*}
{\mathbb{E}} \left[ \widetilde{\widehat{\mathbb{IF}}}_{3, 3, k} \right] = \frac{1}{n} {\mathbb{E}} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} Q_{1, 2} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right]
\end{align*}
exactly cancels \eqref{bias_j1}, the leading-order part of the kernel estimation bias of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ corresponding to $j = 1$.
Next, observe that the expectation of the oracle version of $\widetilde{\widehat{\mathbb{IF}}}_{4, 4, k} (\widehat{\Omega})$ is
\begin{align*}
{\mathbb{E}} \left[ \widetilde{\widehat{\mathbb{IF}}}_{4, 4, k} \right] = & \ \frac{1}{n^{2}} {\mathbb{E}} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} Q_{1, 2}^{2} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right] \\
& - \frac{1}{n} {\mathbb{E}} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} [Q_{3}^{2}] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right]
\end{align*}
which again cancels the kernel estimation bias of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (3, 3), k} (\widehat{\Omega})$ truncated at $j = 2$, dominated by
\begin{equation}
\label{toif_eb}
- \frac{1}{n^{2}} {\mathbb{E}} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} Q_{1, 2}^{2} \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right] + \frac{1}{n} {\mathbb{E}} \left[ {\mathcal{E}}_{a} (\widehat{a}; O_{1}) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] {\mathbb{E}} [Q_{3}^{2}] {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{2}) {\mathcal{E}}_{b} (\widehat{b}; O_{2}) \right]
\end{equation}
which can be derived from \eqref{bias_j2_1}, \eqref{bias_j2_2} and the kernel estimation bias of $\widehat{\mathbb{IF}}_{3, 3, k} (\widehat{\Omega})$ truncated at level $j = 1$; see Appendix \ref{app:toif_eb} for a more detailed calculation. Hence $\widehat{\mathbb{IF}}_{(3, 3) \rightarrow (4, 4), k} (\widehat{\Omega})$ further reduces the kernel estimation bias of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$.
\subsection{Characterization of the bias and variance of the sHOIF estimators}
\label{sec:properties}
We now state the main theoretical result of this paper.
\begin{theorem}
\label{thm:properties}
Under Assumptions \ref{cond:nuis} -- \ref{cond:b}, with $k \lesssim \frac{n}{\log^{2} n}$ and $m \gtrsim \sqrt{\log n}$, one has the following:
\begin{enumerate}[label = {\normalfont(\roman*)}]
\item The kernel estimation bias of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ satisfies
\begin{equation} \label{EBm}
\begin{split}
& \mathsf{kern\mbox{-}bias}_{m, k} (\widehat{\psi}_{1}) \coloneqq {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) \right] - \mathsf{bias}_{\theta, k} (\widehat{\psi}_{1}) \equiv {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \\
& \lesssim \left( \frac{k m}{n} \right)^{\lceil \frac{\lceil \frac{m - 1}{2} \rceil - 1}{2} \rceil \vee 1} \left\{ \begin{array}{c}
\left\Vert \dfrac{\widehat{a} - 1}{a} \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2} + \left\Vert \dfrac{\widehat{a} - a}{a} \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2} \\
+ \left( \left\Vert \dfrac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \dfrac{\widehat{a} - 1}{a} \right\Vert_{\infty} \left\Vert \widehat{b} - b \right\Vert_{2} \right)
\end{array} \right\}.
\end{split}
\end{equation}
\item For $m \geq 2$, the variance of $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ satisfies
\begin{equation} \label{var_if_m}
\begin{split}
\mathsf{var} \left[ \widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega}) \right] \lesssim \frac{1}{n} \left\{ \frac{k}{n} + \left( \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \right) \right\}.
\end{split}
\end{equation}
And thus the variance of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ satisfies
\begin{equation} \label{var_if_2_m}
\begin{split}
\mathsf{var} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) \right] \lesssim \frac{1}{n} \left\{ \frac{k}{n} + \left( \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \right) \right\}.
\end{split}
\end{equation}
\end{enumerate}
\end{theorem}
The proof of the above theorem can be found in Appendix \ref{app:main} (for kernel estimation bias bound) and the online supplements (for variance bound).
\begin{remark}[Asymptotic normality and the bootstrap approximation]
As shown in \citet{liu2020nearly}, the asymptotic normality of the oracle statistic $\frac{\widehat{\mathbb{IF}}_{2, 2, k} - \mathsf{bias}_{\theta, k} (\widehat{\psi}_{1})}{\mathsf{se}_{\theta} (\widehat{\mathbb{IF}}_{2, 2, k})}$ follows from Theorem 1 of \citet{bhattacharya1992class} whence $1 \ll k \ll n^{2}$. Thus to show CLT of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ for any $m \geq 2$, it is sufficient to demonstrate under what conditions $\mathsf{kern\mbox{-}bias}_{m, k} (\widehat{\Omega}) \ll \mathsf{se}_{\theta} (\widehat{\mathbb{IF}}_{2, 2, k})$. Bootstrap approximation (and its rate) of the distribution of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$, or even of $\widehat{\mathbb{IF}}_{2, 2, k}$ is still an important open problem, though \citet{liu2021assumption} have made some partial progress. A more thorough study of the conditions under which central limit theorem (CLT) or bootstrap approximation holds is beyond the scope of this paper.
\end{remark}
\subsection{Numerical stability and time complexity of sHOIF estimators}
\label{sec:numerical}
In what follows we consider the numerical and computational properties of sHOIF estimators, which extends the results in Section \ref{sec:numerical_soif} to higher-order. The first result in this section, Theorem \ref{thm:stability}, earmarks the ``stability'' of sHOIF estimators in terms of their independence of the eigenvalues of $\widehat{\Omega}$, the root cause of the instability of eHOIF estimators.
\begin{theorem} \label{thm:stability}
$\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ does not depend on the eigenvalues of $\widehat{\Omega}$.
\end{theorem}
\begin{proof}
The proof resembles the proof of Proposition \ref{prop:stability} closely by realizing that, for any $i, j \in [n]$,
\begin{align*}
A_{i} \bar{{\mathsf{z}}}_{k} (X_{i})^{\top} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{j}) = U^{A}_{i, \bullet} U^{A \top} U U_{j, \bullet},
\end{align*}
which is completely independent of the eigenvalues of $\widehat{\Omega}$.
\end{proof}
Hence sHOIF estimators do not suffer from any numerical instability resulted from the large condition number of the sample Gram matrix when we let $k$ near $n$ in practice.
\begin{remark}
\label{rem:algorithm}
Theorem \ref{thm:stability} also suggests a better way to compute sHOIF estimators. Instead of computing the sample Gram matrix $\widehat{\Sigma}$ and its inverse $\widehat{\Omega}$ using numerical methods, we should instead perform SVD on the basis matrices $\bar{{\mathsf{Z}}}_{n, k}$ and $\bar{{\mathsf{Z}}}_{n, k}^{A}$ and then compute $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$. In fact, the upcoming R package \citep{wanis2023machine} for computing HOIF related statistics exactly uses this strategy.
\end{remark}
Since sHOIF estimators are numerically stable and thus are potentially useful tools for statistical practice \citep{liu2020nearly, wanis2023machine}, it is worth discussing the computational complexity of sHOIF estimators for general order $m$ as well.
\begin{theorem} \label{thm:time}
The time complexity of computing $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ is $O (\max \{(n k)^{\lceil (m - 1) / 2 \rceil}, n k^{2}\})$.
\end{theorem}
\begin{proof}
Similar to the proof of Proposition \ref{prop:stability}, we need to rewrite $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ in the form of a linear combination of $V$-statistics. Without loss of generality, we take ${\mathcal{E}}_{a} (\widehat{a}; O) \equiv {\mathcal{E}}_{b} (\widehat{b}; O) \equiv 1$. But let us first represent $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ as the following series:
\begin{align*}
\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) \equiv \sum_{j = 1}^{m} (-1)^{j} \binom{m - 1}{j - 1} {\mathbb{U}}_{n, j} \left[ \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} \cdot \prod_{s = 3}^{m} \left( Q_{s} \widehat{\Omega} \right) \cdot \bar{{\mathsf{z}}}_{k} (X_{2}) \right].
\end{align*}
Note that the number of summations in $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ is
\begin{align*}
n (n - 1) \cdots (n - m + 1) = \sum_{j = 1}^{m} (-1)^{m - j} s (m, j) n^{j}
\end{align*}
where $s (m, j)$ are unsigned Stirling numbers of the first kind, or the number of permutations on $m$ elements with $j$ cycles. Accordingly one can write an $m$-th order $U$-statistic into a linear combination of $V$-statistics from order $1$ to order $m$, with the number of $j$-th order $V$-statistics, for $j = 1, \cdots, m$, equal to $s (m, j)$.
The proof is completed by leveraging the special structure of the $U$-statistic kernel for sHOIF estimators.
\end{proof}
\begin{remark}
\label{rem:comp-stat}
Considering Theorem \ref{thm:properties} and Theorem \ref{thm:time} in tandem, there is a clear statistical-computational trade-off.
\noindent However, whether or not such statistical-computational trade-off is an emanation of possibly intrinsic computational hardness of estimating certain smooth statistical functionals is still an open problem
Finally, we briefly comment on our philosophical stance on the usefulness of sHOIF estimators. sHOIF estimators are effectively infinite-order $U$-statistics, so given the current computing devices, there is no doubt that practitioners are not using sHOIF estimators in practice in near term. This is ``conditional'' on the availability of hardware. The numerical stability or lack thereof, however, is an issue regardless of the availability of more powerful computing resources.
\end{remark}
\begin{remark}
\label{rem:ehoif-comp}
Theorem \ref{thm:time} also applies to eHOIF estimators \citep{liu2017semiparametric} and the original HOIF estimators of \citet{robins2016technical}, that needs an estimate of the density of the covariates $X$, if the time for density estimation is not counted.
\end{remark}
\section{Applications of the statistical properties of sHOIF estimators}
\label{sec:applications}
\subsection{Semiparametric efficiency under minimal \text{H\"{o}lder}{} assumptions on the nuisance functions}
\label{sec:semi}
In nonparametric statistics, the optimality of a statistical procedure is often evaluated under the \text{H\"{o}lder}{} nuisance models.
The above calculations culminate into the following theorem, which is the second main result of this paper.
\begin{theorem}
\label{thm:semi}
If ${\mathcal{A}} \times {\mathcal{B}} \subseteq {\mathcal{H}} (s_{a}, {\mathcal{X}}) \times {\mathcal{H}} (s_{b}, {\mathcal{X}})$ with $(s_{a} + s_{b}) / 2 \geq d / 4$, and choosing $m \asymp \sqrt{\log n}$ and $k \lesssim n / \log (n)^{2}$,
\begin{equation}
\sqrt{n} \left( \widehat{\psi}_{m, k} - \psi (\theta) \right) \overset{{\mathcal{L}}}{\rightarrow} {\mathcal{N}} (0, {\mathbb{E}} [\mathsf{IF}_{1} (\theta)^{2}])
\end{equation}
where ${\mathbb{E}} [\mathsf{IF}_{1} (\theta)^{2}]$ is the semiparametric efficiency bound of $\psi (\theta)$.
\end{theorem}
\begin{remark}
According to the lower bound of \citet{robins2009semiparametric} under the \text{H\"{o}lder}{} nuisance model, $(s_{a} + s_{b}) / 2 \geq d / 4$ is the minimal condition for the existence of a semiparametric efficient estimator of $\psi (\theta)$. It is not unreasonable to expect that this minimal condition also holds for most, if not all, of the DRFs.
\end{remark}
\subsection{Implications on the assumption-free bias testing procedure of \citet{liu2020nearly} and \citet{liu2021assumption}}
\label{sec:bias_test}
In light of the growing interest in understanding the performance of deep-learning-based causal inference \citep{farrell2021deep, chen2020causal} and the gap between these theoretical results and empirical performance \citep{xu2022deepmed}, \citet{liu2020nearly} proposed the following oracle assumption-free valid nominal $\alpha$-level test statistic:
\begin{equation}
\label{oracle_test}
\widehat{\chi} \coloneqq \mathbbm{1} \left\{ \frac{\widehat{\mathbb{IF}}_{2, 2, k}}{\widehat{\mathsf{se}} [\widehat{\psi}_{1}]} - z_{\alpha / 2} \frac{\widehat{\mathsf{se}} [\widehat{\mathbb{IF}}_{2, 2, k}]}{\widehat{\mathsf{se}} [\widehat{\psi}_{1}]} > \delta \right\}
\end{equation}
for the following null hypothesis:
\begin{equation}
\label{h0}
{\mathsf{H}}_{0} (\delta): \frac{\mathsf{cs\mbox{-}bias} (\widehat{\psi}_{1})}{\mathsf{se} [\widehat{\psi}_{1}]} \leq \delta
\end{equation}
where
\begin{equation}
\mathsf{cs\mbox{-}bias} (\widehat{\psi}_{1}) \coloneqq \left\{ {\mathbb{E}} \left[ \lambda (X) (\widehat{a} (X) - a (X))^{2} \right] {\mathbb{E}} \left[ \lambda (X) (\widehat{b} (X) - b (X))^{2} \right] \right\}^{1 / 2}.
\end{equation}
\citet{liu2021assumption} in turn constructed a feasible assumption-lean valid nominal $\alpha$-level test statistic
\begin{equation}
\widehat{\chi}_{3, k} (\widehat{\Omega}_{\mathrm{nuis}}) \coloneqq \mathbbm{1} \left\{ \frac{\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (3, 3), k} (\widehat{\Omega}_{\mathrm{nuis}})}{\widehat{\mathsf{se}} [\widehat{\psi}_{1}]} - z_{\alpha / 2} \frac{\widehat{\mathsf{se}} [\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (3, 3), k} (\widehat{\Omega}_{\mathrm{nuis}})]}{\widehat{\mathsf{se}} [\widehat{\psi}_{1}]} > \delta \right\}
\end{equation}
and the following higher-order test statistic based on eHOIF estimators:
\begin{equation}
\widehat{\chi}_{m, k} (\widehat{\Omega}_{\mathrm{nuis}}) \coloneqq \mathbbm{1} \left\{ \frac{\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}_{\mathrm{nuis}})}{\widehat{\mathsf{se}} [\widehat{\psi}_{1}]} - z_{\alpha / 2} \frac{\widehat{\mathsf{se}} [\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}_{\mathrm{nuis}})]}{\widehat{\mathsf{se}} [\widehat{\psi}_{1}]} > \delta \right\}.
\end{equation}
\citet{liu2021assumption} showed that all the standard errors in the above test statistics can be estimated consistently by certain bootstrapping procedure. More importantly, they proved the following.
\begin{proposition}
Let
\begin{align*}
\mathsf{cs\mbox{-}bias}_{k} (\widehat{\psi}_{1}) \coloneqq \left\{ {\mathbb{E}} \left[ \Pi [(A \widehat{a} - 1) | \bar{{\mathsf{z}}}_{k}] (X)^{2} \right] {\mathbb{E}} \left[ \Pi [A (\widehat{b} - b) | \bar{{\mathsf{z}}}_{k}] (X)^{2} \right] \right\}^{1 / 2}.
\end{align*}
Under the assumptions of Theorem \ref{thm:semi}, with $k \lesssim n / (\log n)^{2}$ and the following extra condition:
\begin{equation}
\label{extra}
\vert \mathsf{kern\mbox{-}bias}_{3, k} (\widehat{\psi}_{1}) \vert \not\ll \mathsf{cs\mbox{-}bias}_{k} (\widehat{\psi}_{1})
\end{equation}
then $\widehat{\chi}_{3, k} (\widehat{\Omega}_{tr})$ is a valid nominal $\alpha$-level test of ${\mathsf{H}}_{0} (\delta)$ \eqref{h0}. The extra condition \eqref{extra} can be relaxed to
\begin{equation}
\label{extra_relax}
\vert \mathsf{kern\mbox{-}bias}_{m, k} (\widehat{\psi}_{1}) \vert \not\ll \mathsf{cs\mbox{-}bias}_{k} (\widehat{\psi}_{1}) \left( \frac{k \log k}{n} \right)^{\frac{m - 1}{2}}
\end{equation}
if one uses $\widehat{\chi}_{m, k} (\widehat{\Omega}_{tr})$ instead of $\widehat{\chi}_{3, k} (\widehat{\Omega}_{tr})$.
\end{proposition}
We can similarly define the following sHOIF-based test statistics: for $m \geq 2$,
\begin{align*}
\widehat{\chi}_{m, k} (\widehat{\Omega}) \coloneqq \mathbbm{1} \left\{ \frac{\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})}{\widehat{\mathsf{se}} [\widehat{\psi}_{1}]} - z_{\alpha / 2} \frac{\widehat{\mathsf{se}} [\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})]}{\widehat{\mathsf{se}} [\widehat{\psi}_{1}]} > \delta \right\}.
\end{align*}
Then as an immediate corollary of Theorem \ref{thm:properties}, we have
\begin{theorem}
Under the assumptions of Theorem \ref{thm:semi}, with $k \lesssim n / (\log n)^{2}$ and a different relaxed extra condition from \eqref{extra_relax}:
\begin{equation}
\label{extra_relax_stable}
\vert \mathsf{kern\mbox{-}bias}_{m, k} (\widehat{\psi}_{1}) \vert \not\ll \mathsf{cs\mbox{-}bias}_{k} (\widehat{\psi}_{1}) \left( \frac{k}{n} \right)^{m - 1}
\end{equation}
then $\widehat{\chi}_{m, k} (\widehat{\Omega})$ is a valid nominal $\alpha$-level test of ${\mathsf{H}}_{0} (\delta)$ \eqref{h0}.
\end{theorem}
Given the above theoretical guarantees, and further considering that the sHOIF estimators and tests have better finite-sample performance than the corresponding eHOIF estimators and tests, we recommend using $\widehat{\chi}_{m, k} (\widehat{\Omega})$ in practice. For more examples of its application, see \citet{wanis2023machine}.
\section{Further extensions of sHOIF estimators}
\label{sec:extensions}
\subsection{Generalization to the entire class of DRFs}
\label{sec:drf}
In this subsection, we briefly comment on how our results can be generalized to the entire class of DRFs characterized in \citet{rotnitzky2021characterization}. The class of DRFs includes many other functionals that arise in substantive studies in (bio)statistics, epidemiology, economics, and social sciences, including:
\begin{itemize}
\item the expected conditional variance, which is useful for constructing confidence/predictive sets \citep{robins2006adaptive};
\item the expected conditional covariance, which is useful for both causal inference and conditional independence testing \citep{shah2020hardness};
\item average causal effect of continuous treatment, which is important for treatment allocations \citep{ai2021unified, bonvini2022fast}.
\end{itemize}
\citet{rotnitzky2021characterization} defined the class of DRFs as follows:
\begin{definition}[Doubly Robust Functionals; Definition 1 of \citet{rotnitzky2021characterization}]
$\psi (\theta)$ is a doubly robust functional if, for each $\theta \in \Theta$ there exists $a: x \mapsto a (x) \in {\mathcal{A}}$ and $b: x \mapsto b (x) \in {\mathcal{B}}$ such that (i) $\theta = (a, b, \theta \setminus \{b, p\})$ and $\Theta = {\mathcal{A}} \times {\mathcal{B}} \times \Theta \setminus \{{\mathcal{A}}, {\mathcal{B}}\}$ and (ii) for any $\theta, \theta' \in \Theta$
\begin{equation}
\psi (\theta) - \psi (\theta^{\prime}) + {\mathbb{E}} \left[ \mathsf{IF}_{1} (\theta^{\prime}) \right] = {\mathbb{E}} \left[ S (a (X) - a^{\prime} (X)) (b (X) - b^{\prime} (X)) \right] \label{eq:drbias}
\end{equation}
where $S \equiv s (O)$ with $o \mapsto s (o)$ a known function that does not depend on $\theta$ or $\theta'$ satisfying either ${\mathbb{P}}_{\theta} (S \geq 0) = 1$ or ${\mathbb{P}}_{\theta} (S \leq 0) = 1$. We also denote $\lambda (x) \coloneqq {\mathbb{E}} [S | X = x]$. Then the first-order influence function of $\psi (\theta)$ has the following form: given $\theta' \equiv (a', b', \theta' \setminus \{a', b'\})^{\top} \in \Theta$,
\begin{equation}
\label{drf_if1}
\mathsf{IF}_{1} (\theta') \equiv S a' (X) b' (X) + m_{a} (O, a') + m_{b} (O, b') + S_{0}
\end{equation}
where $S_{0}$ is some known statistic that does not depend on $a'$ and $b'$, and $h \mapsto m_{a} (o, h)$ for $h \in {\mathcal{A}}$ and $h \mapsto m_{b} (o, h)$ for $h \in {\mathcal{B}}$ are two known linear maps satisfying
\begin{align*}
& {\mathbb{E}} \left[ S h (X) b (X) + m_{a} (O, h) \right] \equiv 0, \ \forall \ h \in {\mathcal{A}}, \\
& {\mathbb{E}} \left[ S a (X) h (X) + m_{b} (O, h) \right] \equiv 0, \ \forall \ h \in {\mathcal{B}}.
\end{align*}
As a result, $\psi (\theta) \equiv {\mathbb{E}} [m_{a} (O, a) + S_{0}] \equiv {\mathbb{E}} [m_{b} (O, b) + S_{0}]$.
\end{definition}
\begin{remark}
For $\psi (\theta) \equiv {\mathbb{E}} [Y (1)]$ under strong ignorability, $S$, $a$, $b$, $m_{a} (O, a)$, $m_{b} (O, b)$, and $S_{0}$ correspond to $- A$, $\{{\mathbb{E}} [A | X]\}^{-1}$, ${\mathbb{E}} [Y | X, A = 1]$, $A Y a (X)$, $b (X)$, and $0$, respectively. Thus $\lambda (X) = {\mathbb{E}} [- A | X] = - \frac{1}{a (X)}$. We also have
\begin{align*}
& {\mathbb{E}} [{\mathcal{E}}_{a} (\widehat{a}; O) \bar{{\mathsf{z}}}_{k} (X)] = {\mathbb{E}} \left[ \frac{1}{a (X)} (\widehat{a} (X) - a (X)) \bar{{\mathsf{z}}}_{k} (X) \right], \\
& {\mathbb{E}} [{\mathcal{E}}_{b} (\widehat{b}; O) \bar{{\mathsf{z}}}_{k} (X)] = {\mathbb{E}} \left[ \frac{1}{a (X)} (\widehat{b} (X) - b (X)) \bar{{\mathsf{z}}}_{k} (X) \right].
\end{align*}
\end{remark}
We have the following notation correspondence that maps the results for $\psi (\theta) \equiv {\mathbb{E}} [Y (1)]$ under strong ignorability to any DRF $\psi (\theta)$:
\begin{itemize}
\item ${\mathcal{E}}_{a} (\widehat{a}; O) \bar{{\mathsf{z}}}_{k} (X) \Rightarrow {\mathcal{E}}_{a} (\widehat{a}, \bar{{\mathsf{z}}}_{k}; O)$ and ${\mathcal{E}}_{b} (\widehat{b}; O) \bar{{\mathsf{z}}}_{k} (X) \Rightarrow {\mathcal{E}}_{b} (\widehat{b}, \bar{{\mathsf{z}}}_{k}; O)$ where ${\mathcal{E}}_{a} (\widehat{a}, \bar{{\mathsf{z}}}_{k}; O)$ and ${\mathcal{E}}_{b} (\widehat{b}, \bar{{\mathsf{z}}}_{k}; O)$ satisfy
\begin{align*}
& {\mathbb{E}} [{\mathcal{E}}_{a} (\widehat{a}, \bar{{\mathsf{z}}}_{k}; O)] = {\mathbb{E}} \left[ \lambda (X) (\widehat{a} (X) - a (X)) \bar{{\mathsf{z}}}_{k} (X) \right], \\
& {\mathbb{E}} [{\mathcal{E}}_{b} (\widehat{b}, \bar{{\mathsf{z}}}_{k}; O)] = {\mathbb{E}} \left[ \lambda (X) (\widehat{b} (X) - b (X)) \bar{{\mathsf{z}}}_{k} (X) \right].
\end{align*}
\item $\Sigma = {\mathbb{E}} [A \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]$ $\Rightarrow$ $\Sigma = {\mathbb{E}} [S \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]$ and $\widehat{\Sigma} = {\mathbb{P}}_{n} [A \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]$ $\Rightarrow$ $\widehat{\Sigma} = {\mathbb{P}}_{n} [S \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]$.
\end{itemize}
With the above mappings, all the theoretical results for $\psi (\theta) \equiv {\mathbb{E}} [Y (1)]$ developed herein can be applied to those for an arbitrary DRF $\psi (\theta)$ {\it mutatis mutandis}.
\subsubsection{A special case: the expected conditional covariance}
\label{sec:special}
Before concluding our paper, we further study the implications of the sHOIF theory developed so far for a special cases of DRFs: the expected conditional covariance between two random variables $A$ and $Y$ given a third random variable $X$, $\psi \coloneqq {\mathbb{E}} [\mathsf{cov} (A, Y | X)]$. When $A = Y$ almost surely, $\psi$ reduces to the expected conditional variance of $A$ given $X$, $\psi \coloneqq {\mathbb{E}} [\mathsf{cov} (A | X)]$. For differences between these two parameters, see an extended discussion in \citet{liu2020nearly}.
The main feature that distinguishes $\psi$ from many other DRFs is $S = 1$, which leads to its SOIF:
\begin{align*}
\widehat{\mathbb{IF}}_{2, 2, k} = \frac{1}{n (n - 1)} \sum_{1 \leq i_{1} \neq i_{2} \leq n} (A_{i_{1}} - \widehat{a} (X_{i_{1}})) \bar{{\mathsf{z}}}_{k} (X_{i_{1}})^{\top} \bar{\Omega} \bar{{\mathsf{z}}}_{k} (X_{i_{2}}) (Y_{i_{2}} - \widehat{b} (X_{i_{2}})),
\end{align*}
in which $\bar{\Omega} \equiv \{{\mathbb{E}} [S \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]\}^{-1} \equiv \{{\mathbb{E}} [\bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]\}^{-1}$ only depends on the distribution of $X$. Also, $\widehat{a}$ and $\widehat{b}$ are nuisance estimates of $a (x) = {\mathbb{E}} [A | X = x]$ and $b (x) = {\mathbb{E}} [Y | X = x]$ in this context. This leads to the following improved kernel estimation bias bound:
\begin{corollary}
Under Assumptions \ref{cond:nuis} -- \ref{cond:b}, with $k \lesssim \frac{n}{\log^{2} n}$ and $m \gtrsim \sqrt{\log n}$, one has the following:
\noindent The kernel estimation bias of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ satisfies
\begin{equation} \label{special EBm}
\begin{split}
& \mathsf{kern\mbox{-}bias}_{m, k} (\widehat{\psi}_{1}) \coloneqq {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) \right] - \mathsf{bias}_{\theta, k} (\widehat{\psi}_{1}) \equiv {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \\
& \lesssim \left( \frac{k m}{n} \right)^{\lceil \frac{\lceil \frac{m - 1}{2} \rceil - 1}{2} \rceil \vee 1} \left\{ \left\Vert \widehat{a} - a \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2} + \left( \left\Vert \widehat{a} - a \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \widehat{a} - a \right\Vert_{\infty} \left\Vert \widehat{b} - b \right\Vert_{2} \right) \right\}.
\end{split}
\end{equation}
\end{corollary}
Note that the variance bound is improved in a similar manner and is omitted here.
\section{Discussion}
\label{sec:discussion}
In this paper, we propose a novel class of HOIF estimators, stable HOIF (sHOIF) estimators, for the doubly robust functionals (DRFs) characterized in \citet{rotnitzky2021characterization}. They are semiparametric efficient under the minimal \text{H\"{o}lder}{}-smoothness condition $\frac{s_{a} + s_{b}}{2} > \frac{d}{4}$ of \citet{robins2009semiparametric}, allowing the dimension $k$ of the basis function diverging at a rate just slower than the sample size $n$. As can be seen from Theorem \ref{thm:properties}, sHOIF estimators have improved rate of convergence than eHOIF developed in \citet{liu2017semiparametric}. More importantly, as well documented in the simulation studies of \citet{liu2020nearly} and \citet{wanis2023machine}, the sHOIF estimators also have significantly better finite-sample performance over existing higher-order estimators in practice, making them more amenable for tasks such as testing if the bias of a first-order DML estimator $\widehat{\psi}_{1}$ of a causal effect is dominated by its standard error \citep{liu2020nearly, liu2021assumption, wanis2023machine}. Finally, we end our paper by mentioning several future research directions:
\begin{enumerate}[label = (\arabic*)]
\item It will be interesting to study if one can extend the idea of sHOIF estimators to the non-$\sqrt{n}$-estimable regimes by, for instance, estimating $\Omega$ via some shrinkage or regularized algorithms. As conjectured in \citet{robins2016technical}, the minimax convergence rate of the functionals studied in this paper may depend on the regularity of the density of the covariates $X$. Hence it is expected that the shrinkage or regularization also depends on the density of $X$. Simulation studies in \citet{liu2020nearly} suggest the nonlinear shrinkage covariance matrix estimators of \citet{ledoit2012nonlinear} could be a viable option. Preliminary simulation studies in \citet{liu2020nearly} and \citet{wanis2023machine} suggest that the performance of these shrinkage covariance matrix estimators does degrade with the smoothness of the design density.
\item As pointed out in \citet{kennedy2022minimax}, their Second-Order R-Learner (SORL) for CATE also involves inverting large Gram matrices of certain basis functions (in which they use the Legendre polynomials) under additional complexity-reducing assumptions on the covariates $X$. It will be interesting to investigate if the sHOIF estimators can be generalized to the CATE estimation problems and stabilize their SORL or even HORL estimators.
\item Another important open problem was also mentioned in \citet{van2014higher, liu2020nearly, liu2021assumption}. To define HOIFs for DRFs, one needs to choose a set of $k$-dimensional basis functions $\bar{{\mathsf{z}}}_{k}$ or an approximation kernel $K$ of the Kronecker delta function, ideally in prior to the data analysis. However, such a strategy seems to go against the current data analytic paradigm, which strongly advocates learning representations (e.g. in the form of bases or kernels) adaptively from data rather than choosing some fixed bases/frames {\it a priori}. Prominent examples include DNNs, autoencoders, and GANs. It is thus interesting to construct HOIF estimators along different basis directions and then select one or aggregate all, guided by certain optimality criterion. We leave this important and difficult problem to future endeavor.
\item It will be interesting to also derive HOIFs and sHOIFs for identifiable causal effect functionals in graphical models with latent variables \citep{bhattacharya2022semiparametric} and implicitly defined functionals \citep{robins2016technical, ai2021unified} in general semiparametric regression problems for improved quality of estimation and statistical inference, which however requires extension of the current work to $U$-processes, a much more difficult research problem that we are working on in a separate paper.
\end{enumerate}
\bibliographystyle{plainnat}
\bibliography{Master.bib}
\newpage