Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
77,346 characters · 39 sections · 87 citation commands
Direct Debiased Machine Learning via Bregman Divergence Minimization
This study considers parameters of interest that depend on regression functions, such as treatment effects and policy effects. Empirical analysis often employs machine learning methods to estimate regression functions and then plugs the estimates into the estimating equations for the parameters of interest. However, while machine learning methods are effective in important situations, such as regression with high-dimensional covariates and complex regression functions, they often yield bias in the estimating equations, which prevents us from guaranteeing root-n convergence for the resulting estimator of the parameter of interest.
To reduce this bias, debiased machine learning methods have been investigated Chernozhukov2018doubledebiased. In debiased machine learning, we utilize Neyman orthogonal estimating equations, under which we can asymptotically eliminate the bias caused by regression function estimation with cross-fitting and mild convergence conditions for the regression function estimators. Chernozhukov2022automaticdebiased develops the automatic debiased machine learning framework, which allows us to define and estimate the Riesz representer without specifying a particular functional form.
The Neyman orthogonal scores usually require estimation of the Riesz representer to debias initial estimates of the regression functions. Therefore, we need to estimate the nuisance parameters, the regression functions, and the Riesz representer to construct an estimator of the parameter of interest. Although the Neyman orthogonal score mitigates the bias problem, high-quality estimation of the nuisance parameters remains important. First, even with the Neyman orthogonal estimating equation, to obtain asymptotic efficiency, we must construct nuisance parameter estimators that satisfy certain convergence rate conditions. Second, even if these conditions hold asymptotically, there is room to improve finite-sample performance.
For example, in the treatment effect estimation literature, various methods have been proposed to improve the accuracy of treatment effect estimation. In treatment effect estimation, the Riesz representer usually depends on the propensity score, the probability of assigning treatment. Such a probability is often estimated via a logistic regression model with maximum likelihood estimation. Imai2013estimatingtreatment proposes covariate balancing propensity scores, which fit the logistic regression model to match the first moment of the covariates between treated and control groups. Hainmueller2012entropybalancing instead proposes not using such an explicit model and estimates weight parameters under covariate balancing conditions. These covariate balancing methods have been further extended by subsequent studies Zhao2019covariatebalancing. From the debiased machine learning literature, Chernozhukov2024automaticdebiased establishes Riesz regression, which aims to directly estimate the Riesz representer. Note that BrunsSmith2025augmentedbalancing shows that Riesz regression can be derived as the dual of stable matching by Zubizarreta2015stableweights, and we also show that the formulation is identical to least-squares importance fitting (LSIF) in the density-ratio estimation literature Kanamori2009aleastsquares. For regression function estimation, vanderLaan2011targetedlearning proposes Super Learner and targeted maximum likelihood estimation (TMLE). Super Learner estimates the regression function, and TMLE debiases its initial regression function estimates by using the Riesz representer.
We propose direct debiased machine learning (DDML), a novel framework that unifies existing debiased machine learning frameworks and nuisance parameter estimation methods. Our proposed framework consists of Neyman targeted estimation and generalized Riesz regression, with the covariate balancing property. We first formulate estimation of the nuisance parameters as targeted estimation of the Neyman orthogonal score with the true nuisance parameters. We refer to this as Neyman targeted estimation, which involves estimation of the regression function and the Riesz representer. We then propose generalized Riesz regression, which estimates the Riesz representer by minimizing the Bregman divergence between the true Riesz representer and its model. The Bregman divergence is defined via a convex function chosen by us, and we can derive various estimation methods as special cases. For example, if we use the squared loss in the Bregman divergence, we obtain Riesz regression proposed in Chernozhukov2024automaticdebiased. If we choose a well-designed convex function (see Section (ref)) in treatment effect estimation, we can derive covariate balancing propensity scores, which are proven to be equal to empirical balancing Hainmueller2012entropybalancing by Zhao2019covariatebalancing. Through this argument, in treatment effect estimation, we point out that specific choices of the Bregman divergence guarantee the covariate balancing property, under which the Riesz representer balances covariates between the treated and control groups.
We summarize these important elements below:
Our proposed Riesz representer estimation method is novel, not merely a generalization of existing methods. We define the objective as weighted risk minimization to improve performance. Existing methods such as Riesz regression and covariate balancing propensity scores are special cases where we use identity weights. Using the proposed method, we demonstrate concrete implementations in several applications, such as Average Treatment Effect (ATE) estimation, ATT estimation, average marginal effect estimation, and covariate shift adaptation.
Furthermore, in the course of our arguments, we find the following points:
Thus, our DDML framework includes various existing methods, such as Riesz regression, covariate balancing propensity scores, targeted maximum likelihood estimation, density-ratio estimation, and nearest neighbor matching, as special cases. Also see Figure (ref).
In Section (ref), we first propose a method for ATE estimation based on MSE minimization targeting the oracle score functions. In Section (ref), we formulate the problem more generally and establish the direct debiased machine learning framework. Direct debiased machine learning consists of two important components, the Neyman targeted estimation and generalized Riesz regression, which are introduced in Sections (ref) and (ref), respectively. We demonstrate applications of our framework to ATE estimation, Average Treatment Effect on the Treated (ATT) estimation, Average Marginal Effect (AME) estimation, and covariate shift adaptation in Sections (ref)--(ref). In Section (ref), we define automatic covariate balancing, which bridges Riesz regression and covariate balancing methods.
This work is partially based on our earlier paper, Kato2025directbias. For several theoretical and simulation study results, also see Kato2025directbias.
Our work is mainly based on the literature on debiased machine learning, covariate balancing, TMLE, and direct density-ratio estimation. As briefly explained above and shown in the main text, these strands of literature essentially discuss the same topics from different perspectives. In this section, we introduce existing studies in these literatures.
Efficient estimation of semiparametric models has been intensively studied in various fields, including statistics, economics, machine learning, and epidemiology. A powerful tool for discussing the optimality of semiparametric models is the asymptotic efficiency bound for regular estimators Bickel1998efficientadaptive,Vaart1998asymptoticstatistics. It is known that efficient estimators, whose asymptotic variance matches the efficiency bound, are regular and asymptotically linear (RAL) for the semiparametric efficient influence function. To construct such efficient estimators, existing studies suggest one-step bias correction Vaart2002semiparametricstatistics, estimating equations, or TMLE, given estimated nuisance parameters. These estimators are constructed via the efficient influence function or the Riesz representer, defined later, and under certain conditions they are shown to be efficient. The conditions usually require that the estimators of the nuisance parameters satisfy complexity and convergence rate conditions. The complexity condition is often represented by the Donsker condition, and if it is not satisfied, we typically apply sample splitting to construct the estimators Klaassen1987consistentestimation,Zheng2011crossvalidatedtargeted.
These arguments about constructing efficient estimators have been systematized as debiased machine learning Chernozhukov2018doubledebiased,Chernozhukov2022automaticdebiased. Chernozhukov2018doubledebiased packages the estimating equation approach with sample splitting as Neyman orthogonal estimating equations with cross-fitting. In this framework, we first estimate the nuisance parameters using cross-fitting and plug these estimators into the Neyman orthogonal score. Then, we estimate the parameter of interest by finding a parameter such that the score is zero. The Neyman orthogonal score depends on the Riesz representer, and Chernozhukov2022automaticdebiased shows that we can define the Neyman orthogonal estimating equation without explicitly knowing the Riesz representer. Chernozhukov2024automaticdebiased develops Riesz regression to estimate the Riesz representer under the squared loss. Note that the Neyman orthogonal score usually corresponds to the semiparametric efficient score, and the Riesz representer corresponds to the bias-correction term used in one-step bias correction and clever covariates in TMLE.
Semiparametric models are often used for causal inference, where the main task is to estimate treatment effects. For better estimation of treatment effects, estimation of the propensity score, the probability of assigning treatment, has been intensively studied. For example, the propensity score can be used to eliminate selection bias\footnote{The meanings of bias and debias here differ from those in debiased machine learning.}. In addition, the Riesz representer usually includes the inverse of the propensity score. Thus, better estimation of the propensity score is expected to improve treatment effect estimation. The propensity score can also be interpreted as weights for the outcomes.
Randomized controlled trials are the gold standard for causal inference, where treatments are assigned while maintaining balance between treatment groups. However, they are not always feasible, and we aim to estimate causal effects from observational data, where there often exists imbalance between treatment groups. To correct the imbalance, propensity scores or balancing weights have been proposed. Covariate balancing is a popular approach for propensity score or balancing weight estimation. The propensity score has been proven to be a balancing score, and based on this property, existing studies propose estimating the propensity score or the weights so that the weighted covariate moments between treated and control groups match. Imai2013estimatingtreatment proposes estimating the propensity score by matching first moments, and Hazlett2020kernelbalancing extends the method to higher-moment matching by mapping the covariates into a high-dimensional space with basis functions. On the other hand, estimation methods that do not directly specify the propensity score model have also been proposed. Such methods are called empirical balancing and include entropy balancing Hainmueller2012entropybalancing and stable matching Zubizarreta2015stableweights. These two approaches appear different but are essentially the same via a duality relationship, as proven by Zhao2019covariatebalancing and BrunsSmith2025augmentedbalancing.
Our proposed framework is also inspired by the direct density-ratio estimation framework Sugiyama2012densityratio. We refer to the ratio between two densities as the density ratio. Density-ratio estimation has attracted considerable attention as an essential task in various machine learning problems, such as regression under covariate shift Shimodaira2000improvingpredictive,Reddi2015doublyrobust,Kato2024doubledebiasedcovariateshift, learning with noisy labels Liu2014classificationwith,Fang2020rethinkingimportance, anomaly detection Smola2009relativenovelty,Hido2008inlierbased,Abe2019anomalydetection, two-sample testing Keziou2005testof,Kanamori2010fdivergence,Sugiyama2011leastsquarestwosample, change-point detection Kawahara2009changepointdetection, causal inference Uehara2020offpolicy, recommendation systems Togashi2021densityratiobased. While the density ratio can be estimated by separately estimating each density, such an approach may magnify estimation errors by compounding two independent estimations. To address this issue, end-to-end, direct density-ratio estimation methods have been studied, including moment matching Huang2007correctingsample,Gretton2009covariateshift, probabilistic classification Qin1998inferencesfor,Cheng2004semiparametricdensity, density matching Nguyen2010estimatingdivergence, density-ratio fitting Kanamori2009aleastsquares, and positive–unlabeled learning Kato2019learningfrom. It is also known that when complicated models such as neural networks are used for this task, the loss function can diverge in finite samples Kiryo2017positiveunlabeledlearning. Therefore, density-ratio estimation methods with neural networks have been investigated Kato2021nonnegativebregman,Rhodes2020.
To explain a basic idea of our DDML framework, we introduce an example in ATE estimation. Let $X = (D, Z)$, where $D \in \{1, 0\}$ denotes a binary treatment indicator, and $Z \in {\mathcal{Z}}$ denotes covariates with covariate space ${\mathcal{Z}}$. Let $Y = D \cdot Y(1) + (1 - D) \cdot Y(0)$, where $(Y(1), Y(0))$ are conditionally independent of treatment $D$ given covariates $Z$. We denote the pair of $Y$ and $X$ by $W = (X, Y)$.
The parameter of interest is the ATE, defined as \[\tau_0 = \mathbb{E}\left[Y(1) - Y(0)\right].\] Our goal is to estimate $\tau_0$ by using observations $\{W_i\}^n_{i=1}$, where $W_i$ is an i.i.d. copy of $W$. Let $\gamma(X)\coloneqq \mathbb{E}\left[Y\mid X = (D,Z)\right]$ be the conditional expected outcome $Y$ given $X$. Let $\pi_0(Z) \coloneqq \Pr(D=1\mid Z)$ be the propensity score, the probability of assigning treatment $1$. We assume that $0 < \pi_0(Z) < 1$ almost surely.
We can estimate the ATE by using the Augmented Inverse Probability Weighting (AIPW) estimator, defined as follows: \[\widehat{\tau}^{\text{AIPW}} \coloneqq \frac{1}{n}\sum^n_{i=1}h^{\mathrm{AIPW}}(W; \widehat{\eta}),\] where for $\eta = (\gamma, \pi)$, $h^{\mathrm{AIPW}}(W; \eta)$ is defined as
and $\widehat{\eta} = (\widehat{\gamma}, \widehat{\pi})$ is a pair of estimators of $\gamma_0$ and $\pi_0$. This estimator has various desirable theoretical properties. In particular, its asymptotic efficiency and double robustness play important roles in ATE estimation.
In ATE estimation, the efficient score function is given as \[\psi(W; \eta_0) = h^{\text{AIPW}}(W; \eta_0) - \tau_0,\] where $\eta_0 \coloneqq (\gamma_0, \pi_0)$ is the pair of the true functions $\gamma_0$ and $\pi_0$. Regular and asymptotically linear (RAL) estimators based on the efficient score are known to be asymptotically efficient, that is, among regular estimators, there is no estimator whose asymptotic variance is lower than that of the regular and asymptotically linear (RAL) estimator, $\widehat{\tau}$.
In order to use the AIPW estimator in practice, we need to estimate $\eta_0$ and replace it with its estimator $\widehat{\eta}$ in the estimator. Such functions $\gamma_0$ and $\pi_0$ are called nuisance parameters.
Define an oracle as \[\widetilde{\tau} \coloneqq \frac{1}{n}\sum^n_{i=1}h^{\text{AIPW}}(W_i; \eta_0).\] Then, we decompose the error between the oracle $\tau_0$ and $\widehat{\tau}^{\text{AIPW}}$ as
where $\tau_0 - h^{\text{AIPW}}(W_i; \eta_0)$ corresponds to the efficient score, and $h^{\text{AIPW}}(W_i; \eta_0) - h^{\text{AIPW}}(W_i; \widehat{\eta})$ is the error incurred when using estimated nuisance parameters instead of the true ones.
A straightforward approach is to minimize the empirical mean squared error (MSE): \[\widehat{\eta} \coloneqq \operatorname*{arg\,min}_{\eta \in {\mathcal{M}}}\frac{1}{n}\sum^n_{i=1}\Big(h^{\text{AIPW}}(W_i; \eta_0) - h^{\text{AIPW}}(W_i; \eta)\Big)^2,\] where ${\mathcal{M}} = {\mathcal{G}} \times \Pi$ is a model of $\eta$, and ${\mathcal{G}}$ and $\Pi$ are models of $\gamma_0$ and $\pi_0$. Since this optimization problem is infeasible because $\eta_0$ is unknown, we consider an alternative objective function.
To construct an alternative, we estimate $\eta_0$ by minimizing the population MSE:
Surprisingly, for this population MSE minimization, we can obtain an equivalent objective function that does not include the unknown $\eta_0$.
\paragraph{Estimation of the propensity score.} Given $\gamma$, we estimate $\pi_0$ by the following mean squared error minimization:
where recall that $\Pi$ is a model of $\pi_0$. If $\pi_0 \in \Pi$, then $\pi^* = \pi_0$ holds and $\eta_0 = (\gamma_0, \pi_0)$.
This MSE minimization problem can be written in a form that does not include the unknown $\pi_0$:
The proof is shown in Section (ref). If $\pi_0 \in \Pi$, we have $\pi^\dagger(\gamma) = \pi_0$ for any $\gamma$.
\paragraph{Estimation of the regression function.} Given $\pi$, we aim to estimate $\gamma_0$ as
where recall that ${\mathcal{G}}$ is a model of $\gamma_0$. If $\gamma_0 \in {\mathcal{G}}$, then $\gamma^* = \gamma_0$ holds. This is equivalent to \[\gamma^*(\pi) \coloneqq \operatorname*{arg\,min}_{\gamma \in {\mathcal{G}}}\mathbb{E}\Bigg[\Bigg(\left(1 - \frac{D}{\pi(Z)}\right)\Big(\gamma_0((1, Z)) - \gamma((1, Z))\Big) - \left(1 - \frac{1 - D}{1 - \pi(Z)}\right)\Big(\gamma_0((0, Z)) - \gamma((0, Z))\Big)\Bigg)^2\Bigg].\] Unfortunately, unlike the propensity score estimation, we cannot eliminate $\gamma_0$ in the optimization. Hence, we consider minimizing an upper bound of the objective function:
Here, we derive the above upper bound as follows:
\paragraph{Joint estimation of the nuisance parameters.} Based on the above arguments, we present the first basic two-step algorithm. In the subsequent subsection, we provide an iterative algorithm as an extension of this algorithm.
First, we estimate the nuisance parameters using the following two-step procedure:
We can also estimate the nuisance parameters as
In this section, we derive Theorem (ref); that is, the optimization problem ((ref)) is equivalent to ((ref)). While ((ref)) includes the unknown $\pi_0$, ((ref)) does not include $\pi_0$.
First, we have
We compute $\mathbb{E}\Big[\Big(h^{\mathrm{AIPW}}(W; (\gamma_0, \pi_0)) - h^{\mathrm{AIPW}}(W; (\gamma_0, \pi))\Big)^2\Big]$ as
We compute $\mathbb{E}\Big[\Big(h^{\mathrm{AIPW}}(W; (\gamma_0, \pi)) - h^{\mathrm{AIPW}}(W; (\gamma, \pi))\Big)^2\Big]$ as
Therefore, we have
By generalizing the procedure in Section (ref), we formulate a direct debiased machine learning framework.
We consider data that consist of observations $\{W_i\}_{i=1}^n$, where each $W_i$ is an i.i.d. copy of $W$ following a distribution $F_0$. The variable $W$ includes an outcome variable $Y$ and regressors $X$; that is $W = (X, Y)$. We focus on parameters that depend on the regression function (conditional mean) of $Y$ given $X$. We denote the true regression function under $F_0$ by $\gamma_0(x) \coloneqq \mathbb{E}\left[Y\mid X = x\right]$, while we denote a model of $\gamma_0$ by $\gamma$.
We aim to estimate a parameter of interest of the form \[\theta_0 \coloneqq \mathbb{E}\big[m(W, \gamma_0)\big],\] where $m(w, \gamma)$ is a functional that depends on a data observation $w$ and a candidate regression function $\gamma$.
For simplicity, we assume that the expected functional $\gamma \mapsto \mathbb{E}\left[m(W, \gamma)\right]$ is linear and continuous in $\gamma$, which implies that there is a constant $C > 0$ such that $\mathbb{E}\left[m(W, \gamma)\right]^2 \leq C \mathbb{E}\left[\gamma(X)^2\right]$ holds for all $\gamma$ with $\mathbb{E}\left[\gamma(X)^2\right] < \infty$. Then, from the Riesz representation theorem, there exists a function $v_m$ with $\mathbb{E}\left[v_m(X)^2\right] < \infty$ such that \[\mathbb{E}\left[m(W, \gamma)\right] = \mathbb{E}\left[v_m(X)\gamma(X)\right]\] for all $\gamma$ with $\mathbb{E}\left[\gamma(X)^2\right] < \infty$. This formulation follows Chernozhukov2022automaticdebiased and can be generalized to non-linear maps $\gamma \mapsto \mathbb{E}\left[m(W, \gamma)\right]$. The function $v_m$ is often referred to as the Riesz representer.
Let $\eta_0=(\gamma_0,\alpha_0)$, where $\alpha_0$ is the Riesz representer associated with $m$. The Neyman orthogonal score (often equal to the semiparametric efficient influence function) is \[ \psi(W;\eta,\theta) \coloneqq m(W,\gamma) + \alpha(X)\{Y-\gamma(X)\} - \theta. \] It obeys $\mathbb{E}\left[\psi(W;\eta_0,\theta_0)\right]=0$ and the Neyman orthogonality (Gateaux derivative) with respect to $\eta$ at $\eta_0$ vanishes: \[ \partial_{\eta}\mathbb{E}\left[\psi(W;\eta,\theta_0)\right]\big|_{\eta=\eta_0}=0. \] Orthogonality ensures that first-order errors from estimating $\eta_0$ do not inflate the asymptotic distribution of the final $\theta$-estimator, provided cross-fitting and mild rate conditions hold; see Chernozhukov2018doubledebiased,Chernozhukov2022automaticdebiased.
We propose DDML, which estimates $\eta_0=(\gamma_0,\alpha_0)$ end-to-end by directly targeting the oracle score $\psi(W;\eta_0,\theta_0)$. DDML has two pillars:
We obtain $\widehat{\theta}$ as the value satisfying \[\frac{1}{n}\sum^n_{i=1}\psi(W_i;\widehat{\eta},\widehat{\theta})=0,\] where $\widehat{\eta}$ denotes estimates of $\eta_0$ plugged into the Neyman orthogonal estimating equations.
We can show minimax-optimal convergence rates for the nuisance parameter estimators. For details, see Kato2025directbias and related studies such as Chernozhukov2024automaticdebiased, Kanamori2012statisticalanalysis, and Kato2021nonnegativebregman.
If the estimators satisfy suitable convergence rate conditions and either the Donsker condition holds or cross-fitting is employed, the resulting estimator $\widehat{\theta}$ is asymptotically normal with asymptotic variance equal to the efficiency bound.
Under specific choices of Riesz representer models, we can attain covariate balancing without explicitly solving a covariate balancing objective. For example, in ATE estimation, consider the logistic model for the propensity score \[\pi_\beta(Z) \coloneqq \frac{1}{1 + \exp\big(\beta^\top \Phi(Z)\big)},\] and define the Riesz representer model as \[\alpha_\beta(X) = \frac{\mathbbm{1}[D_i = 1]}{\pi_\beta(X)} - \frac{\mathbbm{1}[D_i = 0]}{1 - \pi_\beta(X)}.\] In this case, if we use the KL divergence, we obtain the following covariate balancing: \[\frac{1}{n}\sum^n_{i=1}\alpha_\beta(X_i) \Phi(Z_i) = \frac{1}{n}\sum^n_{i=1}\left(\frac{\mathbbm{1}[D_i = 1]}{\pi_{\widehat{\beta}}(Z_i)} \Phi(Z_i) - \frac{\mathbbm{1}[D_i = 0]}{1 - \pi_{\widehat{\beta}}(Z_i)} \Phi(Z_i)\right) = 0.\] We call this property automatic covariate balancing.
Given any candidates $\eta = (\gamma,\alpha)$ for $\eta_0$ and $\theta^*$ for $\theta_0$, we have
This is derived from the following computation:
Thus, targeting $\eta$ amounts to reducing two interpretable components: (i) estimation error for $\alpha_0$, and (ii) estimation error for $\gamma$. Our algorithm alternates or jointly optimizes a Bregman objective (for $\alpha$) and a weighted-risk objective (for $\gamma$).
We estimate the Riesz representer by minimizing the error \[\left|\big(\alpha_0(X) - \alpha(X)\big) \big(Y - \gamma_0(X)\big) - \alpha(X)\big(\gamma_0(X) - \gamma(X)\big)\right|.\] Since $\alpha(X)\big(\gamma_0(X) - \gamma(X)\big) = 0$ if $\gamma = \gamma_0$, we focus on $\big(\alpha_0(X) - \alpha(X)\big) \big(Y - \gamma_0(X)\big)$.
\paragraph{Generalized Riesz Regression.} There exist several measures for evaluating $\big(\alpha_0(X) - \alpha(X)\big) \big(Y - \gamma_0(X)\big)$. In this study, we propose using the Bregman divergence to measure the discrepancy between $\alpha_0(X)$ and $\alpha(X)$, weighted by $\mathbb{E}\left[Y - \gamma_0(X)\mid X\right]^2$, the conditional variance of $Y$ given $X$. We explain the details of the Riesz representer estimation via the Bregman divergence in Section (ref). We refer to this estimation method as Generalized Riesz regression. Here, we briefly introduce the idea of Bregman divergence estimation.
The Bregman divergence depends on a convex function, and by specifying this convex function, we obtain various divergence metrics, including loss functions such as the squared loss and log loss. For example, under the squared loss at the population level, we estimate the Riesz representer as \[\alpha^* \coloneqq \operatorname*{arg\,min}_{\alpha \in {\mathcal{A}}}\mathbb{E}\left[\big(\alpha_0(X) - \alpha(X)\big)^2 \big(Y - \gamma_0(X)\big)^2\right],\] where ${\mathcal{A}}$ is a set of candidates of $\alpha_0$. If $\alpha_0 \in {\mathcal{A}}$, then $\alpha^* = \alpha_0$ holds. This objective includes the unknown $\alpha_0(X)$ (we replace $\gamma_0$ with a consistent estimator, as it only weights the loss). Although the presence of $\alpha_0(X)$ seems to make the optimization infeasible, we can obtain an equivalent objective that does not include $\alpha_0(X)$: \[\alpha^* \coloneqq \operatorname*{arg\,min}_{\alpha \in {\mathcal{A}}}\mathbb{E}\Big[ - g(\alpha(X)) + \partial g(\alpha(X)) \alpha(X) - m\big(\partial g(\alpha(X))\big)\Big],\] where $g$ differentiable and strictly convex function. Estimation of $\alpha_0$ based on this optimization is referred to as Riesz regression in automatic debiased machine learning Chernozhukov2024automaticdebiased; when $\alpha_0(X)$ is related to a density ratio, it is a special case of LSIF in the density-ratio estimation literature Kanamori2009aleastsquares.
The Bregman divergence also includes other loss functions. In particular, under specific choices, we obtain the covariate balancing propensity score Imai2013estimatingtreatment and empirical balancing Hainmueller2012entropybalancing by bypassing the tailored loss argument of Zhao2019covariatebalancing, which also shows the equivalence between covariate balancing propensity scores and empirical balancing.
\paragraph{Automatic Covariate Balancing.} By selecting different convex functions in the Bregman divergence, Generalized Riesz regression induces different losses. In treatment effect estimation, certain choices yield a desirable property we call automatic covariate balancing. We provide details in Section (ref).
There are two main approaches to estimating the regression functions.
\paragraph{Case 1: Double Estimation} Estimate the regression function $\gamma_0$ and the Riesz representer separately, using the two-step or iterative procedures in Section (ref).
\paragraph{Case 2: TMLE} After double estimation, given estimates $\widehat{\gamma}^{(0)}$ and $\widehat{\alpha}$, update the regression function as \[\widehat{\gamma}^{(1)} \coloneqq \widehat{\gamma}^{(0)} + \frac{\sum_{i=1}^n \widehat{\alpha}(X_i)\big(Y_i - \widehat{\gamma}(X_i)\big)}{\sum_{i=1}^n \widehat{\alpha}(X_i)^2}\widehat{\alpha}(X_i).\] Then estimate the parameter of interest as \[\widehat{\theta}^{\text{TMLE}} \coloneqq \frac{1}{n}\sum^n_{i=1}\left(\widehat{\gamma}^{(1)}((1, Z_i)) - \widehat{\gamma}^{(1)}((0, Z_i))\right) = \frac{1}{n}\sum^n_{i=1}m(X_i, \widehat{\gamma}^{(1)}(X_i)).\] Note that this update is derived as the solution in $\epsilon$ of \[\sum^n_{i=1}\left(Y_i - \left(\widehat{\gamma}^{(0)}(X_i) + \epsilon\widehat{\alpha}(X_i)\right)\right))\right)} = 0,\] which is given as $\widehat{\epsilon} \coloneqq \frac{\sum_{i=1}^n \widehat{\alpha}(X_i)\{Y_i-\widehat{\gamma}^{(0)}(X_i)\}} {\sum_{i=1}^n \widehat{\alpha}(X_i)^2},\qquad \widehat{\gamma}^{(1)}(x)=\widehat{\gamma}^{(0)}(x)+\widehat{\epsilon}\widehat{\alpha}(x)$.
Under this update, for $\widehat{\eta} = (\widehat{\gamma}, \widehat{\pi})$ and $\widehat{\theta}^{\text{TMLE}} = m(X_i, \widehat{\gamma}^{(1)})$, we have
Therefore,
Thus, the estimation error between the oracle and plug-in Neyman orthogonal scores boils down to the estimation error of the Riesz representer. Hence, Neyman targeted estimation reduces to estimating the Riesz representer.
The essential idea in Section (ref) is to estimate $\gamma_0$ and $\pi_0$ by minimizing the error between $h^{\mathrm{AIPW}}(W; (\gamma_0, \pi_0)) - h^{\mathrm{AIPW}}(W; (\gamma, \pi))$. While we measure the error using the squared loss in Section (ref), here we measure the error using the Bregman divergence. Bregman divergence is a general measure of deviation between two quantities. It not only includes the MSE but also includes the KL divergence as a special case, where the KL divergence can be interpreted as a likelihood. We then formulate the nuisance parameter estimation problem as a Bregman divergence minimization problem.
Through the lens of Bregman divergence minimization, we find several useful connections. As noted above, when using the squared loss, the target estimation method includes Riesz regression in Chernozhukov2024automaticdebiased as a special case. In addition, when using the KL divergence, we can interpret the method as the tailored loss in covariate balancing Zhao2019covariatebalancing, which includes the covariate balancing propensity scores and entropy balancing as special cases. Furthermore, techniques used in our arguments are based on the density-ratio estimation literature. For example, Riesz regression is essentially the same formulation as LSIF in Kanamori2009aleastsquares, and the generalization by Bregman divergence appears in Sugiyama2011densityratio. Thus, Riesz regression, covariate balancing, and density-ratio estimation essentially share the same formulation.
Let $g\colon {\mathbb{R}} \to {\mathbb{R}}$ be a differentiable and strictly convex function. Given $x \in {\mathcal{X}}$, define the Bregman divergence between $\alpha_0(x), \alpha(x) \colon {\mathcal{X}} \to {\mathbb{R}}$ as \[\text{BR}^\dagger_g\big(\alpha_0(x)\mid \alpha(x)\big) \coloneqq g(\alpha_0(x)) - g(\alpha(x)) - \partial g(\alpha(x)) \big(\alpha_0(x) - \alpha(x)\big),\] where $\partial g$ denotes the derivative of $g$. Then define the average Bregman divergence as \[\text{BR}^\dagger_g\big(\alpha_0\mid \alpha\big) \coloneqq \mathbb{E}\Big[g(\alpha_0(X)) - g(\alpha(X)) - \partial g(\alpha(X)) \big(\alpha_0(X) - \alpha(X)\big)\Big].\] We estimate $\alpha_0$ by \[\alpha^* = \operatorname*{arg\,min}_{\alpha\in {\mathcal{A}}} \text{BR}^\dagger_g\big(\alpha_0\mid \alpha\big).\] where ${\mathcal{A}}$ is a set of candidates of $\alpha_0$. If $\alpha_0 \in {\mathcal{A}}$, then $\alpha^* = \alpha_0$ holds.
Although $\alpha_0$ is unknown, we can define an equivalent optimization problem without using $\alpha_0$: \[\alpha^* = \operatorname*{arg\,min}_{r\in {\mathcal{A}}} \text{BR}_g\big(\alpha\big),\] where \[\text{BR}_g\big(\alpha\big) \coloneqq \mathbb{E}\Big[ - g(\alpha(X)) + \partial g(\alpha(X)) \alpha(X) - m\big(\partial g(\alpha(X))\big)\Big].\] Here, we used the linearity of $m$ \[ \mathbb{E}\big[\partial g(\alpha(Z))\alpha_0(Z)\big]=\mathbb{E}\big[m\{\partial g(\alpha(Z))\}\big]. \] By choosing different $g$, we derive various objectives for Riesz representer estimation, including Riesz regression.
We estimate the Riesz representer $\alpha_0$ by minimizing an empirical Bregman divergence: \[ \widehat{\alpha}_n \coloneqq \operatorname*{arg\,min}_{\alpha \in {\mathcal{H}}}\widehat{\text{BR}}_g\big(\alpha\big) + \lambda J(\alpha), \] where $J(\alpha)$ is some regularization function, and \[ \widehat{\text{BR}}_g(\alpha) \coloneqq \frac{1}{n}\sum^n_{i=1}\Big( - g(\alpha(X_i)) + \partial g(\alpha(X_i)) \alpha(X_i) - m\big(\partial g(\alpha(X_i))\big)\Big). \]
We can also define a weighted version of the Bregman divergence minimization: \[ \widehat{\alpha}_n \coloneqq \operatorname*{arg\,min}_{\alpha \in {\mathcal{H}}}\widehat{\text{BR}}^W_g\big(\alpha\big) + \lambda J(\alpha), \] where \[ \widehat{\text{BR}}^W_g(\alpha) \coloneqq \frac{1}{n}\sum^n_{i=1}\Big( - g(\alpha(X_i)) + \partial g(\alpha(X_i)) \alpha(X_i) - m\big(\partial g(\alpha(X_i))\big)\Big)\Big(Y_i - \widehat{\gamma}(X_i)\Big)^2, \] and $\widehat{\gamma}$ is some estimate of $\gamma_0$. This weighted version multiplies each summand by $\{Y-\widehat{\gamma}(D,Z)\}^2$ to stabilize variance and improve small-sample behavior.
We revisit the ATE estimation. The parameter of interest is the ATE, defined as $\tau_0 = \mathbb{E}\left[m^{\text{ATE}}(W, \gamma_0)\right]$, where \[m^{\text{ATE}}(W, \gamma_0) \coloneqq \gamma_0((1, Z)) - \gamma_0((0, Z)),\] where $\gamma_0((d, Z)) = \mathbb{E}\left[Y\mid D = d, Z\right]$. In this setting, the Riesz representer is \[\alpha^{\text{ATE}}_0(X) \coloneqq \frac{D}{\pi_0(Z)} - \frac{1 - D}{1 - \pi_0(Z)},\] where $\pi_0(Z) = \Pr(D = 1\mid Z)$ is the propensity score. Then the Neyman orthogonal score is \[\psi(W; \eta_0, \tau_0) = h^{\text{AIPW}}(W; \eta_0) - \tau_0,\] where recall that $h^{\mathrm{AIPW}}(W; \eta) = \left(\frac{D}{\pi(Z)} - \frac{1 - D}{1 - \pi(Z)}\right)\Big(Y - \gamma(X)\Big) + \gamma((1, Z)) - \gamma((0, Z))$. Then, the estimation error in the Neyman targeted step is given as
We first introduce the squared function for the Bregman divergence, given as \[g^{\text{LS}}(\alpha) = (\alpha - 1)^2.\] Under this choice of $g$, we have \[\text{BR}_{g^{\text{LS}}}\big(\alpha\big) = \mathbb{E}\left[- 2\big(\alpha((1, Z)) + \alpha((0, Z))\big) + \mathbbm{1}[D = 1]\alpha((1, Z))^2 + \mathbbm{1}[D = 0]\alpha((0, Z))^2 \right].\] Then, we estimate $\alpha_0$ by minimizing the empirical objective: \[ \widehat{\alpha} \coloneqq \operatorname*{arg\,min}_{\alpha \in {\mathcal{A}}}\widehat{\text{BR}}_{g^{\mathrm{LS}}}\big(\alpha\big), \] where \[\widehat{\text{BR}}_{g^{\mathrm{LS}}}\big(\alpha\big) \coloneqq \frac{1}{n}\sum^n_{i=1}\left(- 2\big(\alpha((1, Z_i)) + \alpha((0, Z_i))\big) + \mathbbm{1}[D_i = 1]\alpha((1, Z_i))^2 + \mathbbm{1}[D_i = 0]\alpha((0, Z_i))^2 \right).\] This objective function matches LSIF for density-ratio estimation proposed in Kanamori2009aleastsquares,Kanamori2012statisticalanalysis and the Riesz regression for automatic debiased machine learning proposed in Chernozhukov2024automaticdebiased. Furthermore, under specific choice of ${\mathcal{A}}$, a model of $\alpha_0$, this objective function also yields nearest neighbor matching, as Kato2025nearestneighbor shows the equivalence between LSIF (equivalently, Riesz regression) and density-ratio estimation method proposed in Lin2023estimationbased. Kato2022learningcausal also applies LSIF for nonparametric instrumental variable regression.
For $\alpha$, we can use various models, such as random forests and neural networks. Note that there can be occurred train-loss hacking problems and appropriate correction methods are also required Rhodes2020telescopingdensityratio,Kato2021nonnegativebregman,Kiryo2017positiveunlabeledlearning.
We can also derive (constrained) maximum likelihood estimation. Let us consider the convex function \[g^{\mathrm{KL}}(\alpha) = |\alpha|\log |\alpha| - |\alpha|.\] This choice connects the Bregman divergence to the KL divergence. Here, we have \[\partial g^{\mathrm{KL}}(\alpha(X)) =
.\]
Substituting this $g^{\mathrm{KL}}$ into the Bregman divergence gives \[\text{BR}_{g^{\mathrm{KL}}}\big(\alpha\big) \coloneqq \mathbb{E}\big[ \operatorname{sign}(\alpha(X))\log|\alpha(X)| - \log \alpha((1, Z)) - \log \alpha((0, Z))\big].\] Then, we estimate $\alpha_0$ by minimizing the empirical objective: \[ \widehat{\alpha} \coloneqq \operatorname*{arg\,min}_{\alpha \in {\mathcal{A}}}\widehat{\text{BR}}_{g^{\mathrm{KL}}}\big(\alpha\big) + \lambda J(\alpha), \] where \[\widehat{\text{BR}}_{g^{\mathrm{KL}}}\big(\alpha\big) \coloneqq \frac{1}{n}\sum^n_{i=1}\left( \big|\alpha(X_i)\big| - \log \big|\alpha((1, Z_i))\big| - \log \big|\alpha((0, Z_i))\big|\right).\] This corresponds to unnormalized KL (UKL) divergence minimization in density-ratio estimation Sugiyama2011densityratio, which generalizes Kullback-Leibler Importance Estimation. Procedure Sugiyama2008directimportance,Nguyen2007estimatingdivergence.
\paragraph{Inverse propensity score modeling} Consider modeling the Riesz representer via a model of the inverse propensity score. Let $r((d, Z))$ model the inverse propensity score, that is, $r((1, Z))$ models $r_0((1, Z)) \coloneqq 1/\pi_0(Z)$ and $r((0, Z))$ models $r_0((0, Z)) \coloneqq 1/(1 - \pi_0(Z))$. Using this model, define \[\alpha_r((D, Z)) = \mathbbm{1}[D = 1]r((1, Z)) - \mathbbm{1}[D = 0]r((0, Z)).\]
By definition, we use $r$ such that $r(X) \in (1, \infty)$. Then, we have $\alpha((1, Z)) = r((1, Z)) \in (1, \infty)$ and $\alpha((0, Z)) = -r((0, Z)) \in (-\infty, -1)$. Therefore, we have
Then we obtain the objective \[\text{BR}_{g^{\mathrm{KL}}}\big(\alpha_r\big) \coloneqq \mathbb{E}\big[ -\log(r((1, Z))) - \log(r((0, Z))) + \mathbbm{1}[D_i = 1]r((1, Z_i)) + \mathbbm{1}[D_i = 0]r((0, Z_i))\big].\] Let ${\mathcal{R}}$ be a set of $r$. Then, we estimate
Since we do not know the expected value, by replacing it with an empirical estimate, we estimate $r_0$ in the Riesz representer by
where
Solving ((ref)) is equivalent to
This technique is known as Silverman's trick Silverman1982onestimation. For details, see Theorem 3.3 in Kato2023unifiedperspective. Replacing expectations with sample means yields the estimation problem
Next, we derive empirical balancing as a special case of Bregman divergence minimization. Empirical balancing is a specific form of covariate balancing and can be derived from a tailored loss function Zhao2019covariatebalancing.
Consider a convex function defined as \[g^{\mathrm{E}}(\alpha) = (|\alpha|-1)\log\left(\left|\alpha\right| - 1\right) - |\alpha|.\] This choice also yields another KL-type divergence. Note that $\alpha < 0$ or $\alpha > 1$ always holds, and for $\alpha \in (-\infty, 0)$ and $\alpha \in (1, \infty)$, $g^{\mathrm{E}}$ is convex with derivative \[\partial g^{\mathrm{E}}(\alpha) =
.\] Substituting this $g^{\mathrm{E}}$ and $\alpha(X) = \mathbbm{1}[D = 1]r((1, Z)) + \mathbbm{1}[D = 0]r((0, Z))$, we obtain
Then, we estimate $\alpha_0$ by minimizing the empirical objective: \[ \widehat{\alpha} \coloneqq \operatorname*{arg\,min}_{\alpha \in {\mathcal{A}}}\widehat{\text{BR}}_{g^{\mathrm{E}}}\big(\alpha\big) + \lambda J(\alpha), \] where \[\widehat{\text{BR}}_{g^{\mathrm{E}}}\big(\alpha\big) \coloneqq \frac{1}{n}\sum^n_{i=1}\Big(\log\left(|\alpha(X_i)| - 1\right) + |\alpha(X_i)| - \log\left(\alpha((1, Z_i)) - 1\right) - \log\left(-\alpha((0, Z_i)) - 1\right)\Big).\]
\paragraph{Inverse propensity score modeling} As in Section (ref), model the Riesz representer via the inverse propensity score. Let $r((d, Z))$ model the inverse propensity score, that is, $r((1, Z))$ models $r_0((1, Z)) \coloneqq 1/\pi_0(Z)$ and $r((0, Z))$ models $r_0((0, Z)) \coloneqq 1/(1 - \pi_0(Z))$. Based on this model, define \[\alpha_r(X) = \mathbbm{1}[D_i = 1]r((1, Z)) - \mathbbm{1}[D_i = 0]r((0, Z)).\] Under this model, we write the Bregman divergence as \[\text{BR}_{g^{\mathrm{E}}}\big(\alpha_r\big) \coloneqq \sum_{d \in \{1, 0\}}\mathbb{E}\Big[\mathbbm{1}[D = d]\Big(\log\left(r((d, Z)) - 1\right) + r((d, Z))\Big) - \log\left(r((d, Z)) - 1\right)\Big].\] We can simplify this Bregman divergence as
Then, we estimate $r_0$ in the Riesz representer by \[ \widehat{r}_n \coloneqq \operatorname*{arg\,min}_{r \in {\mathcal{R}}}\widehat{\text{BR}}_{g^{\mathrm{E}}}\big(\alpha_r\big), \] where the empirical Bregman divergence is
\paragraph{Connection to the tailored loss.} Since $r((1, Z)) = 1/\pi(X)$ and $r((0, Z)) = 1/(1 - \pi(X))$, we have $r((1, Z)) - 1 = 1/\big(r((0, Z)) - 1\big)$. Therefore, it holds that
This objective is equivalent to the tailored loss in Zhao2019covariatebalancing. From this objective, we can also derive empirical balancing Chan2015globallyefficient.
We consider the ATT, \[ \tau^{\text{ATT}}_0 \coloneqq \mathbb{E}\left[Y(1)-Y(0)\mid D=1\right] = \frac{1}{p}\mathbb{E}\left[ \gamma_0((1,Z))-\gamma_0((0,Z)) \mid D=1\right], \] where $p\coloneqq \Pr(D=1)$ and, as in Section (ref), $\gamma_0((d,Z))=\mathbb{E}\left[Y\mid (D, Z) = (d, Z)\right]$. Define the linear functional \[ m^{\text{ATT}}(W,\gamma)\coloneqq \frac{D}{p}\Big(\gamma((1,Z))-\gamma((0,Z))\Big),\qquad \tau^{\text{ATT}}_0 = \mathbb{E}\left[m^{\text{ATT}}(W,\gamma_0)\right]. \] Its Riesz representer $\alpha^{\text{ATT}}_0$ is characterized by $\mathbb{E}\big[m^{\text{ATT}}(W,\gamma)\big] = \mathbb{E}\big[\alpha^{\text{ATT}}_0(X)\gamma(X)\big]$ for all square–integrable $\gamma$. A standard calculation yields \[ \alpha^{\text{ATT}}_0(X)=\frac{D}{p}-\frac{1-D}{p}\cdot\frac{\pi_0(Z)}{1-\pi_0(Z)}, \] where $\pi_0(Z)=\Pr(D=1\mid Z)$.
The Neyman orthogonal score is \[ \psi^{\text{ATT}}\left(W;\eta,\tau^{\text{ATT}}\right)\coloneqq m^{\text{ATT}}(W,\gamma)+\alpha^{\text{ATT}}(X)\Big(Y-\gamma(X)\Big)-\tau^{\text{ATT}}, \] with $\eta=(\gamma,\alpha^{\text{ATT}})$. Plugging in estimates gives a DDML estimator solving \[ \frac{1}{n}\sum_{i=1}^n\psi^{\text{ATT}}\left(W_i;\widehat\eta,\widehat{\tau}^{\text{ATT}}\right)=0, \] which is root-$n$ and efficient under mild convergence–rate and complexity conditions (via the Donsker condition or cross–fitting), without requiring an explicit model for $\alpha^{\text{ATT}}_0$.
As in Section (ref), estimate $\alpha^{\text{ATT}}_0$ by minimizing an empirical Bregman risk \[ \widehat{\alpha}^{\text{ATT}}\in\operatorname*{arg\,min}_{\alpha\in{\mathcal{A}}}\ \frac1n\sum_{i=1}^n \Big\{-g(\alpha(X_i))+\partial g(\alpha(X_i))\alpha(X_i)-m^{\text{ATT}}\big(W_i,\partial g(\alpha(\cdot))\big)\Big\} +\lambda J(\alpha), \] optionally with the variance–stabilizing weight $\left(Y_i-\widehat\gamma(X_i)\right)^2$ as in Section (ref). Squared loss $g(\alpha)=(\alpha-1)^2$ recovers Riesz regression, while KL-type $g$ induces entropy–style balancing losses.
\paragraph{Squared loss (LSIF/Riesz regression).} With $g^{\text{LS}}(\alpha)=(\alpha-1)^2$, the empirical objective reduces to Riesz regression, which also coincides with LSIF. Furthermore, the dual matches the stable weights in Zubizarreta2015stableweights when linear models are used for $\alpha$.
\paragraph{KL-type loss (entropy balancing).} With $g^{\text{KL}}(\alpha)=(|\alpha| - 1)\log(|\alpha| - 1)-|\alpha|$, the dual problem enforces moment balance between treated units and reweighted controls (automatic covariate balancing), aligning with entropy–balancing–style weights but learned through the Riesz objective.
Let $X=(D,Z)$ with a (scalar) continuous treatment $D$, and define the AME as \[ \theta^{\text{AME}}_0 \coloneqq \mathbb{E}\left[\partial_d \gamma_0((D,Z))\right]. \] Here, linear functional is given as \[m^{\text{AME}}(W,\gamma)=\partial_d \gamma((D, Z)).\] The Riesz representer that satisfies $\mathbb{E}\left[m^{\text{AME}}(W,\gamma)\right]=\mathbb{E}\left[\alpha^{\text{AME}}_0(X)\gamma(X)\right]$ is the (negative) score of the joint density of $X = (D, Z)$ with respect to $d$: \[ \alpha^{\text{AME}}_0(X) = - \partial_d \log f_0((D,Z)), \] where $f_0(X)$ is the joint probability density of $X$.
The Neyman orthogonal score is given as \[ \psi^{\text{AME}}(W; \eta, \theta) = m^{\text{AME}}(W, \gamma)+\alpha^{\text{AME}}(X)\big(Y-\gamma(X)\big)-\theta. \] Then, we define the Neyman targeted estimation for estimating $\alpha^{\text{AME}}_0$ and $\gamma_0$ with generalized Riesz regression.
Generalized Riesz regression estimates the AME Riesz representer $\alpha^{\mathrm{AME}}_0$ directly, without explicitly modeling the density (or its score) $\partial_d\log f_{D\mid Z}(D\mid Z)$. Recall the population Bregman objective for a differentiable, strictly convex $g$: \[ \mathrm{BR}_g(\alpha) \coloneqq \mathbb{E}\Big[-g\big(\alpha(X)\big)+\partial g\big(\alpha(X)\big)\alpha(X)-m\big(\partial g(\alpha)\big)\Big], \] where, for AME, the linear functional is $m(\alpha)=\mathbb{E}\left[\partial_d \alpha(X)\right]$. Hence $m\big(\partial g(\alpha)\big)=\mathbb{E}\big[\partial_d\{\partial g(\alpha(X))\}\big]$.
\paragraph{Squared loss.} Let $g^{\text{LS}}(u)=(u-1)^2$ with $\partial g^{\text{LS}}(u)=2(u-1)$. Then \[ \mathrm{BR}_{g^{\text{LS}}}(\alpha) = \mathbb{E}\Big[ -\big(\alpha(X)-1\big)^2 +2\big(\alpha(X)-1\big)\alpha(X) -\partial_d\big\{2\big(\alpha(X)-1\big)\big\}\Big]. \] Here, we have the equivalent form (up to an additive constant independent of $\alpha$) \[ \mathrm{BR}_{g^{\text{LS}}}(\alpha) = \mathbb{E}\big[\alpha(X)^2-2\partial_d \alpha(X)\big] +\text{const}. \] Thus the squared-loss Bregman objective targets $\alpha^{\text{AME}}_0$ in $L_2$. This method corresponds to Riesz regression, shown in Section 2.2 in Chernozhukov2024automaticdebiased.
\paragraph{KL-type divergence.} For the signed KL-type convex function \[g^{\text{KL}}(\alpha)=|\alpha|\log|\alpha|-|\alpha|,\] the AME objective is given as
Minimizing the empirical version $\widehat{\mathrm{BR}}_g$ over a chosen class ${\mathcal{A}}$ gives us an estimator of $\alpha_0$.
Let $X\sim F_0$ denote a source covariate distribution that generates labeled data $(X,Y)$, and let $\widetilde X\sim G_0$ denote a target covariate distribution. Let $\{(X_i,Y_i)\}_{i\in {\mathcal{I}}_S}$ be i.i.d. from $F_0$ (source) and $\{\widetilde{X}_j\}_{j\in {\mathcal{I}}_T}$ i.i.d.\ from $G_0$ (target), independent.
Suppose we train $\gamma_0(x)=\mathbb{E}\left[Y\mid X=x\right]$ on a source population with covariates $X\sim F_0$, but the target parameter averages $\gamma_0$ over a shifted covariate distribution $\widetilde{X}\sim G_0$: \[ \theta^{\text{CS}}_0 \coloneqq \mathbb{E}\big[\gamma_0(\widetilde{X})\big] = \mathbb{E}\big[m^{\text{CS}}(\widetilde{X},\gamma_0)\big],\qquad m^{\text{CS}}(x,\gamma)\coloneqq \gamma(x). \] When $G_0$ is absolutely continuous with respect to $F_0$ with density ratio $r_0(x)\coloneqq \frac{dG_0}{dF_0}(x)$, the Riesz representer is \[ \alpha^{\text{CS}}_0(X)=r_0(X)=\frac{g_0(X)}{f_0(X)}. \]
The orthogonal score \[ \psi^{\text{CS}}(W;\eta,\theta)=\gamma(X)-\theta+\alpha^{\text{CS}}(X)\{Y-\gamma(X)\} \] delivers a debiased estimator by combining source residuals with target averaging, accommodating independent source/target samples via cross–fitting or data fusion.
We estimate $\alpha^{\text{CS}}_0$ directly by density–ratio matching under a Bregman divergence. Let $g:{\mathbb{R}}_+\to{\mathbb{R}}$ be differentiable and strictly convex. The population Bregman risk for a ratio model $\alpha$ is \[ \text{BR}_g(\alpha) \coloneqq {\mathbb{E}}_{G_0}\big[\partial g(\alpha(X))\alpha(X)-g(\alpha(X))\big] - {\mathbb{E}}_{F_0}\big[\partial g(\alpha(X))\big], \] which is equivalent to the Bregman divergence $\text{BR}'_g(r_0\mid r)$ up to a constant independent of $r$; the same decomposition and its empirical version are standard in the density-ratio estimation framework.
Given target samples $\{\widetilde{X}_j\}$ and source samples $\{X_i\}$, the empirical counterpart is \[ \widehat{\text{BR}}_g(\alpha) \coloneqq \frac{1}{|{\mathcal{I}}_T|}\sum_{j\in{\mathcal{I}}_T}\big\{\partial g(\alpha(\widetilde{X}_j))\alpha(\widetilde{X}_j)-g(r(\widetilde{X}_j))\big\} - \frac{1}{|{\mathcal{I}}_S|}\sum_{i\in{\mathcal{I}}_S}\partial g\big(\alpha(X_i)\big), \] and we set \[ \widehat{\alpha} \coloneqq \arg\min_{\alpha\in{\mathcal{A}}} \widehat{\text{BR}}_g(\alpha)+\lambda J(\alpha), \] with a model class ${\mathcal{R}}$ and regularizer $J$.
\paragraph{Squared loss.} For $g^{\text{LS}}(\alpha)=(\alpha-1)^2$, we have \[ \widehat{\text{BR}}_{g^{\text{LS}}}(\alpha) = \frac{1}{|{\mathcal{I}}_T|}\sum_{j} \alpha(\widetilde{X}_j)^2 - \frac{2}{|{\mathcal{I}}_S|}\sum_{i} \alpha(X_i), \] which yields the classical LSIF used in Kanamori2009aleastsquares. While Chernozhukov2025automaticdebiased proposes Riesz regression for covariate shift adaptation, their Riesz regression problem is the same as LSIF, and their covariate shift adaptation method basically coincides with Kanamori2009aleastsquares except for the use of regression functions. In parallel, Kato2024doubledebiasedcovariateshift also proposes a doubly robust form for covariate shift adaptation by extending Kanamori2009aleastsquares.
\paragraph{KL divergence.} For $g^{\text{KL}}(t)=t\log t - t$, we have \[ \widehat{\text{BR}}_{g^{\text{KL}}}(\alpha) = \frac{1}{|{\mathcal{I}}_T|}\sum_{j}\big\{\alpha(\widetilde{X}_j)\log \alpha(\widetilde{X}_j) - \alpha(\widetilde{X}_j)\big\} - \frac{1}{|{\mathcal{I}}_S|}\sum_{i}\log \alpha(X_i), \] which is the standard UKL risk used by KLIEP–style procedures.
\paragraph{Power–divergence family (robust Bregman).} For $g^{\text{PD}}_b(\alpha)=\big(\alpha^{1+b}-\alpha\big) / b$ with $b>0$, we obtain a continuum between Pearson–type ($b=1$) and KL–type ($b\to 0$) risks. Power–divergence choices trade efficiency and robustness: larger $b$ damp numerator outliers and can stabilize small–sample ratio learning, with convexity for $0 < b\le 1$.
Under specific choices of Riesz regression models and Bregman divergence, we can automatically guarantee the covariate balancing property. The key tool is the duality between Bregman divergence and covariate balancing methods.
Consider the linear model \[\alpha_\beta(X) = \Phi(X)^\top \beta,\] where $\Phi \colon \{1, 0\} \times {\mathcal{Z}} \to {\mathbb{R}}^p$ is a basis function. For this model, using the squared loss (Riesz regression) automatically attains covariate balancing, as discussed in BrunsSmith2025augmentedbalancing.
Specifically, under linear models, from the duality, this MSE minimization problem is equivalent to solving
where $\bm{0}_p$ is the $p$-dimensional zero vector. This optimization problem matches that used to obtain stable weights Zubizarreta2015stableweights.
It enforces the covariate balancing condition
where $\widehat{\alpha}_i = \Phi(X_i)^\top \widehat{\beta}$.
The advantage of the use of linear models is that we can write the whole ATE estimation with a single linear models, as shown by BrunsSmith2025augmentedbalancing.
We can model the Riesz representer via modeling the propensity score as \[\alpha_\beta(X) = \mathbbm{1}[D = 1]r_\beta((1, Z)) - \mathbbm{1}[D = 0]r_\beta((0, Z)),\] where
and $\Phi \colon {\mathcal{Z}} \to {\mathbb{R}}^p$ is a basis function. Note that we do not include $D$ unlike the basis function in linear models. For this model, if we use the KL-divergence–flavored convex function defined in Section (ref), we automatically attain covariate balancing, as discussed in Zhao2019covariatebalancing.
Define
and denote $r_{\widehat{\beta}}$ by $\widehat{r}$.
Specifically, under logistic models, from the duality, the KL divergence-flavored loss is equivalent to solving
This optimization problem matches that used in entropy balancing Hainmueller2012entropybalancing.
As a result, we obtain \[\sum^n_{i= 1}\Big(\mathbbm{1}[D_i = 1]\widehat{w}_i \Phi((1, Z_i)) - \mathbbm{1}[D_i = 0]\widehat{w}_i \Phi((0, Z_i))\Big) = \bm{0}_p,\] where $\widehat{w}_i =
.$
This model has the advantage that we can use a basis function $\Phi(Z)$ independent of $D$. Moreover, it naturally achieves covariate balance in the sense that the covariate distributions match between the treated and control groups. Additionally, we can automatically impose nonnegativity on $\alpha(1,Z)$ and $\alpha(0,Z)$, which can be violated in linear models.
This study proposed Direct Debiased Machine Learning (DDML), a unified framework that targets the oracle Neyman score and estimates the Riesz representer via Bregman divergence minimization. DDML integrates and generalizes several lines of work in semiparametric efficiency and modern causal inference. In particular, Neyman targeted estimation provides a systematic route to estimate nuisance parameters by directly minimizing the estimation error between the oracle score and its plug-in counterpart, while generalized Riesz regression offers a flexible estimator of the Riesz representer that encompasses squared-loss Riesz regression, tailored-loss covariate balancing, and direct density-ratio estimation as special cases.
Conceptually, the framework connects strands that are often treated separately. We show that Riesz regression, covariate balancing propensity scores and empirical balancing, and least squares and KL type density-ratio procedures all arise from the same Bregman divergence minimization once the relevant score and representer are specified. This perspective explains when ostensibly different estimators coincide and when they differ, and it guides principled choices of losses and model classes for a given target parameter. Methodologically, the Bregman divergence shows that, by changing the convex function, we obtain objective functions that include Riesz regression and KL divergence minimization, enabling objectives that enforce automatic covariate balancing.
In summary, DDML provides a single, coherent blueprint for constructing efficient, robust, and practically effective estimators for parameters defined through regression functionals. By casting nuisance estimation as Neyman targeted estimation and Riesz learning as Bregman divergence minimization, the framework unifies existing approaches and offers new tools that improve both theoretical guarantees and empirical performance.