{'fraction_non_alphanumeric': 0.04979991107158737, 'fraction_numerical': 0.01477424309794252, 'mean_word_length': 4.289043292357029, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Information Newton's flow: second-order optimization method in probability space

We introduce a framework for Newton's flows in probability space with information metrics, named information Newton's flows. Here two information metrics are considered, including both the Fisher-Rao metric and the Wasserstein-2 metric. Several examples of information Newton's flows for learning objective/loss functions are provided, such as Kullback-Leibler (KL) divergence, Maximum mean discrepancy (MMD), and cross entropy. The asymptotic convergence results of proposed Newton's methods are provided. A known fact is that overdamped Langevin dynamics correspond to Wasserstein gradient flows of KL divergence. Extending this fact to Wasserstein Newton's flows of KL divergence, we derive Newton's Langevin dynamics. We provide examples of Newton's Langevin dynamics in both one-dimensional space and Gaussian families. For the numerical implementation, we design sampling efficient variational methods to approximate Wasserstein Newton's directions. Several numerical examples in Gaussian families and Bayesian logistic regression are shown to demonstrate the effectiveness of the proposed method.


Introduction
Optimization problems in probability space are of great interest in inverse problems, information science, physics, and scientific computing, with applications in machine learning (Amari, 2016;Stuart, 2010;Liu, 2017;Amari, 1998;Villani, 2003). One typical problem here comes from Bayesian inference, which provides an optimal probability formulation for learning models from observed data. Given a prior distribution, the problem is to generate samples from a (target) posterior distribution (Stuart, 2010). From an optimization perspective, such a problem often refers to minimizing an objective function, such as the Kullback-Leibler (KL) divergence, in the probability space. The update relates to finding a sampling representation for the evolution of the probability.
In practice, one often needs to transfer probability optimization problems into samplingbased formulations, and then design efficient updates in the form of samples. Here firstorder methods, such as gradient descent methods, play essential roles. We notice that gradient directions for samples rely on the metric over the probability space, which reflects the change of objective/loss functions. In practice, there are several important metrics, often named information metrics from information geometry and optimal transport, including the Fisher-Rao metric (Amari, 1998) and the Wasserstein-2 metric (in short, Wasserstein metric) (Otto, 2001;Lafferty, 1988). In literature, along with a given information metric, the probability space can be viewed as a Riemannian manifold, named density manifold (Lafferty, 1988).
For the Fisher-Rao metric, its gradient flow, known as birth-death dynamics, are important in modeling population games and designing evolutionary dynamics (Amari, 2016). It is also important for optimization problems in discrete probability (Malagò and Pistone, 2014) and machine learning (Ollivier et al., 2017). Recently, the Fisher-Rao gradient has also been applied for accelerating Bayesian sampling problems in continuous sample space (Lu et al., 2019). The Fisher-Rao gradient direction also inspires the design of learning algorithms for probability models. Several optimization methods in machine learning approximate the Fisher-Rao gradient direction, including the Kronecker-factored Approximate Curvature (K-FAC) (Martens and Grosse, 2015) method and adaptive estimates of lower-order moments (Adam) method (Kingma and Ba, 2014).
For the Wasserstein metric, its gradient direction deeply connects with stochastic differential equations and the associated Markov chain Monte Carlo methods (MCMC). An important fact is that the Wasserstein gradient of KL divergence forms the Kolmogorov forward generator of overdamped Langevin dynamics (Jordan et al., 1998). Hence, many MCMC methods can be viewed as Wasserstein gradient descent methods. In recent years, there are also several generalized Wasserstein metrics, such as Stein metric (Liu and Wang, 2016;Liu, 2017), Hessian transport (mobility) metrics (Carrillo et al., 2010;Dolbeault et al., 2009;Li and Ying, 2019) and Kalman-Wasserstein metric (Garbuno-Inigo et al., 2019). These metrics introduce various first-order methods with sampling efficient properties. For instance, the Stein variational gradient descent (Liu and Wang, 2016, SVGD) introduces a kernelized interacting Langevin dynamics. The Kalman-Wasserstein metric introduces a particular mean-field interacting Langevin dynamics (Garbuno-Inigo et al., 2019), known as ensemble Kalman sampling. On the other hand, many approaches design fast algorithms on modified Langevin dynamics. These methods can also be viewed and analyzed by the modified Wasserstein gradient descent, see details in (Ma et al., 2019;Simsekli et al., 2016;Li, 2019). By viewing sampling as optimization problems in the probability space, many efficient sampling algorithms are inspired by classical optimization methods. E.g., Bernton (2018); Wibisono (2019) apply the operator splitting technique to improve the unadjusted Langevin algorithm. Liu et al. (2018); Taghvaei and Mehta (2019); Wang and Li (2019) study Nesterov's accelerated gradient methods in probability space.
In optimization, the Newton's method is a fundamental second-order method to accelerate optimization computations. For optimization problems in probability space, several natural questions arise: Can we systematically design Newton's methods to accelerate sampling related optimization problems? What is the Newton's flow in probability space under information metrics? Focusing on the Wasserstein metric, can we extend the relation between Wasserstein gradient flow of KL divergence and Langevin dynamics? In other words, what is the Wasserstein Newton's flow of KL divergence and which Langevin dynamics does it corresponds to?
In this paper, following (Li, 2018;Wang and Li, 2019), we complete these questions. We derive Newton's flows in probability space with general information metrics. By studying these Newton's flows, we provide the convergence analysis.Focusing on Wasserstein Newton's flows of KL divergence, we derive several analytical examples in one-dimensional space and Gaussian families. Besides, we design two algorithms as particle implementations of Wasserstein Newton's flows in high dimensional sample space. This is to restrict the dual variable (cotangent vector) associated with Newton's direction into a finitedimensional affine function space. A hybrid update of Newton's direction and gradient direction is also introduced. For the concreteness of presentation, we demonstrate the Wasserstein Newton's flow of KL divergence in Theorem 1.
Theorem 1 (Wasserstein Newton's flow of KL divergence). Given a target density ρ * (x) ∝ exp (−f (x)) , where f is a given function, denote the KL divergence between ρ and ρ * by where Z = exp(−f (x))dx. Then the Wasserstein Newton's flow of KL divergence follows where Φ Newton t satisfies the following equation Here we notice that Φ Newton t is the solution to the Wasserstein Newton's direction equation (3). In Figure 1, we provide a sampling (particle) formulation of Wasserstein Newton's flows. We compare formulations among Wasserstein Newton's flows, Wasserstein gradient flows and overdamped Langevin dynamics.
In literature, second-order methods are developed for optimization problems on Riemannian manifold, see (Smith, 1994;Yang, 2007). Here we are interested in density manifold, i.e., probability space with information metrics. Compared to known results in Riemannian optimization, we not only develop methods in probability space but also find efficient sampling representations of the algorithms. In discrete probability simplex with Fisher-Rao metric and exponential family models, the Newton's method has also been studied by (Malagò and Pistone, 2014), known as the second order method in information geometry. Also, Detommaso et al. (2018); Chen et al. (2019) design second-order methods for the Stein variational gradient descent direction. Our approach generalizes these results to information metrics, especially for the Wasserstein metric. On the other hand, the Newton-type MCMC method has been studied in (Simsekli et al., 2016), known as Hessian Approximated MCMC (HAMCMC) method. The differences between HAMCMC and our proposed Newton's Langevin dynamics can be observed from evolutions in probability space. HAMCMC utilizes the Hessian matrix of logarithm of target density function and derives the associated drift-diffusion process. In density space, it is still a linear local partial differential equation (PDE). Newton's Langevin dynamics apply the Hessian operator of KL divergence based on the Wasserstein metric. In density space, the Wasserstein Newton's flow is a nonlocal PDE. A careful comparison of all related Langevin dynamics in analytical (subsection 4.2) and numerical examples (subsection 8.1 8.2) are provided.
We organize this paper as follows. In section 2, we briefly review information metrics and corresponding gradient operators in probability space. We introduce properties of Hessian operators and derive various forms of Newton's flows for different objective functions in section 3. Focusing on Wasserstein Newton's flows of KL divergence, we formulate the Newton's direction equation in one-dimensional sample space (see section 4) and Gaussian families (see section 5). In section 6, the asymptotic convergence rate of Newton's method is proved. Two sampling efficient numerical algorithms are presented in section 7. Several numerical examples for sampling problems are provided in section 8.

Review on Newton's flows and information metrics
In this section, we briefly review Newton's methods and Newton's flows in Euclidean spaces and Riemannian manifolds. Then, we focus on a probability space, in which we introduce information metrics with the associated gradient and Hessian operators. Based on them, we will derive the Newton's flow under information metrics later on. Throughout this paper, we use ·, · and · to denote the Euclidean inner product and norm in R n .
2.1. Finite dimensional Newton's flow. We first briefly review Newton's methods and Newton's flows in Euclidean spaces. Given an objective function f : R n → R, consider an optimization problem: min x∈R n f (x). The update rule of the (damped) Newton's method follows Here α k > 0 is a step size and p k is called the Newton's direction. With α k = 1, we recover the classical Newton's methods. By taking a limit α k → 0, the Newton's method in continuous-time, namely Newton's flow, writeṡ (Euclidean Newton's flow) We next consider an optimization problem on a Riemannian manifold M ⊂ R n . Given an objective function f : M → R, consider min x∈M f (x).
The tangent space T x M and the cotangent space T * x M at x are identical to a linear subspace of R n . For p, q ∈ T x M, let p, q x = p T G(x)q denote an inner product in tangent space T x M at x. Here G(x) is called the metric tensor, which corresponds to a symmetric semi-positive definite matrix in R n×n . For the Euclidean case, we can view x M = R n and G(x) = I, where I is an identity matrix. The Riemannian gradient of f at x is the unique tangent vector v such that the following equality holds for all p ∈ T x M.
The Riemannian Hessian of f at x is a linear mapping from T x M to T x M defined by Here ∇ p grad f (x) is the covariant derivative of grad f (x) w.r.t. the tangent vector p. Detailed definitions of gradient and Hessian operators on a Riemannian manifold can be found in (Huang, 2013, Chapter 1). The update rule of the Newton's method writes Here R x k can be the exponential mapping or the retraction (first-order approximation of the exponential mapping) at x k . Based on the Riemannian metric of M, the exponential mapping uniquely maps a tangent vector to a point in M along the geodesic curve. Different from the Euclidean case, the update of x k+1 is based on the (approximated) geodesic curve of M. In continuous time, the Newton's flow followṡ (Riemannian Newton's flow) Example 1 (Rayleigh's quotient on a sphere). Let M = S n−1 ⊂ R n be a unit sphere in R n and f (x) = x T Ax, where A is a symmetric matrix. Minimizing f (x) is equivalent to finding the eigenvector corresponding to A's smallest eigenvalue. The tangent space and the cotangent space of M follow The inner product p, q x = p T q is identical to the Euclidean inner product in R n . Based on previous definitions, we derive If (A − f (x)I) is invertible, then the inverse of the Hessian operator writes The exponential map at x with direction p satisfies Hence, the update rule of the Newton's method in S n−1 follows From now on, we consider optimization problems in probability space. Suppose that sample space Ω is a region in R n . Let F(Ω) represent the set of smooth functions on Ω. Denote the set of probability density Here E(ρ) is the objective or loss functional. It evaluates certain divergence or metric functional between ρ and a target density ρ * ∈ P(Ω). In machine learning problems, typical examples of E(ρ) include the KL divergence, Maximum mean discrepancy (MMD), cross entropy, etc. Similar to (Euclidean Newton's flow) and (Riemannian Newton's flow), the Newton's flow in probability space (density manifold) takes the form (Information Newton's flow) Here grad and Hess represent the gradient and the Hessian operator with respect to certain information metric, respectively. To understand (Information Newton's flow), we briefly review the information metrics with the associated gradient operators.
2.2. Information metrics. We first define the tangent space and the cotangent space in probability space. The tangent space at ρ ∈ P(Ω) is defined by The cotangent space T * ρ P(Ω) is equivalent to F(Ω)/R, which represents the set of functions in F(Ω) defined up to addition of constants.
2.3. Gradient operators. The gradient operator for the objective functional E(ρ) in (P(Ω), G(ρ)) satisfies Here δE δρ is the L 2 first variation w.r.t. ρ. The gradient flow follows We present gradient operators under either Fisher-Rao metric or Wasserstein metric.
Example 4 (Fisher-Rao gradient operator). The Fisher-Rao gradient operator satisfies The Fisher-Rao gradient flow follows Example 5 (Wasserstein gradient operator). The Wasserstein gradient operator writes The Wasserstein gradient flow satisfies

Information Newton's flow
In this section, we introduce and discuss properties of Hessian operators in probability space. Then, we formulate Newton's flows under information metrics. This is based on the previous definition of gradient operators and the inverse of Hessian operators.
3.1. Information Hessian operators. In this subsection, we review the definition of Hessian operators in probability space and provide the exact formulations of Hessian operators.
For σ ∈ T ρ P(Ω), there exists a unique geodesic curveρ s , which satisfiesρ s | s=0 = ρ and∂ s ρ s | s=0 = σ. The Hessian operator of E(ρ) w.r.t. metric tensor G(ρ) is a mapping Hess E(ρ) : T ρ P(Ω) → T ρ P(Ω), which is defined by We briefly review the concept of self-adjoint operator. Definition 2 (Self-adjoint). Suppose that V is a Hilbert space and let H : V → V * be a linear operator. V * is the adjoint space of V , which consists of all linear functionals on V .
Remark 1. If V = R n is the Euclidean space, then the linear operator H can be viewed as a matrix in R n×n . Then, to say that H is self-adjoint operator is equivalent to say that H is a symmetric matrix.
Combining with the metric tensor, the Hessian operator uniquely defines a self-adjoint mapping H E (ρ) : T * ρ P(Ω) → T ρ P(Ω), which satisfies In Proposition 1, we give an exact formulation of ΦH E (ρ)Φdx and a relationship between H E (ρ) and Hess E(ρ).
Proposition 1. The quantity g ρ (σ, Hess E(ρ)σ) is a bi-linear form of Φ: Here δ 2 E δρ 2 (x, y) is defined by where δ(x) is the Dirac delta function. Here A(ρ) : T * ρ P(Ω) × T * ρ P(Ω) → T * ρ P(Ω) is a bi-linear operator which satisfies Moreover, the operator H E (ρ) satisfies Proof. The geodesic curveρ s satisfies geodesic equation with initial valuesρ s | s=0 = ρ and Φ s | s=0 = Φ. For the first-order derivative, it follows where we utilize the fact that G(ρ s ) is self-adjoint. For the second-order derivative, Based on the definition of H E (ρ), (4) is proved by setting s = 0 in the above formula. To prove (5), we introduce Lemma 1.

Proof.
Because H is self-adjoint and linear, for any Φ ∈ T * ρ P(Ω), it follows which completes the proof.
For simplicity, we define positive definite operators as follows.
Definition 3. Suppose that V is a Hilbert space and let H : V → V * be a self-adjoint linear operator. We say that Proposition 2 provides a sufficient condition to ensure that the Hessian operator is injective (invertible).
The inverse of the Hessian operator Hess E(ρ) follows Now, we are ready to present the Newton's flow in probability space, named information Newton's flow.
We next show several examples of Newton's flows under either Fisher-Rao metric or Wasserstein metric.

Newton's flows under Fisher-Rao metric. For Fisher-Rao metric, the geodesic curveρ
And the bi-linear operator A F (ρ) follows For simplicity, we let E ρ [Φ] = Φρdx, where Φ ∈ T * ρ P(Ω). Proposition 4 (Fisher-Rao Newton's flow). For an objective function E : P(Ω) → R, the Fisher-Rao Newton's flow follows Proof. Based on Proposition 1, we only need to prove that

The left hand side follows
The right hand side satisfies We also observe that Hence, the left hand side is equal to the right hand side.
Example 6 (Fisher-Rao Newton's flow of KL divergence). Suppose that E(ρ) evaluates the KL divergence from ρ to ρ * ∼ exp(−f ). This objective functional also writes We derive that Based on Proposition 4, we can compute that (4) is equivalent to Example 7 (Fisher-Rao Newton's flow of interaction energy). Consider an interaction energy is a kernel function. The interaction energy also formulates the MMD, see details in (Gretton et al., 2012). We can compute that Example 8 (Fisher-Rao Newton's flow of cross entropy). Suppose that E(ρ) is the cross entropy, i.e., reverse KL divergence. It evaluates the KL divergence from a given density

Proposition 4 indicates that
Hence, the operator H F E (ρ) follows 3.3. Newton's flows under Wasserstein metric. For Wasserstein metric, the geodesic curveρ s satisfies Proposition 5 (Wasserstein Newton's flow). For an objective functional E : P(Ω) → R, the Wasserstein Newton's flow follows Here H E (ρ) : T * ρ P(Ω) → T ρ P(Ω) defines a bi-linear form: for Φ ∈ T * ρ P(Ω), Proof. Based on integration by parts, we observe that Combining above two observations with Proposition 1, we derive This proves Proposition 5.
Example 9 (Wasserstein Newton's flow of KL divergence). In this example we prove Theorem 1. As a known fact in (Otto and Villani, 2000) and Gamma calculus (Bakry and Emery, 1985;Li, 2018), the Hessian operator of KL divergence under the Wasserstein metric follows where σ = −∇ · (ρ∇Φ) and · F is the Frobenius norm of a matrix in R n×n . Via integration by parts, we validate that the operator H W E (ρ) follows Example 10 (Wasserstein Newton's flow of interaction energy). Consider an interaction energy Combining with previous computations, Proposition 5 yields that Based on integration by parts, the operator H W E (ρ) is given by Example 11 (Wasserstein Newton's flow of cross entropy). Suppose that E(ρ) evaluates the KL divergence from a given density ρ * to ρ Remark 2. For simplicity of presentations, we only present the Hessian formulas for Fisher-Rao and Wasserstein information metrics. In fact, there are many interesting generalized Hessian formulas in Li (2019) from Hessian transport metrics. We leave systematic studies of Newton's flows for general metrics in future works.

Newton's Langevin dynamics
In this section, we primarily focus on the Wasserstein Newton's flow of KL divergence. We formulate it into the Newton's Langevin dynamics for Bayesian sampling problems. The connection and difference with Newton's flows in Euclidean space and Hessian Approximated MCMC (HAMCMC) method (Simsekli et al., 2016) are discussed.
Let the objective functional E(ρ) = D KL (ρ ρ * ) evaluate the KL divergence from ρ to a target density ρ * (x) ∝ exp(−f (x)) with exp(−f (x))dx < ∞. This specific optimization problem is important since it corresponds to sampling from the target density ρ * . Classical Langevin MCMC algorithms evolves samples following overdamped Langevin dynamics (OLD), which satisfies where B t is the standard Brownian motion. Denote ρ t as the density function of the distribution of X t . The evolution of ρ t satisfies the Fokker-Planck equation A known fact is that the Fokker-Planck equation is the Wasserstein gradient flow (WGF) of KL divergence, i.e.
where we use the fact that δ δρ D KL (ρ t ρ * ) = log ρ + t + f + 1 and ρ∇ log ρ = ∇ρ. It is worth mentioning that OLD can be viewed as particle implementations of WGF (14). From the viewpoint of fluid dynamics, WGF also has a Lagrangian formulation We name above dynamics by the Lagrangian Langevin Dynamics (LLD). Here 'Lagrangian' refers to the Lagrangian coordinates (flow map) in fluid dynamics (Villani, 2008).
Overall, many sampling algorithms follow OLD or LLD. The evolution of corresponding density follows the Wasserstein gradient flow (14). E.g. the classical Langevin MCMC (unadjusted Langevin algorithm) is the time discretization of OLD. The Particle-based Variational Inference methods (ParVI), (Liu et al., 2019) can be viewed as the discretetime approximation of LLD.
In short, we notice that the Langevein dynamics can be viewed as first-order methods for Bayesian sampling problems. Analogously, we derive the Wasserstein Newton's flow of KL divergence in Example 9. It corresponds to certain Langevin dynamics of particle systems, named Newton's Langevin dynamics.
Theorem 2. Consider the Newton's Langevin dynamics where Φ Newton t (x) is the solution to Wasserstein Newton's direction equation (3): Here X 0 follows an initial distribution ρ 0 and ρ t is the distribution of X t . Then, ρ t is the solution to Wasserstein Newton's flow with an initial value ρ 0 = ρ 0 .
Proof. Note that ρ t is the distribution of X t . The dynamics of X t implies Because Φ t satisfies the Wasserstein Newton's direction equation (3), ρ t is the solution to Wasserstein Newton's flow.
Remark 3. We notice that the Newton's Langevien dynamics is different from HAMCMC (Simsekli et al., 2016). HAMCMC approximates the dynamics of Here Γ(x) is a correction term to ensure that ρ t converges to ρ * . The evolution of ρ t follows We formulate the above equation as On the other hand, replacing ∇Φ Newton 4.1. Connections and differences with Newton's flows in Euclidean space. We recall that the density evolution of particle's gradient flow in Euclidean space corresponds to the Wasserstein gradient flow (Villani, 2008). We notice that this relationship does not hold for the Wasserstein Newton's flow.
Consider an objective function: where f (x) is a given smooth function. Here we notice that minimize ρ for E(ρ) in probability space is equivalent to minimize x for f (x) in Euclidean space. Namely, the support of the optimal solution ρ contains all global minimizers of f (x). The gradient flow in Euclidean space of each particle follows A known fact is that the density evolution of particles satisfies the following continuity equation , which is the Wasserstein gradient flow of E(ρ) in probability space.
We next show that Newton's flow in Euclidean space of each particle does not coincide with the Wasserstein Newton's flow in probability space. For simplicity, we assume that f (x) is strictly convex so ∇ 2 f (x) is invertible for all x. Here, the Euclidean Newton's flow of each particle follows dX t = −(∇ 2 f (X t )) −1 ∇f (X t )dt. The density evolution of particles satisfies the continuity equation On the other hand, the Wasserstein Newton's flow writes where Φ Newton t is the unique solution to We note that in general equation (17) can be different from equation (18). Later on in Lemma 2, we formulate the following Hodge decomposition of the Euclidean Newton's direction Here, the constraint on ξ t does not necessarily ensure that ∇ · (ρ t ξ t ) = 0. Hence, equation (17) can be different from equation (18).
Hence, ξ also exists. We next prove the uniqueness. Suppose that Because H is positive definite, this yields that Φ 1 − Φ 2 = 0 (up to a spatial constant).

4.2.
Newton's Langevin dynamics in one dimensional sample space. In this subsection, we provide examples of Newton's Langevin dynamics in one dimensional sample space. In particular, similar to the OrnsteinUhlenbeck (OU) process in classical Langevin dynamics, we derive a closed form solution to Newton's OU process.
Here we assume that Ω = R and f is strictly convex. The essence of Newton's Langevin dynamics is to compute Φ Newton t from the Wasserstein Newton's direction equation (3). Proposition 2 ensures the uniqueness of the solution to (3). For the simplicity of notations, we neglect the subscript t.
We consider the case where f (x) and (log ρ) (x) are affine functions. Then, ODE (20) has a closed-form solution. Applying ODE (20), we present the exact formulation of Newton's Langevin dynamics.
Proposition 7. Assume that f (x) = (2Σ * ) −1 (x − µ * ) 2 , where Σ * > 0 and µ * are given. Suppose that the particle system X 0 follows the Gaussian distribution. Then X t follows a Gaussian distribution with mean µ t and variance Σ t . The corresponding NLD satisfies And the evolution of µ t and Σ t satisfies The explicit solutions of µ t and Σ t satisfy Proof. In section 5 Proposition 10, we show that if the evolution of X t follows NLD, then X t follows the Gaussian distribution. We first solve the Newton's direction from ODE (20). Suppose that (log ρ) (x) = Σ −1 (x − µ). The ODE turns to be We can examine that the following u is a solution to the above ODE.
The dynamics of µ t satisfies This indicates that µ t = µ * + e −t (µ 0 − µ * ). The dynamics of Σ t follows We can rewrite that Integrating both sides of the above equation yields Hence, the solution Σ t follows Now, we are ready to compare the NLD with OLD, LLD and HAMCMC. Here we consider f (x) = (2Σ * ) −1 (x − µ * ) 2 , where Σ * > 0 and µ * are given. The OLD satisfies which is also known as the Ornstein-Uhlenbeck process. And LLD writes The mean µ t and variance Σ t of OLD and LLD both satisfy On the other hand, HAMCMC follows the dynamics For HAMCMC, the evolution of mean µ t follows and the evolution of variance Σ t satisfies The mean µ t and variance Σ t of HAMCMC follows We summarized our results in Table 1.

Dynamics
Particle Mean and variance Compared to OLD and LLD, the exponential convergence rate of µ t and Σ t in NLD does not depend on Σ * . This fact shows that the NLD is the Newton's flow for both the evolution of mean and variance in Gaussian process. We also note that the convergence rates of mean and variance are different in HAMCMC, while they are same in NLD. Later on in section 8, we use numerical examples to further demonstrate the differences between NLD and HAMCMC.

Newton's flows in Gaussian families
In this section, we study information Newton's flows in Gaussian families with respect to Wasserstein metric. Let P n and S n represent the space of symmetric positive definite matrices and symmetric matrices with size n × n respectively.
We let N 0 n denote multivariate Gaussian densities with zero means. Each ρ ∈ N 0 n is uniquely determined by its covariance matrix Σ ∈ P n . So we can view N 0 n P n . The Wasserstein metric G W (ρ) on P(R n ) induces the Wasserstein metric G W (Σ) on P n , see (Takatsu, 2008;Modin, 2016;Malagò et al., 2018). For Σ ∈ P n , tangent space and cotangent space follow T Σ P n T * Σ P n S n . Definition 4 (Wasserstein metric in Gaussian families). Given Σ ∈ P n , the Wasserstein metric tensor G W (Σ) : S n → S n is defined by It defines an inner product on the tangent space T Σ P n . Namely, for A 1 , A 2 ∈ T Σ P n S n Here S i ∈ T * Σ P n S n is the solution to discrete Lyapunov equation For Σ ∈ P n , there exits a unique solution to discrete Lyapunov equation. Again, we focus on the case where the objective functional E(Σ) evaluates the KL divergence from ρ with covariance matrix Σ to a target Gaussian density ρ * with covariance matrix Σ * . Then, E(Σ) satisfies Proposition 8 (Gradient and Hessian operators in P n ). The gradient operator follows And the Hessian operator satisfies that for all A ∈ S n , g W Σ (A, Hess W E(Σ)A) = 4 tr(SΣS(Σ * ) −1 ) + 4 tr(S 2 ), where S is the unique solution to A = 2(ΣS + SΣ).
This completes the proof.
Example 12. In one dimensional case, the second equation in (22) has an explicit solution where Y t > 0. Then, the first equation in (22) turns to Hence, the Newton's flow (23) coincides with Newton's flow of f (X) in Euclidean space. We also note that (23) is identical to the evolution of Σ t in Proposition 7 by substituting Σ t = Y 2 t .
Proposition 10 ensures the existence of information Newton's flows in Gaussian families.
Proof. We first prove that Σ t is positive definite. We formulate that As a result, E(Σ t ) is non-increasing. Applying the idea of proof in (Wang and Li, 2019, Theorem 1), we can establish that Σ t is positive definite. Then, we examine that Φ t satisfies (3). We observe that Hence, we derive all four terms in the above equation as follows. First, it is easy to observe that We can also compute that Taking S t = I into the above equation yields Because (Σ t , S t ) satisfies (22), we have This completes the proof.
Remark 5. Proposition 10 also indicates the fact that if the evolution of X t follows NLD, then X t follows the Gaussian distribution.

Information Newton's method and convergence analysis
In this section, we introduce a general update rule of information Newton's method for probability densities and analyze its convergence rate.
The general update rule of the Newton's method follows Here α k > 0 is a step size. Here Exp α k ρ k (Φ k ) is the exponential map, which we elaborate in subsection 6.1. 6.1. Riemannian structure of probability space. We first provide some background knowledge for the Riemannian structure of probability space. Given a metric tensor G(ρ) and two probability densities ρ 0 , ρ 1 ∈ P(Ω), we define the distance D(ρ 0 , ρ 1 ) as follows The minimizerρ s is the geodesic curve (6) connecting ρ 0 and ρ 1 . For the Wasserstein metric G(ρ) = G W (ρ), D(ρ 0 , ρ 1 ) is the Wasserstein-2 distance between ρ 0 and ρ 1 .
For simplicity, we define the exponential map and other Riemannian operators on cotangent space.
Definition 6 (Parallelism). Suppose thatρ s is the geodesic curve connecting ρ 0 and ρ 1 . We say that τ : We define high-order derivatives on the cotangent-space in Proposition 11.
Proof. We first show that From the definition, it follows that From (25), we can recursively compute that The Taylor expansion of E(Exp s ρ (Φ)) w.r.t. s completes the proof.
Remark 6. For the first-order derivative of E(ρ), it follows For the second-order derivative of E(ρ), it satisfies 6.2. Convergence analysis. To ensure the convergence of Newton's method in Euclidean space, it is assumed that ∇ 2 f (x) is positive definite around a small neighbour of the optimal solution x * . To prove the convergence rate of Newton's method in the probability space, we assume the following fact analogously.
Assumption 1. Assume that there exists > 0, such that for all ρ satisfying D(ρ, ρ * ) < , it follows Relying on Assumption 1, Theorem 3 shows the quadratic convergence rate of the Newton's method in the probability space.
Because the arbitrary choice of s, it follows that This completes the proof.

Proof. Consider the auxiliary function
A(s) = ∇E(Exp s ρ k (T k ))(τ s Φ). Directly from the definition of high-order derivatives, it follows ∂ ∂s . Hence, we can compute the Taylor expansion . On the other hand, we notice that Let τ be the parallelism from ρ k to ρ k+1 . There exists a unique R k ∈ T * ρ k P(Ω) such that For arbitrary Ψ ∈ T * ρ k P(Ω), we have Here the last equality comes from Lemma 4. Based on the definition of parallelism, we notice the fact (26), applying Assumption 1 and utilizing Lemma 3 yields In order to provide an estimation on R k ρ k , we introduce Lemma 5.
Lemma 5. We have following estimations Proof. From Assumption 1 and Lemma 3, it follows that We also notice that from Lemma 4, As a result, we have Φ k ρ k = O δE δρ k ρ k = O ( T k ρ k ). We also note the triangle This yields T k+1 ρ k+1 = O( T k ρ k ).
We finally show the estimation of R k ρ k .
Lemma 6. For all Ψ ∈ T * ρ k P(Ω), it follows Proof. Based on the first-order approximation of the exponential map and the parallelsim, we have the following estimations This completes the proof.
. We also note that T k 2 ρ k = D(ρ k , ρ * ) 2 . As a result, we have Remark 7. We follow the proof technique in (Smith, 1994, Section 4, Theorem 4.4). Nevertheless, our estimation of R k ρ k is much coarser. With future developments in the theories of curvature estimations and parallelism over density manifold, we expect to prove that R k ρ k = O(D(ρ k , ρ * ) 3 ). The other issue is that our proof is based on the formulation in probability space. In practice, our algorithm is often formulated into samples. In a future work, we shall conduct related proofs in term of samples.

Particle implementation of Wasserstein Newton's method
In this section, we design sampling efficient implementations of Wasserstein Newton's method. Focusing on Wasserstein Newton's flow of KL divergence, we introduce a variational formulation for computing the Wasserstein Newton's direction. By restricting the Newton's direction in a linear subspace, we derive a sampling efficient algorithm. We also provide a modified algorithm to simplify the calculation of the Newton's direction. Besides, a hybrid method between Newton's Langevin dynamics and overdamped Langevin dynamics is provided.
We breifly review update rules of Newton's methods and hybrid methods in Euclidean space. In each iteration, the update rule of Newton's method follows Suppose that f (x) is strictly convex. Namely, ∇ 2 f (x) is positive definite for all x ∈ R n . To compute the Newton's direction p k , it is equivalent to solve the following variational problem min p∈R n p T ∇ 2 f (x k )p + 2∇f (x k ) T p. In practice, the Newton's direction may not lead to the decrease in the objective function, especially when f (x) is non-convex. Nevertheless, the Newton's method often converges when the update is close to the minimizer. One way to overcome this problem is the hybrid method. Consider a hybrid update of the Newton's direction and the gradient's direction Following above ideas in Euclidean space, we present a particle implementation of information Newton's method. Here we connect density ρ k ∈ P(Ω) with a particle system Here ρ k can be approximated by the average of delta measures 1 N N i=1 δ(x−x i k ). We update each particle by Here Φ k is the solution to Wasserstein Newton's direction equation (3). We next formulate a variational formulation for estimating the Newton's direction as follows.
Proposition 12. Suppose that H : T * ρ P(Ω) → T ρ P is a linear self-adjoint operator and H is positive definite. Let u ∈ T ρ P. Then the minimizer of variational problem min Φ∈T * ρ P(Ω) Proof. Since H is linear and self-adjoint, then the optimal solution of satisfies Hence, Φ satisfies HΦ = u. On the other hand, let Φ satisfy HΦ = u. Then, for any Ψ ∈ T * ρ P(Ω), it follows The last inequality is based on the fact that H is positive definite. Hence, Φ is the optimal solution to the proposed variational problem. This completes the proof.
Suppose that f is strongly convex, or equivalent, ∇ 2 f (x) is positive definite for x ∈ Ω. Then, the operator H E (ρ) defined in (13) is positive definite. Proposition 12 indicates that solving Wasserstein Newton's direction equation (3) is equivalent to optimizing the following variational problem.
Here we denote v 2 A = v T Av. For general f , suppose that there exists ≥ 0 such that ∇ 2 f (x) + I is strictly convex for x ∈ Ω. Then, we consider a regularized problem Remark 8. The regularized problem (27) is equivalent to Namely, we penalize the squared norm of Φ, which is induced by the Wasserstein metric.
Based on (27), we are able to present a particle implementation of the Newton's direction. For the case where ρ k is Gaussian and ρ * is Gaussian, we have wherex k andΣ k are mean and empirical covariance matrix of {x i k } N i=1 , respectively. For general ρ k and ρ * , we can approximate ∇ log ρ k via kernel density estimation (KDE) (Gretton et al., 2012). Namely, we approximate ∇ log ρ k by .
Here k(x, y) is a given positive kernel. A typical choice of k(x, y) is a Gaussian kernel with bandwidth h > 0, such that The bandwidth h is critical to the estimation of ∇ log ρ k .
Remark 9. We can also approximate ∇ log ρ k via a kernel function by the blob method (Carrillo et al., 2019) and the diffusion map (Taghvaei and Mehta, 2019).
7.1. Affine Wasserstein Newton's method. In high dimensional sample space, directly solving (27) for Φ ∈ T * ρ k P(Ω) can be difficult. In this subsection, we present a practical method to approximate the Newton's direction Φ Newton t . To deal with this, we restrict Φ in infinite dimensional function space F(Ω) to a finite-dimensional subspace.
Here the subspace is designed as the following linear vector space spanned by basis functions.
Proposition 13. Suppose that Φ(x) = a, ψ(x) . Then, the optimization problem (27) constrained in Φ ∈ Ψ turns to be where B k , D k ∈ R m×m and c k ∈ R m . The detailed formulations of B k , D k and c k are provided as follows.
If B k + D k is positive definite, the optimal solution follows a = −(B k + D k ) −1 c k . The optimal solution Φ Newton follows Φ Newton (x) = a, ψ(x) .
Proof. We denote the Jacobian ∇ψ(x) ∈ R n×m . As a result, J(a) turns to be We can further compute that This completes the proof.
This affine approximation technique has been used in approximating natural gradient direction in . Hence, we call our method affine information Newton's method.
We further consider a special case where ψ i are linear and quadratic polynomials. Namely, we set m = n(n + 3)/2 and In other words, we assume that Φ(x) takes the form Φ(x) = 1 2 xSx + b T x, where S ∈ R n×n is a symmetric matrix. For simplicity, we denote v i k = ∇f (x i k ) + ξ k (x i k ). Then, the variational problem (27) turns to be min S∈S n ,b∈R n J(S, b) = tr(S 2 ) + We further restrict S to be a diagonal matrix S = diag(s), where s ∈ R n . Here we use diag(s) to represent the diagonal matrix in R n×n with elements in s as its diagonal components. Then, we can write the objective function to be where we denote H k ∈ R 2n×2n via , and u k ∈ R 2n via Hence, the optimal solution for minimizing J(s, b) follows The overall algorithm is summarized in 1.

Algorithm 1 Particle implementation of affine Wasserstein Newton's method
Require: initial positions {x i 0 } N i=1 , ≥ 0, step sizes α k , maximum iteration K. 1: Set k = 0. 2: while k < K and the convergence criterion is not met do 3: Here ξ k is an approximation of ∇ log ρ k .

4:
Calculate H k by , and formulate u k by

5:
Compute s k and b k by s k b k = −(H k ) −1 u k .

6:
Update particle positions by

7:
Set k = k + 1. 8: end while 7.2. Modified affine Wasserstein Newton's method. Optimizing J(S, b) defined in (29) for S ∈ S n can be challenging and time-consuming. Here we consider a modification of Newton's method. In other words, we consider the following problem. Here where ≥ 0 is a parameter. Namely, we use the averaged Hessian F k to replace the Hessian ∇ 2 f (x i k ) + I of individual particle. We notice that this modification does not change the original variational problem (29) if ρ * is a Gaussian distribution.
With this modification, we can compute the direction by solving a quadratic matrix optimization problem as follows.
Proposition 14. The optimal solution to (30) follows and S k is the solution to the quadratic matrix optimization problem min S∈S n 2 tr(S 2 ) + tr(SΣ k SF k ) + tr(SF k SΣ k ) − 2 tr(T k S).
HereΣ k is an empirical covariance matrix of x i k and barT k is an empirical covariance matrix between x i k and v i k .
Proof. For a fixed S, the optimal b for minimizingĴ(S, b) defined in (30) For simplicity, letx k = 1 Then, the above problem is equivalent to Then, problem (30) is equivalent to problem (31).
Because (31) is a quadratic optimization problem on S n , it can be efficiently solved by the conjugate gradient (CG) method. Detailed description of the CG method can be found in Appendix B. Denote the solution by S k . The update rule for x i k+1 follows The overall algorithm is summarized in Algorithm 2 Algorithm 2 Particle implementation of modified affine Wasserstein Newton's method , ≥ 0, step sizes α k , maximum iteration K. 1: Set k = 1 and S 0 = I. 2: while k ≤ K and the convergence criterion is not met do 3: Here ξ k is an approximation of ∇ log ρ k .

4:
Compute the averaged Hessian with regularization by F k = I + 1

8:
Solve S k from the problem (30) via the CG method using the initial guess S k−1 .

9:
Update particle positions by . 10: end while 7.3. The hybrid method. In practice, the Newton's direction may not be a descent direction if the update is far away from the minimizer. To overcome this issue, we propose a hybrid update of Newton's direction and gradient direction.
Let λ k ≥ 0 be a parameter. Here we recall that there are two choices for using the gradient direction. Namely, if we use overdamped Langevin dynamics as gradient direction, the hybrid update rule follows where B k ∼ N (0, I). If we use Lagrangian Langevin dynamics as gradient direction, the hybrid update rule satisfies Remark 10. It worths mentioning that our algorithm corresponds to the following hybrid Langvien dynamics where B t is the standrad Brownian motion, λ t , > 0 are given parameters and Φ −Newton t is the minimizer of variational problem (27).

Numerical experiments
In this section, we present numerical experiments to demonstrate the strength of information Newton's flows.
8.1. Toy examples. We compare particle implementations among affine Wasserstein Newton's method (WNewton), modified affined Wasserstein Newton's method (mWNewton), Wasserstein gradient flow (WGF) and the Hessian Approximated Lagrangian Langevin dynamics (HALLD). We notice that the update rule of WGF satisfies ). The update rule of HALLD follows We note that the density evolutions of HALLD and HAMCMC are identical to each other. In other words, we replace the Brownian motion in HAMCMC by ξ k in HALLD. Here ξ k is an approximation of ∇ log ρ k . For all compared methods, we use constant step sizes. For the calculation of ξ k , we apply KDE with Gaussian kernels and the kernel bandwidth is selected by Brownian Motion method (Wang and Li, 2019)[section 5.1]. This method adaptively learns the bandwidth from samples generated by Brownian motions.
We first consider a 1D target density ρ * (x) ∝ exp (−f (x)), where f (x) = 1 2 (x 2 − 1) 2 . For WGF, we set α k = 0.01. For WNewton and mWNewton, we let α k = 1, = 0 and λ k = 0. Namely, we do not apply the hybrid update. For HALLD, we set α = 0.01. The sample number follows N = 100. The initial distribution follows N (0, 0.01). We plot the distribution after 2, 5, 10, 20 iterations in Figure 2. Although we use affine approximations to compute the Newton's direction, WNewton tends to converge to the target density and it is faster than WGF. Even with a small step size, some samples generated by HALLD tends to go to infinity. This is because f (x) is not a convex function. Then, we let the target density ρ * to be a toy bimodal distribution (Rezende and Mohamed, 2015) in R 2 . For WGF, we set α k = 0.1. For WNewton and mWNewton, we apply the hybrid update and set α k = 0.2, = 0 and λ k = 0.5. The sample number follows N = 200. The initial distribution follows N ([0, 10] , I). We plot the distribution after 5, 10, 20, 40 iterations in Figure 3. HALLD fails because ∇ 2 f becomes singular on certain sample points. We observe that the Wasserstein Newton's direction helps samples converge faster towards the target density. 8.2. Gaussian families. The target density ρ * is a Gaussian distribution with zero mean on R 100 . The covariance matrix of ρ * is Σ * and W * = (Σ * ) −1 . Let L and β be the largest/smallest eigenvalue of W * . The condition number of W * is defined as κ = L/β. The initial distribution follows N (0, I) and the number of particle N = 600. For a Gaussian target density, mWNewton is exactly the Wasserstein Newton's method. We compare mWNewton with WGF, the particle implementation of Accelerated Information Gradient flow (Wang and Li, 2019, AIG) and HALLD. For WGF and AIG, we set the step size to be α k = 1/(2L). For AIG, we use the restart technique. For mWNewton, we set α k = 1, = 0 and λ k = 0. Namely, we add no regularization and do not apply the hybrid update. For HALLD, we choose two step sizes α k = 1 and α k = 0.5. For a particle system {X i k } N i=1 , we record the KL divergence E(Σ k ) (21) using the empirical covariance matrixΣ k . The results are collected in Figure 4. We observe that mWNewton converges to the optimal distribution in less than 30 iterates, while AIG takes nearly 2000 iterations. In CPU-time, mWNewton also has competitive performance with AIG. With α k = 1, HALLD does not converge while with α k = 0.5, HALLD converges rapidly.
8.3. Bayesian logistic regression. We perform the standard Bayesian logistic regression experiment on the Covertype dataset, following the settings in (Liu and Wang, 2016). We compare WNewton and mWNewton with MCMC, SVGD (Liu and Wang, 2016), WNes (Liu et al., 2019), WGF and AIG. For the calculation of ξ k in WGF, AIG, WNes, WNewton and mWNewton we use KDE with Gaussian kernel and the bandwidth is selected by the median method, which is the same as (Liu and Wang, 2016). The sample number follows We first elaborate on the choice of step sizes. The initial step sizes for the compared methods are given in Table 2. The step size of SVGD is adjusted by Adagrad[cite] , which is same as (Liu and Wang, 2016). For WRes, the step size is give by τ l = τ 0 /l 0.9 for l ≥ 1. For MCMC, WGF and AIG, the step size is multiplied by 0.9 every 100 iterations. For WNewton and mWNewton, the step size is multiplied by 0.82 every 100 iterations.
We then elarborate on the implementation details of compared methods. The parameters for Wnes are identical to (Liu et al., 2019). For AIG, we apply the adaptive restart technique. For WNewton and mWNewton, we apply the hybrid update and set = 0, λ k = 0.005.
From Figure 5, we observe that AIG and mWNewton has competitive performance in test accuracy and they outperform other methods. Although the performance of WNewton is not stable in first 500 iterations, WNewton achieves larger test log-likelihood value than other methods (Larger is better). And mWNewton has slightly smaller test log-likelihood value than WNewton.

Conclusion and discussion
In this paper, we introduce information Newton's flows (second-order optimization methods) for optimization problems in probability space arising from Bayesian statistics, inverse problems, and machine learning. Here two information metrics, such as Fisher-Rao metric and Wasserstein-2 metric, are considered. Several examples and convergence analysis of the proposed second-order methods are provided. Following the fact that the Wasserstein gradient flow of KL divergence formulates the Langevin dynamics, we derive the Wasserstein Newton's flow of KL divergence as Newton's Langevin dynamics. Focusing on Newton's Langevin dynamics, we study analytical examples in one-dimensional sample space and Gaussian families. We further propose practical sampling efficient algorithms to Our work opens a door to design high-order optimization methods for Bayesian sampling, machine learning and inverse problems. We observe that high-order derivatives from information metrics are very useful in designing sampling efficient optimization methods. In future works, we shall study the following questions. For example, what is the sampling efficient Quasi-Newton's method for information metrics? Several natural choices of quasi-Newton's method include the Rank one formulation (SR1) and BFGS. What is the other efficient method to approximate the Wasserstein Newton's direction? In our paper, we apply the affine method to restrict the Newton's direction in a finite dimension subspace. We expect to apply kernelized methods to improve the approximation. Lastly, we shall further propose the information Newton's flow not only for probability space but also for probability models, especially generative models. Objective functional E(ρ) H W E (ρ)Φ KL divergence: (ρ log ρ + f ρ)dx.
Algorithm 3 The CG method in the space of symmetric matrices.

7:
Set i = i + 1. 8: end while For the subproblem in Algorithm 2, the operator H follows HS = 2S + Σ k SF k + F k SΣ k , and we have T =T k .