{'fraction_non_alphanumeric': 0.04949732218359485, 'fraction_numerical': 0.013887061918631965, 'mean_word_length': 4.451341938127433, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimizing Variational Representations of Divergences and Accelerating their Statistical Estimation

Variational representations of divergences and distances between high-dimensional probability distributions offer significant theoretical insights and practical advantages in numerous research areas. Recently, they have gained popularity in machine learning as a tractable and scalable approach for training probabilistic models and for statistically differentiating between data distributions. Their advantages include: 1) They can be estimated from data as statistical averages. 2) Such representations can leverage the ability of neural networks to efficiently approximate optimal solutions in function spaces. However, a systematic and practical approach to improving the tightness of such variational formulas, and accordingly accelerate statistical learning and estimation from data, is currently lacking. Here we develop such a methodology for building new, tighter variational representations of divergences. Our approach relies on improved objective functionals constructed via an auxiliary optimization problem. Furthermore, the calculation of the functional Hessian of objective functionals unveils the local curvature differences around the common optimal variational solution; this quantifies and orders the tightness gains between different variational representations. Finally, numerical simulations utilizing neural network optimization demonstrate that tighter representations can result in significantly faster learning and more accurate estimation of divergences in both synthetic and real datasets (of more than 1000 dimensions), often accelerated by nearly an order of magnitude.


I. INTRODUCTION
Divergences and distances between multivariate probability distributions play a central role in many mathematical, engineering, and scientific fields ranging from statistical physics, large deviations theory, and uncertainty quantification to information theory, statistics, and machine learning. Variational representation formulas for divergences, also referred to as dual formulations, convert divergence calculation into an optimization problem over a function space and offer a valuable mathematical tool to build, train, and analyze probabilistic models and measure the similarity between data collections. Typical examples of variational representations are, among others, the Legendre transformation (LT) of an f -divergence [1], [2], the Donsker-Varadhan (DV) formula for the Kullback-Leibler (KL) divergence [3], [4] and the Rubinstein-Kantorovich duality formula for Wasserstein distance [5]. Variational representations have been used in statistical mechanics and interacting particles systems [6], large deviations [4], divergence estimation [7], [8], [9], determining variable independence through mutual information estimation [10], adversarial learning of generative models [11], [12], [13], uncertainty quantification (UQ) of stochastic processes [14], [15], bounding risk in probably approximately correct (PAC) learning [16], [17], [18], as well as in parameter estimation [19].
Two main mathematical ingredients are involved in the construction of a variational formula.
First, the function space where the optimal solution will be searched for and, second, the representation expression, called here the 'objective functional', whose optimization leads to the value of the divergence. Crucial practical advantages of variational formulas in statistics and machine learning include: a) they do not require an explicit form of the probability distributions (or their density ratio); related probabilistic quantities can be approximated by statistical estimators over available data; b) they can exploit the capacity of rich regression models such as neural networks to efficiently search the function space for optimal solutions; the optimal solution is typically related to the density ratio.
A single divergence can be derived from several different objective functionals. The key contribution of this paper is a systematic methodology that uses families of transformations (e.g., shifts, affine, and powers) to build new, tighter variational representations for divergences by creating improved objective functionals, as described in our main Theorem 1. This idea is both simple and powerful; it provides a general framework that unifies many of the previous variational formulas in the literature, reveals new connections between them, drives the derivation of new variational formulas, and has practical implications in terms of accelerated statistical training, learning, and estimation from data.
Striking consequences of the proposed framework include: (i) the connection between LTbased KL, the DV representation formula, and a new, improved DV-type formula, (ii) a concrete representation of the abstract objective functional in [9], and (iii) a derivation of new representation formulas for α-divergence and connections with a recently derived, DV-type variational representation of Rényi divergences. Moreover, when the optimization over the transformation family is not analytically tractable, a second-order approximation is employed resulting in new variational representations.
The improved objective functionals constructed via our framework have the same optimal solution, but they are tighter in the sense that the same approximation of the optimum will provide a better approximation of the divergence, i.e., they are flatter around the optimal solution. We employ (functional) Hessians of the objective functionals to quantify and order relative tightness gains between different variational representations of divergences, in terms of the local curvature around the optimal solution. Although our goal here is not variance reduction, we do also study the asymptotic variance of the tighter variational representations. We obtain theoretical evidence that our optimization procedures do not increase the variance and provide some numerical evidence that the variance as well as the averaged error can be decreased.
Finally, we demonstrate that these tighter representation formulas can accelerate numerical optimization and estimation of divergences in a series of synthetic and real examples, such as the statistical estimation of f -divergences and mutual information, including cases with real and/or high-dimensional data (in excess of 1000 dimensions). Similarly to [10], we parameterize the function space using neural networks, hence the (parameter) optimization is efficiently performed with back-propagation algorithms. As an example of our method, we develop neural-based estimators of controlled sample complexity for Hellinger-based mutual information (Hellinger-MINE). Overall, we find that the improved, tighter representation formulas converge several times faster than the initial representation formula, often by nearly an order of magnitude in high-dimensional problems.

Background
Define F 1 (a, b) to be the set of all convex functions f : (a, b) → R with f (1) = 0. If a (resp. b) is finite, we extend f to a (resp. b) by continuity and set f (x) = ∞ for x ∈ [a, b]. Such functions are appropriate for defining f -divergences, D f , which have the variational characterization where M b (Ω) denotes the set of all bounded measurable functions and f * is the Legendre transform of f [20], [8]. Under appropriate assumptions (see Theorem 4.4 in [20]), the maximum is achieved at There are cases, such as α-divergence with α ∈ (0, 1), where D f (Q P ) can be given a meaningful finite value even if Q ≪ P (see [21]). However, the right hand side of (1) is always +∞ if Q ≪ P and so here we use the convention that D f (Q P ) ≡ ∞ when Q ≪ P . The special case of Kullback-Leibler (KL) divergence (i.e., D f with f (x) = x log(x)) has a well-known alternative variational representation, the Donsker-Varadhan (DV) variational formula It is known that the objective functional in (3) is tighter than that of (1), in the sense that (Ω), [9]. In this paper, we present a general procedure for obtaining tighter variational representations of any f -divergence, for which the transition from (1) to (3) is just one special case.

Theoretical Results
Our method for deriving tighter variational representations for f -divergences is described in the following Theorem; a proof can be found in Section VI.
Theorem 1: Let f ∈ F 1 (a, b) and suppose D f (Q P ) < ∞. With φ * defined by Eq. (2) (when it exists), let Φ be a family of functions (the test functions) with Consider any family of transformations that includes the identity map. Then where and the maximum in (6) is achieved at φ * . Furthermore, the objective functional, H T ,f , in the variational representation (6) is tighter than the objective functional in (1), in the sense that for all φ ∈ Φ.

Remark 1:
When the choice of f and/or T is unimportant or clear from the context, we will omit the corresponding subscript on H T ,f .
(i) The main new insight and primary mathematical tool in this paper are formulas (6) and (7) alternatively and due to its low dimensionality, the corresponding optimization can be easily incorporated into the gradient descent framework without significant additional computational cost.

Families of Transformations
Next, we present several useful families of transformations that, in conjunction with Theorem This result was first obtained in [23], and in this sense Theorem 1 is a broad and systematic generalization of Theorem 4.4 in [23], where only shift transformations where considered.
c. Scaling Transformation: T scale η (φ) = ηφ, η ∈ R, which lead to the scaling or η-improved variational formula d. Affine Transformation: The above two cases can be combined into a two parameter family e. Power Transformation: T power β (φ) = φ β , β ∈ R, which lead to the power or β-improved variational formula As with the affine transformations, the power transformations can be combined with the shift and/or scaling transformations to form a multiparameter family. These are related to the wellknown Box-Cox transformation [24], used in statistics to transform non-normal data sets to approximately normal.

New Derivations of Existing Variational Representations
In this and next subsections, we explore several specific cases of the above framework, focusing primarily on examples where the optimization over the shift and/or scaling parameter can be done analytically. Here, we uncover connections with previously known variational formulas.
The maximum over the shift transformations T shif t in (6) occurs at ν * = log E P [e φ ] − 1, and hence The result is the objective functional in the well-known Donsker-Varadhan variational formula (3) and so this framework provides the connection between (3) and (1). We also note that in [25] a connection is derived between (3) and (1) by using a logarithmic change of variables in the function space, based on (2) for the KL case.
2) Connection with the results of [9]: In Theorem 1 of [9] the following improved variational formula was derived: This is another special case of our framework, as can be seen by first rewriting the Legendre-Fenchel transform and then using Theorem 4.2 in [23]: Hence the variational formula (13) is in fact the same as the ν-improved variational formula (9).
3) χ 2 -Divergence: We use Theorem 1 to provide a variational perspective on the classical Hammersley-Chapman-Robbins bound for the χ 2 -divergence. To our knowledge, the tightness result of the bound is novel. The χ 2 -divergence is a special case of the α-divergences, Here we optimize (6) over T af f ine to obtain where Φ = {φ ∈ L 1 (Q) ∩ L 2 (P ) : φ ≥ 0, Var P [φ] > 0} and the maximum is achieved at φ * = dQ/dP ≥ 0. Eq. (15) implies the Hammersley-Chapman-Robbins bound for the χ 2 divergence (see, e.g., Eq. 4.13 in [26]) and shows tightness over the set Φ. The objective functional in (15) was proposed as loss function for χ 2 -GANs, [27]; thus, (15) provides a complete and rigorous justification for this choice. Finally, if we instead optimize (6) over T shif t , we obtain the objective functional for χ 2 derived in [25], which is thus less tight than (15); see Appendix A.3.

Deriving New Variational Representations and Further Connections
In this subsection, we produce new variational representations and reveal further connections between divergences. We do not claim these examples are exhaustive. Nevertheless, they cover many important cases and illustrate the power and flexibility of Theorem 1.
The inclusions T shif t ⊂ T af f ine implies that (16) is tighter than DV (3). Calculations and numerical results that quantify this improved tightness are found in Section III. The optimization over η in Eq. (16) cannot be evaluated analytically in general, but it can be done numerically (as discussed in Section V).
2) α-Divergences (scaling transformations): The α-divergences are the family of f -divergences corresponding to f α (t) = t α −1 α(α−1) , α = 0, 1. See [21] for properties, related families, and further references. It includes the KL, Hellinger and χ 2 divergences as special cases [28], is closely related to the Tsallis entropies [29], and appears also in the context of information geometry [30]. By optimizing (6) over the family of scaling transformations, T scale (restricted to η > 0), we obtain a new variational representation of the α-divergences: if 0 < α < 1. Eq. (17) has the exact optimizers φ * α = (dQ/dP ) α−1 . Theorem 1 guarantees that the objective functionals in the new variational representations (17) are tighter than that of the LT f -divergence objective functional from (1). See Appendix A for details on the calculations that lead to Eq. (17), as well as for connections to the KL divergence in the limits as α → 0, 1.
Moreover, UQ bounds for risk-sensitive functionals in terms of Rényi divergences, which were obtained recently in [33], readily follow from (18) after appropriate manipulations and an optimization over α > 1.
The optimization over scalings was evaluated as in Eq. (17) but the optimization over the power transformations, β, cannot be done analytically. In practice, β can be included as an additional parameter in a numerical optimization procedure; see Section V for further discussion.

5) Exponential
Families: If P and Q are members of a parametric family then the set of test functions, Φ, in (6) can be reduced to a finite dimensional manifold. For instance, if P = P θp and Q = P θq are members of the same exponential family dP θ = h(x)e κ(θ)·T (x)−β(θ) , θ ∈ Θ with T : Ω → R n the vector of sufficient statistics then the explicit optimizer φ * = f ′ (dQ/dP ) lies on an (n+1)-dimensional manifold of functions, parameterized by the sufficient statistics and constants, φ (β,κ) = f ′ (exp(κ · T + β)), (β, κ) ∈ R × R n , and computation of the fdivergence reduces to the following finite-dimensional optimization problem: This variational representation can be further combined with any appropriate family of transformations T . We refer to Appendix E for further details.
where the optimal ∆η * is obtained by maximizing the second-order Taylor approximation and it is given by with dP φ = e φ dP/E P [e φ ] being the tilted measure. Since we quadratically approximate the improved Donsker-Varadhan, there is no guarantee that the objective functional in Eq. (21) is tighter than that of Donsker-Varadhan's representation. Despite no general guarantee, it is expected to be tighter when ∆η * (φ) is sufficiently small. This example demonstrates that even when there is no explicit formula for the transformation's optimization, one is still able to derive a closed-form formula for an approximate version of it and thereby obtain a new, rigorous variational formula. This same methodological approach can be used to analytically approximate the optimization over other families of transformations. As another demonstration, we refer to Theorem 7 in Appendix C where a new approximate variational formula is derived for the power-improved Rényi variational representation.

7)
Connections with Uncertainty Quantification: The improved DV representation (16) provides an alternative and arguably more general derivation of model uncertainty bounds derived recently in [14], [15]. These results quantify the effects of model uncertainty by bounding expectations of an observable φ under an alternative model, Q, in terms of the behavior under a baseline model P and the model discrepancy, measured by D KL (Q P ). Specifically, (16) implies after straightforward manipulation that 14]. More generally, we can obtain similar UQ bounds when model discrepancy is measured by an f -divergence by using T af f ine η,ν (φ) = ηφ − ν in (6) and performing the analogous manipulations: The Hammersley-Chapman-Robbins bound can also be viewed as a special case of (23) in this UQ context.

III. VARIATIONAL DERIVATIVES AND TIGHTNESS GAINS
In Theorem 1, we established the general methodology for building tighter variational representations of f -divergences, by constructing suitable objective functionals H T . Here we will quantify relative tightness gains corresponding to different transformation families T : for all such families the maximizer in (6) is always φ * given by (2). Therefore, our approach relies on building quadratic variational approximations of each objective functional (7) around the common maximizer φ * , and subsequently comparing the corresponding (variational) curvatures; see Figure   1 for a demonstration. Specifically, using that the maximum occurs at φ * , an asymptotic expansion yields where we formally define and ψ is any functional perturbation of the maximizer φ * . Formally, the second order term ∇ 2 H T [φ * ], i.e., a variational Hessian, is necessarily non-positive and determines the behavior in a neighborhood of the maximizer. By for different families T , we will quantify the 'tightness gains' provided by different transformation families. All these calculations can be made rigorous under appropriate assumptions as demonstrated next in Theorem 2 and in Theorem 8 in Appendix D.
Here we focus our analysis on affine transformations, T af f ine η,ν (φ) = ηφ − ν, but a similar analysis can be performed for any family T with a smooth, finite-dimensional parameterization.
The LT f -divergence variational formula (1) corresponds to T containing only the identity, and we write the corresponding objective functional as H id [φ]. Specializing (7) to the affine case, we define the functional which leads to four different objective functionals and variational representations of the f - and the corresponding Hessians  Here, φǫ = φ * +ǫψ where ψ = x 2 (top panel) and ψ(x) = x (bottom panel). P and Q are 1-dimensional Gaussians. Black curves: LT-based f -divergence objective functional, Blue curves: shift-improved (i.e., Donkser-Varadhan), Magenta curves: η-improved, Red curves: shift-scaling-improved. Note that ψ = x 2 is related to φ * by a shift and scaling, hence the shift-scaling-improved objective functional in the top plot has zero curvature in this direction, a manifestation of its shift and scale invariance.
Motivated by the formal calculation (24), we define The second derivative, when it exists, d 2 in the direction ψ and describes the curvature of the objective functional H T [φ] along any given direction ψ ∈ M b (Ω) at the maximizer φ * . Therefore we will define a corresponding variational Hessian in the direction ψ as We remark that relation (28) is simply a notation that only intends to draw a parallel to finite dimensional optimization calculations using conventional Hessians at maxima or minima of finite dimensional functions. All rigorous results are stated in terms of the 2nd Gateaux derivative x log(x). In addition to the assumptions of Theorem 1, we assume that the maximizer φ * = log(dQ/dP ) + 1 in (2) is bounded (i.e., φ * ∈ M b (Ω)) and select the function space Φ = (27) is twice differentiable at ǫ = 0. Furthermore, using the notation (28), the corresponding 2nd Gateaux T af f ine satisfy the following: corresponding to (1) Eq. (31) as computed by these 2nd Gateaux derivatives; note that they are progressively smaller in magnitude. These curvature computations demonstrate how one can rigorously and precisely quantify heuristics such as those presented in Figure 1 from [9]. Furthermore, our Hessian computations in the form of Theorem 8 also quantify and extend to f -divergences the accuracy gains observed in the neural estimation of mutual information in [10].

IV. f -DIVERGENCE ESTIMATOR BIAS AND VARIANCE
The variational formula Eq. (6) suggests the following natural f -divergence estimator where Q n and P n are the n-sample empirical measures and Φ ⊂ Φ is a function space that can be optimized over numerically (e.g., a family of neural networks). A natural question is therefore the bias and variance of this estimator. In practice, the optimizations are performed via some stochastic gradient descent (SGD) algorithm, and one is actually interested in the bias and variance after a finite number of training steps. Addressing this complicated problem, which depends heavily on the choice of space Φ and the SGD algorithm, is outside the scope of the present work. However, in this section we will follow the prior work in [35] for KL divergences and discuss the bias and variance of the objective functional in (32). Finally, we emphasize that our goal of this paper is to develop tighter objective functionals and study the impact of improved curvature on the speed of convergence of numerical estimators. This is a separate question from that of variance reduction, which we do not pursue here. However, we will show that, for α-divergences, optimizing over scalings does not worsen the variance at the optimizer.

Objective Functional Bias
The estimator (32) can be viewed either as a single-stage optimization problem over both φ and T , in which case the objective functional is unbiased, or as a two-stage optimization, with a biased objective functional. If the optimization over T must be performed numerically then we take the former view, but if part of the optimization can be performed analytically (such as in (12) and (17)) then we take the latter view. In any case, the divergence estimator as a whole, Eq. (32), is not unbiased. This can be seen in the simple case of discrete measures on a finite sample space, where one can numerically optimize over all of M b (Ω), which is biased. This fact renders the biased/unbiased objective functional question less relevant in practice, as one's goal is generally to estimate the divergence and not just the objective functional at a fixed φ.

Objective Functional Variance
For general f , the objective functional estimator for the LT variational formula is The two terms are independent, therefore In particular, for α divergences (36) reduces to where for α ∈ (0, 1) we made a change of variables φ → −φ.
To compute the asymptotic variance of an optimized objective functional, we rely on the delta method. This method can be applied to any objective functional that can be expressed as a function of expectations. We provide the details in the case of scaling-optimized α-divergences Eq. (17).
Then we have: (a) The scaling-optimized α-divergence objective functional (17), has asymptotic variance See Appendix F for the proof.
Theorem 3 implies that we need (dQ/dP ) 2α ∈ L 1 (P ) to achieve a finite variance; this fact can motivate an appropriate choice of α; we also refer to the Hellinger-MINE discussion and Corollary 1 in Section V. The asymptotic variance (40) agrees with the variance of the LT objective functional, Eq. (37), at its optimizer φ * α = (dQ/dP ) α−1 /|α − 1|. Hence, from the perspective of the variance at the optimizer, neither method has an advantage. Away from the optimizer there is no consistent relationship between the two variances. In practice, both methods will take different paths to the optimizer and this further complicates any variance comparison away from the optimizer. However, empirically we found that the estimators constructed via the improved variational formulas have smaller variance at the estimated optimal which is different in general from the theoretical optimizer φ * (i.e., after a finite number of SGD steps); see Figure   3.

V. NUMERICAL EXAMPLES: FASTER STATISTICAL ESTIMATION AND LEARNING
Next we discuss practical implications of using tighter variational representations developed in Theorem 1, focusing on accelerating neural-based statistical learning and estimation. In recent works, variational representations such as (1) or (3) were used to estimate f -divergences and likelihood ratios based solely on available data [8]. This variational perspective proved also to be a crucial mathematical step in training generative adversarial networks (GAN) [11], [12], [13] and towards developing neural-based estimators for mutual information, [10], taking advantage of the ability of neural networks to search efficiently through function spaces.
Improved variational formulas for statistical estimation and learning were previously studied in: i) [9], using Eq. (13) and assuming a Hilbert space (RKHS) function space, ii) [10], where the DV and LT formulas for the KL divergence were used to estimate mutual information with improved accuracy. Both of these implicitly rely on the shift-improved variational formula (see Eq. (12) and Eq. (14)). Our Theorem 1 provides a broad generalization of these ideas to other transformation families, allows for practical implementation of the method in [9] to more general functions space parametrizations (e.g., neural networks), and generalizes the ideas in [10] to other f -divergences beyond KL, where it can provide improved mutual information estimators based on (6).
In the following, we employ the outcomes from Sections II & III and build several variational neural network estimators, in the general spirit of [8], [10]. We demonstrate the performance improvements that result from representations such as Eq. (6). We start with the heuristic observation, illustrated in Figure 1, that tighter representations can improve the accuracy of statistical estimators for f -divergences, in the sense that the same approximation of the optimal φ * will provide a better approximation of the divergence. Moreover, tighter variational formulas can lead to faster convergence of the search algorithm, as we now motivate: suppose one minimizes a convex function f (x) by the simple gradient descent algorithm x n+1 = x n − γ∇f (x n ). If ∇f is L-Lipschitz (i.e., the Hessian is bounded by L) then this algorithm converges if 0 < γ < 2/L, and the analysis suggests the optimal learning rate of γ = 1/L and leads to the error bound [36]). If f has the same optimizer and optimal value, but has a smaller Hessian bound, L, then the optimal learning rate, γ = 1/ L is larger, and the error bound after an equal number of steps is smaller, i.e., the use of f in place of f can lead to faster convergence.
The above argument is only heuristic; the constant learning rate algorithm is far from optimal in most cases and the above analysis does not capture the complexity of the current setting.
Nonetheless, it does provide important insight into the numerical results, presented below, which demonstrate that, in practice, the improved variational formulas do generally lead to faster convergence of the estimators, letting all other factors be equal.
The examples below use the new variational formulas (17) and (19) for α-divergence, largely focusing on the well-known Hellinger divergence (α = 1/2), defined as where dQ = qdµ, dP = pdµ. Computations were done in TensorFlow using the AdamOptimizer [37], an adaptive learning-rate SGD optimizer, with all methods given the same initial learning rate. When working with neural-network based estimators of α-divergence, we enforce positivity of the test functions (see (17) and (19)) via the parameterization φ θ = exp(g θ ) where g θ , θ ∈ Θ, is a neural network family with ReLU activation functions. We compare the LT method, (1) with the improved estimators based on our formulas (17) and (19). For instance, in the case of Hellinger divergences we compare the following estimators.
LT Hellinger Estimator: Scaling-Improved Hellinger Estimator: Scaling-Power-Improved Hellinger Estimator: In the above, Q N and P N denote the expectation under the empirical distributions using N iid samples from Q and P respectively.
If the optimization over a parameterized family of transformations, T δ , cannot be performed analytically then we solve the minimization problem (6) -(7) by performing stochastic gradient descent (SGD) on the full collection of parameters, (θ, δ). In such cases, our two-step formulation can be thought of as parametric enhancement of the neural network architecture. The nested nature of the minimization over φ and T also allows for more sophisticated methods (not explored here), e.g., for each φ one can perform several SGD steps for T , thus solving the (generally low dimensional) problem (7) to high accuracy, before performing another SGD step for φ in (6); this is reminiscent of multiscale numerical methods [38]. The parameterization, δ, of the families of transformations considered here is at most two dimensional. For example in the case of T af f ine we have δ = (ν, η), where ν, η ∈ R; see also the remaining examples in Section II. Including this small number of additional parameters in the stochastic gradient descent iterations is expected to add a negligible additional computational cost, as compared to the (generally) much larger number of neural-network parameters, θ. In practice, we do find the additional computational cost to be negligible.

Hellinger-MINE
Here we consider the problem of computing Hellinger mutual information (Hellinger-MI), is chosen to be KL. However, one can consider a whole array of different f -divergences for this purpose, see for instance [39]. A motivation for choosing Hellinger over KL is rigorously based on the variance calculations in Section IV. In particular, as a direct consequence of Theorem 3 we obtain the following: Corollary 1: Under the assumptions of Theorem 3, the relative variance for the Hellinger Therefore, the sample complexity of the estimator H scale,f 1/2 [φ * ; Q n , P n ] at the optimizer φ * is Comparing Corollary 1 with the corresponding result for the KL divergence from Theorem 2 of Ref. [35], we see that in practice the KL divergence requires n = O(e D KL (Q,P ) ) samples, while the Hellinger divergence requires n = O(1) samples due to (46).
In Figure 2, we present the computation of Hellinger mutual information (Hellinger-MI), formulas (red and blue) yield more accurate results, i.e., they converge faster than the LT fdivergence method (1) (black). Moreover, optimizing over both scalings and powers (19) (red) provides a non-trivial improvement over the scaling-improved method (17) (blue). This is a generalization of the findings in [10], which compared the DV variational formula (3) with (1) for the KL divergence. We emphasize that despite the lack of an analytical formula for the optimization over β, the inclusion of this single additional parameter in the variational formula (a negligible addition to the computational cost) leads to a clear performance gain.
As a followup, we demonstrate the effectiveness of our method in estimating Hellinger-MI for high-dimensional problems with low-dimensional structure. Specifically, in Figure 3 we compare 20-dimensional Gaussians embedded in high dimensional space via a nonlinear map. The left panels demonstrate the performance gain when using the curvature-improved objective functionals; we find that the optimized methods significantly outperform the LT method, especially in higher dimensions. The right panels are a (partial) demonstration that the variances are comparable, if not improved, when using the optimized objective functionals. Specifically, we find that the optimized methods approach the exact asymptotic variance faster while the LT method, which has not yet converged, has a larger variance.

Submanifold Parameterization for Exponential Families
Our method allows for a great deal of flexibility in the choice of function space parameterization. In a 'small-data' setting, the assumption of an exponential family structure can serve as an effective regularization. We illustrate this with Figure 4, which shows the estimation of the αdivergence with α = 0.25 between 10-dimensional Gaussians using a data set of 5000 samples from each distribution for SGD (minibatch size of 100) and using another 5000 samples for Monte Carlo estimation of the value of the objective functional. Using the submanifold estimation formula (20) and its scaling-improved variant (see Appendix E) we obtain the magenta and red curves, respectively. In blue, we show the result from the scaling-improved variational formula (17) and in black we show the result using the LT f -divergence objective functional (1); both use  Fig. 4. Estimation of the α-divergence with α = 0.25 between two 10-dimensional Gaussians with randomly generated variances and one of the means randomly perturbed from zero. We compare the convergence performance between the neural network and submanifold parameterizations. The relative error was averaged over 50 runs.
neural network families with one fully connected hidden layer (5 nodes). The number of nodes was chosen so that all methods use approximately the same number of parameters. The neural network parameterization converges faster, but ends up with a larger bias than the submanifold parameterization. The scaling-improved variational formulas lead to faster convergence than the LT variational formula in both cases as expected by our theory.

MNIST Dataset Examples
Next we illustrate the accelerated speed of convergence on high-dimensional (28 × 28 = 784 dimensional) realistic data by estimating the Hellinger divergence between two distributions obtained by (iid) randomly translating the MNIST handwritten digits image dataset [40]. This provides an effective test case wherein we know the exact answer (D f 1/2 (Q P ) = 0). Figure 5 shows the error, as a function of the number of SGD iterations, and once again demonstrates that the improved variational formulas lead to faster convergence; in this case, nearly one order of magnitude fewer SGD iterations are required to reach an accuracy of 10 −2 when using the shifts, rounded to the nearest pixel and with periodic boundary conditions). Each step of SGD uses two independent minibatches of 100 such samples (one minibatch for P and one for Q). Monte Carlo estimation of the value of the Hellinger divergence was done using the corresponding objective functional and with samples coming from two separate datasets of 10000 randomly shifted MNIST images (one collection of images for P and an independent collection for Q). The function space was parameterized via fully-connected neural networks with one hidden layer of 128 nodes. The results were averaged over 50 runs. tighter objective functionals. In practice, this means that one can more quickly detect whether or not the two data streams are in fact coming from the same distribution.
To further illustrate that our estimators are behaving appropriately, we perform a pair of consistency checks using the MNIST dataset, similar to the tests in [35]. While not a perfect substitute for computing the relative error, tests such as these are very useful in situations where the exact value of the divergence is nonzero and unknown. In Figure 6(a) we test the data processing inequality for f -divergences: where κ is a probability kernel. Here we let Q be the MNIST dataset, P be the MNIST dataset of digits 0 through N P where N P ranges from 0 to 8, and we let κ x be the distribution of random translations of the image x (specifically, N(0, 1) translations, with components rounded to the nearest integer). The plot shows the ratio of the estimators for D f 1/2 (Q ⊗ κ P ⊗ κ) and D f 1/2 (Q P ), using various objective functionals. In Figure 6(b) we test the product measure identity for α-divergences: Here we let Q i be copies of the MNIST dataset and P i be copies of the MNIST dataset of digits 0 through N P , where N P again ranges from 0 to 8. The plot shows the ratio of the righthand-side of Eq. (49) to the left-hand-side. The black horizontal lines in both panels show the exact value of 1. We find that all methods perform well on these consistency checks, though the best performing method is different for the two tests. This is unsurprising, as our methods are designed to accelerate the convergence to the optimum value, which is a different goal than preserving these the above two properties. Note that as N P increases the distributions Q and P become more similar and so both the numerator and denominator approach zero, making the task of estimating the ratio more difficult; this is reflected in Figure 6. To conclude these examples, we note that although the proposed optimization framework was applied on statistical learning and estimation, it can be of broader interest, among others, in epistemic uncertainty quantification [14], in coarse-graining and model reduction [41], [42], [43], as well as in PAC learning [18] and adversarial learning [11], [13]. In particular, we intend to explore the use of variational formulas derived via the quadratic approximation method (e.g., Eq. (21) and Eq. (100)) for uncertainty quantification, along with the UQ bound (23).

VI. PROOF OF THEOREM 1
In this section we provide a detailed proof of Theorem 1. For the convenience of the reader, we will recall the relevant definitions and notation below.
Let P, Q be probability measures on a measurable space (Ω, M) and, for any −∞ ≤ a < with respect to P is defined by Our starting point is the the following variational characterization [20], [8]: where M b (Ω) denotes the set of all bounded and measurable functions and is the Legendre transform of f . As a first step, we show that Eq. (51) can be extended to certain unbounded φ. This is similar to results in [20], [8] but we will prove explicit conditions for which the expectations exist. To do this, we will need the following lemmas: a, b). Then one of the following holds 1) f * is bounded below.
2 To show f * is non-decreasing, suppose that we have x 1 < x 2 with f * (x 1 ) > f * (x 2 ). Taking y n as above, find an n such that y n < x 1 and f * (y n ) < f * (x 2 ). f * is convex and so, letting This is a contradiction, hence f * is non-decreasing. Then Remark 5: We use the notation g + ≡ g1 g≥0 and g − ≡ −g1 g≤0 for positive and negative parts of a function g = g + − g − .
Proof: Fix φ ∈ L 1 (Q). If f * is bounded below then (f * ) − is bounded above and the result is trivial, so suppose not. Lemma 1 then implies that f * is non-decreasing. General properties of Legendre transforms on the real line imply that f * is convex, lower semicontinuous, and f * Hence is bounded above and so φ b,n ∈ M b (Ω) and we can use Eq. (51) to find We have φ b,n → φ b pointwise, |φ b,n | ≤ |φ b |, and φ b ∈ L 1 (Q), so we can use the dominated convergence theorem to obtain (here it was important that we are in the case where D f (Q P ) < ∞). We also have φ b,n+1 ≤ φ b,n , Therefore the monotone convergence theorem implies We therefore conclude that We can now prove that Eq. (51) can be extended to φ ∈ L 1 (Q).
Theorem 4: Let f ∈ F 1 (a, b) and suppose either f * is bounded below or D f (Q P ) < ∞. Then where the objective functional is valued in [−∞, ∞).
then the claimed result will follow by using Eq. (51).
suppose not. Then f * (φ) < ∞ P -a.s. We are in the case where D f (Q P ) < ∞, and so Q ≪ P and f * (φ) < ∞ Q-a.s. as well.
We now prove Theorem 1 from the main text, which we restate below. For completeness, we also provide a derivation of the formula for the optimizer, φ * , which was obtained in [20].
Theorem 5: Let f ∈ F 1 (a, b) and suppose either f * is bounded below or D f (Q P ) < ∞. Then: 1) Suppose Q ≪ P , f is C 1 , f ′ is strictly increasing, and one of the following holds: exists and is finite.
latter case, we also assume that the conditions from part (1) hold, so that φ * is the optimizer).

Consider any family of transformations T ⊂ {T
includes the identity map. Then 3) If Q ≪ P and (Ω, M) is a metric space with the Borel σ-algebra then one can replace M b (Ω) with C b (Ω) (bounded continuous functions) and L 1 (Q) with L 1 c (Q) (continuous L 1 (Q) functions) in the above.

Remark 6:
Recall that −∞ ≤ a < 1 < b ≤ ∞. In particular, one could have b = ∞ and so the assumptions of part (1) do not necessary imply dQ/dP is bounded.
Proof: Suppose that f satisfies the additional conditions from item (1). Then the Legendre transform can be computed: for all y ∈ range(f ′ ). This implies that for x ∈ (a, b) we have f (x) = xf ′ (x) − f * (f ′ (x)). By taking limits and using the assumptions on dQ/dP , we find f (dQ/dP ) = dQ/dP f ′ (dQ/dP ) − f * (f ′ (dQ/dP )). Therefore, assuming φ * ≡ f ′ (dQ/dP ) ∈ L 1 (Q), we have which proves the claim. Now, let Φ and T be as in item (2). Since every T ∈ T maps Φ into L 1 (Q) and T contains the identity we can use Theorem 4 to obtain for all φ ∈ Φ. Maximizing over Φ and using the fact that either along with Eq. (59) gives which proves the claim.
Finally, on a metric space, one can approximate measurable functions with continuous functions via Lusin's theorem (see, e.g., Appendix D in [22]). Using this fact it is straightforward to show that the same results are obtained if one replaces M b (Ω) with C b (Ω) and L 1 (Q) with L 1 c (Q) (see the Proof of Theorem 1 in [31] for details on the use of this technique in a similar context). This proves item (3).

APPENDIX A DERIVATION OF VARIATIONAL FORMULAS FOR α AND RÉNYI DIVERGENCES
Here we provide additional details regarding the derivation of the new variational formulas for α-divergences, as well as their connection to the Rényi divergences.
For a, c ≥ 0, b > 0 we have the following general solution to the optimization problem where the maximum occurs at η * = (a(α − 1)/(bα)) α−1 . Hence we can use this formula to .
Using this, for Q ≪ P with (dQ/dP ) α ∈ L 1 (P ) we find and the maximum is achieved at In Eq. (74) one should interpret 0/0 ≡ 1 and c/∞ ≡ 0, or else restrict to the subset of functions with 0 < E P [φ α/(α−1) ] < ∞ for which these cases do not occur. Also note that the objective functional (74) is now invariant under scaling and so we were able to drop the factor of 1/(α−1) in Eq. (75).
If α ∈ (0, 1) and y < 0 then and, if Q ≪ P , a similar computation to the above yields If dQ/dP > 0 then the exact optimizer is φ * α = (dQ/dP ) α−1 . Note that we reparameterized φ → −φ in Eq. (77) so that we optimize over strictly positive functions, rather than strictly negative functions. Also, note that it is not necessary to restrict the optimization in (77) to φ ∈ L 1 (Q); if E Q [φ] = ∞ then the objective functional equals −∞ and hence including such φ's does not change the value of the supremum. Theorem 1 guarantees that the objective functionals in the variational representations Eq. (74), and Eq. (77) are tighter than that of the LT f -divergence objective functional from (51).

2.
In particular, for the Hellinger distance (α = 1/2) we find the scaling-improved formula as compared to the LT f -divergence variational formula (after changing variables φ → −φ) The improved tightness, 2 ), is guaranteed by the general result of Theorem 5, but it can also be seen to be a consequence of the inequality 3. For the χ 2 -divergence (equal to 2D f 2 ) one can evaluate the optimization over all affine tranformations. We have f 2 (x) = (x 2 − 1)/2 and f * 2 (y) = (y 2 + 1)/2. Suppose D f 2 (Q P ) < ∞. Optimizing first over shifts we find, for φ ∈ L 1 (Q) ∩ L 2 (P ), with the maximum occurring at ν * = E P [φ]−1. Therefore we obtain the variational representation Note that the objective functional in the last line is the same as that obtained in Eq. (48) of [25].
Further optimizing the objective functional over the scaling parameter η ∈ R we find then the supremum equals zero, since Q ≪ P ). Therefore Equality is achieved at φ * = dQ/dP and so one can further restrict the optimization to φ ≥ 0.
This provides a rigorous justification of the loss function for χ 2 -GANs proposed in, [27]. We emphasize that the objective functional in (83) is tighter than the one in (81).
and so which is the the DV objective functional for D KL (Q P ).
Similarly, reparameterizing via φ α = e −(α−1)g gives and so which is the objective functional for the Donsker-Varadhan representation of D KL (P Q).

5.
Using the same reparametrization as in Eq. (84) we can also derive a connection with the Rényi family of divergences: Fix α ∈ (0, 1), reparametrize φ = e (α−1)g (g is any measurable function), and rewrite Eq. (77) in terms of an optimization over g: Using the connection between α-divergences and Rényi divergences (see [21], but note that our definition of Rényi divergence differs from theirs by a factor of 1/α), we obtain .
A similar calculation applies when α > 1; in either case, we obtain an independent derivation of the Rényi divergence variational formula from [31].

APPENDIX B IMPROVED DONSKER-VARADHAN VARIATIONAL FORMULA
Here we collect some additional properties of the improved Donsker-Varadhan (DV) variational formula (16): One obviously has T id ⊂ T shif t ⊂ T af f ine . As a consequence, the objective functional in (91) is tighter than DV, which in turn is tighter than the LT f -Divergence objective functional from Eq. (51): Although the supremum over η in Eq. (91) cannot in general be evaluated analytically, one can obtain an explicit approximation as follows: and Taylor expand around η = 1 to obtain We approximate the optimal ∆η by maximizing the quadratic approximation (92) to find where dP φ = e φ dP/E P [e φ ] is the tilted measure and ∆η * (φ) is defined to be 0 if Var P φ [φ] = 0.
Eq. (21) is a new variational representation of the KL-divergence, however we make no claim that the objective functional in (21) is tighter than Donsker-Varadhan for every φ; Eq. (21) is only guaranteed to be tighter than DV when ∆η * (φ) is sufficiently small. This is because we only maximized the quadratic approximation in η, and hence only obtained an approximation to the scaling-improved objective functional from Eq. (91). One can of course circumvent this by using the maximum of the two: Although Eq. (96) is less than aesthetically appealing, note that its objective functional is certainly no worse than DV and will be tighter when ∆η * (φ) is sufficiently small. Also note that no additional expectations need to be computed to evaluate the objective functional in Eq. (96) as compared to that of Eq. (21).

APPENDIX C APPROXIMATING THE POWER-IMPROVED RÉNYI VARIATIONAL REPRESENTATION
The same method used to derive the second-order approximation of the improved Donsker-Varadhan variational formula (21) can be applied to the Rényi divergence with power transfor-mations. Recall that the power-improved Rényi variational is given by to be the objective functional as a function of β. A Taylor expansion around β = 1 gives with dQ α = e αg dQ/E Q [e αg ] being the α-tilted measure and similarly for P . The quadraticallyoptimal ∆β * is obtained by maximizing the second-order Taylor approximation and is given by where we assumed 0 < α < 1.
A new Rényi representation formula is then obtained as the following theorem asserts.

APPENDIX D TIGHTNESS GAINS AND VARIATIONAL DERIVATIVES
Here we prove an extension of Theorem 2 for the more general case of f -divergences.

Theorem 8 (Tightness gains for f -divergences):
In addition to the assumptions of Theorem 1 suppose f ′ is strictly increasing and f * ∈ C 3 ((c, d)) with (f * ) ′′ > 0. We select the function space Φ = M b (Ω) and assume that the maximizer φ * in (2) is valued in a compact subset of (c, d). Then, for T = T id , T shif t , T scale , T af f ine we have that J T in (27) is twice differentiable at ǫ = 0 for any ψ ∈ M b (Ω). Furthermore, using the notation (28), the corresponding 2nd T af f ine satisfy the following: corresponding to (1), (9), (10) and the combination of the last two respectively. The tilted probability measure P * is defined as while ρ P * (φ * , ψ) denotes the correlation between φ * and ψ.
The following general computation will facilitate the computation of the variational Hessians that appear in Eq. (111): The maxima in Eq. (110) are achieved at η = 1, ν = 0, φ = φ * ≡ f ′ (dQ/dP ). Given C 2 families φ ǫ = φ * + ǫψ, η ǫ with η 0 = 1, and ν ǫ with ν 0 = 0, we can Here we used the boundedness assumptions on f * , the differentiablity of f * and the dominated convergence theorem [34] (Theorem 2.27); the latter allowed us to rigorously exchange derivatives and expectations when we calculated the 1st and the 2nd Gateaux derivatives of the objective functionals above. The same argument is also used in the evaluation of all Gateaux derivatives below. Furthermore, for any choice of ψ, ν ǫ , η ǫ that have the above properties, the maximum of ǫ → H[φ ǫ , η ǫ , ν ǫ ] is achieved at ǫ = 0. Therefore the first derivative vanishes at ǫ = 0. In particular, by considering the case η ǫ ≡ 1, ν ǫ ≡ 0 we see that for all ψ. The second derivative at the maximum is given by where we used Eq. (113) to cancel several terms. We will specialize this to compute all four Hessians from Eq. (111).

A. LT f -Divergence objective functional
Here we fix η ǫ ≡ 1, ν ǫ ≡ 0: and hence We have assumed (f * ) ′′ > 0 and so we can define the probability measure and thereby write (compare with the KL-case Eq. (29)).

B. Optimization over shifts
Here η ǫ ≡ 1 and ν ǫ is chosen so that By using the convexity of the objective functional in ν along with implicit function theorem, one can see that the assumptions on f are sufficient to ensure that such a smooth ν ǫ exists for ǫ in a neighborhood of 0. We can simplify The derivative of ν ǫ can be computed as follows: Eq. (119) implies for all ǫ. Differentiating with respect to ǫ gives hence where P * is the probability measure defined in Eq. (117). Therefore (compare with the KL-case Eq. (30)).
Recalling the result (118) we arrive at This shows that the magnitude of the second derivative (which is guaranteed to be non-positive) has been reduced by the amount

C. Optimization over scaling transformations
Here ν ǫ ≡ 0 and η ǫ is chosen so that Again, the existence and smoothness of η ǫ is guaranteed by convexity and the implicit function theorem. Next note that for all ǫ. Taking the derivative with respect to ǫ gives Solving for η ′ 0 and using Eq. (113) to simplify yields Therefore, the Hessian is Once again, the magnitude of the second derivative has been reduced. 41

D. Optimization over affine transformations
Finally, let η ǫ and ν ǫ be defined by The existence and smoothness of η ǫ and ν ǫ are again guaranteed by convexity and the implicit function theorem. They satisfy for all ǫ. Simplifying this and differentiating with respect to ǫ we obtain the two equations If we define then these have the solution  (note that a 2 < bc follows from the Cauchy Schwarz inequality together with the assumptions that f ′ is strictly increasing and Q = P ). We can compute We have b + c ≥ 0 and bc − a 2 > 0 hence the matrix in Eq. (136) is positive semi-definite.
Therefore, computing the f -divergence reduces to the following finite-dimensional optimization problem: In some cases one can further reduce the dimension by optimizing over an appropriate family of transformations.
1) For KL divergence, we have φ (β,κ) = log(exp(κ · T + β)) + 1 = κ · T + (β + 1) , and so the optimizer lies on a linear manifold. If one optimizes over the family of shifts (i.e., one uses the Donsker-Varadhan formula) then the β + 1 factor is eliminated and one finds that for any Q, P that are members of the exponential family, i.e., the optimization is over the n-dimensional subspace spanned by the sufficient statistics: φ κ = κ · T, κ ∈ R n . Also, note that κ → E Q [κ · T ] − log E P [e κ·T ] is concave.

FUNCTIONAL
In this section we provide details regarding the computation of asymptotic variance of the objective functional estimators for α-divergences (see Theorem 3), generalizing some of the work on KL divergences from [35]. The main tool is the following lemma, which is based on the delta method.
Convergence of the mean and variance will follow from (148) if we can show uniform integrability of the random variables | √ n(h(E n [X]) − h(µ))| 2 ; see, e.g., Theorem 5.9 in [45]. If sup n E[n 2 E n [X] − µ 4 ] < ∞ then the limit of (151) as c → ∞ equals zero and we are done. Expanding this expression and using the assumption that the X i 's are iid we find, after a somewhat long but straightforward calculation, which is bounded above uniformly in n. This completes the proof.
We now prove the asymptotic variance result for α-divergences, Theorem 3.
Simplifying this gives the claimed result (39) .