<figure><img src="image_1.png" /></figure>

HAL
open science

# Optimal Riemannian metric for Poincaré inequalities and how to ideally precondition Langevin dymanics

Tiangang Cui, Xin Tong, Olivier Zahm

<figure><img src="image_2.png" /></figure>

## To cite this version:

Tiangang Cui, Xin Tong, Olivier Zahm. Optimal Riemannian metric for Poincaré inequalities and how to ideally precondition Langevin dymanics. 2024. hal-04526677

HAL Id: hal-04526677

https://inria.hal.science/hal-04526677

Preprint submitted on 2 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal Riemannian metric for Poincaré inequalities and how to ideally precondition Langevin dymanics

Tiangang Cui, Xin Tong, Olivier Zahm

<figure><img src="image_1.png" /></figure>

## To cite this version:

Tiangang Cui, Xin Tong, Olivier Zahm. Optimal Riemannian metric for Poincaré inequalities and how to ideally precondition Langevin dymanics. 2024. hal-04526677

HAL Id: hal-04526677

https://inria.hal.science/hal-04526677v1

Preprint submitted on 2 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal Riemannian metric for Poincaré inequalities
and how to ideally precondition Langevin dynamics

Tiangang Cui*, Xin Tong† Olivier Zahm‡

April 2, 2024

## Abstract

Poincaré inequality is a fundamental property that rises naturally in different branches of mathematics. The associated Poincaré constant plays a central role in many applications, since it governs the convergence of various practical algorithms. For instance, the convergence rate of the Langevin dynamics is exactly given by the Poincaré constant. This paper investigates a Riemannian version of Poincaré inequality where a positive definite weighting matrix field (i.e. a Riemannian metric) is introduced to improve the Poincaré constant, and therefore the performances of the associated algorithm. Assuming the underlying measure is a *moment measure*, we show that an optimal metric exists and the resulting Poincaré constant is 1. We demonstrate that such optimal metric is necessarily a *Stein kernel*, offering a novel perspective on these complex but central mathematical objects that are hard to obtain in practice. We further discuss how to numerically obtain the optimal metric by deriving an implementable optimization algorithm. The resulting method is illustrated on a few simple but nontrivial examples, where solutions are revealed to be rather sophisticated. We also demonstrate how to design efficient Langevin-based sampling schemes by utilizing the precomputed optimal metric as a preconditioner.

**Keywords:** Poincaré inequalities, Langevin SDE, Stein kernels, Moment measures.

## 1 Introduction

In 1890, Poincaré inequality was first proposed in [46] as a result of calculus. Later, its deep connection with various directions of mathematics, such as functional analysis, geometry, and probability theory have been discovered and thoroughly studied [5, 9, 22, 36, 57]. Recently, the significance of Poincaré in applied mathematics have also been recognized in a wide range of applications. For examples, the convergence analysis of stochastic differential equations [7, 55], MCMC samplers [3, 4, 24], sensitivity analysis [33, 49, 52] and dimension reduction [29, 37, 42, 56], to name just a few. Given a

<footnote>*tiangang.cui@sydney.edu.au</footnote>

<footnote>†xin.t.tong@nus.edu.sg</footnote>

<footnote>‡olivier.zahm@inria.fr</footnote>
Optimal Riemannian metric for Poincaré inequalities
and how to ideally precondition Langevin dynamics

Tiangang Cui,* Xin Tong,† Olivier Zahm‡

April 2, 2024

**Abstract**

Poincaré inequality is a fundamental property that rises naturally in different branches of mathematics. The associated Poincaré constant plays a central role in many applications, since it governs the convergence of various practical algorithms. For instance, the convergence rate of the Langevin dynamics is exactly given by the Poincaré constant. This paper investigates a Riemannian version of Poincaré inequality where a positive definite weighting matrix field (i.e. a Riemannian metric) is introduced to improve the Poincaré constant, and therefore the performances of the associated algorithm. Assuming the underlying measure is a *moment measure*, we show that an optimal metric exists and the resulting Poincaré constant is 1. We demonstrate that such optimal metric is necessarily a *Stein kernel*, offering a novel perspective on these complex but central mathematical objects that are hard to obtain in practice. We further discuss how to numerically obtain the optimal metric by deriving an implementable optimization algorithm. The resulting method is illustrated on a few simple but nontrivial examples, where solutions are revealed to be rather sophisticated. We also demonstrate how to design efficient Langevin-based sampling schemes by utilizing the precomputed optimal metric as a preconditioner.

**Keywords:** Poincaré inequalities, Langevin SDE, Stein kernels, Moment measures.

# 1 Introduction

In 1890, Poincaré inequality was first proposed in [46] as a result of calculus. Later, its deep connection with various directions of mathematics, such as functional analysis, geometry, and probability theory have been discovered and thoroughly studied [5, 9, 22, 36, 57]. Recently, the significance of Poincaré in applied mathematics have also been recognized in a wide range of applications. For examples, the convergence analysis of stochastic differential equations [7, 55], MCMC samplers [3, 4, 24], sensitivity analysis [33, 49, 52] and dimension reduction [29, 37, 42, 56], to name just a few. Given a

<footnote>*tiangang.cui@sydney.edu.au</footnote>

<footnote>†xin.t.tong@nus.edu.sg</footnote>

<footnote>‡olivier.zahm@inria.fr</footnote>
probability measure $\mu$ on $\mathbb{R}^d$, the Poincaré constant $C(\mu) \in [0, \infty)$ is the smallest constant such that the Poincaré inequality

$$ \mathrm{Var}_{\mu}(f) \leq C(\mu) \int \|\nabla f(x)\|^2 d\mu(x) \quad (1) $$

holds for any sufficiently smooth function $f: \mathbb{R}^d \to \mathbb{R}$. Here, $\mathrm{Var}_{\mu}(f) = \int (f - \int f d\mu)^2 d\mu$ denotes the variance of $f$ under $\mu$ and $\| \cdot \|$ is the Euclidean norm of $\mathbb{R}^d$. Showing that the Poincaré constant $C(\mu)$ is finite, and providing a tractable estimate of its optimal value, has long been a challenging and active problem in pure and applied mathematics [1, 44, 45].

The Poincaré constant plays a central role in determining the performance and the convergence rate of numerous stochastic algorithms. A classical example is the (over-damped) Langevin stochastic differential equation (SDE)

$$ dX_t = -\nabla V(X_t)dt + \sqrt{2}dB_t, \quad (2) $$

where $B_t$ is the standard Brownian motion in $\mathbb{R}^d$ and $V: \mathbb{R}^d \to \mathbb{R}$ is a smooth potential function. Denoting by $\mu_t$ the probability measure of the solution $X_t$ to (2), it is well known that $\mu_t$ converges to the equilibrium distribution

$$ d\mu(x) \propto \exp(-V(x))dx, \quad (3) $$

if and only if $\mu$ satisfies the Poincaré inequality with $C(\mu) < \infty$, for instance, see [7] and Proposition 1 below. In particular, $\mu_t$ converges to $\mu$ at a rate of $e^{-2t/C(\mu)}$, which highlights the central role of the Poincaré constant.

The Bakry-Émery theorem [6, 7] provides the baseline estimate for the Poincaré constant. It states that $C(\mu) \le 1/\rho$ for any measure $d\mu(x) \propto \exp(-V(x))dx$ with $\rho$-convex potential $V$, i.e. $V$ is twice-differentiable and its Hessian is uniformly bounded as $Hess(V(x)) \succeq \rho I_d$ for some constant $\rho > 0$. Here, $\succeq$ denotes the Loewner order and $I_d \in \mathbb{R}^{d \times d}$ the identity matrix. However, this convexity assumption may not be satisfied by many probability measures of interest, such as those of multi-modal or heavy-tailed distributions. In these situations, the Poincaré constant can be arbitrarily large, or it may even be infinite.

In this paper, we consider a Riemannian version of the Poincaré inequality:

$$ \mathrm{Var}_{\mu}(f) \leq C(\mu, W) \int \nabla f(x)^{\top} W(x) \nabla f(x) d\mu(x), \quad (4) $$

where $W: x \mapsto W(x) \in S_+^d$ is a field taking value in the convex cone $S_+^d \subset \mathbb{R}^{d \times d}$ of symmetric positive semi-definite matrices. Again, $C(\mu, W)$ is the smallest constant such that (4) holds for any sufficiently smooth function $f$. The field $W$ provides a Riemannian metric on the space $\mathbb{R}^d$ that assigns to each $x \in \mathbb{R}^d$ an inner product $\langle u, v \rangle_x = u^\top W(x) v$, where $u, v \in \mathbb{R}^d$. The inequality (4) extends existing work in several ways. It encompasses the Brascamp-Lieb inequality [13, 18] that takes the same form as (4) but requires a strongly log-concave $\mu$ to have $W(x) = (\mathrm{Hess} V(x))^{-1}$ and $C(\mu, W) = 1$. It also generalizes the scalar-weighted Poincaré inequalities considered in [14, 16], in which $W(x) = w(x)I_d$ with a positive scalar field $w: x \to w(x) \in \mathbb{R}_+$. Such scalar-weighted Poincaré inequalities have been thoroughly analyzed in [15, 51, 52] in dimension $d=1$.
probability measure $\mu$ on $\mathbb{R}^d$, the Poincaré constant $C(\mu) \in [0, \infty)$ is the smallest constant such that the Poincaré inequality

$$ \mathrm{Var}_{\mu}(f) \leq C(\mu) \int \|\nabla f(x)\|^2 d\mu(x) \quad (1) $$

holds for any sufficiently smooth function $f: \mathbb{R}^d \to \mathbb{R}$. Here, $\mathrm{Var}_{\mu}(f) = \int (f - \int f d\mu)^2 d\mu$ denotes the variance of $f$ under $\mu$ and $\| \cdot \|$ is the Euclidean norm of $\mathbb{R}^d$. Showing that the Poincaré constant $C(\mu)$ is finite, and providing a tractable estimate of its optimal value, has long been a challenging and active problem in pure and applied mathematics [1, 44, 45].

The Poincaré constant plays a central role in determining the performance and the convergence rate of numerous stochastic algorithms. A classical example is the (over-damped) Langevin stochastic differential equation (SDE)

$$ dX_t = -\nabla V(X_t)dt + \sqrt{2}dB_t, \quad (2) $$

where $B_t$ is the standard Brownian motion in $\mathbb{R}^d$ and $V: \mathbb{R}^d \to \mathbb{R}$ is a smooth potential function. Denoting by $\mu_t$ the probability measure of the solution $X_t$ to (2), it is well known that $\mu_t$ converges to the equilibrium distribution

$$ d\mu(x) \propto \exp(-V(x))dx, \quad (3) $$

if and only if $\mu$ satisfies the Poincaré inequality with $C(\mu) < \infty$, for instance, see [7] and Proposition 1 below. In particular, $\mu_t$ converges to $\mu$ at a rate of $e^{-2t/C(\mu)}$, which highlights the central role of the Poincaré constant.

The Bakry-Émery theorem [6, 7] provides the baseline estimate for the Poincaré constant. It states that $C(\mu) \le 1/\rho$ for any measure $d\mu(x) \propto \exp(-V(x))dx$ with $\rho$-convex potential $V$, i.e. $V$ is twice-differentiable and its Hessian is uniformly bounded as $Hess(V(x)) \succeq \rho I_d$ for some constant $\rho > 0$. Here, $\succeq$ denotes the Loewner order and $I_d \in \mathbb{R}^{d \times d}$ the identity matrix. However, this convexity assumption may not be satisfied by many probability measures of interest, such as those of multi-modal or heavy-tailed distributions. In these situations, the Poincaré constant can be arbitrarily large, or it may even be infinite.

In this paper, we consider a Riemannian version of the Poincaré inequality:

$$ \mathrm{Var}_{\mu}(f) \leq C(\mu, W) \int \nabla f(x)^{\top} W(x) \nabla f(x) d\mu(x), \quad (4) $$

where $W: x \mapsto W(x) \in S_+^d$ is a field taking value in the convex cone $S_+^d \subset \mathbb{R}^{d \times d}$ of symmetric positive semi-definite matrices. Again, $C(\mu, W)$ is the smallest constant such that (4) holds for any sufficiently smooth function $f$. The field $W$ provides a Riemannian metric on the space $\mathbb{R}^d$ that assigns to each $x \in \mathbb{R}^d$ an inner product $\langle u, v \rangle_x = u^\top W(x) v$, where $u, v \in \mathbb{R}^d$. The inequality (4) extends existing work in several ways. It encompasses the Brascamp-Lieb inequality [13, 18] that takes the same form as (4) but requires a strongly log-concave $\mu$ to have $W(x) = (\mathrm{Hess} V(x))^{-1}$ and $C(\mu, W) = 1$. It also generalizes the scalar-weighted Poincaré inequalities considered in [14, 16], in which $W(x) = w(x)I_d$ with a positive scalar field $w: x \to w(x) \in \mathbb{R}_+$. Such scalar-weighted Poincaré inequalities have been thoroughly analyzed in [15, 51, 52] in dimension $d = 1$.
The inequality (4) can be used in many ways. We motivate the work presented here by considering the case where the metric W is used to precondition the Langevin dynamic (2):

$$dX_t = (\mathrm{div}(W) - W\nabla V)dt + \sqrt{2W}dB_t, \quad (5)$$

where $\mathrm{div}(W) = (\sum_{j=1}^d \partial_t W_{i,j})_{1\le i \le d}$ is the divergence field of W and $\sqrt{W}$ is any square root of W such that $\sqrt{W}(x)\sqrt{W}(x)^\top = W(x)$. Following the discussion in [7, Section 1.15.2], such preconditioning strategy can be interpreted as a local and anisotropic time change of the SDE (2). It is local because $W(X_t)$ may vary in space and anisotropic as the matrix $W(X_t)$ may not be proportional to the identity matrix. As a result, the diffusion term $\sqrt{2W}dB_t$ can locally accelerate or decelerate the particle $X_t$ in different directions. The preconditioning in the form of (5) is referred to as the Riemannian Langevin dynamics, see [34, 38, 43] for applications. The following proposition, whose proof is given in Appendix A, establishes the exponential convergence rate of (5) toward the equilibrium measure $\mu$, provided $\mu$ satisfies the Riemannian Poincaré inequality (4).

**Proposition 1.** Let $\mu$ be a probability measure on $\mathbb{R}^d$ such that $d\mu(x) \propto \exp(-V(x))dx$ and let $W: \mathbb{R}^d \to S_+^d$. Then for any $C \ge 0$ the following assertions are equivalent.

1. $\mu$ satisfies the Riemannian Poincaré inequality (4) with $C(\mu, W) \le C$.

2. For any probability measure $\mu_0$, the solution $X_t \sim \mu_t$ to the Riemannian Langevin dynamics (5) with $X_0 \sim \mu_0$ satisfies

$$\chi^2(\mu_t, \mu) \le e^{-2t/C} \chi^2(\mu_0, \mu), \quad (6)$$

for all $t \ge 0$, where $\chi^2(\mu_t, \mu) = \mathrm{Var}_\mu(d\mu_t/d\mu)$ denotes the chi-square divergence.

In this setting, the Riemannian metric W offers great flexibility in improving the convergence rate of the SDE (5) by directly minimizing the Poincaré constant $C(\mu, W)$ over W. In this paper, we analyze the optimization problem

$$\min_{W: \mathbb{R}^d \to S_+^d} C(\mu, W) \qquad (7) \\ \int \mathrm{tr}(W) d\mu = \mathrm{tr}(\mathrm{Cov}_\mu)$$

for dimensions $d > 1$, where $\mathrm{tr}(\cdot)$ denotes the trace operator and $\mathrm{Cov}_\mu = \int (x-m)(x-m)^\top d\mu(x)$ and $m = \int x d\mu$ are respectively the covariance matrix and the mean of $\mu$, assuming they both exist. The normalization $\int \mathrm{tr}(W)d\mu = \mathrm{tr}(\mathrm{Cov}_\mu)$ in (7) removes the trivial invariance of the Poincaré constant $C(\mu, \alpha W) = C(\mu, W)/\alpha$ for any positive scalar $\alpha > 0$. In addition, this normalization is convenient because it yields<fnref n="1" /> the lower bound

$$C(\mu, W) \ge 1. \qquad (8)$$

Our contributions are the following. First, in Section 2, we use the notion of moment measure [27] to establish the existence of solutions to the optimization problem (7) with $C(\mu, W) = 1$. Then, in Section 3, we provide characterizations and properties of the optimal metric $W(x)$, showing that they are necessarily symmetric positive-semidefinite Stein kernels and relating them to the central limit theorem, see [21, 53]. Next, we show in Section 4 that (7) can be formulated as a concave optimization problem on the

<footnote n="1">¹Evaluating the Poincaré inequality (4) with the affine function $f(x) = x_i$ and summing over $i$ yields $\mathrm{tr}(\mathrm{Cov}_\mu) \le C(\mu, W) \mathrm{tr}(\int W d\mu)$ and then (8).</footnote>
The inequality (4) can be used in many ways. We motivate the work presented here by considering the case where the metric $W$ is used to precondition the Langevin dynamic (2):

$$dX_t = (\mathrm{div}(W) - W\nabla V)dt + \sqrt{2W}dB_t, \quad (5)$$

where $\mathrm{div}(W) = (\sum_{j=1}^d \partial_t W_{i,j})_{1\le i \le d}$ is the divergence field of $W$ and $\sqrt{W}$ is any square root of $W$ such that $\sqrt{W}(x)\sqrt{W}(x)^\top = W(x)$. Following the discussion in [7, Section 1.15.2], such preconditioning strategy can be interpreted as a local and anisotropic time change of the SDE (2). It is local because $W(X_t)$ may vary in space and anisotropic as the matrix $W(X_t)$ may not be proportional to the identity matrix. As a result, the diffusion term $\sqrt{2W}dB_t$ can locally accelerate or decelerate the particle $X_t$ in different directions. The preconditioning in the form of (5) is referred to as the Riemannian Langevin dynamics, see [34, 38, 43] for applications. The following proposition, whose proof is given in Appendix A, establishes the exponential convergence rate of (5) toward the equilibrium measure $\mu$, provided $\mu$ satisfies the Riemannian Poincaré inequality (4).

**Proposition 1.** Let $\mu$ be a probability measure on $\mathbb{R}^d$ such that $d\mu(x) \propto \exp(-V(x))dx$ and let $W: \mathbb{R}^d \to S_+^d$. Then for any $C \ge 0$ the following assertions are equivalent.

1. $\mu$ satisfies the Riemannian Poincaré inequality (4) with $C(\mu, W) \le C$.

2. For any probability measure $\mu_0$, the solution $X_t \sim \mu_t$ to the Riemannian Langevin dynamics (5) with $X_0 \sim \mu_0$ satisfies

$$\chi^2(\mu_t, \mu) \le e^{-2t/C} \chi^2(\mu_0, \mu), \quad (6)$$

for all $t \ge 0$, where $\chi^2(\mu_t, \mu) = \mathrm{Var}_\mu(d\mu_t/d\mu)$ denotes the chi-square divergence.

In this setting, the Riemannian metric $W$ offers great flexibility in improving the convergence rate of the SDE (5) by directly minimizing the Poincaré constant $C(\mu, W)$ over $W$. In this paper, we analyze the optimization problem

$$\min_{W: \mathbb{R}^d \to S_+^d} C(\mu, W) \qquad (7) \\ \int \mathrm{tr}(W) d\mu = \mathrm{tr}(\mathrm{Cov}_\mu)$$

for dimensions $d > 1$, where $\mathrm{tr}(\cdot)$ denotes the trace operator and $\mathrm{Cov}_\mu = \int (x-m)(x-m)^\top d\mu(x)$ and $m = \int x d\mu$ are respectively the covariance matrix and the mean of $\mu$, assuming they both exist. The normalization $\int \mathrm{tr}(W)d\mu = \mathrm{tr}(\mathrm{Cov}_\mu)$ in (7) removes the trivial invariance of the Poincaré constant $C(\mu, \alpha W) = C(\mu, W)/\alpha$ for any positive scalar $\alpha > 0$. In addition, this normalization is convenient because it yields<fnref n="1" /> the lower bound

$$C(\mu, W) \ge 1. \quad (8)$$

Our contributions are the following. First, in Section 2, we use the notion of moment measure [27] to establish the existence of solutions to the optimization problem (7) with $C(\mu, W) = 1$. Then, in Section 3, we provide characterizations and properties of the optimal metric $W(x)$, showing that they are necessarily symmetric positive-semidefinite Stein kernels and relating them to the central limit theorem, see [21, 53]. Next, we show in Section 4 that (7) can be formulated as a concave optimization problem on the

<footnote n="1">¹Evaluating the Poincaré inequality (4) with the affine function $f(x) = x_i$ and summing over $i$ yields $\mathrm{tr}(\mathrm{Cov}_\mu) \le C(\mu, W) \mathrm{tr}(\int W d\mu)$ and then (8).</footnote>
spectrum of the diffusion operator associated with the Riemannian Langevin dynamics (5). In Section 5, we propose a gradient-based algorithm to numerically solve (7) using the finite element method. In Section 6, we demonstrate the proposed algorithms on three benchmarks of dimension $d = 2$. In particular, we show the benefit of using the preconditioned Langevin dynamic (5) in the numerical discretization.

## 2 Existence of optimal metric

The existence of a solution to (7) can be established using moment measures [27, 50]. A probability measure $\mu$ is said to be the moment measure of a locally Lipschitz convex function $\varphi: \mathbb{R}^d \to \mathbb{R} \cup \{+\infty\}$ if $\mu$ is the pushforward measure of $e^{-\varphi(z)}dz$ under the transformation $z \mapsto \nabla\varphi(z)$, meaning

$$ \int f(x)d\mu(x) = \int f(\nabla\varphi(z))e^{-\varphi(z)}dz, \qquad (9) $$

for any $\mu$-integrable function $f$. The convexity of $\varphi$ ensures that the map $z \mapsto \nabla\varphi(z)$ is invertible and the measure $e^{-\varphi(z)}dz$ is log-concave. From the perspective of the optimal transport theory, the invertible map $\nabla\varphi$ is precisely the Brenier map that pushes forward the source measure $e^{-\varphi(z)}dz$ to $\mu$ and uniquely minimizes the quadratic cost $\int \|z - \nabla\varphi(z)\|^2 e^{-\varphi(z)}dz$.

Theorem 2 in [27] shows that any probability measure $\mu$ is a moment measure if it is not supported on a hyperplane, and satisfies $\int |x|d\mu < \infty$ and $\int xd\mu = 0$. Moreover, the corresponding moment map $\varphi$ is unique up to translation. The centred assumption on $\mu$ is a mild one: if $m = \int xd\mu \neq 0$, it is sufficient to replace $\mu$ with the translated measure $d\mu'(x) = d\mu(x - m)$ to satisfy this assumption.

**Theorem 2.1 (Existence of a solution to (7)).** Assume that $\mu$ is the moment measure of a convex function $\varphi$ which is twice continuously differentiable. Then the field of symmetric positive-semidefinite matrices

$$ W(x) = \mathrm{Hess} \varphi((\nabla\varphi)^{-1}(x)), \qquad (10) $$

is a solution to (7) with $C(\mu, W) = 1$.

*Proof.* We recall that the Brascamp-Lieb inequality [18] states that any log-concave measure $d\nu(z) = e^{-\varphi(z)}dz$ satisfies $\mathrm{Var}_{\nu}(f) \le \int \nabla f^{\top}(\mathrm{Hess} \varphi)^{-1} \nabla f d\nu$ for any sufficiently smooth function $f$. Following Section 4 in [31], we can write

$$ \begin{aligned} \mathrm{Var}_{\mu}(f) &= \mathrm{Var}_{e^{-\varphi}}(f \circ \nabla\varphi) \\ &\le \int (\nabla(f \circ \nabla\varphi))^{\top} (\mathrm{Hess} \varphi)^{-1} \nabla(f \circ \nabla\varphi) e^{-\varphi} dz \\ &= \int (\nabla f \circ \nabla\varphi)^{\top} (\mathrm{Hess} \varphi) (\nabla f \circ \nabla\varphi) e^{-\varphi} dz \\ &= \int \nabla f^{\top} (\mathrm{Hess} \varphi \circ (\nabla\varphi)^{-1}) \nabla f d\mu. \end{aligned} $$

Thus, $C(\mu, W) = 1$ for $W$ defined in (10). Next, we show that $W$ satisfies $\int \mathrm{tr}(W)d\mu = \mathrm{tr}(\mathrm{Cov}_{\mu})$. Theorem 2.3 in [31] ensures that $W$ defined in (10) is a *Stein kernel*, meaning that $\int xfd\mu = \int W\nabla f d\mu$ holds for any $f$. Letting $f$ going over all affine functions,
spectrum of the diffusion operator associated with the Riemannian Langevin dynamics (5). In Section 5, we propose a gradient-based algorithm to numerically solve (7) using the finite element method. In Section 6, we demonstrate the proposed algorithms on three benchmarks of dimension $d = 2$. In particular, we show the benefit of using the preconditioned Langevin dynamic (5) in the numerical discretization.

## 2 Existence of optimal metric

The existence of a solution to (7) can be established using moment measures [27, 50]. A probability measure $\mu$ is said to be the moment measure of a locally Lipschitz convex function $\varphi: \mathbb{R}^d \to \mathbb{R} \cup \{+\infty\}$ if $\mu$ is the pushforward measure of $e^{-\varphi(z)}dz$ under the transformation $z \mapsto \nabla\varphi(z)$, meaning

$$ \int f(x)d\mu(x) = \int f(\nabla\varphi(z))e^{-\varphi(z)}dz, \qquad (9) $$

for any $\mu$-integrable function $f$. The convexity of $\varphi$ ensures that the map $z \mapsto \nabla\varphi(z)$ is invertible and the measure $e^{-\varphi(z)}dz$ is log-concave. From the perspective of the optimal transport theory, the invertible map $\nabla\varphi$ is precisely the Brenier map that pushes forward the source measure $e^{-\varphi(z)}dz$ to $\mu$ and uniquely minimizes the quadratic cost $\int \|z - \nabla\varphi(z)\|^2e^{-\varphi(z)}dz$.

Theorem 2 in [27] shows that any probability measure $\mu$ is a moment measure if it is not supported on a hyperplane, and satisfies $\int |x|d\mu < \infty$ and $\int xd\mu = 0$. Moreover, the corresponding moment map $\varphi$ is unique up to translation. The centred assumption on $\mu$ is a mild one: if $m = \int xd\mu \neq 0$, it is sufficient to replace $\mu$ with the translated measure $d\mu'(x) = d\mu(x - m)$ to satisfy this assumption.

**Theorem 2.1 (Existence of a solution to (7)).** Assume that $\mu$ is the moment measure of a convex function $\varphi$ which is twice continuously differentiable. Then the field of symmetric positive-semidefinite matrices

$$ W(x) = \mathrm{Hess} \varphi((\nabla\varphi)^{-1}(x)), \qquad (10) $$

is a solution to (7) with $C(\mu, W) = 1$.

*Proof.* We recall that the Brascamp-Lieb inequality [18] states that any log-concave measure $d\nu(z) = e^{-\varphi(z)}dz$ satisfies $\mathrm{Var}_{\nu}(f) \le \int \nabla f^{\top}(\mathrm{Hess} \varphi)^{-1}\nabla f d\nu$ for any sufficiently smooth function $f$. Following Section 4 in [31], we can write

$$ \begin{aligned} \mathrm{Var}_{\mu}(f) &= \mathrm{Var}_{e^{-\varphi}}(f \circ \nabla\varphi) \\ &\le \int (\nabla(f \circ \nabla\varphi))^{\top} (\mathrm{Hess} \varphi)^{-1} \nabla(f \circ \nabla\varphi) e^{-\varphi} dz \\ &= \int (\nabla f \circ \nabla\varphi)^{\top} (\mathrm{Hess} \varphi) (\nabla f \circ \nabla\varphi) e^{-\varphi} dz \\ &= \int \nabla f^{\top} (\mathrm{Hess} \varphi \circ (\nabla\varphi)^{-1}) \nabla f d\mu. \end{aligned} $$

Thus, $C(\mu, W) = 1$ for $W$ defined in (10). Next, we show that $W$ satisfies $\int \mathrm{tr}(W)d\mu = \mathrm{tr}(\mathrm{Cov}_{\mu})$. Theorem 2.3 in [31] ensures that $W$ defined in (10) is a *Stein kernel*, meaning that $\int xfd\mu = \int W\nabla f d\mu$ holds for any $f$. Letting $f$ going over all affine functions,
we deduce that $\int W d\mu = \int xx^\top d\mu = \text{Cov}_\mu$, where we used the fact that, because $\mu$ is a moment measure, it is necessarily centered (see Proposition 1 in [27]). This yields $\int \text{tr}(W)d\mu = \text{tr}(\text{Cov}_\mu)$ so that, by (8) and $C(\mu, W) = 1$, we have that $W$ is a solution to (7). $\square$

Theorem 2.1 ensures the existence (but not the uniqueness) of solutions to the optimization problem (7) provided the moment map $\varphi$ is sufficiently smooth. This moment map is the solution to the toric Kähler-Einstein equation

$$e^{-\varphi} = \mu(\nabla\varphi) \det(\text{Hess } \varphi), \qquad (11)$$

which can be obtained by applying the change of variable $x = \nabla\varphi(z)$ to the left-hand side of (9). Note that (11) is a variant of the Monge-Ampère equation [32]. As shown in [11], if $\mu$ is centered, supported on a compact and convex set and has density bounded from above and below by positive constants, then $\varphi$ is smooth enough so that Theorem 2.1 applies.

In principle, Theorem 2.1 offers a constructive way to build an optimal metric $W$: one can first numerically solve (11) for $\varphi$ and then compute $W$ via (10). However, the numerical solution to (11) is not an easy task in general, mostly because of the convexity constraint of $\varphi$. Later on in Section 5, we propose an alternative way to compute $W$ which does not require solving (11) for $\varphi$. Nonetheless, in the particular case where $\mu$ is a product measure, we have a closed-form expression for the optimal $W$.

**Proposition 2.** Let $\mu = \mu_1 \otimes \dots \otimes \mu_L$ be a product probability measure on $\mathbb{R}^d$ where each $\mu_\ell$ is a marginal probability measure on $\mathbb{R}^{d_\ell}$ with $d = \sum_{\ell=1}^L d_\ell$. We denote $\tau_\ell = (1 + \sum_{i=1}^{\ell-1} d_i, \dots, \sum_{i=1}^\ell d_i)$ a multi-index with cardinality $d_\ell$ and $x_{\tau_\ell}$ the vector containing the coordinates of $x$ associated with $\mu_\ell$. Assume that, for any $1 \le \ell \le L$, there exists a symmetric positive-semidefinite matrix field $W_\ell: \mathbb{R}^{d_\ell} \to \mathcal{S}_+^{d_\ell}$ such that $C(\mu_\ell, W_\ell) = 1$ with $\int \text{tr}(W_\ell)d\mu_\ell = \text{tr}(\text{Cov}_{\mu_\ell})$. Then, the block-diagonal metric

$$W(x) = \begin{pmatrix} W_1(x_{\tau_1}) & 0 \\ \vdots & \vdots \\ 0 & W_L(x_{\tau_L}) \end{pmatrix}, \qquad (12)$$

is a solution to (7). For any $\tau_\ell$ with cardinality one, i.e., $d_\ell = 1$, suppose the corresponding $\mu_\ell$ has connected support on $\mathbb{R}$ with $\text{Cov}_{\mu_\ell} = \int (x_{\tau_\ell} - m_{\tau_\ell})^2 d\mu < \infty$, then we have the closed-form expression

$$W_\ell(x_{\tau_\ell}) = \frac{1}{\mu_{\tau_\ell}(x_{\tau_\ell})} \int_{x_{\tau_\ell}}^{\infty} (t - m_{\tau_\ell}) d\mu_\ell(t). \qquad (13)$$

*Proof.* Because $\int \text{tr}(W_\ell)d\mu_\ell = \text{tr}(\text{Cov}_{\mu_\ell})$, we have $\int \text{tr}(W)d\mu = \sum_{\ell=1}^L \text{tr}(\text{Cov}_{\mu_\ell}) = \text{tr}(\text{Cov}_\mu)$. By the stability of Poincaré inequalities under tensorization, e.g., see Proposition 4.3.1 in [7], we have $C(\mu, W) = \max_\ell\{C(\mu_\ell, W_\ell)\} = 1$, because $C(\mu_\ell, W_\ell) = 1$. This way, $W$ is a solution to (7). Finally, for the case $d_\ell = 1$, Theorem 6 in [51] ensures that $W_\ell$ defined in (13) satisfies $C(\mu_\ell, W_\ell) = 1$. It remains to check that $\int \text{tr}(W_\ell)d\mu_\ell = \text{tr}(\text{Cov}_{\mu_\ell})$. Applying integration by part, we have

$$\int W_\ell d\mu_\ell = \int_{-\infty}^{+\infty} \left( \int_{x_{\tau_\ell}}^{+\infty} (t - m_{\tau_\ell}) d\mu_\ell(t) \right) dx_{\tau_\ell} = \int_{-\infty}^{+\infty} (x_{\tau_\ell} - m_{\tau_\ell})^2 d\mu_\ell(x_{\tau_\ell}) = \text{Cov}_{\mu_\ell},$$

which concludes the proof. $\square$
we deduce that $\int W d\mu = \int xx^\top d\mu = \text{Cov}_\mu$, where we used the fact that, because $\mu$ is a moment measure, it is necessarily centered (see Proposition 1 in [27]). This yields $\int \text{tr}(W)d\mu = \text{tr}(\text{Cov}_\mu)$ so that, by (8) and $C(\mu, W) = 1$, we have that $W$ is a solution to (7). $\square$

Theorem 2.1 ensures the existence (but not the uniqueness) of solutions to the optimization problem (7) provided the moment map $\varphi$ is sufficiently smooth. This moment map is the solution to the toric Kähler-Einstein equation

$$e^{-\varphi} = \mu(\nabla\varphi) \det(\text{Hess } \varphi), \qquad (11)$$

which can be obtained by applying the change of variable $x = \nabla\varphi(z)$ to the left-hand side of (9). Note that (11) is a variant of the Monge-Ampère equation [32]. As shown in [11], if $\mu$ is centered, supported on a compact and convex set and has density bounded from above and below by positive constants, then $\varphi$ is smooth enough so that Theorem 2.1 applies.

In principle, Theorem 2.1 offers a constructive way to build an optimal metric $W$: one can first numerically solve (11) for $\varphi$ and then compute $W$ via (10). However, the numerical solution to (11) is not an easy task in general, mostly because of the convexity constraint of $\varphi$. Later on in Section 5, we propose an alternative way to compute $W$ which does not require solving (11) for $\varphi$. Nonetheless, in the particular case where $\mu$ is a product measure, we have a closed-form expression for the optimal $W$.

**Proposition 2.** Let $\mu = \mu_1 \otimes \dots \otimes \mu_L$ be a product probability measure on $\mathbb{R}^d$ where each $\mu_\ell$ is a marginal probability measure on $\mathbb{R}^{d_\ell}$ with $d = \sum_{\ell=1}^L d_\ell$. We denote $\tau_\ell = (1 + \sum_{i=1}^{\ell-1} d_i, \dots, \sum_{i=1}^\ell d_i)$ a multi-index with cardinality $d_\ell$ and $x_{\tau_\ell}$ the vector containing the coordinates of $x$ associated with $\mu_\ell$. Assume that, for any $1 \le \ell \le L$, there exists a symmetric positive-semidefinite matrix field $W_\ell: \mathbb{R}^{d_\ell} \to \mathcal{S}_+^{d_\ell}$ such that $C(\mu_\ell, W_\ell) = 1$ with $\int \text{tr}(W_\ell)d\mu_\ell = \text{tr}(\text{Cov}_{\mu_\ell})$. Then, the block-diagonal metric

$$W(x) = \begin{pmatrix} W_1(x_{\tau_1}) & 0 \\ \vdots & \vdots \\ 0 & W_L(x_{\tau_L}) \end{pmatrix}, \qquad (12)$$

is a solution to (7). For any $\tau_\ell$ with cardinality one, i.e., $d_\ell = 1$, suppose the corresponding $\mu_\ell$ has connected support on $\mathbb{R}$ with $\text{Cov}_{\mu_\ell} = \int (x_{\tau_\ell} - m_{\tau_\ell})^2 d\mu < \infty$, then we have the closed-form expression

$$W_\ell(x_{\tau_\ell}) = \frac{1}{\mu_{\tau_\ell}(x_{\tau_\ell})} \int_{x_{\tau_\ell}}^{\infty} (t - m_{\tau_\ell}) d\mu_\ell(t). \qquad (13)$$

*Proof.* Because $\int \text{tr}(W_\ell)d\mu_\ell = \text{tr}(\text{Cov}_{\mu_\ell})$, we have $\int \text{tr}(W)d\mu = \sum_{\ell=1}^L \text{tr}(\text{Cov}_{\mu_\ell}) = \text{tr}(\text{Cov}_\mu)$. By the stability of Poincaré inequalities under tensorization, e.g., see Proposition 4.3.1 in [7], we have $C(\mu, W) = \max_\ell\{C(\mu_\ell, W_\ell)\} = 1$, because $C(\mu_\ell, W_\ell) = 1$. This way, $W$ is a solution to (7). Finally, for the case $d_\ell = 1$, Theorem 6 in [51] ensures that $W_\ell$ defined in (13) satisfies $C(\mu_\ell, W_\ell) = 1$. It remains to check that $\int \text{tr}(W_\ell)d\mu_\ell = \text{tr}(\text{Cov}_{\mu_\ell})$. Applying integration by part, we have

$$\int W_\ell d\mu_\ell = \int_{-\infty}^{+\infty} \left( \int_{x_{\tau_\ell}}^{+\infty} (t - m_{\tau_\ell}) d\mu_\ell(t) \right) dx_{\tau_\ell} = \int_{-\infty}^{+\infty} (x_{\tau_\ell} - m_{\tau_\ell})^2 d\mu_\ell(x_{\tau_\ell}) = \text{Cov}_{\mu_\ell},$$

which concludes the proof. $\square$
We give a few examples of measures where (12) and (13) have closed-form expressions. We refer to [51, Section 3.1] for in-depth discussion.

**Example 1.** For the standard Gaussian measure $d\mu_N(x) \propto \prod_{i=1}^d \exp(-|x_i|^2/2)dx$ on $\mathbb{R}^d$, we have

$$W_N(x) = \begin{pmatrix} 1 & & & 0 \\ & \ddots & & \\ 0 & & & 1 \end{pmatrix}.$$

For the Laplace measure $d\mu_L(x) \propto \prod_{i=1}^d \exp(-|x_i|)dx$ on $\mathbb{R}^d$, we have

$$W_L(x) = \begin{pmatrix} 1 + |x_1| & & & 0 \\ & \ddots & & \\ 0 & & 1 + |x_d| & \end{pmatrix}.$$

The generalized Cauchy distribution $d\mu_C \propto \prod_{i=1}^d (1 + |x_i|^2)^{-\beta}dx$ on $\mathbb{R}^d$ with parameter $\beta > 1$ yields

$$W_C(x) = \frac{1}{2(\beta - 1)} \begin{pmatrix} 1 + |x_1|^2 & & & 0 \\ & \ddots & & \\ 0 & & 1 + |x_d|^2 & \end{pmatrix}.$$

The following corollary characterizes the optimal metric for sums of independent random variables.

**Corollary 1.** Let $X, Y$ taking values in $\mathbb{R}^d$ be two independent random vectors and consider $Z = \alpha X + \beta Y$ for some $\alpha, \beta \in \mathbb{R}$. Assume there exists optimal metrics $W_X$ and $W_Y$ such that $C(\mu_X, W_X) = C(\mu_Y, W_Y) = 1$, as well as $\int \mathrm{tr}(W_X)d\mu_X = \mathrm{tr}(\mathrm{Cov}(X))$ and $\int \mathrm{tr}(W_Y)d\mu_Y = \mathrm{tr}(\mathrm{Cov}(Y))$. Then, the conditional expectation

$$W_Z(z) = \mathbb{E}[\alpha^2 W_X(X) + \beta^2 W_Y(Y) | Z = z],$$

is also an optimal metric for $Z$. In particular, for $\bar{X}_N = \frac{1}{\sqrt{N}} \sum_{i=1}^N X_i$ where $X_1, \dots, X_N$ are $N$ independent copies of $X$, we have

$$W_{\bar{X}_N}(x) = \mathbb{E} \left[ \frac{W_X(X_1) + \dots + W_X(X_N)}{N} \middle| \bar{X}_N = x \right]. \quad (14)$$

*Proof.* Let $\mu = \mu_X \otimes \mu_Y$ be the measure of the random vector $(X, Y)$. By Proposition 2 we have an optimal metric $\mathrm{diag}(W_X(X), W_Y(Y))$ for $(X, Y)$. Thus, applying the Riemannian Poincaré inequality (4) with $f(x, y) = h(\alpha x + \beta y)$ yields

$$\begin{align*} \mathrm{Var}(h(Z)) &= \mathrm{Var}(f(X, Y)) \\ &\le \mathbb{E}[\nabla f(X, Y)^{\top} \mathrm{diag}(W_X(X), W_Y(Y)) \nabla f(X, Y)] \\ &= \mathbb{E}[\alpha^2 \nabla h(Z)^{\top} W_X(X) \nabla h(Z)^{\top} + \beta^2 \nabla h(Z)^{\top} W_Y(Y) \nabla h(Z)^{\top}] \\ &= \mathbb{E}[\nabla h(Z)^{\top} (\alpha^2 W_X(X) + \beta^2 W_Y(Y)) \nabla h(Z)^{\top}]. \end{align*}$$

We conclude that $C(\mu_Z, W_Z) = 1$. Furthermore, using the independence of $X$ and $Y$, we have $\int \mathrm{tr}(W_Z)d\mu_Z = \alpha^2 \mathrm{tr} \mathrm{Cov}(X) + \beta^2 \mathrm{Cov}(Y) = \mathrm{Cov}(Z)$. We conclude that $W_Z$ is an optimal metric for $Z$. Then, (14) follows by recursion and concludes the proof. $\square$
We give a few examples of measures where (12) and (13) have closed-form expressions. We refer to [51, Section 3.1] for in-depth discussion.

**Example 1.** For the standard Gaussian measure $d\mu_N(x) \propto \prod_{i=1}^d \exp(-|x_i|^2/2)dx$ on $\mathbb{R}^d$, we have

$$W_N(x) = \begin{pmatrix} 1 & & & 0 \\ & \ddots & & \\ & & \ddots & \\ 0 & & & 1 \end{pmatrix}.$$

For the Laplace measure $d\mu_L(x) \propto \prod_{i=1}^d \exp(-|x_i|)dx$ on $\mathbb{R}^d$, we have

$$W_L(x) = \begin{pmatrix} 1 + |x_1| & & & 0 \\ & \ddots & & \\ & & \ddots & \\ 0 & & 1 + |x_d| \end{pmatrix}.$$

The generalized Cauchy distribution $d\mu_C \propto \prod_{i=1}^d (1 + |x_i|^2)^{-\beta}dx$ on $\mathbb{R}^d$ with parameter $\beta > 1$ yields

$$W_C(x) = \frac{1}{2(\beta - 1)} \begin{pmatrix} 1 + |x_1|^2 & & & 0 \\ & \ddots & & \\ & & \ddots & \\ 0 & & 1 + |x_d|^2 \end{pmatrix}.$$

The following corollary characterizes the optimal metric for sums of independent random variables.

**Corollary 1.** Let $X, Y$ taking values in $\mathbb{R}^d$ be two independent random vectors and consider $Z = \alpha X + \beta Y$ for some $\alpha, \beta \in \mathbb{R}$. Assume there exists optimal metrics $W_X$ and $W_Y$ such that $C(\mu_X, W_X) = C(\mu_Y, W_Y) = 1$, as well as $\int \text{tr}(W_X)d\mu_X = \text{tr}(\text{Cov}(X))$ and $\int \text{tr}(W_Y)d\mu_Y = \text{tr}(\text{Cov}(Y))$. Then, the conditional expectation

$$W_Z(z) = \mathbb{E} [\alpha^2 W_X(X) + \beta^2 W_Y(Y) | Z = z],$$

is also an optimal metric for $Z$. In particular, for $\bar{X}_N = \frac{1}{\sqrt{N}} \sum_{i=1}^N X_i$ where $X_1, \dots, X_N$ are $N$ independent copies of $X$, we have

$$W_{\bar{X}_N}(x) = \mathbb{E} \left[ \frac{W_X(X_1) + \dots + W_X(X_N)}{N} \middle| \bar{X}_N = x \right]. \quad (14)$$

*Proof.* Let $\mu = \mu_X \otimes \mu_Y$ be the measure of the random vector $(X, Y)$. By Proposition 2 we have an optimal metric $\text{diag}(W_X(X), W_Y(Y))$ for $(X, Y)$. Thus, applying the Riemannian Poincaré inequality (4) with $f(x, y) = h(\alpha x + \beta y)$ yields

$$\begin{align*}
\text{Var}(h(Z)) &= \text{Var}(f(X, Y)) \\
&\le \mathbb{E} [\nabla f(X, Y)^\top \text{diag}(W_X(X), W_Y(Y)) \nabla f(X, Y)] \\
&= \mathbb{E} [\alpha^2 \nabla h(Z)^\top W_X(X) \nabla h(Z)^\top + \beta^2 \nabla h(Z)^\top W_Y(Y) \nabla h(Z)^\top] \\
&= \mathbb{E} [\nabla h(Z)^\top (\alpha^2 W_X(X) + \beta^2 W_Y(Y)) \nabla h(Z)^\top].
\end{align*}$$

We conclude that $C(\mu_Z, W_Z) = 1$. Furthermore, using the independence of $X$ and $Y$, we have $\int \text{tr}(W_Z)d\mu_Z = \alpha^2 \text{tr} \text{Cov}(X) + \beta^2 \text{Cov}(Y) = \text{Cov}(Z)$. We conclude that $W_Z$ is an optimal metric for $Z$. Then, (14) follows by recursion and concludes the proof. $\square$
# 3 Characterization and properties of optimal metrics

We start with a simple characterization of the solutions to the optimization problem (7) with $C(\mu, W) = 1$.

**Theorem 3.1.** Let $W: \mathbb{R}^d \to S_+^d$ satisfy $\int \mathrm{tr}(W)d\mu = \mathrm{tr}(\mathrm{Cov}_\mu)$. Then, $C(\mu, W) = 1$ if and only if the Poincaré inequality (4) becomes an equality for all affine functions, that is

$$ \mathrm{Var}_\mu(f) = C(\mu, W) \int \nabla f^\top W \nabla f d\mu, \quad \forall f: \text{affine function.} \quad (15) $$

In this case we have

$$ \int W d\mu = \mathrm{Cov}_\mu. \quad (16) $$

*Proof.* First, we assume that the identity in (15) holds. Choosing $f(x) = x_i$ in (15) and summing over $1 \le i \le d$ yield $\mathrm{tr}(\mathrm{Cov}_\mu) = C(\mu, W) \mathrm{tr}(\int W d\mu)$, and thus $C(\mu, W) = 1$. Next, we assume $C(\mu, W) = 1$ and show (15) by contradiction. Assume there exists an affine function $f(x) = \alpha^\top x + \beta$ with $\alpha \in \mathbb{R}^d, \beta \in \mathbb{R}$ such that $\mathrm{Var}_\mu(f) < \int \nabla f^\top W \nabla f d\mu$, which is equivalent to $\alpha^\top \mathrm{Cov}_\mu \alpha < \alpha^\top (\int W d\mu) \alpha$. Define $u_1 = \alpha/\|\alpha\|$ and choose vectors $u_2, \dots, u_d$ to form a unitary matrix $U = [u_1, \dots, u_d] \in \mathbb{R}^{d \times d}$. We have

$$ \mathrm{tr}(\mathrm{Cov}_\mu) = \sum_{i=1}^{d} u_i^\top \mathrm{Cov}_\mu u_i < \sum_{i=1}^{d} u_i^\top \left( \int W d\mu \right) u_i = \mathrm{tr} \left( \int W d\mu \right), $$

which contradicts $\int \mathrm{tr}(W)d\mu = \mathrm{tr}(\mathrm{Cov}_\mu)$. Finally, the identity in (16) is a direct consequence of (15) with $C(\mu, W) = 1$, which concludes the proof. □

Next, we show that any solutions to (7) is necessarily a *Stein kernel*. A *Stein kernel* of a probability measure $\mu$ is a matrix field $W: \mathbb{R}^d \to \mathbb{R}^{d \times d}$ (not necessarily symmetric or positive-semidefinite) such that

$$ \int (x-m)fd\mu = \int W\nabla f d\mu, \quad (17) $$

for any smooth function $f: \mathbb{R}^d \to \mathbb{R}$.

**Theorem 3.2.** If a symmetric positive-semidefinite matrix field $W$ is a solution to the optimization problem (7) with $C(\mu, W) = 1$, then $W$ is a *Stein kernel* of $\mu$.

*Proof.* By Theorem 3.1 we have $\mathrm{Cov}_\mu = \int W d\mu$. Then, for any $f: \mathbb{R}^d \to \mathbb{R}$, $\alpha \in \mathbb{R}^d$ and $h \in \mathbb{R}$, the function $g = \alpha^\top x + hf$ satisfies

$$ \mathrm{Var}_\mu(g) = \alpha^\top \mathrm{Cov}_\mu \alpha + 2h\alpha^\top \int (x-m)fd\mu + h^2 \mathrm{Var}_\mu(f), \\ \int \nabla g^\top W \nabla g d\mu = \alpha^\top \mathrm{Cov}_\mu \alpha + 2h\alpha^\top \int W \nabla f d\mu + h^2 \int \nabla f^\top W \nabla f d\mu. $$

Applying the Riemannian Poincaré inequality (4) to $g$ yields

$$ 0 \le h\alpha^\top \left( \int W \nabla f d\mu - \int (x-m) f d\mu \right) + \mathcal{O}(h^2). $$

For this to hold for any $h \in \mathbb{R}$ and $\alpha \in \mathbb{R}^d$, it is necessary that $\int W \nabla f d\mu - \int (x-m) f d\mu = 0$, which defines $W$ as a *Stein kernel*. □
# 3 Characterization and properties of optimal metrics

We start with a simple characterization of the solutions to the optimization problem (7) with $C(\mu, W) = 1$.

**Theorem 3.1.** Let $W: \mathbb{R}^d \to S_+^d$ satisfy $\int \mathrm{tr}(W)d\mu = \mathrm{tr}(\mathrm{Cov}_\mu)$. Then, $C(\mu, W) = 1$ if and only if the Poincaré inequality (4) becomes an equality for all affine functions, that is

$$ \mathrm{Var}_\mu(f) = C(\mu, W) \int \nabla f^\top W \nabla f d\mu, \quad \forall f: \text{affine function.} \quad (15) $$

In this case we have

$$ \int W d\mu = \mathrm{Cov}_\mu. \quad (16) $$

*Proof.* First, we assume that the identity in (15) holds. Choosing $f(x) = x_i$ in (15) and summing over $1 \le i \le d$ yield $\mathrm{tr}(\mathrm{Cov}_\mu) = C(\mu, W) \mathrm{tr}(\int W d\mu)$, and thus $C(\mu, W) = 1$. Next, we assume $C(\mu, W) = 1$ and show (15) by contradiction. Assume there exists an affine function $f(x) = \alpha^\top x + \beta$ with $\alpha \in \mathbb{R}^d, \beta \in \mathbb{R}$ such that $\mathrm{Var}_\mu(f) < \int \nabla f^\top W \nabla f d\mu$, which is equivalent to $\alpha^\top \mathrm{Cov}_\mu \alpha < \alpha^\top (\int W d\mu) \alpha$. Define $u_1 = \alpha/\|\alpha\|$ and choose vectors $u_2, \dots, u_d$ to form a unitary matrix $U = [u_1, \dots, u_d] \in \mathbb{R}^{d \times d}$. We have

$$ \mathrm{tr}(\mathrm{Cov}_\mu) = \sum_{i=1}^{d} u_i^\top \mathrm{Cov}_\mu u_i < \sum_{i=1}^{d} u_i^\top \left( \int W d\mu \right) u_i = \mathrm{tr} \left( \int W d\mu \right), $$

which contradicts $\int \mathrm{tr}(W)d\mu = \mathrm{tr}(\mathrm{Cov}_\mu)$. Finally, the identity in (16) is a direct consequence of (15) with $C(\mu, W) = 1$, which concludes the proof. □

Next, we show that any solutions to (7) is necessarily a *Stein kernel*. A Stein kernel of a probability measure $\mu$ is a matrix field $W: \mathbb{R}^d \to \mathbb{R}^{d \times d}$ (not necessarily symmetric or positive-semidefinite) such that

$$ \int (x-m)fd\mu = \int W\nabla f d\mu, \quad (17) $$

for any smooth function $f: \mathbb{R}^d \to \mathbb{R}$.

**Theorem 3.2.** If a symmetric positive-semidefinite matrix field $W$ is a solution to the optimization problem (7) with $C(\mu, W) = 1$, then $W$ is a Stein kernel of $\mu$.

*Proof.* By Theorem 3.1 we have $\mathrm{Cov}_\mu = \int W d\mu$. Then, for any $f: \mathbb{R}^d \to \mathbb{R}, \alpha \in \mathbb{R}^d$ and $h \in \mathbb{R}$, the function $g = \alpha^\top x + hf$ satisfies

$$ \mathrm{Var}_\mu(g) = \alpha^\top \mathrm{Cov}_\mu \alpha + 2h\alpha^\top \int (x-m)fd\mu + h^2 \mathrm{Var}_\mu(f), \\ \int \nabla g^\top W \nabla g d\mu = \alpha^\top \mathrm{Cov}_\mu \alpha + 2h\alpha^\top \int W \nabla f d\mu + h^2 \int \nabla f^\top W \nabla f d\mu. $$

Applying the Riemannian Poincaré inequality (4) to $g$ yields

$$ 0 \le h\alpha^\top \left( \int W \nabla f d\mu - \int (x-m) f d\mu \right) + \mathcal{O}(h^2). $$

For this to hold for any $h \in \mathbb{R}$ and $\alpha \in \mathbb{R}^d$, it is necessary that $\int W \nabla f d\mu - \int (x-m) f d\mu = 0$, which defines $W$ as a Stein kernel. □
Since the seminal work [53], Stein kernels have been implicitly used in numerous works on Stein's method, and they have received growing attention in several recent investigations, see [25, 28, 30] for instance. In particular, it is shown in [31, Proposition 3.1] that the so-called *Stein discrepancy*

$$S_p(\mu|\gamma) = \inf_{W: \text{Stein kernel as in (17)}} \left( \int \|W - I_d\|_F^p d\mu \right)^{1/p},$$

measures the divergence from $\mu$ to the Gaussian measure $\gamma = \mathcal{N}(0, I_d)$, which remarkably bounds the Wasserstein-$p$ distance between $\mu$ and $\gamma$ for $p \ge 2$, see also [36] for the original proof for $p = 2$. A notable application of Stein kernels is that they often lead to non-asymptotic central limit theorems (CLT) [8, 23]. The following corollary provides such a CLT based on the optimal metric. The proof is given in Appendix B.

**Corollary 2.** Let $X_1, \dots, X_N$ be independent copies of $X \sim \mu$ such that $\mathbb{E}[X] = 0$ and $\text{Cov}(X) = I_d$, and let $\mu_N$ be the measure of $\bar{X}_N = \frac{1}{\sqrt{N}} \sum X_i$. Consider the Wasserstein-$p$ distance between $\mu_N$ and $\gamma = \mathcal{N}(0, I_d)$ given by

$$\mathcal{W}_p(\mu_N, \gamma) = \left( \inf_{\pi \in \Pi(\mu_N, \gamma)} \int \|x - y\|^p \pi(dx, dy) \right)^{1/p},$$

where $\Pi(\mu_N, \gamma)$ denotes the set of couplings between $\mu_N$ and $\gamma$. Then, for any $p \ge 2$, we have

$$\mathcal{W}_p(\mu_N, \gamma) \le \frac{C_p d^{1-2/p}}{\sqrt{N}} \left( \int \|W(X) - I_d\|_p^p d\mu \right)^{1/p},$$

where $W$ is a solution to (7) with $C(\mu, W) = 1$ and $C_p$ is a constant depending only on $p$ and where $\| \cdot \|_p$ is the matrix norm defined by $\|A\|_p^p = \sum_{ij} A_{ij}^p$.

Compared to Theorem 3.3 in [31], the main improvement of Corollary 2 is that it does not require $\mu$ to be log-concave. It applies to more general probability measures, for example, those heavy-tail distributions in Example 1, as long as the optimal metric $W$ for the underlying moment measure $\mu$ satisfies $\mathbb{E}[\|W(X) - I_d\|_F^p] < \infty$ for the chosen $p$.

Another immediate property of Stein kernels is that they permit to derive a lower bound on the variance of a function. The following corollary is a generalization of the results of [20] for dimension $d > 1$.

**Corollary 3.** Let a symmetric positive-semidefinite matrix field $W$ be a solution to (7) with $C(\mu, W) = 1$. Then for any smooth function $f: \mathbb{R}^d \to \mathbb{R}$ we have

$$\left\| \int W \nabla f d\mu \right\|_{\text{Cov}_{\mu}^{-1}}^2 \le \text{Var}_{\mu}(f) \le \int \nabla f^{\top} W \nabla f d\mu, \quad (18)$$

where $\|v\|_{\text{Cov}_{\mu}^{-1}}^2 = v^{\top} \text{Cov}_{\mu}^{-1} v$.

*Proof.* By Theorem 3.2, $W$ is a Stein kernel of $\mu$. The Cauchy-Schwarz inequality yields

$$\left| \alpha^{\top} \int W \nabla f d\mu \right|^{(17)} \stackrel{(17)}{=} \left| \alpha^{\top} \int (x-m) \left( f - \int f d\mu \right) d\mu \right|^2 \le (\alpha^{\top} \text{Cov}_{\mu} \alpha) \text{Var}_{\mu}(f),$$

for any $\alpha \in \mathbb{R}^d$. Taking the supremum over $\alpha \in \mathbb{R}^d$ such that $\alpha^{\top} \text{Cov}_{\mu} \alpha = 1$ yields the left-hand side of (18). The right-hand side of (18) is just the Riemmanian Poincaré inequality (4) with $C(\mu, W) = 1$. □
Since the seminal work [53], Stein kernels have been implicitly used in numerous works on Stein's method, and they have received growing attention in several recent investigations, see [25, 28, 30] for instance. In particular, it is shown in [31, Proposition 3.1] that the so-called *Stein discrepancy*

$$S_p(\mu|\gamma) = \inf_{W: \text{Stein kernel as in (17)}} \left( \int \|W - I_d\|_F^p d\mu \right)^{1/p},$$

measures the divergence from $\mu$ to the Gaussian measure $\gamma = \mathcal{N}(0, I_d)$, which remarkably bounds the Wasserstein-$p$ distance between $\mu$ and $\gamma$ for $p \ge 2$, see also [36] for the original proof for $p = 2$. A notable application of Stein kernels is that they often lead to non-asymptotic central limit theorems (CLT) [8, 23]. The following corollary provides such a CLT based on the optimal metric. The proof is given in Appendix B.

**Corollary 2.** Let $X_1, \dots, X_N$ be independent copies of $X \sim \mu$ such that $\mathbb{E}[X] = 0$ and $\text{Cov}(X) = I_d$, and let $\mu_N$ be the measure of $\bar{X}_N = \frac{1}{\sqrt{N}} \sum X_i$. Consider the Wasserstein-$p$ distance between $\mu_N$ and $\gamma = \mathcal{N}(0, I_d)$ given by

$$\mathcal{W}_p(\mu_N, \gamma) = \left( \inf_{\pi \in \Pi(\mu_N, \gamma)} \int \|x - y\|^p \pi(dx, dy) \right)^{1/p},$$

where $\Pi(\mu_N, \gamma)$ denotes the set of couplings between $\mu_N$ and $\gamma$. Then, for any $p \ge 2$, we have

$$\mathcal{W}_p(\mu_N, \gamma) \le \frac{C_p d^{1-2/p}}{\sqrt{N}} \left( \int \|W(X) - I_d\|_p^p d\mu \right)^{1/p},$$

where $W$ is a solution to (7) with $C(\mu, W) = 1$ and $C_p$ is a constant depending only on $p$ and where $\| \cdot \|_p$ is the matrix norm defined by $\|A\|_p^p = \sum_{ij} A_{ij}^p$.

Compared to Theorem 3.3 in [31], the main improvement of Corollary 2 is that it does not require $\mu$ to be log-concave. It applies to more general probability measures, for example, those heavy-tail distributions in Example 1, as long as the optimal metric $W$ for the underlying moment measure $\mu$ satisfies $\mathbb{E}[\|W(X) - I_d\|_F^p] < \infty$ for the chosen $p$.

Another immediate property of Stein kernels is that they permit to derive a lower bound on the variance of a function. The following corollary is a generalization of the results of [20] for dimension $d > 1$.

**Corollary 3.** Let a symmetric positive-semidefinite matrix field $W$ be a solution to (7) with $C(\mu, W) = 1$. Then for any smooth function $f: \mathbb{R}^d \to \mathbb{R}$ we have

$$\left\| \int W \nabla f d\mu \right\|_{\text{Cov}_{\mu}^{-1}}^2 \le \text{Var}_{\mu}(f) \le \int \nabla f^{\top} W \nabla f d\mu, \quad (18)$$

where $\|v\|_{\text{Cov}_{\mu}^{-1}}^2 = v^{\top} \text{Cov}_{\mu}^{-1} v$.

*Proof.* By Theorem 3.2, $W$ is a Stein kernel of $\mu$. The Cauchy-Schwarz inequality yields

$$\left| \alpha^{\top} \int W \nabla f d\mu \right|^{(17)} \stackrel{(17)}{=} \left| \alpha^{\top} \int (x-m) \left( f - \int f d\mu \right) d\mu \right|^2 \le (\alpha^{\top} \text{Cov}_{\mu} \alpha) \text{Var}_{\mu}(f),$$

for any $\alpha \in \mathbb{R}^d$. Taking the supremum over $\alpha \in \mathbb{R}^d$ such that $\alpha^{\top} \text{Cov}_{\mu} \alpha = 1$ yields the left-hand side of (18). The right-hand side of (18) is just the Riemmanian Poincaré inequality (4) with $C(\mu, W) = 1$. □
We conclude this section with an excursion. Let $f$ be a 1-Lipschitz function so that $\|\nabla f(x)\| \le 1$. The right inequality of (18) yields $\mathrm{Var}_\mu(f) \le \int \lambda_{\max}(W) \|\nabla f\|^2 d\mu \le \int \lambda_{\max}(W) d\mu$. Maximizing over the set of 1-Lipschitz function then yields

$$ \sup_{f: \mathbb{R}^d \to \mathbb{R} \atop \mathrm{1-Lipschitz}} \mathrm{Var}_\mu(f) \le \int \lambda_{\max}(W) d\mu. \quad (19) $$

This gives an upper bound on the maximal variance of 1-Lipschitz function, which is known to be a difficult question related to the isoperimetric problem [12].

# 4 Convex analysis and spectral gap

While Theorem 3.1 gives a necessary and sufficient condition for $C(\mu, W) = 1$, it does not provide any constructive way to build such an optimal metric $W$. In the following, we discuss how to numerically obtain $W$ using optimization methods. Toward this goal, we introduce the following Sobolev spaces

$$ \begin{aligned} L^2(\mu) &= \left\{ f : \mathbb{R}^d \to \mathbb{R} : \|f\|_{L^2(\mu)}^2 = \int f^2 d\mu < \infty \right\}, \\ H^1(\mu, W) &= \left\{ f : \mathbb{R}^d \to \mathbb{R} : \|f\|_{H^1(\mu, W)}^2 = \int f^2 d\mu + \int \nabla f^\top W \nabla f d\mu < \infty \right\}. \end{aligned} $$

The space $H^1(\mu, W)$ corresponds to the largest function space for $f$ for which the terms $\mathrm{Var}_\mu(f)$ and $\int \nabla f^\top W \nabla f d\mu$ of (4) are both finite. The Poincaré constant $C(\mu, W)$ can thus be expressed as

$$ \frac{1}{C(\mu, W)} = \inf_{\substack{f \in H^1(\mu, W) \\ \int f d\mu = 0}} \frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu}. \quad (20) $$

The next proposition shows that the variational problem (7) can be solved by a convex optimization problem.

**Proposition 3.** The function $W \mapsto C(\mu, W)^{-1}$ is concave over the convex cone of symmetric positive-semidefinite matrix field $W$. Furthermore, for any symmetric positive-semidefinite matrix fields $W$ and $W'$ such that $W' \le \Omega W$ for some constant $\Omega$, we have

$$ \frac{1}{C(\mu, W')} \le \frac{1}{C(\mu, W)} + \int \mathrm{tr}((W' - W)G) d\mu, \quad (21) $$

for any $G \in \mathcal{G}(W)$ where

$$ \mathcal{G}(W) = \mathrm{conv} \left\{ \frac{\nabla u \nabla u^\top}{\int u^2 d\mu} : u \in \arg \min_{\substack{f \in H^1(\mu, W) \\ \int f d\mu = 0}} \frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu} \right\}. \quad (22) $$

Here $\mathrm{conv}\{\cdot\}$ denotes the convex hull of a set.

*Proof.* For any symmetric positive-semidefinite matrix field $W_0, W_1$ and any $0 < t < 1$, we have that $W_t = tW_1 + (1-t)W_0$ is a symmetric positive-semidefinite matrix field. By construction, for any $f \in H^1(\mu, W_t)$ we have

$$ \|f\|_{H^1(\mu, W_t)}^2 = t \|f\|_{H^1(\mu, W_1)}^2 + (1-t) \|f\|_{H^1(\mu, W_0)}^2, $$
We conclude this section with an excursion. Let $f$ be a 1-Lipschitz function so that $\|\nabla f(x)\| \le 1$. The right inequality of (18) yields $\mathrm{Var}_\mu(f) \le \int \lambda_{\max}(W) \|\nabla f\|^2 d\mu \le \int \lambda_{\max}(W) d\mu$. Maximizing over the set of 1-Lipschitz function then yields

$$ \sup_{f:\mathbb{R}^d \to \mathbb{R} \atop \mathrm{1-Lipschitz}} \mathrm{Var}_\mu(f) \le \int \lambda_{\max}(W) d\mu. \quad (19) $$

This gives an upper bound on the maximal variance of 1-Lipschitz function, which is known to be a difficult question related to the isoperimetric problem [12].

# 4 Convex analysis and spectral gap

While Theorem 3.1 gives a necessary and sufficient condition for $C(\mu, W) = 1$, it does not provide any constructive way to build such an optimal metric $W$. In the following, we discuss how to numerically obtain $W$ using optimization methods. Toward this goal, we introduce the following Sobolev spaces

$$ \begin{aligned} L^2(\mu) &= \left\{ f : \mathbb{R}^d \to \mathbb{R} : \|f\|_{L^2(\mu)}^2 = \int f^2 d\mu < \infty \right\}, \\ H^1(\mu, W) &= \left\{ f : \mathbb{R}^d \to \mathbb{R} : \|f\|_{H^1(\mu, W)}^2 = \int f^2 d\mu + \int \nabla f^\top W \nabla f d\mu < \infty \right\}. \end{aligned} $$

The space $H^1(\mu, W)$ corresponds to the largest function space for $f$ for which the terms $\mathrm{Var}_\mu(f)$ and $\int \nabla f^\top W \nabla f d\mu$ of (4) are both finite. The Poincaré constant $C(\mu, W)$ can thus be expressed as

$$ \frac{1}{C(\mu, W)} = \inf_{\substack{f \in H^1(\mu, W) \\ \int f d\mu = 0}} \frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu}. \quad (20) $$

The next proposition shows that the variational problem (7) can be solved by a convex optimization problem.

**Proposition 3.** The function $W \mapsto C(\mu, W)^{-1}$ is concave over the convex cone of symmetric positive-semidefinite matrix field $W$. Furthermore, for any symmetric positive-semidefinite matrix fields $W$ and $W'$ such that $W' \le \Omega W$ for some constant $\Omega$, we have

$$ \frac{1}{C(\mu, W')} \le \frac{1}{C(\mu, W)} + \int \mathrm{tr}((W' - W)G) d\mu, \quad (21) $$

for any $G \in \mathcal{G}(W)$ where

$$ \mathcal{G}(W) = \mathrm{conv} \left\{ \frac{\nabla u \nabla u^\top}{\int u^2 d\mu} : u \in \arg \min_{\substack{f \in H^1(\mu, W) \\ \int f d\mu = 0}} \frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu} \right\}. \quad (22) $$

Here $\mathrm{conv}\{\cdot\}$ denotes the convex hull of a set.

*Proof.* For any symmetric positive-semidefinite matrix field $W_0, W_1$ and any $0 < t < 1$, we have that $W_t = tW_1 + (1-t)W_0$ is a symmetric positive-semidefinite matrix field. By construction, for any $f \in H^1(\mu, W_t)$ we have

$$ \|f\|_{H^1(\mu, W_t)}^2 = t \|f\|_{H^1(\mu, W_1)}^2 + (1-t) \|f\|_{H^1(\mu, W_0)}^2, $$
so that $\|f\|_{H^1(\mu, W_1)}^2 \le t^{-1}\|f\|_{H^1(\mu, W_1)}^2$ and $\|f\|_{H^1(\mu, W_0)}^2 \le (1-t)^{-1}\|f\|_{H^1(\mu, W_1)}^2$. We deduce that $H(\mu, W_t) \subseteq H(\mu, W_0)$ and $H(\mu, W_t) \subseteq H(\mu, W_1)$, and thus we can write

$$ \begin{aligned} \frac{1}{C(\mu, W_t)} &= \inf_{\substack{f \in H^1(\mu, W_t) \\ \int fd\mu = 0}} \frac{\int \nabla f^\top W_t \nabla f d\mu}{\int f^2 d\mu} \\ &\ge t \inf_{\substack{f \in H^1(\mu, W_t) \\ \int fd\mu = 0}} \frac{\int \nabla f^\top W_1 \nabla f d\mu}{\int f^2 d\mu} + (1-t) \inf_{\substack{f \in H^1(\mu, W_t) \\ \int fd\mu = 0}} \frac{\int \nabla f^\top W_0 \nabla f d\mu}{\int f^2 d\mu} \\ &\ge t \inf_{\substack{f \in H^1(\mu, W_1) \\ \int fd\mu = 0}} \frac{\int \nabla f^\top W_1 \nabla f d\mu}{\int f^2 d\mu} + (1-t) \inf_{\substack{f \in H^1(\mu, W_0) \\ \int fd\mu = 0}} \frac{\int \nabla f^\top W_0 \nabla f d\mu}{\int f^2 d\mu} \\ &= \frac{t}{C(\mu, W_1)} + \frac{1-t}{C(\mu, W_0)}, \end{aligned} $$

which shows that $W \mapsto C(\mu, W)^{-1}$ is concave. To show (21), we first notice that assumption $W' \subseteq \Omega W$ permits to write

$$ \|f\|_{H^1(\mu, W')}^2 \le \int f^2 d\mu + \Omega \int \nabla f^\top W \nabla f d\mu \le (1+\Omega) \|f\|_{H^1(\mu, W)}^2, $$

so that $H^1(\mu, W) \subseteq H^1(\mu, W')$. Next, we let $G \in \mathcal{G}(W)$ be such that

$$ G = \sum_{i=1}^{m} \omega_i \frac{\nabla u_i \nabla u_i^\top}{\int u_i^2 d\mu}, $$

for some $m \ge 1$, where $\omega_i \ge 0$ are such that $\sum_{i=1}^m \omega_i = 1$ and where $u_i \in \arg \min \left\{ \frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu} : f \in H^1(\mu, W), \int fd\mu = 0 \right\}$. Because $u_i \in H^1(\mu, W) \subseteq H^1(\mu, W')$ we can write

$$ \begin{aligned} \frac{1}{C(\mu, W')} &= \inf_{\substack{f \in H^1(\mu, W') \\ \int fd\mu = 0}} \frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu} + \frac{\int \nabla f^\top (W' - W) \nabla f d\mu}{\int f^2 d\mu} \\ &\le \frac{\int \nabla u_i^\top W \nabla u_i d\mu}{\int u_i^2 d\mu} + \frac{\int \nabla u_i^\top (W' - W) \nabla u_i d\mu}{\int u_i^2 d\mu} \\ &= \frac{1}{C(\mu, W)} + \int \operatorname{tr} \left( (W' - W) \frac{\nabla u_i \nabla u_i^\top}{\int u_i^2 d\mu} \right) d\mu. \end{aligned} $$

Multiplying the above inequality with $\omega_i \ge 0$ and summing over $i = 1, \dots, m$ yields (21) and concludes the proof. $\square$

Let us comment on Proposition 3. First, because $W \mapsto C(\mu, W)^{-1}$ is concave, it admits a superdifferential $\partial_W(C(\mu, W)^{-1})$ at any $W$. The superdifferential is the set of all $G$ such that (21) holds for any $W'$, see e.g. [17, 48]. Here, the superdifferential is identified using the scalar product $\langle W, G \rangle = \int \operatorname{tr}(WG)d\mu$. The set $\mathcal{G}(W)$ as in (22) is thus a subset of $\partial_W(C(\mu, W)^{-1})$. This will be used later on in Section 5.1 to propose a gradient algorithm for solving (7). Second, the definition (22) of $\mathcal{G}(W)$ requires that a minimizer of $\frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu}$ over $f \in H^1(\mu, W)$ exists. To analyze the existence of these minimizers, let us introduce the linear operator $\mathcal{L}_\mu^W$ from $H^1(\mu, W)$ to its dual defined by

$$ \int f \mathcal{L}_\mu^W g d\mu = - \int \nabla f^\top W \nabla g d\mu, \quad (23) $$
so that $\|f\|_{H^1(\mu, W_1)}^2 \le t^{-1}\|f\|_{H^1(\mu, W_1)}^2$ and $\|f\|_{H^1(\mu, W_0)}^2 \le (1-t)^{-1}\|f\|_{H^1(\mu, W_1)}^2$. We deduce that $H(\mu, W_t) \subseteq H(\mu, W_0)$ and $H(\mu, W_t) \subseteq H(\mu, W_1)$, and thus we can write

$$ 
\begin{aligned}
\frac{1}{C(\mu, W_t)} &= \inf_{\substack{f \in H^1(\mu, W_t) \\ \int f d\mu = 0}} \frac{\int \nabla f^\top W_t \nabla f d\mu}{\int f^2 d\mu} \\
&\ge t \inf_{\substack{f \in H^1(\mu, W_t) \\ \int f d\mu = 0}} \frac{\int \nabla f^\top W_1 \nabla f d\mu}{\int f^2 d\mu} + (1-t) \inf_{\substack{f \in H^1(\mu, W_t) \\ \int f d\mu = 0}} \frac{\int \nabla f^\top W_0 \nabla f d\mu}{\int f^2 d\mu} \\
&\ge t \inf_{\substack{f \in H^1(\mu, W_1) \\ \int f d\mu = 0}} \frac{\int \nabla f^\top W_1 \nabla f d\mu}{\int f^2 d\mu} + (1-t) \inf_{\substack{f \in H^1(\mu, W_0) \\ \int f d\mu = 0}} \frac{\int \nabla f^\top W_0 \nabla f d\mu}{\int f^2 d\mu} \\
&= \frac{t}{C(\mu, W_1)} + \frac{1-t}{C(\mu, W_0)},
\end{aligned}
$$

which shows that $W \mapsto C(\mu, W)^{-1}$ is concave. To show (21), we first notice that assumption $W' \subseteq \Omega W$ permits to write

$$ \|f\|_{H^1(\mu, W')}^2 \le \int f^2 d\mu + \Omega \int \nabla f^\top W \nabla f d\mu \le (1+\Omega) \|f\|_{H^1(\mu, W)}^2, $$

so that $H^1(\mu, W) \subseteq H^1(\mu, W')$. Next, we let $G \in \mathcal{G}(W)$ be such that

$$ G = \sum_{i=1}^{m} \omega_i \frac{\nabla u_i \nabla u_i^\top}{\int u_i^2 d\mu}, $$

for some $m \ge 1$, where $\omega_i \ge 0$ are such that $\sum_{i=1}^m \omega_i = 1$ and where $u_i \in \arg \min\{\frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu} : f \in H^1(\mu, W), \int f d\mu = 0\}$. Because $u_i \in H^1(\mu, W) \subseteq H^1(\mu, W')$ we can write

$$ 
\begin{aligned}
\frac{1}{C(\mu, W')} &= \inf_{\substack{f \in H^1(\mu, W') \\ \int f d\mu = 0}} \frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu} + \frac{\int \nabla f^\top (W' - W) \nabla f d\mu}{\int f^2 d\mu} \\
&\le \frac{\int \nabla u_i^\top W \nabla u_i d\mu}{\int u_i^2 d\mu} + \frac{\int \nabla u_i^\top (W' - W) \nabla u_i d\mu}{\int u_i^2 d\mu} \\
&= \frac{1}{C(\mu, W)} + \int \operatorname{tr} \left( (W' - W) \frac{\nabla u_i \nabla u_i^\top}{\int u_i^2 d\mu} \right) d\mu.
\end{aligned}
$$

Multiplying the above inequality with $\omega_i \ge 0$ and summing over $i = 1, \dots, m$ yields (21) and concludes the proof. $\square$

Let us comment on Proposition 3. First, because $W \mapsto C(\mu, W)^{-1}$ is concave, it admits a superdifferential $\partial_W(C(\mu, W)^{-1})$ at any $W$. The superdifferential is the set of all $G$ such that (21) holds for any $W'$, see e.g. [17, 48]. Here, the superdifferential is identified using the scalar product $\langle W, G \rangle = \int \operatorname{tr}(WG)d\mu$. The set $\mathcal{G}(W)$ as in (22) is thus a subset of $\partial_W(C(\mu, W)^{-1})$. This will be used later on in Section 5.1 to propose a gradient algorithm for solving (7). Second, the definition (22) of $\mathcal{G}(W)$ requires that a minimizer of $\frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu}$ over $f \in H^1(\mu, W)$ exists. To analyze the existence of these minimizers, let us introduce the linear operator $\mathcal{L}_\mu^W$ from $H^1(\mu, W)$ to its dual defined by

$$ \int f \mathcal{L}_\mu^W g d\mu = - \int \nabla f^\top W \nabla g d\mu, \qquad (23) $$
for any $f, g \in H^1(\mu, W)$. The operator $\mathcal{L}_\mu^W$ is the infinitesimal generator associated with the SDE (5). Using one integration by part, $\mathcal{L}_\mu^W$ can be expressed as

$$ \mathcal{L}_\mu^W f = W : \nabla^2 f + (\mathrm{div}(W) + W \nabla \log \mu)^\top \nabla f, $$

where $W : \nabla^2 f = \sum_{i,j=1}^d W_{i,j} \partial_{i,j}^2 f$. The following proposition relates $C(\mu, W)$ with the spectrum of $\mathcal{L}_\mu^W$, and ensures that minimizers of $\frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu}$ over $f \in H^1(\mu, W)$ exists.

**Proposition 4.** Assume that $W$ is such that the injection of $H^1(\mu, W)$ in $L^2(\mu)$ is compact and that $H^1(\mu, W)$ is dense in $L^2(\mu)$. Then, if $C(\mu, W) < \infty$, there exists an orthonormal basis $\{u_i^W\}_{i=1}^\infty$ of $L^2(\mu)$ which forms an eigenbasis for $-\mathcal{L}_\mu^W$ such that

$$ -\mathcal{L}_\mu^W u_i^W = \lambda_i^W u_i^W, \qquad (24) $$

where $\lambda_i^W \ge 0$ is the $i$-th smallest eigenvalue of $-\mathcal{L}_\mu^W$. In addition we have that:

1. the eigenvalues $\lambda_i^W$ tends to $+\infty$ when $i \to +\infty$,

2. the first eigenvalue is trivial $\lambda_1^W = 0$ and the associated eigenfunction is the constant function $u_1^W = 1$,

3. the second eigenvalue is positive $\lambda_2^W > 0$ and satisfies

$$ \lambda_2^W = \frac{1}{C(\mu, W)}. \qquad (25) $$

4. the set of minimizers of $\frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu}$ over $f \in H^1(\mu, W)$, $\int f d\mu = 0$, is the eigenspace span$\{u_2^W, \dots, u_{m+1}^W\}$ associated with the eigenvalue $\lambda_2^W$, and its dimension $m$ is necessarily finite. Moreover, the superdifferential of the concave function $W \mapsto \lambda_2^W$ is given by

$$ \partial_W \lambda_2^W = \overline{\mathrm{conv}} \left\{ \frac{\nabla u \nabla u^\top}{\int u^2 d\mu} : u \in \mathrm{span}\{u_2^W, \dots, u_{m+1}^W\} \right\}. \qquad (26) $$

where $\overline{\cdot}$ denotes the closure with respect to the topology induced by the $L^2(\mu; \mathbb{R}^{d \times d})$-norm on the matrix field $\|A\|_{L^2(\mu; \mathbb{R}^{d \times d})} = (\int \|A\|_F^2 d\mu)^{1/2}$. In particular if $m = 1$ then $\lambda_2^W$ is differentiable and

$$ \nabla \lambda_2^W \stackrel{m=1}{=} \frac{(\nabla u_2^W)(\nabla u_2^W)^\top}{\int (u_2^W)^2 d\mu}. \qquad (27) $$

*Proof.* To show that $-\mathcal{L}_\mu^W$ admits the eigendecomposition (24), it is sufficient to show that the bi-linear form $(f, g) \mapsto \int f(-\mathcal{L}_\mu^W)g d\mu = \int \nabla f^\top W \nabla g d\mu$ is continuous and coercive on $H^1(\mu, W) \subset L^2(\mu)$, see e.g. Theorem 7.3.2 in [2]. To show the continuity, a Cauchy-Schwarz inequality permits to write

$$ \begin{aligned} & \left| \int f(-\mathcal{L}_\mu^W)g d\mu - \int f'(-\mathcal{L}_\mu^W)g' d\mu \right| \\ &= \left| \int (\sqrt{W} \nabla f)^\top (\sqrt{W} (\nabla g - \nabla g')) d\mu - \int (\sqrt{W} (\nabla f' - \nabla f))^\top \sqrt{W} \nabla g' d\mu \right| \\ &\le \|g - g'\|_{H^1(\mu, W)} \left( \int \nabla f^\top W \nabla f d\mu \right)^{1/2} + \|f - f'\|_{H^1(\mu, W)} \left( \int \nabla g'^\top W \nabla g' d\mu \right)^{1/2}. \end{aligned} $$
for any $f, g \in H^1(\mu, W)$. The operator $\mathcal{L}_\mu^W$ is the infinitesimal generator associated with the SDE (5). Using one integration by part, $\mathcal{L}_\mu^W$ can be expressed as

$$ \mathcal{L}_\mu^W f = W : \nabla^2 f + (\mathrm{div}(W) + W \nabla \log \mu)^\top \nabla f, $$

where $W : \nabla^2 f = \sum_{i,j=1}^d W_{i,j} \partial_{i,j}^2 f$. The following proposition relates $C(\mu, W)$ with the spectrum of $\mathcal{L}_\mu^W$, and ensures that minimizers of $\frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu}$ over $f \in H^1(\mu, W)$ exists.

**Proposition 4.** Assume that $W$ is such that the injection of $H^1(\mu, W)$ in $L^2(\mu)$ is compact and that $H^1(\mu, W)$ is dense in $L^2(\mu)$. Then, if $C(\mu, W) < \infty$, there exists an orthonormal basis $\{u_i^W\}_{i=1}^\infty$ of $L^2(\mu)$ which forms an eigenbasis for $-\mathcal{L}_\mu^W$ such that

$$ -\mathcal{L}_\mu^W u_i^W = \lambda_i^W u_i^W, \qquad (24) $$

where $\lambda_i^W \ge 0$ is the $i$-th smallest eigenvalue of $-\mathcal{L}_\mu^W$. In addition we have that:

1. the eigenvalues $\lambda_i^W$ tends to $+\infty$ when $i \to +\infty$,

2. the first eigenvalue is trivial $\lambda_1^W = 0$ and the associated eigenfunction is the constant function $u_1^W = 1$,

3. the second eigenvalue is positive $\lambda_2^W > 0$ and satisfies

$$ \lambda_2^W = \frac{1}{C(\mu, W)}. \qquad (25) $$

4. the set of minimizers of $\frac{\int \nabla f^\top W \nabla f d\mu}{\int f^2 d\mu}$ over $f \in H^1(\mu, W)$, $\int f d\mu = 0$, is the eigenspace span$\{u_2^W, \dots, u_{m+1}^W\}$ associated with the eigenvalue $\lambda_2^W$, and its dimension $m$ is necessarily finite. Moreover, the superdifferential of the concave function $W \mapsto \lambda_2^W$ is given by

$$ \partial_W \lambda_2^W = \overline{\mathrm{conv}} \left\{ \frac{\nabla u \nabla u^\top}{\int u^2 d\mu} : u \in \mathrm{span}\{u_2^W, \dots, u_{m+1}^W\} \right\}. \qquad (26) $$

where $\overline{\cdot}$ denotes the closure with respect to the topology induced by the $L^2(\mu; \mathbb{R}^{d \times d})$-norm on the matrix field $\|A\|_{L^2(\mu; \mathbb{R}^{d \times d})} = (\int \|A\|_F^2 d\mu)^{1/2}$. In particular if $m = 1$ then $\lambda_2^W$ is differentiable and

$$ \nabla \lambda_2^W \stackrel{m=1}{=} \frac{(\nabla u_2^W)(\nabla u_2^W)^\top}{\int (u_2^W)^2 d\mu}. \qquad (27) $$

*Proof.* To show that $-\mathcal{L}_\mu^W$ admits the eigendecomposition (24), it is sufficient to show that the bi-linear form $(f, g) \mapsto \int f(-\mathcal{L}_\mu^W)g d\mu = \int \nabla f^\top W \nabla g d\mu$ is continuous and coercive on $H^1(\mu, W) \subset L^2(\mu)$, see e.g. Theorem 7.3.2 in [2]. To show the continuity, a Cauchy-Schwarz inequality permits to write

$$ \begin{aligned} & \left| \int f(-\mathcal{L}_\mu^W)g d\mu - \int f'(-\mathcal{L}_\mu^W)g' d\mu \right| \\ &= \left| \int (\sqrt{W} \nabla f)^\top (\sqrt{W} (\nabla g - \nabla g')) d\mu - \int (\sqrt{W} (\nabla f' - \nabla f))^\top \sqrt{W} \nabla g' d\mu \right| \\ &\le \|g - g'\|_{H^1(\mu, W)} \left( \int \nabla f^\top W \nabla f d\mu \right)^{1/2} + \|f - f'\|_{H^1(\mu, W)} \left( \int \nabla g'^\top W \nabla g' d\mu \right)^{1/2}. \end{aligned} $$
Then $\int f(-\mathcal{L}_{\mu}^{W})g d\mu \rightarrow \int f'(-\mathcal{L}_{\mu}^{W})g' d\mu$ whenever $f \rightarrow f'$ and $g \rightarrow g'$ in $H^{1}(\mu, W)$. The coercivity is a direct consequence of $C(\mu, W) < \infty$ since

$$ \|f\|_{H^1(\mu, W)}^2 = \int f^2 d\mu + \int \nabla f^\top W \nabla f d\mu \le (C(\mu, W) + 1) \int f(-\mathcal{L}_\mu^W) f d\mu, $$

holds for any $f \in H^1(\mu, W)$. Next, because the injection of $H^1(\mu, W)$ in $L^2(\mu)$ is compact and that $H^1(\mu, W)$ is dense in $L^2(\mu)$, Theorem 7.3.2 in [2] ensures that $-\mathcal{L}_{\mu}^{W}$ admits the eigendecomposition (24) and that $\lambda_{i}^{W} \rightarrow +\infty$ when $i \rightarrow +\infty$. By (23), the first eigenpaire is trivial with $\lambda_{1}^{W} = 0$ and $u_{1}^{W} = 1$. Using the Courant-Fischer-Weyl min-max principle, we have

$$ \lambda_2^W = \inf_{\substack{f \in H^1(\mu, W) \\ \int f u_1^W d\mu = 0}} \frac{\int f(-\mathcal{L}_\mu^W) f d\mu}{\int f^2 d\mu}, $$

which, by (20), gives (25). Because $\lambda_{i}^{W} \rightarrow +\infty$, we have that the multiplicity $m$ of $\lambda_{2}^{W}$ is necessarily finite. As a consequence, the set of the minizers of

$$ \frac{\int f(-\mathcal{L}_{\mu}^{W}) f d\mu}{\int f^2 d\mu} = \int \operatorname{tr} \left( W \frac{\nabla f \nabla f^{\top}}{\int f^2 d\mu} \right) d\mu, $$

over $f \in H^1(\mu, W)$ with $\int f d\mu = 0$ is exactly span{$u_2^W, \dots, u_{m+1}^W$}. Thus, Theorem 2.4.18 from [58] establishes that (26) holds. Finally, (27) is a direct consequence of (26), which concludes the proof. $\square$

Let us comment on the assumptions of Proposition 4. First, the density of $H^1(\mu, W)$ in $L^2(\mu)$ is classically obtained by showing that $C_c^{\infty}(\mathbb{R}^d)$, the space of functions of class $C^{\infty}$ with compact support in $\mathbb{R}^d$, is dense in both $H^1(\mu, W)$ and $L^2(\mu)$. In the case of classical (unweighted) Sobolev space, it is well known that $C_c^{\infty}(\mathbb{R}^d)$ is dense in $L^2(\mathbb{R}^d)$ and in $H^1(\mathbb{R}^d)$, see e.g. Theorem 4.12 and Corollary 9.8 from [19]. For weighted Sobolev spaces, the result still holds under the assumption that the Lebesgue density of $\mu$ is in the Muckenhoupt's $A_2$-class, see Theorem 1.1 in [40]. In particular, any measure $\mu$ with bounded support $\Omega = \text{supp}(\mu)$ and whose Lebesgue density is bounded away from 0 and from $+\infty$ is in the Muckenhoupt's $A_2$-class. If we further assume that $\alpha I_d \le W \le \beta I_d$ for some constant $0 < \alpha \le \beta < \infty$, then $H^1(\mu, W) = H^1(\mu)$ is dense in $L^2(\mu)$.

Second, the compact injection of $H^1(\mu, W)$ in $L^2(\mu)$ can be obtained via the Rellich-Kondrachov theorem, see Theorem 9.16 [19]. This theorem states that if $\Omega \subset \mathbb{R}^d$ is a bounded subset of class $C^1$, then the injection of the (unweighted) Sobolev space $H^1(\Omega)$ in $L^2(\Omega)$ is compact. Therefore, if $\Omega = \text{supp}(\mu)$ and $\alpha I_d \le W \le \beta I_d$ and $\alpha \le \frac{d\mu}{dx} \le \beta$ for some constant $0 < \alpha \le \beta < \infty$, then the $H^1(\Omega) = H^1(\mu, W)$ is dense in $L^2(\Omega) = L^2(\mu)$.

# 5 Solution algorithm

In order to account for the two constraints $W(x) \in S_+^d$ and $\int \operatorname{tr}(W) d\mu = \operatorname{tr}(\operatorname{Cov}_\mu)$ in Problem (7), we propose to use the following parametrization of $W$

$$ W(x) = V(x)^2 \frac{\operatorname{tr}(\operatorname{Cov}_\mu)}{\int \|V\|_F^2 d\mu}, \qquad (28) $$

where $V: \mathbb{R}^d \rightarrow S^d$ is a matrix field to be determined. Here, $S^d \subset \mathbb{R}^{d \times d}$ denotes the set of symmetric matrices and $\| \cdot \|_F$ denotes the Frobenius matrix norm such that
Then $\int f(-\mathcal{L}_{\mu}^{W})g d\mu \rightarrow \int f'(-\mathcal{L}_{\mu}^{W})g' d\mu$ whenever $f \rightarrow f'$ and $g \rightarrow g'$ in $H^{1}(\mu, W)$. The coercivity is a direct consequence of $C(\mu, W) < \infty$ since

$$ \|f\|_{H^1(\mu, W)}^2 = \int f^2 d\mu + \int \nabla f^\top W \nabla f d\mu \le (C(\mu, W) + 1) \int f(-\mathcal{L}_\mu^W) f d\mu, $$

holds for any $f \in H^1(\mu, W)$. Next, because the injection of $H^1(\mu, W)$ in $L^2(\mu)$ is compact and that $H^1(\mu, W)$ is dense in $L^2(\mu)$, Theorem 7.3.2 in [2] ensures that $-\mathcal{L}_{\mu}^{W}$ admits the eigendecomposition (24) and that $\lambda_{i}^{W} \rightarrow +\infty$ when $i \rightarrow +\infty$. By (23), the first eigenpaire is trivial with $\lambda_{1}^{W} = 0$ and $u_{1}^{W} = 1$. Using the Courant-Fischer-Weyl min-max principle, we have

$$ \lambda_2^W = \inf_{\substack{f \in H^1(\mu, W) \\ \int f u_1^W d\mu = 0}} \frac{\int f(-\mathcal{L}_\mu^W) f d\mu}{\int f^2 d\mu}, $$

which, by (20), gives (25). Because $\lambda_{i}^{W} \rightarrow +\infty$, we have that the multiplicity $m$ of $\lambda_{2}^{W}$ is necessarily finite. As a consequence, the set of the minizers of

$$ \frac{\int f(-\mathcal{L}_{\mu}^{W}) f d\mu}{\int f^2 d\mu} = \int \mathrm{tr} \left( W \frac{\nabla f \nabla f^{\top}}{\int f^2 d\mu} \right) d\mu, $$

over $f \in H^1(\mu, W)$ with $\int f d\mu = 0$ is exactly $\mathrm{span}\{u_2^W, \dots, u_{m+1}^W\}$. Thus, Theorem 2.4.18 from [58] establishes that (26) holds. Finally, (27) is a direct consequence of (26), which concludes the proof. $\square$

Let us comment on the assumptions of Proposition 4. First, the density of $H^1(\mu, W)$ in $L^2(\mu)$ is classically obtained by showing that $C_c^{\infty}(\mathbb{R}^d)$, the space of functions of class $C^{\infty}$ with compact support in $\mathbb{R}^d$, is dense in both $H^1(\mu, W)$ and $L^2(\mu)$. In the case of classical (unweighted) Sobolev space, it is well known that $C_c^{\infty}(\mathbb{R}^d)$ is dense in $L^2(\mathbb{R}^d)$ and in $H^1(\mathbb{R}^d)$, see e.g. Theorem 4.12 and Corollary 9.8 from [19]. For weighted Sobolev spaces, the result still holds under the assumption that the Lebesgue density of $\mu$ is in the Muckenhoupt's $A_2$-class, see Theorem 1.1 in [40]. In particular, any measure $\mu$ with bounded support $\Omega = \mathrm{supp}(\mu)$ and whose Lebesgue density is bounded away from 0 and from $+\infty$ is in the Muckenhoupt's $A_2$-class. If we further assume that $\alpha I_d \le W \le \beta I_d$ for some constant $0 < \alpha \le \beta < \infty$, then $H^1(\mu, W) = H^1(\mu)$ is dense in $L^2(\mu)$.

Second, the compact injection of $H^1(\mu, W)$ in $L^2(\mu)$ can be obtained via the Rellich-Kondrachov theorem, see Theorem 9.16 [19]. This theorem states that if $\Omega \subset \mathbb{R}^d$ is a bounded subset of class $C^1$, then the injection of the (unweighted) Sobolev space $H^1(\Omega)$ in $L^2(\Omega)$ is compact. Therefore, if $\Omega = \mathrm{supp}(\mu)$ and $\alpha I_d \le W \le \beta I_d$ and $\alpha \le \frac{d\mu}{dx} \le \beta$ for some constant $0 < \alpha \le \beta < \infty$, then the $H^1(\Omega) = H^1(\mu, W)$ is dense in $L^2(\Omega) = L^2(\mu)$.

# 5 Solution algorithm

In order to account for the two constraints $W(x) \in S_+^d$ and $\int \mathrm{tr}(W) d\mu = \mathrm{tr}(\mathrm{Cov}_\mu)$ in Problem (7), we propose to use the following parametrization of $W$

$$ W(x) = V(x)^2 \frac{\mathrm{tr}(\mathrm{Cov}_\mu)}{\int \|V\|_F^2 d\mu}, \qquad (28) $$

where $V: \mathbb{R}^d \rightarrow S^d$ is a matrix field to be determined. Here, $S^d \subset \mathbb{R}^{d \times d}$ denotes the set of symmetric matrices and $\| \cdot \|_F$ denotes the Frobenius matrix norm such that
$\|A\|_F^2 = \text{tr}(A^\top A)$. Because any $W: \mathbb{R}^d \to S_+^d$ with $\int \text{tr}(W)d\mu = \text{tr}(\text{Cov}_\mu)$ can be written as in (28) for some $V$, we have that any maximizer of

$$J(V) = \inf_{\substack{f \in H^1(\mu, V^2) \\ \int fd\mu = 0}} \frac{\int \nabla f^\top V^2 \nabla f d\mu}{(\int \|V\|_F^2 d\mu)(\int f^2 d\mu)} \stackrel{(20)\& (28)}{=} \frac{1}{\text{tr}(\text{Cov}_\mu)C(\mu, W)}, \quad (29)$$

yields a solution $W$ to (7). By parametrizing $W$ with $V$, we lose the concavity of the constrained problem (7) (recall Proposition 3) but instead we obtain an unconstrained maximization problem

$$\max_{V:\mathbb{R}^d \to S^d} J(V), \quad (30)$$

which is much easier to solve with in practice.

## 5.1 Gradient ascent algorithms

Proposition 4 permits to write $J(V) = \lambda_2^{V^2} / \int \|V\|_F^2 d\mu$, where $\lambda_2^{V^2}$ denotes the first nonzero eigenvalue of the operator $\mathcal{L}_\mu^{V^2}$ as in (23). If we assume that $\lambda_2^{V^2}$ has multiplicity one, then $W \mapsto \lambda_2^W$ is differentiable around $W = V^2$. In this case, $J$ is also differentiable around $V$ as the product and composition of differentiable functions, and its gradient is given by

$$\nabla J(V) = \frac{GV + VG - 2J(V)V}{\int \|V\|_F^2 d\mu}, \quad \text{with } G = \frac{(\nabla u_2^{V^2})(\nabla u_2^{V^2})^\top}{\int (u_2^{V^2})^2 d\mu}. \quad (31)$$

If the multiplicity of $\lambda_2^{V^2}$ is larger than one, then $J$ is no longer differentiable nor admits a superdifferential (because $J$ is not concave). Still, $J$ possesses a generalized gradient in the sense of [26, Chapter 10] which is given with a similar expression at in (31), albeit $G$ now belonging to the superdifferential $\partial_W(\lambda_2^{V^2})$. This is, however, not further analyzed in the present article.

The gradient ascent algorithm with constant step size $\rho \ge 0$ consists in defining a sequence of $\{V^{(k)}\}_{k \ge 0}$ as follow

$$V^{(k+1/2)} = V^{(k)} + \rho \nabla J(V^{(k)}), \quad (32)$$

$$V^{(k+1)} = \frac{V^{(k+1/2)}}{\left(\int \|V^{(k+1/2)}\|_F^2 d\mu\right)^{1/2}}. \quad (33)$$

Since $J(\alpha V) = J(V)$ for all $\alpha \ne 0$, the last step above permits to normalize the iterates without modifying the value of the objective $J(V^{(k+1)}) = J(V^{(k+1/2)})$. While gradient ascent is universally popular, alternative methods such as momentum [47] and Nesterov's Accelerated Gradient [41] can result in significantly faster convergence to the optimum. The momentum algorithm is

$$D^{(k+1)} = \alpha D^{(k)} + \rho \nabla J(V^{(k)}), \quad (34)$$

$$V^{(k+1/2)} = V^{(k)} + D^{(k+1)}, \quad (35)$$

$$V^{(k+1)} = \frac{V^{(k+1/2)}}{\left(\int \|V^{(k+1/2)}\|_F^2 d\mu\right)^{1/2}}. \quad (36)$$

It consists in continuing updating along the previous update direction $D^{(k)}$, scaled by a factor $0 \le \alpha < 1$. The Nesterov's Accelerated Gradient algorithm replaces (34) with

$$D^{(k+1)} = \alpha_k D^{(k)} + \rho \nabla J(V^{(k)} + \alpha_k D^{(k)}), \quad (37)$$
$\|A\|_F^2 = \text{tr}(A^\top A)$. Because any $W: \mathbb{R}^d \to S_+^d$ with $\int \text{tr}(W)d\mu = \text{tr}(\text{Cov}_\mu)$ can be written as in (28) for some $V$, we have that any maximizer of

$$J(V) = \inf_{\substack{f \in H^1(\mu, V^2) \\ \int fd\mu = 0}} \frac{\int \nabla f^\top V^2 \nabla f d\mu}{\int \|V\|_F^2 d\mu} \stackrel{(20)\& (28)}{=} \frac{1}{\text{tr}(\text{Cov}_\mu)C(\mu, W)}, \quad (29)$$

yields a solution $W$ to (7). By parametrizing $W$ with $V$, we lose the concavity of the constrained problem (7) (recall Proposition 3) but instead we obtain an unconstrained maximization problem

$$\max_{V:\mathbb{R}^d \to S^d} J(V), \quad (30)$$

which is much easier to solve with in practice.

## 5.1 Gradient ascent algorithms

Proposition 4 permits to write $J(V) = \lambda_2^{V^2} / \int \|V\|_F^2 d\mu$, where $\lambda_2^{V^2}$ denotes the first nonzero eigenvalue of the operator $\mathcal{L}_\mu^{V^2}$ as in (23). If we assume that $\lambda_2^{V^2}$ has multiplicity one, then $W \mapsto \lambda_2^W$ is differentiable around $W = V^2$. In this case, $J$ is also differentiable around $V$ as the product and composition of differentiable functions, and its gradient is given by

$$\nabla J(V) = \frac{GV + VG - 2J(V)V}{\int \|V\|_F^2 d\mu}, \quad \text{with } G = \frac{(\nabla u_2^{V^2})(\nabla u_2^{V^2})^\top}{\int (u_2^{V^2})^2 d\mu}. \quad (31)$$

If the multiplicity of $\lambda_2^{V^2}$ is larger than one, then $J$ is no longer differentiable nor admits a superdifferential (because $J$ is not concave). Still, $J$ possesses a generalized gradient in the sense of [26, Chapter 10] which is given with a similar expression at in (31), albeit $G$ now belonging to the superdifferential $\partial_W(\lambda_2^{V^2})$. This is, however, not further analyzed in the present article.

The gradient ascent algorithm with constant step size $\rho \ge 0$ consists in defining a sequence of $\{V^{(k)}\}_{k \ge 0}$ as follow

$$V^{(k+1/2)} = V^{(k)} + \rho \nabla J(V^{(k)}), \quad (32)$$

$$V^{(k+1)} = \frac{V^{(k+1/2)}}{\int \|V^{(k+1/2)}\|_F^2 d\mu} \quad (33)$$

Since $J(\alpha V) = J(V)$ for all $\alpha \ne 0$, the last step above permits to normalize the iterates without modifying the value of the objective $J(V^{(k+1)}) = J(V^{(k+1/2)})$. While gradient ascent is universally popular, alternative methods such as momentum [47] and Nesterov's Accelerated Gradient [41] can result in significantly faster convergence to the optimum. The momentum algorithm is

$$D^{(k+1)} = \alpha D^{(k)} + \rho \nabla J(V^{(k)}), \quad (34)$$

$$V^{(k+1/2)} = V^{(k)} + D^{(k+1)}, \quad (35)$$

$$V^{(k+1)} = \frac{V^{(k+1/2)}}{\int \|V^{(k+1/2)}\|_F^2 d\mu} \quad (36)$$

It consists in continuing updating along the previous update direction $D^{(k)}$, scaled by a factor $0 \le \alpha < 1$. The Nesterov's Accelerated Gradient algorithm replaces (34) with

$$D^{(k+1)} = \alpha_k D^{(k)} + \rho \nabla J(V^{(k)} + \alpha_k D^{(k)}), \quad (37)$$
where $0 \le \alpha_k < 1$ is now depending on the iteration number, see [54]. Nesterov proposed $\alpha_k = 1 - 3/(5+k)$, which we adopt in our experiments.

In the previous equations, we compute the gradient $\nabla J(V^{(k)})$ as in (31), which assumes $\lambda_2^{(V^{(k)})^2}$ has multiplicity one. As shown in the numerical results, the first two nonzero eigenvalues cross each other several times during the iteration process and, at optimality, they tend to be equal. Even though the first two nonzero eigenvalues are always numerically different, this indicates that the objective function $J$ possesses several edges (discontinuities in the gradient). The use of either the momentum algorithm or NAG turns out to be very efficient for solving (30) when the iterates approach optimality.

## 5.2 Finite element discretization

In order to numerically compute the gradient $\nabla J(V^{(k)})$ in (31), we employ the finite element method [2, 59]. We consider a triangulation $\mathcal{T}_M$ of the support of $\mu$ made with $M$ simplices (i.e. triangles when $d=2$) and $N$ nodes such that $\text{supp}(\mu) = \omega_1 \cup \dots \cup \omega_M$. Here, we assume that the support of $\mu$ is bounded. We define the finite element space $\mathcal{V}_N \subset H^1(\mu, W)$ as

$$ \mathcal{V}_N = \operatorname{span}\{\phi_1, \dots, \phi_N\}, \quad \phi_i: \text{continuous and piecewise affine,} $$

such that $\phi_i(x_j) = \delta_{i,j}$, where $x_j$ is the $j$-th node of $\mathcal{T}_M$. The function

$$ J_N(V) = \inf_{\substack{f \in \mathcal{V}_N \\ \int f d\mu = 0}} \frac{\int \nabla f^\top V^2 \nabla f d\mu}{(\int \|V\|_F^2 d\mu)(\int f^2 d\mu)}, $$

yields an approximation to $J(V)$ such that $J_N(V) \to J(V)$ when $N \to \infty$, see e.g. Theorem 7.4.4 from [2]. We propose to replace $J$ by $J_N$ in the gradient algorithms (32), (34) and (37). Next, we show how to compute $\nabla J_N(V)$.

For a given matrix field $W$ as in (28), we let $K^W \in \mathbb{R}^{N \times N}$ and $M \in \mathbb{R}^{N \times N}$ be the matrices defined by

$$ K_{i,j}^W = \int \nabla \phi_i^\top W \nabla \phi_j d\mu, \quad \text{and} \quad M_{i,j} = \int \phi_i \phi_j d\mu, \qquad (38) $$

for all $1 \le i, j \le N$. Consider the generalized eigendecomposition of the matrix pair $(K^W, M)$

$$ K^W U_i^W = \lambda_{N,i}^W M U_i^W, $$

where $(U_j^W)^\top M U_i^W = \delta_{i,j}$ for any $1 \le i, j \le N$ and where the eigenvalues $\lambda_{N,1}^W, \dots, \lambda_{N,N}^W$ are listed in the increasing order. Because the finite element space $\mathcal{V}_N$ contains the constant function, the first eigenvalue is trivial $\lambda_{N,1}^W = 0$ and the corresponding eigenvector $U_1^W = U_1 = (1, \dots, 1)^\top$ represents the constant function $1 = \sum_{i=1}^N U_{1,i} \phi_i$ onto the finite element basis. Also, because $\int f d\mu = U^\top M U_1$ for any $f = \sum_{i=1}^N U_i \phi_i$, we can express $J_N(V)$ as follow

$$ J_N(V) = \frac{1}{\int \|V\|_2^2 d\mu} \inf_{\substack{U \in \mathbb{R}^N \\ U^\top M U_1 = 0}} \frac{U^\top K^2 U}{U^\top M U} = \frac{\lambda_{N,2}^{V^2}}{\int \|V\|_2^2 d\mu}. $$
where $0 \le \alpha_k < 1$ is now depending on the iteration number, see [54]. Nesterov proposed $\alpha_k = 1 - 3/(5+k)$, which we adopt in our experiments.

In the previous equations, we compute the gradient $\nabla J(V^{(k)})$ as in (31), which assumes $\lambda_2^{(V^{(k)})^2}$ has multiplicity one. As shown in the numerical results, the first two nonzero eigenvalues cross each other several times during the iteration process and, at optimality, they tend to be equal. Even though the first two nonzero eigenvalues are always numerically different, this indicates that the objective function $J$ possesses several edges (discontinuities in the gradient). The use of either the momentum algorithm or NAG turns out to be very efficient for solving (30) when the iterates approach optimality.

## 5.2 Finite element discretization

In order to numerically compute the gradient $\nabla J(V^{(k)})$ in (31), we employ the finite element method [2, 59]. We consider a triangulation $\mathcal{T}_M$ of the support of $\mu$ made with $M$ simplices (i.e. triangles when $d=2$) and $N$ nodes such that $\text{supp}(\mu) = \omega_1 \cup \dots \cup \omega_M$. Here, we assume that the support of $\mu$ is bounded. We define the finite element space $\mathcal{V}_N \subset H^1(\mu, W)$ as

$$ \mathcal{V}_N = \text{span}\{\phi_1, \dots, \phi_N\}, \quad \phi_i: \text{continuous and piecewise affine,} $$

such that $\phi_i(x_j) = \delta_{i,j}$, where $x_j$ is the $j$-th node of $\mathcal{T}_M$. The function

$$ J_N(V) = \inf_{\substack{f \in \mathcal{V}_N \\ \int f d\mu = 0}} \frac{\int \nabla f^\top V^2 \nabla f d\mu}{(\int \|V\|_F^2 d\mu)(\int f^2 d\mu)}, $$

yields an approximation to $J(V)$ such that $J_N(V) \to J(V)$ when $N \to \infty$, see e.g. Theorem 7.4.4 from [2]. We propose to replace $J$ by $J_N$ in the gradient algorithms (32), (34) and (37). Next, we show how to compute $\nabla J_N(V)$.

For a given matrix field $W$ as in (28), we let $K^W \in \mathbb{R}^{N \times N}$ and $M \in \mathbb{R}^{N \times N}$ be the matrices defined by

$$ K_{i,j}^W = \int \nabla \phi_i^\top W \nabla \phi_j d\mu, \quad \text{and} \quad M_{i,j} = \int \phi_i \phi_j d\mu, \qquad (38) $$

for all $1 \le i, j \le N$. Consider the generalized eigendecomposition of the matrix pair $(K^W, M)$

$$ K^W U_i^W = \lambda_{N,i}^W M U_i^W, $$

where $(U_j^W)^\top M U_i^W = \delta_{i,j}$ for any $1 \le i, j \le N$ and where the eigenvalues $\lambda_{N,1}^W, \dots, \lambda_{N,N}^W$ are listed in the increasing order. Because the finite element space $\mathcal{V}_N$ contains the constant function, the first eigenvalue is trivial $\lambda_{N,1}^W = 0$ and the corresponding eigenvector $U_1^W = U_1 = (1, \dots, 1)^\top$ represents the constant function $1 = \sum_{i=1}^N U_{1,i} \phi_i$ onto the finite element basis. Also, because $\int f d\mu = U^\top M U_1$ for any $f = \sum_{i=1}^N U_i \phi_i$, we can express $J_N(V)$ as follow

$$ J_N(V) = \frac{1}{\int \|V\|_2^2 d\mu} \inf_{\substack{U \in \mathbb{R}^N \\ U^\top M U_1 = 0}} \frac{U^\top K^{V^2} U}{U^\top M U} = \frac{\lambda_{N,2}^{V^2}}{\int \|V\|_2^2 d\mu}. $$
As in Equation (31), if the multiplicity of $\lambda_{N,2}^{V^2}$ is one, then $J_N$ admits a gradient at $V$ which is given by

$$ 
\begin{gathered}
\nabla J_N(V) = \frac{G_N V + V G_N - 2J_N(V)V}{\int \|V\|_F^2 d\mu}, \\
\text{with } G_N = \frac{\left(\sum_{i=1}^N (U_2^{V^2})_i \nabla \phi_i\right) \left(\sum_{i=1}^N (U_2^{V^2})_i \nabla \phi_i\right)^\top}{(U_2^{V^2})^\top M (U_2^{V^2})}.
\end{gathered}
\quad (39) $$

Thus, to compute $\nabla J_N(V)$, one needs first to assemble the finite element matrices $K^{V^2}$ and $M$ and to compute the second smallest generalized eigenpair $(\lambda_2^{V^2}, U_2^{V^2})$ of $(K^{V^2}, M)$. Then, $\nabla J_N(V)$ is assembled according to (39). Notice that, because $\phi_i$ are piecewise affine, the matrix field $G_N$ is piecewise constant on each element $\omega_i$ of the mesh. In addition, if the field $V$ is also piecewise constant of the form of

$$ V(x) = \sum_{m=1}^{M} V_m \mathbb{1}_{\omega_m}(x), \quad \text{where } V_m \in \mathbb{R}^{d \times d}, \quad (40) $$

then $\nabla J_N(V)$ as in (39) is also piecewise constant. Here, $\mathbb{1}_{\omega_m}$ is the indicator function of the $m$-th element $\omega_m$. Therefore, any of the iterates of the gradient algorithms in Equations (32)–(37) remains piecewise constant, and so is the resulting symmetric positive field $W$ given by (28). The resulting algorithm is summarized in Algorithm 1.

# 6 Numerical experiments

We demonstrate the proposed method on three benchmarks, each of them in dimension $d=2$. Code for reproducing all of these numerical results is available online<fnref n="2" />. The first is a Gaussian mixture of the form of

$$ d\mu_1(x) \propto \sum_{i=1}^{3} \exp \left( - \frac{\|x - x_i\|^2}{\sigma^2} \right) dx, $$

where $\sigma^2 = 0.025$, $x_1 = (0, 0.5)$, $x_2 = (\sqrt{3}/4, \sqrt{3}/4)$ and $x_3 = (-\sqrt{3}/4, \sqrt{3}/4)$. The second distribution is concentrated on a ring of radius $r$ and is defined by

$$ d\mu_2(x) \propto \exp \left( - \frac{(\|x\| - r)^2}{\sigma^2} \right) dx, $$

where $\sigma^2 = 0.0032$ and $r = 0.65$. The third is a uniform measure on a $H$-shaped domain $\Omega \subset \mathbb{R}^2$, meaning

$$ d\mu_3(x) \propto \mathbb{1}_{\Omega}(x)dx. $$

The fourth is a slight modification of $\mu_3$ which consists in enlarging the support of $\mu_3$ to the convex hull $\text{conv}(\Omega)$ of the $H$-shaped $\Omega$ as follow

$$ d\mu_4(x) \propto d\mu_3(x) + \varepsilon \mathbb{1}_{\text{conv}(\Omega)}(x)dx, $$

<footnote n="2">²https://gitlab.inria.fr/ozahm/optimal-riemannian-metric-for-poincare-inequalities. git</footnote>
As in Equation (31), if the multiplicity of $\lambda_{N,2}^{V^2}$ is one, then $J_N$ admits a gradient at $V$ which is given by

$$ \nabla J_N(V) = \frac{G_N V + V G_N - 2 J_N(V) V}{\int \|V\|_F^2 d\mu}, $$ 
with $G_N = \frac{\left(\sum_{i=1}^N (U_2^{V^2})_i \nabla \phi_i\right) \left(\sum_{i=1}^N (U_2^{V^2})_i \nabla \phi_i\right)^\top}{(U_2^{V^2})^\top M (U_2^{V^2})} \quad (39) $$

Thus, to compute $\nabla J_N(V)$, one needs first to assemble the finite element matrices $K^{V^2}$ and $M$ and to compute the second smallest generalized eigenpair $(\lambda_2^{V^2}, U_2^{V^2})$ of $(K^{V^2}, M)$. Then, $\nabla J_N(V)$ is assembled according to (39). Notice that, because $\phi_i$ are piecewise affine, the matrix field $G_N$ is piecewise constant on each element $\omega_i$ of the mesh. In addition, if the field $V$ is also piecewise constant of the form of

$$ V(x) = \sum_{m=1}^{M} V_m 1_{\omega_m}(x), \quad \text{where } V_m \in \mathbb{R}^{d \times d}, \quad (40) $$

then $\nabla J_N(V)$ as in (39) is also piecewise constant. Here, $1_{\omega_m}$ is the indicator function of the $m$-th element $\omega_m$. Therefore, any of the iterates of the gradient algorithms in Equations (32)–(37) remains piecewise constant, and so is the resulting symmetric positive field $W$ given by (28). The resulting algorithm is summarized in Algorithm 1.

# 6 Numerical experiments

We demonstrate the proposed method on three benchmarks, each of them in dimension $d=2$. Code for reproducing all of these numerical results is available online<fnref n="2" />. The first is a Gaussian mixture of the form of

$$ d\mu_1(x) \propto \sum_{i=1}^{3} \exp \left( - \frac{\|x - x_i\|^2}{\sigma^2} \right) dx, $$

where $\sigma^2 = 0.025$, $x_1 = (0, 0.5)$, $x_2 = (\sqrt{3}/4, \sqrt{3}/4)$ and $x_3 = (-\sqrt{3}/4, \sqrt{3}/4)$. The second distribution is concentrated on a ring of radius $r$ and is defined by

$$ d\mu_2(x) \propto \exp \left( - \frac{(\|x\| - r)^2}{\sigma^2} \right) dx, $$

where $\sigma^2 = 0.0032$ and $r = 0.65$. The third is a uniform measure on a $H$-shaped domain $\Omega \subset \mathbb{R}^2$, meaning

$$ d\mu_3(x) \propto 1_{\Omega}(x) dx. $$

The fourth is a slight modification of $\mu_3$ which consists in enlarging the support of $\mu_3$ to the convex hull $\text{conv}(\Omega)$ of the $H$-shaped $\Omega$ as follow

$$ d\mu_4(x) \propto d\mu_3(x) + \varepsilon 1_{\text{conv}(\Omega)}(x) dx, $$

<footnote n="2">²https://gitlab.inria.fr/ozahm/optimal-riemannian-metric-for-poincare-inequalities. git</footnote>
**Algorithm 1** Optimal metric using finite element method and momentum algorithm

1: **Input:** measure $\mu$, finite element partition supp($\mu$) = $\omega_1 \cup \dots \cup \omega_M$, step size $\rho > 0$ and momentum parameter $0 \le \alpha < 1$.

2: Define the finite element space $V_N = \text{span}\{\phi_1, \dots, \phi_N\}$ of piecewise affine functions.

3: For all $i, j \le N$ and $m \le M$, precompute

$$ \mu_m = \int_{\omega_m} d\mu \quad \text{and} \quad M_{i,j} = \int \phi_i \phi_j d\mu \quad \text{and} \quad K_{i,j,m} = \int_{\omega_m} \nabla \phi_i \nabla \phi_j^\top d\mu \in \mathbb{R}^{d \times d}. $$

4: For all $m \le M$, initialize $V_m^{(0)} = I_d/d$ and $D_m^{(0)} = 0$.

5: Set $k = 0$.

6: **while** not converged **do**

7: For all $m \le M$, compute the metric (see Equation (28))

$$ W_m^{(k)} = (V_m^{(k)})^2 \frac{\text{tr}(\text{Cov}_\mu)}{\sum_{m=1}^{M} \|V_m^{(k)}\|_F^2 \mu_m} $$

8: Assemble the matrix $K^{(k)} \in \mathbb{R}^{N \times N}$ given by (see Equation (38))

$$ K_{i,j}^{(k)} = \sum_{i=1}^{M} \text{tr}(K_{i,j,m} W_m^{(k)}) $$

9: Compute the smallest non-zero generalized eigenvalue $(\lambda^{(k)}, U^{(k)}) \in \mathbb{R}_{>0} \times \mathbb{R}^N$ solution to

$$ K^{(k)} U^{(k)} = \lambda^{(k)} M U^{(k)} $$

10: For all $m \le M$, compute (see Equation (39))

$$ G_m^{(k)} = \frac{\sum_{i,j=1}^{N} U_i^{(k)} U_j^{(k)} K_{i,j,m}}{\sum_{i,j=1}^{N} U_i^{(k)} U_j^{(k)} M_{i,j}} $$

and then $\nabla J_m^{(k)} = G_m^{(k)} V_m^{(k)} + V_m^{(k)} G_m^{(k)} - 2\lambda^{(k)} V_m^{(k)}$

11: Compute the increment $D_m^{(k+1)} = \alpha D_m^{(k)} + \rho \nabla J_m^{(k)}$ (see Equation (34))

12: Update $V_m^{(k+1/2)} = V_m^{(k)} + D_m^{(k+1)}$ (see Equation (35))

13: Normalize $V_m^{(k+1)} = V_m^{(k+1/2)}/(\sum_{m'=1}^M \|V_{m'}^{(k)}\|_F^2 \mu_{m'})^{1/2}$ (see Equation (36))

14: $k \leftarrow k + 1$

15: **end while**

16: **Output:** Metric $W^{(k)} = \sum_{m=1}^M W_m^{(k)} 1_{\omega_m}$ and the associated Poincaré constant $C(\mu, W^{(k)}) = 1/\lambda^{(k)}$.

where $\varepsilon = 10^{-7}$. The goal is to create a measure which is close to $\mu_3$ while being supported on a convex domain. That way, $\mu_4$ satisfies the (sufficient) conditions for being a moment measure with a smooth enough convex potential $\varphi$ so that Theorem 2.1 applies, see [11]. This will be illustrated later on Figure 4.

The measures $\mu_1, \dots, \mu_4$ as well as the finite element discretizations of their supports,
**Algorithm 1** Optimal metric using finite element method and momentum algorithm

1: **Input:** measure $\mu$, finite element partition supp($\mu$) = $\omega_1 \cup \dots \cup \omega_M$, step size $\rho > 0$ and momentum parameter $0 \le \alpha < 1$.

2: Define the finite element space $V_N = \text{span}\{\phi_1, \dots, \phi_N\}$ of piecewise affine functions.

3: For all $i, j \le N$ and $m \le M$, precompute

$$ \mu_m = \int_{\omega_m} d\mu \quad \text{and} \quad M_{i,j} = \int \phi_i \phi_j d\mu \quad \text{and} \quad K_{i,j,m} = \int_{\omega_m} \nabla \phi_i \nabla \phi_j^\top d\mu \in \mathbb{R}^{d \times d}. $$

4: For all $m \le M$, initialize $V_m^{(0)} = I_d/d$ and $D_m^{(0)} = 0$.

5: Set $k = 0$.

6: **while** not converged **do**

7: For all $m \le M$, compute the metric (see Equation (28))

$$ W_m^{(k)} = (V_m^{(k)})^2 \frac{\text{tr}(\text{Cov}_\mu)}{\sum_{m=1}^M \|V_m^{(k)}\|_F^2 \mu_m} $$

8: Assemble the matrix $K^{(k)} \in \mathbb{R}^{N \times N}$ given by (see Equation (38))

$$ K_{i,j}^{(k)} = \sum_{i=1}^{M} \text{tr}(K_{i,j,m} W_m^{(k)}) $$

9: Compute the smallest non-zero generalized eigenvalue $(\lambda^{(k)}, U^{(k)}) \in \mathbb{R}_{>0} \times \mathbb{R}^N$ solution to

$$ K^{(k)} U^{(k)} = \lambda^{(k)} M U^{(k)} $$

10: For all $m \le M$, compute (see Equation (39))

$$ G_m^{(k)} = \frac{\sum_{i,j=1}^{N} U_i^{(k)} U_j^{(k)} K_{i,j,m}}{\sum_{i,j=1}^{N} U_i^{(k)} U_j^{(k)} M_{i,j}} $$

and then $\nabla J_m^{(k)} = G_m^{(k)} V_m^{(k)} + V_m^{(k)} G_m^{(k)} - 2\lambda^{(k)} V_m^{(k)}$

11: Compute the increment $D_m^{(k+1)} = \alpha D_m^{(k)} + \rho \nabla J_m^{(k)}$ (see Equation (34))

12: Update $V_m^{(k+1/2)} = V_m^{(k)} + D_m^{(k+1)}$ (see Equation (35))

13: Normalize $V_m^{(k+1)} = V_m^{(k+1/2)} / (\sum_{m'=1}^M \|V_{m'}^{(k)}\|_F^2 \mu_{m'})^{1/2}$ (see Equation (36))

14: $k \leftarrow k + 1$

15: **end while**

16: **Output:** Metric $W^{(k)} = \sum_{m=1}^M W_m^{(k)} 1_{\omega_m}$ and the associated Poincaré constant $C(\mu, W^{(k)}) = 1/\lambda^{(k)}$.

where $\varepsilon = 10^{-7}$. The goal is to create a measure which is close to $\mu_3$ while being supported on a convex domain. That way, $\mu_4$ satisfies the (sufficient) conditions for being a moment measure with a smooth enough convex potential $\varphi$ so that Theorem 2.1 applies, see [11]. This will be illustrated later on Figure 4.

The measures $\mu_1, \dots, \mu_4$ as well as the finite element discretizations of their supports,
<figure><img src="image_4.png" /><figcaption>Figure 1: The three benchmark measures and the finite element discretization of their support (cyan indicates low probability and purple high probability). The number of elements in the mesh are respectively $M_1 = 5718$, $M_2 = 10256$, $M_3 = 4610$ and $M_4 = 5742$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2: (Tri-modal $\mu_1$) Evaluation of the five first nonzero eigenvalues of the diffusion operator $\mathcal{L}_{\mu}^{W}$ during the iteration process of the gradient-ascent methods, with $\rho = 0.01$ and $\alpha = 0$ (Gradient ascent), $\alpha = 0.5$ (Momentum) and $\alpha_k = 1 - 3/(5 + k)$ (NAG).</figcaption></figure>

are represented on Figure 1. The Poincaré constant of these measures are respectively

$$C(\mu_1) = 61.79, \quad C(\mu_2) = 1.971, \quad C(\mu_3) = 15.01, \quad C(\mu_4) = 15.01.$$

## 6.1 Gradient ascent methods on the tri-modal measure

Figure 2 reports the evolution of the spectrum of $\mathcal{L}_{\mu_1}^{W}$ during the optimization process. We observe that all the gradient-ascent methods we considered yield $\lambda_1 \ge 0.9$ (and therefore $C(\mu, W) \le 1.12$) after 15 iterations. Notice the oscillations between of the two first nonzero eigenvalues $\lambda_2$ and $\lambda_3$, which indicates that those eigenvalues are constantly crossing each other during the iterations. When approaching optimality $\lambda_2, \lambda_3 \to 1$, the NAG yields more stable behaviour of the iterates compared to the alternative methods. Running the NAG algorithm for $k = 100$ iterations, we obtain $\lambda_2^{(k)} = 0.9998$, showing that one can bring the Poincaré constant $C(\mu, W^{(k)})$ arbitrarily close to 1.

Figure 3 represents the first 5 nontrivial eigenvectors of $\mathcal{L}_{\mu_1}^{W}$ for the (initial) constant metric $W^{(0)} = I_d \operatorname{tr}(\operatorname{Cov}_{\mu})/d$ and the (final) metric $W^{(k)}$ obtained using NAG with $k =$
<figure><img src="image_3.png" /><figcaption>Figure 1: The three benchmark measures and the finite element discretization of their support (cyan indicates low probability and purple high probability). The number of elements in the mesh are respectively $M_1 = 5718$, $M_2 = 10256$, $M_3 = 4610$ and $M_4 = 5742$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: (Tri-modal $\mu_1$) Evaluation of the five first nonzero eigenvalues of the diffusion operator $\mathcal{L}_{\mu}^W$ during the iteration process of the gradient-ascent methods, with $\rho = 0.01$ and $\alpha = 0$ (Gradient ascent), $\alpha = 0.5$ (Momentum) and $\alpha_k = 1 - 3/(5+k)$ (NAG).</figcaption></figure>

are represented on Figure 1. The Poincaré constant of these measures are respectively

$$C(\mu_1) = 61.79, \quad C(\mu_2) = 1.971, \quad C(\mu_3) = 15.01, \quad C(\mu_4) = 15.01.$$

## 6.1 Gradient ascent methods on the tri-modal measure

Figure 2 reports the evolution of the spectrum of $\mathcal{L}_{\mu_1}^W$ during the optimization process. We observe that all the gradient-ascent methods we considered yield $\lambda_1 \ge 0.9$ (and therefore $C(\mu, W) \le 1.12$) after 15 iterations. Notice the oscillations between of the two first nonzero eigenvalues $\lambda_2$ and $\lambda_3$, which indicates that those eigenvalues are constantly crossing each other during the iterations. When approaching optimality $\lambda_2, \lambda_3 \to 1$, the NAG yields more stable behaviour of the iterates compared to the alternative methods. Running the NAG algorithm for $k = 100$ iterations, we obtain $\lambda_2^{(k)} = 0.9998$, showing that one can bring the Poincaré constant $C(\mu, W^{(k)})$ arbitrarily close to 1.

Figure 3 represents the first 5 nontrivial eigenvectors of $\mathcal{L}_{\mu_1}^W$ for the (initial) constant metric $W^{(0)} = I_d \operatorname{tr}(\operatorname{Cov}_{\mu})/d$ and the (final) metric $W^{(k)}$ obtained using NAG with $k =$
<figure><img src="image_6.png" /><figcaption>Figure 3: (Tri-modal $\mu_1$) Eigenvectors $u_2^W, \dots, u_6^W$ of the diffusion operator $\mathcal{L}_{\mu_1}^W$ for the constant metric $W^{(0)} = I_d \frac{\text{tr}(\text{Cov}_{\mu_1})}{d}$ (top row) and the computed metric $W^{(k)}$ obtained using $k = 100$ iterations of NAG.</figcaption></figure>

100 iterations. We observe that the first two eigenvectors are affine functions which, according to Theorem 3.1, reflects the optimality of $W^{(k)}$.

## 6.2 Optimal metrics

Figure 4 represents the optimal metric $W$ on the four considered test cases. To represent a field of $2 \times 2$ positive matrices, we plot the trace of $W$ (in log-scale to increase visibility) and, to represent the anisotropy in the field $W$, we plot the ellipses whose main directions are the eigenvectors of $W(x)$.

For the tri-modal measure $\mu_1$, we see on Figure 4a that $W$ is low and isotropic around the three modes, and significantly high and anisotripic across the modes. For the ring measure $\mu_2$, we observe on Figure 4b that $W$ admits a singularity at the origin ($\text{tr}(W) \approx 10^{10}$) and is anisotropic where $\mu_2$ is high (the ring of radius $r$).

For the H-shaped measure $\mu_3$, Figure 4c indicates a high value of $W$ in the “bridge” connecting the left and right parts of the domain. It is important to notice that, for this benchmark, the resulting Poincaré constant is $C(\mu_3, W) = 0.9166$ which is significantly far from 1. A credible explanation for this is that $\mu_3$ is not a moment measure, and so Theorem 2.1 does not apply. In contrast, the measure $\mu_4 \approx \mu_3$ has a convex support, and our algorithm yields a Poincaré constant $C(\mu_4, W^{(k)}) = 0.9989$ which can be made arbitrarily close to 1. We notice that in the region of low-probability ($d\mu_4 \approx \varepsilon dx$), the metric $W$ degenerates ($\text{tr}(W) \approx \varepsilon^{-1}$) and is highly anisotropic.

## 6.3 Application to the unadjusted Langevin algorithm (ULA)

We consider now the Euler-Maruyama discretization of the Riemannian Langevin dynamic $dX_t = (\text{div}(W) - W\nabla V)dt + \sqrt{2W}dB_t$ as in (5), which writes

$$X_{n+1} = X_n + (\text{div}W(X_n) - W(X_n)\nabla V(X_n))\Delta t + \sqrt{2\Delta tW(X_n)}Z_n, \quad (41)$$
<figure><img src="image_5.png" /><figcaption>Figure 3: (Tri-modal $\mu_1$) Eigenvectors $u_2^W, \dots, u_6^W$ of the diffusion operator $\mathcal{L}_{\mu_1}^W$ for the constant metric $W^{(0)} = I_d \frac{\mathrm{tr}(\mathrm{Cov}_{\mu_1})}{d}$ (top row) and the computed metric $W^{(k)}$ obtained using $k = 100$ iterations of NAG.</figcaption></figure>

100 iterations. We observe that the first two eigenvectors are affine functions which, according to Theorem 3.1, reflects the optimality of $W^{(k)}$.

## 6.2 Optimal metrics

Figure 4 represents the optimal metric $W$ on the four considered test cases. To represent a field of $2 \times 2$ positive matrices, we plot the trace of $W$ (in log-scale to increase visibility) and, to represent the anisotropy in the field $W$, we plot the ellipses whose main directions are the eigenvectors of $W(x)$.

For the tri-modal measure $\mu_1$, we see on Figure 4a that $W$ is low and isotropic around the three modes, and significantly high and anisotripic across the modes. For the ring measure $\mu_2$, we observe on Figure 4b that $W$ admits a singularity at the origin ($\mathrm{tr}(W) \approx 10^{10}$) and is anisotropic where $\mu_2$ is high (the ring of radius $r$).

For the H-shaped measure $\mu_3$, Figure 4c indicates a high value of $W$ in the “bridge” connecting the left and right parts of the domain. It is important to notice that, for this benchmark, the resulting Poincaré constant is $C(\mu_3, W) = 0.9166$ which is significantly far from 1. A credible explanation for this is that $\mu_3$ is not a moment measure, and so Theorem 2.1 does not apply. In contrast, the measure $\mu_4 \approx \mu_3$ has a convex support, and our algorithm yields a Poincaré constant $C(\mu_4, W^{(k)}) = 0.9989$ which can be made arbitrarily close to 1. We notice that in the region of low-probability ($d\mu_4 \approx \varepsilon dx$), the metric $W$ degenerates ($\mathrm{tr}(W) \approx \varepsilon^{-1}$) and is highly anisotropic.

## 6.3 Application to the unadjusted Langevin algorithm (ULA)

We consider now the Euler-Maruyama discretization of the Riemannian Langevin dynamic $dX_t = (\mathrm{div}(W) - W\nabla V)dt + \sqrt{2W}dB_t$ as in (5), which writes

$$X_{n+1} = X_n + (\mathrm{div}W(X_n) - W(X_n)\nabla V(X_n))\Delta t + \sqrt{2\Delta tW(X_n)}Z_n, \quad (41)$$
<figure><img src="image_7.png" /><figcaption>Figure 4: Metric $W^{(k)}$ computed with $k = 100$ iterations of NAG. The colorbar corresponds to $\log_{10}(\operatorname{tr}(W^{(k)}))$, which is the amplitude of $W^{(k)}$ in log-scale. The ellipses represent the local anisotropy of $W^{(k)}$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Metric $W^{(k)}$ computed with $k = 100$ iterations of NAG. The colorbar corresponds to $\log_{10}(\text{tr}(W^{(k)}))$, which is the amplitude of $W^{(k)}$ in log-scale. The ellipses represent the local anisotropy of $W^{(k)}$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 5: Drift $(\operatorname{div}(W) - W \nabla V_i)$ for the four considered measures $d\mu_i \propto \exp(-V_i)dx$, $i = 1, \dots, 4$, with either the constant metric $W \propto I_d$ (Figures 5a and 5d) or the optimal metric $W$.</figcaption></figure>

for some time step $\Delta t > 0$. Here, $Z_1, Z_2, \dots$ are independent random vector drawn from $N(0, I_d)$. While Figure 4 already represents the diffusion term $\sqrt{2W}$ of the Riemannian Langevin SDE, Figure 5 represents the drift term $\operatorname{div}(W) - W\nabla V$. From it, we can observe how the new drift field connects different high probability region. In particular, from Figure 5a, we observe the non-preconditioned drift $\nabla V$ points towards the centers of the mixture components. While such drift pushes the particle $X_t$ in high probability regions, it also prevent jumps from one mode to another. In comparison, the new drift field on Figure 5b connects the different mixture components by pushing the particles to escape the modes, which intuitively improves the mixing speed of the SDE. Likewise, on Figure 5d, the drift $\nabla V$ only leads towards the high probability region, but not circulating along the ring. The drift field on Figure 5e actually suggests leaving the high probability region to the center of the ring, where it will be directed to other area of the ring. Regarding the uniform measure $\mu_3$ for which $\nabla V = 0$, we see on Figure 5c that the optimal metric $W$ creates a drift $\operatorname{div}W$ which pushes the particles towards the bottleneck of distribution so they can get to the other side. Regarding $\mu_4$, the cyan area on Figure 5f provides additional leeway for the mixing to happen.

Finally, Figure 6 represents the trace of 5000 iterations of the ULA algorithm (41). We observe that the unpreconditioned ULA (top row) has a large bias compared to the preconditioned ULA (middle row). In particular, the preconditioned ULA is more
<figure><img src="image_7.png" /><figcaption>Figure 5: Drift $(\operatorname{div}(W) - W \nabla V_i)$ for the four considered measures $d\mu_i \propto \exp(-V_i)dx$, $i = 1, \dots, 4$, with either the constant metric $W \propto I_d$ (Figures 5a and 5d) or the optimal metric $W$.</figcaption></figure>

for some time step $\Delta t > 0$. Here, $Z_1, Z_2, \dots$ are independent random vector drawn from $N(0, I_d)$. While Figure 4 already represents the diffusion term $\sqrt{2W}$ of the Riemannian Langevin SDE, Figure 5 represents the drift term $\operatorname{div}(W) - W\nabla V$. From it, we can observe how the new drift field connects different high probability region. In particular, from Figure 5a, we observe the non-preconditioned drift $\nabla V$ points towards the centers of the mixture components. While such drift pushes the particle $X_t$ in high probability regions, it also prevent jumps from one mode to another. In comparison, the new drift field on Figure 5b connects the different mixture components by pushing the particles to escape the modes, which intuitively improves the mixing speed of the SDE. Likewise, on Figure 5d, the drift $\nabla V$ only leads towards the high probability region, but not circulating along the ring. The drift field on Figure 5e actually suggests leaving the high probability region to the center of the ring, where it will be directed to other area of the ring. Regarding the uniform measure $\mu_3$ for which $\nabla V = 0$, we see on Figure 5c that the optimal metric $W$ creates a drift $\operatorname{div}W$ which pushes the particles towards the bottleneck of distribution so they can get to the other side. Regarding $\mu_4$, the cyan area on Figure 5f provides additional leeway for the mixing to happen.

Finally, Figure 6 represents the trace of 5000 iterations of the ULA algorithm (41). We observe that the unpreconditioned ULA (top row) has a large bias compared to the preconditioned ULA (middle row). In particular, the preconditioned ULA is more
<figure><img src="image_9.png" /><figcaption>Figure 6: Trace of 5000 iterations of ULA (41) for different $\Delta t$ using either the unpreconditioned ULA ($W \propto I_d$, top row) or the preconditioned ULA (optimal $W$, bottom row).</figcaption></figure>

robust with respect to large time steps $\Delta t$ compared to the non-preconditioned ULA, see for instance the difference between Figure 6b and Figure 6f. For the ring measure $\mu_2$, having an anisotropic diffusion term (remember Figure 4b) permits a better mixing of the particles.

# 7 Conclusion

This paper considers adding a Riemannian metric in the Poincaré inequality, and investigates the notion of optimal metric by minimizing the associated Poincaré constant. Several important properties of the optimal metric has been obtained and a gradient-based numerical algorithm is further proposed to compute it. It is applied in four different 2D distributions, where the optimal matrix field reveals the geometric structure of the distributions. These metrics are then exploited to precondition the Langevin SDE.

While this work has shown various important properties and benefits of the optimal metric, much is left to be done, especially in the computational side. Next we explain a few possible directions for further developments:

* Although Theorem 2.1 establishes the existence of an optimal metric, the issue of its uniqueness remains unresolved. Given that optimal metrics must be positive-semidefinite Stein kernels (as per Theorem 3.2), one approach to investigating uniqueness would be to ascertain whether positive-semidefinite Stein kernels are indeed unique or not.

* Stein kernels have proven effective in constructing statistical tests to discern whether data samples originate from a particular distribution [39]. Assessing whether the positive-semidefinite Stein kernel we obtain via the optimal metric improves upon existing statistical tests is of natural interest.
<figure><img src="image_8.png" /><figcaption>Figure 6: Trace of 5000 iterations of ULA (41) for different $\Delta t$ using either the unpreconditioned ULA ($W \propto I_d$, top row) or the preconditioned ULA (optimal $W$, bottom row).</figcaption></figure>

robust with respect to large time steps $\Delta t$ compared to the non-preconditioned ULA, see for instance the difference between Figure 6b and Figure 6f. For the ring measure $\mu_2$, having an anisotropic diffusion term (remember Figure 4b) permits a better mixing of the particles.

# 7 Conclusion

This paper considers adding a Riemannian metric in the Poincaré inequality, and investigates the notion of optimal metric by minimizing the associated Poincaré constant. Several important properties of the optimal metric has been obtained and a gradient-based numerical algorithm is further proposed to compute it. It is applied in four different 2D distributions, where the optimal matrix field reveals the geometric structure of the distributions. These metrics are then exploited to precondition the Langevin SDE.

While this work has shown various important properties and benefits of the optimal metric, much is left to be done, especially in the computational side. Next we explain a few possible directions for further developments:

* Although Theorem 2.1 establishes the existence of an optimal metric, the issue of its uniqueness remains unresolved. Given that optimal metrics must be positive-semidefinite Stein kernels (as per Theorem 3.2), one approach to investigating uniqueness would be to ascertain whether positive-semidefinite Stein kernels are indeed unique or not.

* Stein kernels have proven effective in constructing statistical tests to discern whether data samples originate from a particular distribution [39]. Assessing whether the positive-semidefinite Stein kernel we obtain via the optimal metric improves upon existing statistical tests is of natural interest.
* While algorithm 1 permits to accurately compute the optimal metric, its application to large scale problems requires further investigations. In particular, the finite element methods we employed becomes prominently expensive if the underlying dimension $d$ is high, which makes it inapplicable for high-dimension problems. Kernel based as in [45] or neural network based methods in the spirit of [10] may be able to face the dimensionality issue.

* An additional concern regarding Algorithm 1 is its dependency on global knowledge of $\mu$, which may not be readily available beforehand. Algorithms capable of estimating $W$ using available data (e.g. few samples or moments from $\mu$) or dynamically adapting to new information (e.g. during an MCMC procedure which samples from $\mu$) would present a more compelling approach.

## Acknowledgments

The authors appreciate discussion of some of the results with Louis Chen Hsiao Yun, Ramon van Handel and Arthur Vidard.

## Funding

T Cui acknowledges support from the Australian Research Council under the grant DP210103092. Xin T. Tong's research is supported by Singapore MOE grant Tier-1-A-8000459-00-00. O. Zahm acknowledges support from the ANR JCJC project MODENA (ANR-21-CE46-0006-01).

## A Proof of Proposition 1

The Fokker-Planck equation satisfied by the density $\mu_t \propto \exp(-V_t)$ of $X_t$ solution to (5) writes

$$ 
\begin{aligned}
\partial_t \mu_t &= -\nabla \cdot ((\mathrm{div}(W) - W \nabla V)\mu_t) + \nabla \cdot (\mathrm{div}(W \mu_t)) \\
&= -\nabla \cdot ((\mathrm{div}(W) - W \nabla V)\mu_t) + \nabla \cdot (\mathrm{div}(W)\mu_t + W \nabla \mu_t) \\
&= -\nabla \cdot (-W \nabla V)\mu_t - W \nabla \mu_t \\
&= -\nabla \cdot (\mu_t W (\nabla V_t - \nabla V)).
\end{aligned}
\quad (42) $$

Because $\chi^2(\mu_t, \mu) = \int (\frac{\mu_t}{\mu})^2 d\mu - 1$ we can write

$$ 
\begin{aligned}
\partial_t \chi^2(\mu_t, \mu) &= 2 \int \partial_t \left(\frac{\mu_t}{\mu}\right) \mu_t dx \\
&\stackrel{(42)}{=} -2 \int \nabla \cdot (\mu_t W (\nabla V_t - \nabla V)) \left(\frac{\mu_t}{\mu}\right) dx \\
&= 2 \int (\mu_t W (\nabla V_t - \nabla V))^\top \nabla \left(\frac{\mu_t}{\mu}\right) dx \\
&= -2 \int \mu_t \nabla \log \left(\frac{\mu_t}{\mu}\right)^\top W \nabla \left(\frac{\mu_t}{\mu}\right) dx \\
&= -2 \int \nabla \left(\frac{\mu_t}{\mu}\right)^\top W \nabla \left(\frac{\mu_t}{\mu}\right) d\mu
\end{aligned}
\quad (43) $$
* While algorithm 1 permits to accurately compute the optimal metric, its application to large scale problems requires further investigations. In particular, the finite element methods we employed becomes prominently expensive if the underlying dimension $d$ is high, which makes it inapplicable for high-dimension problems. Kernel based as in [45] or neural network based methods in the spirit of [10] may be able to face the dimensionality issue.

* An additional concern regarding Algorithm 1 is its dependency on global knowledge of $\mu$, which may not be readily available beforehand. Algorithms capable of estimating $W$ using available data (e.g. few samples or moments from $\mu$) or dynamically adapting to new information (e.g. during an MCMC procedure which samples from $\mu$) would present a more compelling approach.

## Acknowledgments

The authors appreciate discussion of some of the results with Louis Chen Hsiao Yun, Ramon van Handel and Arthur Vidard.

## Funding

T Cui acknowledges support from the Australian Research Council under the grant DP210103092. Xin T. Tong's research is supported by Singapore MOE grant Tier-1-A-8000459-00-00. O. Zahm acknowledges support from the ANR JCJC project MODENA (ANR-21-CE46-0006-01).

## A Proof of Proposition 1

The Fokker-Planck equation satisfied by the density $\mu_t \propto \exp(-V_t)$ of $X_t$ solution to (5) writes

$$ 
\begin{aligned}
\partial_t \mu_t &= -\nabla \cdot ((\mathrm{div}(W) - W \nabla V)\mu_t) + \nabla \cdot (\mathrm{div}(W \mu_t)) \\
&= -\nabla \cdot ((\mathrm{div}(W) - W \nabla V)\mu_t) + \nabla \cdot (\mathrm{div}(W)\mu_t + W \nabla \mu_t) \\
&= -\nabla \cdot (-W \nabla V)\mu_t - W \nabla \mu_t \\
&= -\nabla \cdot (\mu_t W (\nabla V_t - \nabla V)).
\end{aligned}
\quad (42) $$

Because $\chi^2(\mu_t, \mu) = \int (\frac{\mu_t}{\mu})^2 d\mu - 1$ we can write

$$ 
\begin{aligned}
\partial_t \chi^2(\mu_t, \mu) &= 2 \int \partial_t \left( \frac{\mu_t}{\mu} \right) \mu_t dx \\
&\stackrel{(42)}{=} -2 \int \nabla \cdot (\mu_t W (\nabla V_t - \nabla V)) \left( \frac{\mu_t}{\mu} \right) dx \\
&= 2 \int (\mu_t W (\nabla V_t - \nabla V))^\top \nabla \left( \frac{\mu_t}{\mu} \right) dx \\
&= -2 \int \mu_t \nabla \log \left( \frac{\mu_t}{\mu} \right)^\top W \nabla \left( \frac{\mu_t}{\mu} \right) dx \\
&= -2 \int \nabla \left( \frac{\mu_t}{\mu} \right)^\top W \nabla \left( \frac{\mu_t}{\mu} \right) d\mu
\end{aligned}
\quad (43) $$
To show that $C(\mu, W) \le C$ implies (6), we apply the Riemannian Poincaré inequality (4) to (43) to obtain

$$ \partial_t \chi^2(\mu_t, \mu) \stackrel{(4)}{\le} -\frac{2}{C} \mathrm{Var}_\mu(\mu_t/\mu) = -\frac{2}{C} \chi^2(\mu_t, \mu). \quad (44) $$

This shows that the function $t \mapsto e^{2t/C}\chi^2(\mu_t, \mu)$ is decreasing (because its derivative $2/Ce^{2t/C}\chi^2(\mu_t, \mu) + e^{2t/C}\partial_t\chi^2(\mu_t, \mu) \le 0$ is nonpositive) so that $\chi^2(\mu_0, \mu) \ge e^{2t/C}\chi^2(\mu_t, \mu)$, which is (6).

We now show that if (6) holds for any probability measure $\mu_0$ then $C(\mu, W) \le C$. Let $f_0$ be a positive function such that $\int f_0 d\mu = 1$. Letting $d\mu_0 = f_0 d\mu$, inequality (6) becomes $\chi^2(\mu_t, \mu) \le e^{-2t/C} \mathrm{Var}_\mu(f_0)$. A Taylor expansion around $t=0$ yields

$$ \chi^2(\mu_0, \mu) + \partial_t \chi^2(\mu_t, \mu)|_{t=0} t + \mathcal{O}(t^2) \le (1 - 2t/C) \mathrm{Var}_\mu(f_0) + \mathcal{O}(t^2), \quad (45) $$

for any $t \ge 0$. By (43) we have $\partial_t \chi^2(\mu_t, \mu)|_{t=0} = -2 \int \nabla f_0^\top W \nabla f_0 d\mu$ so that (45) yields

$$ \mathrm{Var}_\mu(f_0) \le C \int \nabla f_0^\top W \nabla f_0 d\mu, \quad (46) $$

for any positive function $f_0$ such that $\int f_0 d\mu = 1$. It remains to show that it also holds for any smooth function $f$ with $\int \nabla f^\top W \nabla f d\mu < \infty$. If such a $f$ is lower bounded, then $f_0 = (f - \inf f)/\alpha$ with $\alpha = \int f d\mu - \inf f$ is a positive function with $\int f_0 d\mu = 1$ and so (46) yields $\alpha^{-2} \mathrm{Var}_\mu(f) \le \alpha^{-2} C \int \nabla f^\top W \nabla f d\mu$. We conclude that any smooth lower-bounded $f$ satisfies the Riemannian Poincaré inequality. Finally, if $f$ is not lower bounded, we introduce the sequence of lower-bounded function $f_n = \max\{f; -n\}$ so that

$$ \mathrm{Var}_\mu(f_n) \le C \int \nabla f_n^\top W \nabla f_n d\mu. $$

Taking the limit when $n \to \infty$ we obtain $\mathrm{Var}_\mu(f) \le C \int \nabla f^\top W \nabla f d\mu$.

## B Proof of Corollary 2

The proof is a simple adaptation of the proof of Theorem 3.3 in [31]. By Theorem 3.2, $W_N$ defined in (14) is a Stein kernel so that

$$ \begin{align*} S_p(\mu_N|\gamma)^p &\le \int \left\| \mathbb{E} \left[ \left\| \sum_{i=1}^N \frac{W(X_i) - I_d}{N} \middle| \bar{X}_N \right\| \right\|_F^p d\mu_N \\ \text{(Jensen's inequality)} \quad &\le \int \mathbb{E} \left[ \left\| \sum_{i=1}^N \frac{W(X_i) - I_d}{N} \right\|_F^p \middle| \bar{X}_N \right] d\mu_N \\ &= \mathbb{E} \left[ \left\| \sum_{i=1}^N \frac{W(X_i) - I_d}{N} \right\|_F^p \right], \end{align*} $$
To show that $C(\mu, W) \le C$ implies (6), we apply the Riemannian Poincaré inequality (4) to (43) to obtain

$$ \partial_t \chi^2(\mu_t, \mu) \stackrel{(4)}{\le} -\frac{2}{C} \mathrm{Var}_\mu(\mu_t/\mu) = -\frac{2}{C} \chi^2(\mu_t, \mu). \quad (44) $$

This shows that the function $t \mapsto e^{2t/C}\chi^2(\mu_t, \mu)$ is decreasing (because its derivative $2/Ce^{2t/C}\chi^2(\mu_t, \mu) + e^{2t/C}\partial_t\chi^2(\mu_t, \mu) \le 0$ is nonpositive) so that $\chi^2(\mu_0, \mu) \ge e^{2t/C}\chi^2(\mu_t, \mu)$, which is (6).

We now show that if (6) holds for any probability measure $\mu_0$ then $C(\mu, W) \le C$. Let $f_0$ be a positive function such that $\int f_0 d\mu = 1$. Letting $d\mu_0 = f_0 d\mu$, inequality (6) becomes $\chi^2(\mu_t, \mu) \le e^{-2t/C} \mathrm{Var}_\mu(f_0)$. A Taylor expansion around $t=0$ yields

$$ \chi^2(\mu_0, \mu) + \partial_t \chi^2(\mu_t, \mu)|_{t=0} t + \mathcal{O}(t^2) \le (1 - 2t/C) \mathrm{Var}_\mu(f_0) + \mathcal{O}(t^2), \quad (45) $$

for any $t \ge 0$. By (43) we have $\partial_t \chi^2(\mu_t, \mu)|_{t=0} = -2 \int \nabla f_0^\top W \nabla f_0 d\mu$ so that (45) yields

$$ \mathrm{Var}_\mu(f_0) \le C \int \nabla f_0^\top W \nabla f_0 d\mu, \quad (46) $$

for any positive function $f_0$ such that $\int f_0 d\mu = 1$. It remains to show that it also holds for any smooth function $f$ with $\int \nabla f^\top W \nabla f d\mu < \infty$. If such a $f$ is lower bounded, then $f_0 = (f - \inf f)/\alpha$ with $\alpha = \int f d\mu - \inf f$ is a positive function with $\int f_0 d\mu = 1$ and so (46) yields $\alpha^{-2} \mathrm{Var}_\mu(f) \le \alpha^{-2} C \int \nabla f^\top W \nabla f d\mu$. We conclude that any smooth lower-bounded $f$ satisfies the Riemannian Poincaré inequality. Finally, if $f$ is not lower bounded, we introduce the sequence of lower-bounded function $f_n = \max\{f; -n\}$ so that

$$ \mathrm{Var}_\mu(f_n) \le C \int \nabla f_n^\top W \nabla f_n d\mu. $$

Taking the limit when $n \to \infty$ we obtain $\mathrm{Var}_\mu(f) \le C \int \nabla f^\top W \nabla f d\mu$.

## B Proof of Corollary 2

The proof is a simple adaptation of the proof of Theorem 3.3 in [31]. By Theorem 3.2, $W_N$ defined in (14) is a Stein kernel so that

$$ \begin{align*} S_p(\mu_N|\gamma)^p &\le \int \left\| \mathbb{E} \left[ \left\| \sum_{i=1}^N \frac{W(X_i) - I_d}{N} \middle| \bar{X}_N \right\| \right\|_F^p d\mu_N \\ \text{(Jensen's inequality)} \quad &\le \int \mathbb{E} \left[ \left\| \sum_{i=1}^N \frac{W(X_i) - I_d}{N} \right\|_F^p \middle| \bar{X}_N \right] d\mu_N \\ &= \mathbb{E} \left[ \left\| \sum_{i=1}^N \frac{W(X_i) - I_d}{N} \right\|_F^p \right], \end{align*} $$
Letting $Y_{ikl} = W_{kl}(X_i) - \delta_{kl}$, we deduce

$$
\begin{align*}
S_p(\mu_N|\gamma)^p &\le \mathbb{E}\left[\left(\frac{1}{d^2}\sum_{kl=1}^d d^2\left(\frac{1}{N}\sum_{i=1}^N Y_{ikl}\right)^2\right)^{p/2}\right] \\
\text{(Jensen's inequality)} &\le \frac{1}{d^2}\sum_{kl=1}^d \mathbb{E}\left[\left(d^2\left(\frac{1}{N}\sum_{i=1}^N Y_{ikl}\right)^2\right)^{p/2}\right] \\
&= d^{p-2}N^{-p}\sum_{kl=1}^d \mathbb{E}\left[\left(\sum_{i=1}^N Y_{ikl}\right)^p\right] \\
\text{(Rosenthal)} &\le d^{p-2}N^{-p}\sum_{kl=1}^d K_p \max\left\{\sum_{i=1}^N \mathbb{E}[Y_{ikl}^p]; \left(\sum_{i=1}^N \mathbb{E}[Y_{ikl}^2]\right)^{p/2}\right\}
\end{align*}
$$

where we applied the Rosenthal inequality [35] in the last step. Here, $K_p$ is a constant depending only on $p$. On one side we have $\sum_{i=1}^N \mathbb{E}[Y_{ikl}^p] = N \mathbb{E}[(W_{kl}(X) - \delta_{kl})^p]$ and, on the other side, a Hölder inequality yields

$$ \left( \sum_{i=1}^{N} \mathbb{E}[Y_{ikl}^2] \right)^{p/2} = N^{p/2} \mathbb{E}[(W_{kl}(X) - \delta_{kl})^2]^{p/2} \le N^{p/2} \mathbb{E}[(W_{kl}(X) - \delta_{kl})^p]. $$

Thus we deduce

$$
\begin{align*}
S_p(\mu_N|\gamma)^p &\le d^{p-2}N^{-p} \sum_{kl=1}^d K_p \mathbb{E}[(W_{kl}(X) - \delta_{kl})^p] \max\{N; N^{p/2}\} \\
&= d^{p-2}N^{-p/2} K_p \int \sum_{kl=1}^d (W_{kl}(X) - \delta_{kl})^p d\mu
\end{align*}
$$

Together with Proposition 3.1 in [31], which establishes $W_p(\mu_N, \gamma) \le C'_p S_p(\mu_N|\gamma)$ for another constant $C'_p$ depending only on $p$, we conclude the proof.

## References

[1] G. ACOSTA AND R. DURÁN, *An optimal poincaré inequality in $l^1$ for convex domains*, Proceedings of the american mathematical society, 132 (2004), pp. 195–202.

[2] G. ALLAIRE, *Numerical analysis and optimization: an introduction to mathematical modelling and numerical simulation*, OUP Oxford, 2007.

[3] C. ANDRIEU, A. LEE, S. POWER, AND A. Q. WANG, *Comparison of markov chains via weak poincaré inequalities with application to pseudo-marginal mcmc*, The Annals of Statistics, 50 (2022), pp. 3592–3618.

[4] ——, *Poincaré inequalities for markov chains: a meeting with cheeger, lyapunov and metropolis*, arXiv preprint arXiv:2208.05239, (2022).
Letting $Y_{ikl} = W_{kl}(X_i) - \delta_{kl}$, we deduce

$$
\begin{align*}
S_p(\mu_N | \gamma)^p &\le \mathbb{E} \left[ \left( \frac{1}{d^2} \sum_{kl=1}^d d^2 \left( \frac{1}{N} \sum_{i=1}^N Y_{ikl} \right)^2 \right)^{p/2} \right] \\
\text{(Jensen's inequality)} &\le \frac{1}{d^2} \sum_{kl=1}^d \mathbb{E} \left[ \left( d^2 \left( \frac{1}{N} \sum_{i=1}^N Y_{ikl} \right)^2 \right)^{p/2} \right] \\
&= d^{p-2} N^{-p} \sum_{kl=1}^d \mathbb{E} \left[ \left( \sum_{i=1}^N Y_{ikl} \right)^p \right] \\
\text{(Rosenthal)} &\le d^{p-2} N^{-p} \sum_{kl=1}^d K_p \max \left\{ \sum_{i=1}^N \mathbb{E}[Y_{ikl}^p]; \left( \sum_{i=1}^N \mathbb{E}[Y_{ikl}^2] \right)^{p/2} \right\}
\end{align*}
$$

where we applied the Rosenthal inequality [35] in the last step. Here, $K_p$ is a constant depending only on $p$. On one side we have $\sum_{i=1}^N \mathbb{E}[Y_{ikl}^p] = N \mathbb{E}[(W_{kl}(X) - \delta_{kl})^p]$ and, on the other side, a Hölder inequality yields

$$ \left( \sum_{i=1}^{N} \mathbb{E}[Y_{ikl}^2] \right)^{p/2} = N^{p/2} \mathbb{E}[(W_{kl}(X) - \delta_{kl})^2]^{p/2} \le N^{p/2} \mathbb{E}[(W_{kl}(X) - \delta_{kl})^p]. $$

Thus we deduce

$$
\begin{align*}
S_p(\mu_N | \gamma)^p &\le d^{p-2} N^{-p} \sum_{kl=1}^d K_p \mathbb{E}[(W_{kl}(X) - \delta_{kl})^p] \max \{N; N^{p/2}\} \\
&= d^{p-2} N^{-p/2} K_p \int \sum_{kl=1}^d (W_{kl}(X) - \delta_{kl})^p d\mu
\end{align*}
$$

Together with Proposition 3.1 in [31], which establishes $W_p(\mu_N, \gamma) \le C'_p S_p(\mu_N | \gamma)$ for another constant $C'_p$ depending only on $p$, we conclude the proof.

## References

[1] G. ACOSTA AND R. DURÁN, *An optimal poincaré inequality in $l^1$ for convex domains*, Proceedings of the american mathematical society, 132 (2004), pp. 195–202.

[2] G. ALLAIRE, *Numerical analysis and optimization: an introduction to mathematical modelling and numerical simulation*, OUP Oxford, 2007.

[3] C. ANDRIEU, A. LEE, S. POWER, AND A. Q. WANG, *Comparison of markov chains via weak poincaré inequalities with application to pseudo-marginal mcmc*, The Annals of Statistics, 50 (2022), pp. 3592–3618.

[4] ——, *Poincaré inequalities for markov chains: a meeting with cheeger, lyapunov and metropolis*, arXiv preprint arXiv:2208.05239, (2022).
[5] D. BAKRY, F. BARTHE, P. CATTIAUX, AND A. GUILLIN, *A simple proof of the Poincaré inequality for a large class of probability measures*, Electronic Communications in Probability, 13 (2008), pp. 60 – 66.

[6] D. BAKRY AND M. ÉMERY, *Diffusions hypercontractives*, in Séminaire de Probabilités XIX 1983/84: Proceedings, Springer, 2006, pp. 177–206.

[7] D. BAKRY, I. GENTIL, M. LEDOUX, ET AL., *Analysis and geometry of Markov diffusion operators*, vol. 103, Springer, 2014.

[8] A. D. BARBOUR AND L. H. Y. CHEN, *An introduction to Stein's method*, vol. 4, World Scientific, 2005.

[9] M. BEBENDORF, *A note on the poincaré inequality for convex domains*, Zeitschrift für Analysis und ihre Anwendungen, 22 (2003), pp. 751–756.

[10] Z. BELKACEMI, M. BIANCIOTTO, H. MINOUX, T. LELIÈVRE, G. STOLTZ, AND P. GKEKA, *Autoencoders for dimensionality reduction in molecular dynamics: Collective variable dimension, biasing, and transition states*, The Journal of Chemical Physics, 159 (2023).

[11] R. J. BERMAN AND B. BERNDTSSON, *Real monge-ampère equations and kähler-ricci solitons on toric log fano varieties*, in Annales de la Faculté des sciences de Toulouse: Mathématiques, vol. 22, 2013, pp. 649–711.

[12] S. G. BOBKOV AND C. HOUDRÉ, *Variance of lipschitz functions and an isoperimetric problem for a class of product measures*, Bernoulli, (1996), pp. 249–255.

[13] S. G. BOBKOV AND M. LEDOUX, *From brunn-minkowski to brascamp-lieb and to logarithmic sobolev inequalities*, Geometric and Functional Analysis, 10 (2000), pp. 1028–1052.

[14] ——, *Weighted poincaré-type inequalities for cauchy and other convex measures*, The Annals of Probability, 37 (2009), pp. 403–427.

[15] M. BONNEFONT, A. JOULIN, AND Y. MA, *A note on spectral gap and weighted poincaré inequalities for some one-dimensional diffusions*, ESAIM: Probability and Statistics, 20 (2016), pp. 18–29.

[16] ——, *Spectral gap for spherically symmetric log-concave probability measures, and beyond*, Journal of Functional Analysis, 270 (2016), pp. 2456–2482.

[17] S. BOYD, S. P. BOYD, AND L. VANDENBERGHE, *Convex optimization*, Cambridge university press, 2004.

[18] H. J. BRASCAMP AND E. H. LIEB, *On extensions of the brunn-minkowski and prékopa-leindler theorems, including inequalities for log concave functions, and with an application to the diffusion equation*, Journal of functional analysis, 22 (1976), pp. 366–389.

[19] H. BREZIS AND H. BRÉZIS, *Functional analysis, Sobolev spaces and partial differential equations*, vol. 2, Springer, 2011.
[5] D. BAKRY, F. BARTHE, P. CATTIAUX, AND A. GUILLIN, *A simple proof of the Poincaré inequality for a large class of probability measures*, Electronic Communications in Probability, 13 (2008), pp. 60 – 66.

[6] D. BAKRY AND M. ÉMERY, *Diffusions hypercontractives*, in Séminaire de Probabilités XIX 1983/84: Proceedings, Springer, 2006, pp. 177–206.

[7] D. BAKRY, I. GENTIL, M. LEDOUX, ET AL., *Analysis and geometry of Markov diffusion operators*, vol. 103, Springer, 2014.

[8] A. D. BARBOUR AND L. H. Y. CHEN, *An introduction to Stein's method*, vol. 4, World Scientific, 2005.

[9] M. BEBENDORF, *A note on the poincaré inequality for convex domains*, Zeitschrift für Analysis und ihre Anwendungen, 22 (2003), pp. 751–756.

[10] Z. BELKACEMI, M. BIANCIOTTO, H. MINOUX, T. LELIÈVRE, G. STOLTZ, AND P. GKEKA, *Autoencoders for dimensionality reduction in molecular dynamics: Collective variable dimension, biasing, and transition states*, The Journal of Chemical Physics, 159 (2023).

[11] R. J. BERMAN AND B. BERNDTSSON, *Real monge-ampère equations and kähler-ricci solitons on toric log fano varieties*, in Annales de la Faculté des sciences de Toulouse: Mathématiques, vol. 22, 2013, pp. 649–711.

[12] S. G. BOBKOV AND C. HOUDRÉ, *Variance of lipschitz functions and an isoperimetric problem for a class of product measures*, Bernoulli, (1996), pp. 249–255.

[13] S. G. BOBKOV AND M. LEDOUX, *From brunn-minkowski to brascamp-lieb and to logarithmic sobolev inequalities*, Geometric and Functional Analysis, 10 (2000), pp. 1028–1052.

[14] ——, *Weighted poincaré-type inequalities for cauchy and other convex measures*, The Annals of Probability, 37 (2009), pp. 403–427.

[15] M. BONNEFONT, A. JOULIN, AND Y. MA, *A note on spectral gap and weighted poincaré inequalities for some one-dimensional diffusions*, ESAIM: Probability and Statistics, 20 (2016), pp. 18–29.

[16] ——, *Spectral gap for spherically symmetric log-concave probability measures, and beyond*, Journal of Functional Analysis, 270 (2016), pp. 2456–2482.

[17] S. BOYD, S. P. BOYD, AND L. VANDENBERGHE, *Convex optimization*, Cambridge university press, 2004.

[18] H. J. BRASCAMP AND E. H. LIEB, *On extensions of the brunn-minkowski and prékopa-leindler theorems, including inequalities for log concave functions, and with an application to the diffusion equation*, Journal of functional analysis, 22 (1976), pp. 366–389.

[19] H. BREZIS AND H. BRÉZIS, *Functional analysis, Sobolev spaces and partial differential equations*, vol. 2, Springer, 2011.
[20] T. CACOULLOS, *On upper and lower bounds for the variance of a function of a random variable*, The Annals of Probability, 10 (1982), pp. 799–809.

[21] T. CACOULLOS, V. PAPATHANASIOU, AND S. A. UTEV, *Variational inequalities with examples and an application to the central limit theorem*, The Annals of Probability, (1994), pp. 1607–1618.

[22] J. CHEEGER, *Differentiability of lipschitz functions on metric measure spaces*, Geometric & Functional Analysis GAFA, 9 (1999), pp. 428–517.

[23] L. H. CHEN, *Stein's method of normal approximation: Some recollections and reflections*, The Annals of Statistics, 49 (2021), pp. 1850–1863.

[24] X. CHENG, B. WANG, J. ZHANG, AND Y. ZHU, *Fast conditional mixing of mcmc algorithms for non-log-concave distributions*, Advances in Neural Information Processing Systems, 36 (2024).

[25] V. CHERNOZHUOKOV, D. CHETVERIKOV, K. KATO, AND Y. KOIKE, *Improved central limit theorem and bootstrap approximations in high dimensions*, The Annals of Statistics, 50 (2022), pp. 2562–2586.

[26] F. CLARKE, *Functional analysis, calculus of variations and optimal control*, vol. 264, Springer, 2013.

[27] D. CORDERO-ERAUSQUIN AND B. KLARTAG, *Moment measures*, Journal of Functional Analysis, 268 (2015), pp. 3834–3866.

[28] T. A. COURTADE, M. FATHI, AND A. PANANJADY, *Existence of stein kernels under a spectral gap, and discrepancy bounds*, in Annales de l'Institut Henri Poincaré, Probabilités et Statistiques, vol. 55, Institut Henri Poincaré, 2019, pp. 777–790.

[29] T. CUI, X. T. TONG, AND O. ZAHM, *Prior normalization for certified likelihood-informed subspace detection of bayesian inverse problems*, Inverse Problems, 38 (2022), p. 124002.

[30] X. FANG AND Y. KOIKE, *High-dimensional central limit theorems by stein's method*, The Annals of Applied Probability, 31 (2021), pp. 1660–1686.

[31] M. FATHI, *Stein kernels and moment maps*, The Annals of Probability, 47 (2019), pp. 2172–2185.

[32] A. FIGALLI, *The Monge-Ampère equation and its applications*, 2017.

[33] R. FLOCK, Y. DONG, F. URIBE, AND O. ZAHM, *Certified coordinate selection for high-dimensional bayesian inversion with laplace prior*, (2023).

[34] M. GIROLAMI AND B. CALDERHEAD, *Riemann manifold langevin and hamiltonian monte carlo methods*, Journal of the Royal Statistical Society Series B: Statistical Methodology, 73 (2011), pp. 123–214.

[35] R. IBRAGIMOV AND S. SHARAKHMETOV, *The exact constant in the rosenthal inequality for random variables with mean zero*, Theory of Probability & Its Applications, 46 (2002), pp. 127–132.
[20] T. CACOULLOS, *On upper and lower bounds for the variance of a function of a random variable*, The Annals of Probability, 10 (1982), pp. 799–809.

[21] T. CACOULLOS, V. PAPATHANASIOU, AND S. A. UTEV, *Variational inequalities with examples and an application to the central limit theorem*, The Annals of Probability, (1994), pp. 1607–1618.

[22] J. CHEEGER, *Differentiability of lipschitz functions on metric measure spaces*, Geometric & Functional Analysis GAFA, 9 (1999), pp. 428–517.

[23] L. H. CHEN, *Stein's method of normal approximation: Some recollections and reflections*, The Annals of Statistics, 49 (2021), pp. 1850–1863.

[24] X. CHENG, B. WANG, J. ZHANG, AND Y. ZHU, *Fast conditional mixing of mcmc algorithms for non-log-concave distributions*, Advances in Neural Information Processing Systems, 36 (2024).

[25] V. CHERNOZHUOKOV, D. CHETVERIKOV, K. KATO, AND Y. KOIKE, *Improved central limit theorem and bootstrap approximations in high dimensions*, The Annals of Statistics, 50 (2022), pp. 2562–2586.

[26] F. CLARKE, *Functional analysis, calculus of variations and optimal control*, vol. 264, Springer, 2013.

[27] D. CORDERO-ERAUSQUIN AND B. KLARTAG, *Moment measures*, Journal of Functional Analysis, 268 (2015), pp. 3834–3866.

[28] T. A. COURTADE, M. FATHI, AND A. PANANJADY, *Existence of stein kernels under a spectral gap, and discrepancy bounds*, in Annales de l'Institut Henri Poincaré, Probabilités et Statistiques, vol. 55, Institut Henri Poincaré, 2019, pp. 777–790.

[29] T. CUI, X. T. TONG, AND O. ZAHM, *Prior normalization for certified likelihood-informed subspace detection of bayesian inverse problems*, Inverse Problems, 38 (2022), p. 124002.

[30] X. FANG AND Y. KOIKE, *High-dimensional central limit theorems by stein's method*, The Annals of Applied Probability, 31 (2021), pp. 1660–1686.

[31] M. FATHI, *Stein kernels and moment maps*, The Annals of Probability, 47 (2019), pp. 2172–2185.

[32] A. FIGALLI, *The Monge-Ampère equation and its applications*, 2017.

[33] R. FLOCK, Y. DONG, F. URIBE, AND O. ZAHM, *Certified coordinate selection for high-dimensional bayesian inversion with laplace prior*, (2023).

[34] M. GIROLAMI AND B. CALDERHEAD, *Riemann manifold langevin and hamiltonian monte carlo methods*, Journal of the Royal Statistical Society Series B: Statistical Methodology, 73 (2011), pp. 123–214.

[35] R. IBRAGIMOV AND S. SHARAKHMETOV, *The exact constant in the rosenthal inequality for random variables with mean zero*, Theory of Probability & Its Applications, 46 (2002), pp. 127–132.
[36] M. LEDOUX, I. NOURDIN, AND G. PECCATI, *Stein's method, logarithmic sobolev and transport inequalities*, Geometric and Functional Analysis, 25 (2015), pp. 256–306.

[37] M. T. LI, Y. MARZOUK, AND O. ZAHM, *Principal feature detection via $\phi$-sobolev inequalities*, arXiv preprint arXiv:2305.06172, (2023).

[38] R. LI, M. TAO, S. S. VEMPALA, AND A. WIBISONO, *The mirror langevin algorithm converges with vanishing bias*, in International Conference on Algorithmic Learning Theory, PMLR, 2022, pp. 718–742.

[39] Q. LIU, J. LEE, AND M. JORDAN, *A kernelized stein discrepancy for goodness-of-fit tests*, in International conference on machine learning, PMLR, 2016, pp. 276–284.

[40] E. NAKAI, N. TOMITA, AND K. YABUTA, *Density of the set of all infinitely differentiable functions with compact support in weighted sobolev spaces*, Scientiae Mathematicae Japonicae, 60 (2004), pp. 121–128.

[41] Y. E. NESTEROV, *A method of solving a convex programming problem with convergence rate $O(k^2)$*, in Doklady Akademii Nauk, vol. 269, Russian Academy of Sciences, 1983, pp. 543–547.

[42] M. T. PARENTE, J. WALLIN, AND B. WOHLMUTH, *Generalized bounds for active subspaces*, Electronic Journal of Statistics, 14 (2020), pp. 917–943.

[43] S. PATTERSON AND Y. W. TEH, *Stochastic gradient riemannian langevin dynamics on the probability simplex*, Advances in neural information processing systems, 26 (2013).

[44] L. E. PAYNE AND H. F. WEINBERGER, *An optimal poincaré inequality for convex domains*, Archive for Rational Mechanics and Analysis, 5 (1960), pp. 286–292.

[45] L. PILLAUD-VIVIEN, F. BACH, T. LELIÈVRE, A. RUDI, AND G. STOLTZ, *Statistical estimation of the poincaré constant and application to sampling multimodal distributions*, in International Conference on Artificial Intelligence and Statistics, PMLR, 2020, pp. 2753–2763.

[46] H. POINCARÉ, *Sur les équations aux dérivées partielles de la physique mathématique*, American Journal of Mathematics, (1890), pp. 211–294.

[47] N. QIAN, *On the momentum term in gradient descent learning algorithms*, Neural networks, 12 (1999), pp. 145–151.

[48] R. T. ROCKAFELLAR, *Convex Analysis*, Princeton University Press, Princeton, 1970.

[49] O. ROUSTANT, F. BARTHE, AND B. IOOSS, *Poincaré inequalities on intervals–application to sensitivity analysis*, (2017).

[50] F. SANTAMBROGIO, *Dealing with moment measures via entropy and optimal transport*, Journal of Functional Analysis, 271 (2016), pp. 418–436.
[36] M. LEDOUX, I. NOURDIN, AND G. PECCATI, *Stein's method, logarithmic sobolev and transport inequalities*, Geometric and Functional Analysis, 25 (2015), pp. 256–306.

[37] M. T. LI, Y. MARZOUK, AND O. ZAHM, *Principal feature detection via $\phi$-sobolev inequalities*, arXiv preprint arXiv:2305.06172, (2023).

[38] R. LI, M. TAO, S. S. VEMPALA, AND A. WIBISONO, *The mirror langevin algorithm converges with vanishing bias*, in International Conference on Algorithmic Learning Theory, PMLR, 2022, pp. 718–742.

[39] Q. LIU, J. LEE, AND M. JORDAN, *A kernelized stein discrepancy for goodness-of-fit tests*, in International conference on machine learning, PMLR, 2016, pp. 276–284.

[40] E. NAKAI, N. TOMITA, AND K. YABUTA, *Density of the set of all infinitely differentiable functions with compact support in weighted sobolev spaces*, Scientiae Mathematicae Japonicae, 60 (2004), pp. 121–128.

[41] Y. E. NESTEROV, *A method of solving a convex programming problem with convergence rate $O(k^2)$*, in Doklady Akademii Nauk, vol. 269, Russian Academy of Sciences, 1983, pp. 543–547.

[42] M. T. PARENTE, J. WALLIN, AND B. WOHLMUTH, *Generalized bounds for active subspaces*, Electronic Journal of Statistics, 14 (2020), pp. 917–943.

[43] S. PATTERSON AND Y. W. TEH, *Stochastic gradient riemannian langevin dynamics on the probability simplex*, Advances in neural information processing systems, 26 (2013).

[44] L. E. PAYNE AND H. F. WEINBERGER, *An optimal poincaré inequality for convex domains*, Archive for Rational Mechanics and Analysis, 5 (1960), pp. 286–292.

[45] L. PILLAUD-VIVIEN, F. BACH, T. LELIÈVRE, A. RUDI, AND G. STOLTZ, *Statistical estimation of the poincaré constant and application to sampling multimodal distributions*, in International Conference on Artificial Intelligence and Statistics, PMLR, 2020, pp. 2753–2763.

[46] H. POINCARÉ, *Sur les équations aux dérivées partielles de la physique mathématique*, American Journal of Mathematics, (1890), pp. 211–294.

[47] N. QIAN, *On the momentum term in gradient descent learning algorithms*, Neural networks, 12 (1999), pp. 145–151.

[48] R. T. ROCKAFELLAR, *Convex Analysis*, Princeton University Press, Princeton, 1970.

[49] O. ROUSTANT, F. BARTHE, AND B. IOOSS, *Poincaré inequalities on intervals—application to sensitivity analysis*, (2017).

[50] F. SANTAMBROGIO, *Dealing with moment measures via entropy and optimal transport*, Journal of Functional Analysis, 271 (2016), pp. 418–436.
[51] A. SAUMARD, *Weighted poincaré inequalities, concentration inequalities and tail bounds related to stein kernels in dimension one*, Bernoulli, 25 (2019), pp. 3978–4006.

[52] S. SONG, T. ZHOU, L. WANG, S. KUCHERENKO, AND Z. LU, *Derivative-based new upper bound of sobol'sensitivity measure*, Reliability Engineering & System Safety, 187 (2019), pp. 142–148.

[53] C. STEIN, *Approximate computation of expectations*, IMS, 1986.

[54] I. SUTSKEVER, J. MARTENS, G. DAHL, AND G. HINTON, *On the importance of initialization and momentum in deep learning*, in International conference on machine learning, PMLR, 2013, pp. 1139–1147.

[55] S. VEMPALA AND A. WIBISONO, *Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices*, Advances in neural information processing systems, 32 (2019).

[56] R. VERDIÈRE, C. PRIEUR, AND O. ZAHM, *Diffeomorphism-based feature learning using poincaré inequalities on augmented input space*, (2023).

[57] C. VILLANI, *Topics in optimal transportation*, Graduate Studies in Mathematics, (2003).

[58] C. ZALINESCU, *Convex analysis in general vector spaces*, World scientific, 2002.

[59] O. C. ZIENKIEWICZ, R. L. TAYLOR, AND R. L. TAYLOR, *The finite element method: solid mechanics*, vol. 2, Butterworth-heinemann, 2000.
[51] A. SAUMARD, *Weighted poincaré inequalities, concentration inequalities and tail bounds related to stein kernels in dimension one*, Bernoulli, 25 (2019), pp. 3978–4006.

[52] S. SONG, T. ZHOU, L. WANG, S. KUCHERENKO, AND Z. LU, *Derivative-based new upper bound of sobol'sensitivity measure*, Reliability Engineering & System Safety, 187 (2019), pp. 142–148.

[53] C. STEIN, *Approximate computation of expectations*, IMS, 1986.

[54] I. SUTSKEVER, J. MARTENS, G. DAHL, AND G. HINTON, *On the importance of initialization and momentum in deep learning*, in International conference on machine learning, PMLR, 2013, pp. 1139–1147.

[55] S. VEMPALA AND A. WIBISONO, *Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices*, Advances in neural information processing systems, 32 (2019).

[56] R. VERDIÈRE, C. PRIEUR, AND O. ZAHM, *Diffeomorphism-based feature learning using poincaré inequalities on augmented input space*, (2023).

[57] C. VILLANI, *Topics in optimal transportation*, Graduate Studies in Mathematics, (2003).

[58] C. ZALINESCU, *Convex analysis in general vector spaces*, World scientific, 2002.

[59] O. C. ZIENKIEWICZ, R. L. TAYLOR, AND R. L. TAYLOR, *The finite element method: solid mechanics*, vol. 2, Butterworth-heinemann, 2000.