Guillaume Mijoule
email: guillaume.mijoule@inria.fr

Gesine Reinert

Yvik Swan
email: yvswan@ulb.ac.be

Stein operators, kernels and discrepancies for multivariate continuous distributions

Keywords:

may come

Introduction

Stein's method is a collection of tools permitting to bound quantities of the form

d H (X, W ) = sup h∈H |E[h(W )] -E[h(X)]|
where H is a measure-determining class and X, W are two random quantities of interest with X, say, following the target distribution µ. The method can be summarised as follows. First find a so-called Stein operator A and a wide class of functions F(A) such that (i) X ∼ µ if and only if E[Af (X)] = 0 for all functions f ∈ F(A) and (ii) for each h ∈ H there exists a well-defined and tractable solution f = f h ∈ F(A) of the Stein equation h(x) -E[h(X)] = Af (x). Then, upon noting that E[h(W )] -E[h(X)] = E[Af h (W )] for all h, the problem of bounding d H (X, W ) has been re-expressed in terms of that of bounding sup h∈H E[Af h (W )]. The success of the method lies in the fact that this last quantity is amenable to a wide variety of approaches.

There exist many frameworks in which Stein's method is well understood. We refer to the surveys [START_REF] Ross | Fundamentals of Stein's method[END_REF][START_REF] Barbour | Stein's (magic) method[END_REF][START_REF] Chatterjee | A short survey of Stein's method[END_REF] as well as the papers [START_REF] Ley | Stein's method for comparison of univariate distributions[END_REF][START_REF] Ley | Distances between nested densities and a measure of the impact of the prior in Bayesian statistics[END_REF] for the univariate setting. We also highlight [START_REF] Arras | On Stein's Method for Infinitely Laws with Finite First Moment[END_REF] who provide an up-to-date overview in the context of infinitely divisible distributions. Comprehensive introductions to some of the most important aspects of the theory are available from the monographs [START_REF] Stein | Approximate computation of expectations[END_REF] as well as [START_REF] Nourdin | Normal approximations with Malliavin calculus : from Stein's method to universality[END_REF][START_REF] Chen | Normal approximation by Stein's method[END_REF], with a particular focus on Gaussian approximation. Although the univariate case is the most studied, many references also tackle multivariate distributions. For discrete multivariate distributions, [START_REF] Barbour | Multivariate approximation in total variation, i: equilibrium distributions of markov jump processes[END_REF] and [START_REF] Barbour | Multivariate approximation in total variation, ii: Discrete normal approximation[END_REF] provide a framework which is applicable in many situations. In [START_REF] Reinert | Approximating stationary distributions of fast mixing Glauber dynamics, with applications to exponential random graphs[END_REF], stationary distributions of Glauber Markov chains are characterised. The multivariate Gaussian case has naturally received the most attention, starting with [START_REF] Barbour | Stein's method for diffusion approximations[END_REF] and [START_REF] Götze | On the rate of convergence in the multivariate clt[END_REF], see also [START_REF] Reinert | Multivariate normal approximation with Stein's method of exchangeable pairs under a general linearity condition[END_REF][START_REF] Fang | Multivariate approximations in Wasserstein distance by Stein's method and Bismut's formula[END_REF] and [START_REF] Gallouët | Regularity of solutions of the Stein equation and rates in the multivariate central limit theorem[END_REF]. The Dirichlet distribution has been treated in [START_REF] Gan | Dirichlet approximation of equilibrium distributions in Cannings models with mutation[END_REF] using a coupling approach. Log-concave densities are treated in [START_REF] Mackey | Multivariate Stein factors for a class of strongly logconcave distributions[END_REF], and more general settings have been dealt with in [START_REF] Gorham | Measuring sample quality with diffusions[END_REF][START_REF] Fang | Multivariate approximations in Wasserstein distance by Stein's method and Bismut's formula[END_REF]. Yet, a general approach for the multivariate case has been elusive. This is the task which our paper addresses.

The starting point for a multivariate Stein's method is a Stein characterization for the Gaussian law which states that a random vector Y is a multivariate Gaussian d-random vector with mean ν and covariance Σ (short: Y ∼ N (ν, Σ)) if and only if

E (Y -ν) t ∇f (Y ) = E ∇ t Σ∇f (Y ) , (1)
for all absolutely continuous function f : R d → R for which the expectations exist; here ∇ denotes the gradient operator. Assume that h : R d → R has three bounded derivatives. Then, if Σ ∈ R d×d is symmetric and positive definite, and Z ∼ MVN (0, Σ), there is a solution f : R d → R to the Stein equation

∇ t Σ∇f (w) -w t ∇f (w) = h(w) -Eh(Σ 1/2 Z), (2)
given at every w ∈ R d by the Mehler formula

f (w) = 1 0 1 2t Eh(Z w,t ) -Eh(Σ 1/2 Z) dt (3)
with Z w,t = √ tw + √ 1 -t Σ 1/2 Z. Generalizations to infinitely-dimensional functionals using Malliavin calculus are available in [START_REF] Nourdin | Normal approximations with Malliavin calculus : from Stein's method to universality[END_REF] and a generalization to Gaussian random measures is employed in [START_REF] Holmes | Stein's method for the bootstrap[END_REF].

Pursuing the analogy, classical Markov theory leads, under appropriate regularity assumptions, to Stein equations of the form a(x), ∇ 2 f (x) + b(x), ∇f (x) = h(x)-E[h(Z)] for targets Z ∼ µ the ergodic measure of SDE's of the form dZ t = a(t)dB t + b(X t )dt; see e.g. [START_REF] Bakry | Analysis and geometry of Markov diffusion operators[END_REF][START_REF] Gorham | Measuring sample quality with diffusions[END_REF][START_REF] Fang | Multivariate approximations in Wasserstein distance by Stein's method and Bismut's formula[END_REF]. As in the Gaussian case, a solution f h to the Stein equation is also provided by semigroup approach to Markov generators.

In this paper we shall use first order directional derivatives as basic building blocks of Stein operators. In Definition 3.1 we define the canonical Stein derivative for p in the direction e as T e,p φ = ∂ e (p φ)/p acting on sufficiently smooth test functions φ. Using this building block, we define the Stein gradient operator acting on real valued functions f , vector valued functions f or matrix valued functions F as • → T ∇,p • = ∇ (p•) p = d i=1 e i T e i ,p where {e i , i = 1, . . . , d} is the canonical basis of R d . Similarly we define the Stein divergence operator acting on vector or matrix valued functions with compatible dimensions as • → T div,p • = div (p•)/p. As in the one-dimensional case, this operator satisfies a product rule for all smooth functions f, g, namely T e,p (f g) = (T e,p f )g + f (∂ e g) for all e in the d-dimensional unit sphere S d- 1 .

Our framework provides a mechanism for obtaining a broad family of Stein operators which we call standardizations of the canonical operator T D,p , where D = ∇ or div -see Section 3. We do not claim that all useful operators from the literature are of this form -for instance [START_REF] Bonis | Rates in the central limit theorem and diffusion approximation via Stein's method[END_REF] proposes an alternative construction. We do demonstrate, however, that the operators we construct provide crucial handles on a wide family of multivariable distributions.

For example, setting ρ = ∇ log p = T ∇,p 1 the multivariate score function, the analog to the one-dimensional score-Stein operator is the vector valued operator A p g = ∇g + ρ g acting on differentiable functions g : R d → R, see Definition 3. [START_REF] Arras | On Stein's Method for Infinitely Laws with Finite First Moment[END_REF] A second set of Stein operators are related to Stein kernels. In dimension d = 1, this is the unique bounded solution x → τ p (x) in F(p) to the ODE T p (τ p (x)) = ν -x, given by

τ p (x) = 1 p(x) ∞ x (y -ν)p(y)dy, x ∈ I (4)
(with ν = E p [X]). Properties of the Stein kernel were first studied in [START_REF] Stein | Approximate computation of expectations[END_REF] (although it had long been a known important handle on smooth densities p, see e.g. [START_REF] Cacoullos | Characterizations of distributions by variance bounds[END_REF][START_REF] Cacoullos | On upper and lower bounds for the variance of a function of a random variable[END_REF], who refer to it as a covariance kernel). The Stein kernel has become a powerful tool in Stein's method, see e.g. [START_REF] Cacoullos | Variational inequalities with examples and an application to the central limit theorem[END_REF][START_REF] Nourdin | Density formula and concentration inequalities with Malliavin calculus[END_REF][START_REF] Nourdin | Normal approximations with Malliavin calculus : from Stein's method to universality[END_REF] or the more recent works [START_REF] Fathi | Stein kernels and moment maps[END_REF][START_REF] Courtade | Existence of Stein kernels under a spectral gap, and discrepancy bound[END_REF]. In one dimension, Stein kernels are unique when they exist -the Stein kernel is the zero bias density from [START_REF] Goldstein | Stein's method and the zero bias transformation with application to simple random sampling[END_REF]. They have been studied in detail in [START_REF] Ernst | First order covariance inequalities via Stein's method[END_REF][START_REF] Saumard | Weighted Poincaré inequalities, concentration inequalities and tail bounds related to the behavior of the Stein kernel in dimension one[END_REF]. In higher dimensions even the definition of a Stein kernel is not obvious; see [START_REF] Nourdin | Entropy and the fourth moment phenomenon[END_REF][START_REF] Nourdin | Integration by parts and representation of information functionals[END_REF][START_REF] Courtade | Existence of Stein kernels under a spectral gap, and discrepancy bound[END_REF][START_REF] Ledoux | Stein's method, logarithmic Sobolev and transport inequalities[END_REF]. The zero-bias coupling definition in [START_REF] Goldstein | Zero biasing in one and higher dimensions, and applications[END_REF] uses a collection of random elements in higher dimensions. In analogy, in Definition 3.8 we define a directional Stein kernel for each canonical direction e i ∈ R d , as any function x → τ p,i (x) ∈ R × R d belonging to F(p) such that T div,p τ p,i (x) = E p [X i ]-x i for all x ∈ supp(p). A Stein kernel is then any square matrix τ τ τ such that each line τ τ τ i = (τ i1 , . . . , τ id ) is a kernel in the direction e i , and the kernel-Stein operator is the operator A p g(x) = τ τ τ p (x)∇g(x) -(x -ν)g(x). Stein kernels need not exist and may not be unique when they exist. Aside from discussing this issue in Section 6, we also give an explicit -formula for Stein kernels of elliptical distributions in Section 5. In particular, we obtain new identities even for the standard Gaussian: aside from Stein's classical Gaussian covariance identity

X ∼ N (ν, Σ) ⇒ E[Σ∇g(X)] = E[(X -ν)g(X)] (5)
we also prove that for all β = 2 and d = 1, if X ∼ N (ν, Σ) then

β (β -2)(d -1) E X T Σ -1 XΣ -XX T ∇g(X) = E 2 β -2
Σ∇g(X) + Xg(X) ;

both identities hold for a wide family of test functions g, see Section 5.2. Also, in the case of the multivariate Student t-distribution, we obtain the identity

X ∼ t k (ν, Σ) ⇒ E (XX T + kI d )∇g(X) = (k -1)E [Xg(X)]
which generalizes the corresponding univariate result (see Section 5.4). We also provide, in Section 6, a connection between Stein kernels and transport maps. These results are closely linked to the seminal works [START_REF] Artstein | Solution of Shannon's problem on the monotonicity of entropy[END_REF][START_REF] Ball | Entropy jumps in the presence of a spectral gap[END_REF][START_REF] Artstein | On the rate of convergence in the entropic central limit theorem[END_REF]; we show that the key quantity they introduce to solve e.g. Shannon's conjecture actually is a Stein kernel. In particular, these results give new and explicit expressions for Stein kernels in d = 2 and d = 3 dimensions. A natural question is which Stein operator to choose. Here practical issues may serve as guidance: solutions to multivariate Stein equations may not generally exist. In the case of densities which possess a Poincaré constant, weak solutions of second-order Stein equations exist, and some regularity results for these solutions are available, see [START_REF] Raič | A multivariate CLT for decomposable random vectors with finite second moments[END_REF][START_REF] Reinert | Multivariate normal approximation with Stein's method of exchangeable pairs under a general linearity condition[END_REF][START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF][START_REF] Valiant | A CLT and tight lower bounds for estimating entropy[END_REF][START_REF] Gorham | Measuring sample quality with diffusions[END_REF][START_REF] Fang | Multivariate approximations in Wasserstein distance by Stein's method and Bismut's formula[END_REF][START_REF] Gallouët | Regularity of solutions of the Stein equation and rates in the multivariate central limit theorem[END_REF], see Section 3.3.

We illustrate our framework by assessing the (1-)Wasserstein distance between two absolutely continuous distributions. This result leads to three applications: assessing the difference between copulas, assessing the distance between prior and posterior in a Bayesian setting, and bounding the Wasserstein distance between skew-normal and standard normal distributions.

In practical applications and for obtaining information inequalities, the Stein discrepancy plays a key role. The Stein discrepancy from distribution p to distribution q using the Stein operator A p for p is given in (104);

S • (q, A p , G) = sup g∈G E [A p g(Y )] .
Here Y ∼ q and G is a class of smooth test functions. A particular case is that of kernelized Stein discrepancies, introduced in [START_REF] Liu | A kernelized Stein discrepancy for goodness-of-fit tests[END_REF][START_REF] Chwialkowski | A kernel test of goodness of fit[END_REF]. We assess the performance of such a kernelized Stein discrepancy measure for simulations from a centred Student t-distribution, employing a Stein discrepancy for a goodness-of-fit test.

The paper is structured as follows. Section 2 introduces the notations and the product rule. Section 3.1 gives the general Stein operators, Stein identities, and Stein characterizations. Score-Stein operators and kernel-Stein operators are introduced, and properties of the solution of Stein equations are given for densities which admit a Poincaré constant. Applications of this framework are given in Section 4. Then the paper turns its attention to Stein kernels. In Section 5, Stein kernels for elliptical distributions are derived. Section 6 provides the general discussion on Stein kernels. Finally, Section 7 discusses information metrics, kernelized Stein discrepancies, and ends with a goodness-of-fit test based on Stein discrepancies.

Notations, gradients and product rules

In this paper we use small caps for elements of R m , m ≥ 1 (by convention column vectors), capitalized letters for matrices in R m × R n , small caps for real-valued functions, boldfaced small caps for vector valued functions and boldfaced capitalized letters for matrix valued functions.

Fix d ∈ N 0 and let e 1 , . . . , e d be the canonical basis for Cartesian coordinates in R d . Given x, y ∈ R d and any symmetric positive-definite d × d matrix A we set x, y A = x T Ay (here • T denotes the transpose) with associated norm x A =

x, x A . With Tr(•) the trace operator, the Hilbert-Schmidt scalar product between matrices A, B of compatible dimension is A, B HS = Tr(AB T ), with associated norm A 2 HS := Tr(A T A) = d i,j=1 a 2 i,j . Let S d-1 denote the unit sphere in R d and let e ∈ S d-1 be a unit vector in R d . The directional derivative of a smooth function v : R d → R in the direction e is the function

∂ e v : R d → R given by ∂ e v(x) = lim h→0 v(x + he) -v(x) h , (6)
at every point where this limit exists. For i = 1, . . . , d we write ∂ i v for the derivative in the direction of the unit vector e i . The gradient of a smooth function v :

R d → R is ∇v = grad(v) = (∂ 1 v, . . . , ∂ d v) T = d i=1 (∂ i v) e i (7)
(by convention, a column vector). The gradient of a 1 × m vector field v = (v 1 , . . . , v m ) (a line vector) is

∇v = ∇v 1 • • • ∇v m =    ∂ 1 v 1 • • • ∂ 1 v m . . . . . . . . . ∂ d v 1 • • • ∂ d v m    . (8)
(by convention, a d × m matrix).

The Jacobian matrix of a R m -valued function w = (w 1 , . . . , w m ) T (a column vector) is the m × d matrix Jac(w) = ∂ i w j 1≤i≤m,1≤j≤d . As Jac(v T ) = ∇v we will simply use the generic notation ∇ for both operations from here onwards. The divergence operator is defined for d valued (line or column) vector fields v with components v j , j = 1, . . . , d as

div(v) = Tr (∇v) = d i=1 ∂ i v i .
More generally, given any m, n ∈ N, the directional derivative of a matrix valued function

F : R d → R m × R d : x → F(x) =    f 1 (x)
. . .

f m (x)    =    f 11 (x) . . . f 1d (x) . . . . . . . . . f m1 (x) . . . f md (x)    (9)
is defined componentwise;

(∂ e F(x)) 1≤i≤m,1≤j≤n = (∂ e f ij (x)) 1≤i≤m,1≤j≤n
for all e ∈ S d-1 and all x ∈ R d . The gradient of a m × r matrix valued function F of the form ( 9) is defined as the d × r × m tensor with entry

(∇F ) i,j,k = ∂ i f jk , 1 ≤ i ≤ d, 1 ≤ j ≤ m, 1 ≤ k ≤ r.
The divergence of F as in ( 9) is defined as the m × 1 column vector with components div(f j ), j = 1, . . . , m; so that

div(F) =    div(f 1 ) . . . div(f m )    =        d i=1 ∂f 1i ∂x i . . . d i=1 ∂f mi ∂x i        . Similarly, the divergence of F = (f 1 , . . . , f m ) with values in R d ×R m is div(F) = (div(F T )) T .

The product rules

Given v and w two sufficiently smooth functions from R d → R and e ∈ S d-1 , the directional derivative satisfies the product rule

∂ e (v(x)w(x)) = (∂ e v(x))w(x) + v(x)(∂ e w(x), or short: ∂ e (vw) = (∂ e v)w + v∂ e w, (10)
for all x ∈ R d at which all derivatives are defined. Gradients and divergences therefore also satisfy the corresponding product rules. We shall mainly consider three instances.

1. If v is a m-dimensional vector field and φ : R d → R then

∇(φv) = φ(∇v) + (∇φ)v (11) (a d × m matrix if v is a line, a m × d matrix if v is a column). In particular if m = d and v is a d × 1 column vector, then div(φv) = φ div(v) + ∇φ, v (12)
(a scalar). When v = ∇ψ is a gradient of a sufficiently smooth function ψ, then div(φ∇ψ) = φ ∆ψ + ∇φ, ∇ψ .

2. For v ∈ R m and F a m × d square matrix we have

div(v T F) = v, divF + ∇v, F HS (14)
(a scalar).

3. If v, w ∈ R m and A is a m × m positive definite symmetric matrix then (with a slight abuse of notation)

∇ ( v, w A ) = ∇v, w A + v, ∇w A . ( 15
)

Stein operators and their standardizations

In this paper, all random vectors are assumed to admit a probability density function (pdf) p with respect to the Lebesgue measure. The support of p, denoted by K p , is the complement of the largest open set O such that p ≡ 0 on O. We further make the following assumption:

Assumption A: There exists an open set Ω p such that the Lebesgue measure of K p \Ω p is zero, p > 0 on Ω p and p is C 1 on Ω p .

If p satisfies this assumption, then p is C 1 on Ω p ∪ K C p , which by Assumption A is a set of full Lebesgue measure. For ease of notation the Lebesgue measure is left out in integrals, so that R d f (x)dx is written as R d f . Similarly we abbreviate E p f = Kp f p = Kp f (x)p(x) dx. defined by T e,p φ = ∂e(p φ) p . We set F 1 (p) the collection of all scalar valued functions f which belong to F 1,e i (p) for all unit vectors e i , i = 1, . . . , d, and define the Stein class of p, denoted F(p), as the collection of all scalar, vector and matrix valued functions whose components belong to F 1 (p).

Definition and some general comments

Note in particular that

if φ ∈ F 1 (p), then R d ∂ e (p φ) = 0 for all unit vectors e. Furthermore, any function φ ∈ C 1 (R d ) with compact support K ⊂ Ω p ∪ K C p lies in F 1,e (p). Indeed, if e = (a 1 , . . . , a d ), since pφ ∈ C 1 (R d ), we have ∂ e (pφ) = d i=1 a i ∂ i (pφ), then R d ∂ i (pφ) = R d-1 R ∂ i (pφ)dx i j =i dx j ,
and the inside integral is zero since p(x)φ(x) cancels out for x away enough from 0.

On F(p) we define the Stein gradient operator acting on real valued functions f , vector valued functions f or matrix valued functions F, as

• → T ∇,p • = ∇ (p•) p
(expressed in terms of T e,p 's as T ∇,p = d i=1 e i T e i ,p ); we also define the Stein divergence operator acting on vector or matrix valued functions with compatible dimensions as

• → T div,p • = div (p•) p
(this operator also can be expressed in terms of T e,p 's). Although in the literature, it is not T e,p but T ∇,p and T div,p that are the operators used, their properties follow from those of (T e,p ), e ∈ B with B a unit basis of R d .

Remark 3.2. Let F (0) (p) denote the collection of functions (in any dimension) with mean 0 under p. The definition of F(p) is tailored to ensure T D,p (F) ∈ F (0) (p) for all F ∈ F(p) (here and below D denotes ∂ e , ∇ or div). Hence, the Stein operators are also a machinery for producing mean 0 functions under p. This can be useful in applications, particularly when p is intractable, since by construction the operators are invariant under scalar multiplication, and therefore do not depend on normalizing constants.

Stein identities

Product rule [START_REF] Bakry | Analysis and geometry of Markov diffusion operators[END_REF] directly leads to Stein-type product rules of the type

T e,p (f g) = (T e,p f )g + f (∂ e g) (16)
on a suitable set of functions. The next definition introduces such sets. Using the product rules from Section 2, we deduce similar identities for Stein gradients and divergences. Here we give two instances.

1. For all scalar functions f ∈ F 1 (p) and g ∈ dom(p, f ), from (11)

T ∇,p (f g) = (T ∇,p f )g + f ∇g. (17)
2. For all

F : R d → R m × R d ∈ F 1 (p) and g : R d → R × R m ∈ dom(p, F), using (14),
T div,p (gF) = T div,p F, g + F, ∇g HS .

We deduce a family of Stein operators for p by fixing F = (F ij ) 1≤i≤m,1≤j≤d and considering A p := A F,p defined by

g → A p g = T div,p F, g + F, ∇g HS (19)
with domain F(A p ) := dom(p, F). Note that A p (g) is a scalar function. We also deduce the Stein identities (one for each F)

E p [ T div,p F, g ] = -E p [ F, ∇g HS ] for all g ∈ F(A p ). (20)
More product rules are provided later, see e.g. ( 25), [START_REF] Döbler | An iterative technique for bounding derivatives of solutions of Stein equations[END_REF].

It follows from the definition of the classes of functions dom(p, F) and the Stein operators and the Stein class that the left hand side in each of ( 16), [START_REF] Brascamp | 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[END_REF], and their variations such as [START_REF] Cacoullos | Characterizations of distributions by variance bounds[END_REF], ( 19), ( 23), [START_REF] Chen | Normal approximation by Stein's method[END_REF], integrates to 0. These probabilistic integration by parts formulas lead to identities known as Stein (covariance) identities; they are inspired by Stein's original Gaussian identity (1). Proposition 3.4 (Stein identities). For all f : R d → R ∈ F 1 (p) and all e ∈ S d-1 we have

E p [f (T e,p g)] = -E p [g(∂ e f )] (21)
for all g:

R d → R ∈ dom(p, f ) Proof. Let f ∈ F 1 (p) and g ∈ dom(p, f
) and e ∈ S d-1 . Then, from [START_REF] Bonis | Rates in the central limit theorem and diffusion approximation via Stein's method[END_REF], f (T e,p g) = T e,p (f g) -g∂ e f . Since E p T e,p (f g) = 0 for all f ∈ F 1 (p) and g ∈ dom(p, f ), identity (

Identity [START_REF] Cacoullos | Variational inequalities with examples and an application to the central limit theorem[END_REF] shows that the Stein operators are the skew adjoint operators to the classical directional derivatives, with respect to integration in p. In this sense the operator from Definition 3.1 is "canonical".

The Stein identities [START_REF] Cacoullos | Variational inequalities with examples and an application to the central limit theorem[END_REF] for (e 1 , . . . , e d ) the standard unit basis in R d yield corresponding Stein identities for the gradient Stein operator T ∇,p and for the divergence Stein operator T div,p . Here we give three instances.

1. For all f ∈ F 1 (p) and g ∈ dom(p, f ),

E p [(T ∇,p f )g] = -E p [f (∇g)] ' (22)
2. For all f ∈ F(p) and g such that ∇g ∈ dom(p, f ), from ( 13) we obtain

T div,p (f ∇g) = f ∆g + T ∇,p f, ∇g . (23)
Consequently,

E p [f ∆g] = -E p [ T ∇,p f, ∇g ] (24)
for all such g.

3.

For all m × d square matrices F ∈ F(p) and g : R d → R m ∈ dom(p, F), using (14) gives T div,p (g T F) = g, T div,p (F) + ∇g, F HS .

Consequently,

E p [ g, T div,p (F) ] = -E p [ ∇g, F HS ] . (26)

Stein characterizations

Recall that E p [T e,p F] = 0 for every F ∈ F(p). Under well chosen assumptions, the collection of operators T e i ,p for (e 1 , . . . , e d ) the standard unit basis of R d characterises p in the sense that if E q [T e i ,p F] = 0, i = 1, . . . , d for a sufficiently wide class of F, then necessarily q = p.

Theorem 3.5 (Stein characterizations). Let X ∼ p and Y ∼ q; assume that p and q both satisfy Assumption A. Assume moreover that Ω p = Ω q , and that this set is connected. Let f : R d → R ∈ F 1 (p), and assume f > 0 on Ω p .

It holds that

Y L = X if and only if for all g : R d → R d ∈ dom(p, f ) E q [(T div,p g)f ] = -E q [ g, ∇f ] . (27)
2. It holds that Y L = X if and only if for all g : R d → R ∈ dom(p, f )

E q [(T ∇,p g)f ] = -E q [g(∇f )] . (28)
Proof. First let f : R d → R ∈ F 1 (p) and g : R d → R d ∈ dom(p, f ). By the product rule [START_REF] Barbour | Stein's method for diffusion approximations[END_REF], f div(pg) = div(f pg) -p ∇f, g and hence f T div,p g = T div,p (f g) -∇f, g . Taking expectations gives

E p [f T div,p g] = E p [T div,p (f g)] -E p [ ∇f, g ].
By construction of dom(p, f ), E p [T div,p (f g)] = 0 and the only-if direction follows.

For the converse direction, let g be an infinitely differentiable vector field with compact support K ⊂ Ω p ; recall that in this case g ∈ dom(p, f ). Suppose that ( 27) holds. Rewriting the expectations,

K div(pg)f q p = - K q ∇f, g .
Now, by the divergence (or Ostrogradski) theorem, and since g ≡ 0 on ∂K,

K div(pg)f q p = - K p g, ∇ f q p .
Thus K p g, ∇ f q p = K g, ∇g q. Hence p∇ f q p = q∇g in the weak sense on Ω p ; but it is also true in the strong sense, because all functions are continuous (recall that Ω p = Ω q ). We deduce that f ∇(q/p) = 0, on Ω p , but since f > 0 on this connected domain, q/p is constant on Ω p . Finally, using Ωp p = Ωp q gives the first result.

For Item 2, proceeding as for the first part, we arrive at the differential equation p∇ f q p = q∇f from which we deduce that p/q = 1.

The proof of the claims in Theorem 3.5 are greatly facilitated by the tailored assumptions which hide some difficulties. It still remains to write such characterizations out in detail for specific targets; this can turn out to be quite complicated and will not be the focus of the paper. Also, there are many more general ways to formulate similar characterizations as in Theorem 3.5 be it by changing the starting operator, by relaxing the assumptions on the test functions or indeed by relaxing the assumptions on Y .

Standardizations

Similarly as in the univariate case (see [START_REF] Ley | Stein's method for comparison of univariate distributions[END_REF]Section 4]), we introduce a broad family of Stein operators A p : F(A p ) → F (0) (p) which we call standardizations of the canonical operator T D,p , where D = ∇ or div. These are operators such that there exists a transformation T : F(A p ) → F(p) : u → T(u) such that A p (u) = T p (T(u)) for all u ∈ F(A p ). Such standardizations were studied in [START_REF] Ley | Stein's method for comparison of univariate distributions[END_REF]; [START_REF] Oates | Posterior Integration on an Embedded Riemannian Manifold[END_REF]Appendix A.2] suggests some standardizations for densities on a manifold. Here we concentrate on standardizations which are connected to score functions and Stein kernels.

Gradient based operators and the score function

In this section we focus on the Stein identity [START_REF] Chatterjee | A short survey of Stein's method[END_REF], which arises from the product rule ( 17) extended as

T ∇,p (f g) = (T ∇,p f )g + f ∇g (29)
for f ∈ F 1 (p) and g : R d → R m ∈ dom(p, f ). From this, we deduce a family of Stein operators for p obtained by fixing some f ∈ F(p) and considering the operator

A p := A f,p given in (29) acting through g → A p g := T ∇,p (f g) on g : R d → ×R m ∈ F(A p ) = dom(p, f ). Note that A p (g) is a d × m matrix, and E p [A p g] = 0 for all g ∈ dom(p, f ).
Each particular f ∈ F 1 (p) thus gives rise to a Stein identity

E p [(T ∇,p f )g] = -E p [f ∇g] for all g ∈ F(A p ). (30)
One particular choice for f stands out: f = 1. The following definition is classical.

Definition 3.6 (Score function). Let p be differentiable. The score function of p is the function defined on Ω p by

ρ p (x) = ∇ log p(x) = 1 p(x) ∇p(x). (31)
The score-Stein operator is the vector valued operator

A p g = ∇g + ρ p g (32)
acting on differentiable functions g : R d → R.

∂ i p is integrable for all i = 1, . . . , d and R d ∂ i p = 0. It is easy to see that C ∞ 0 (R d ) ⊂ dom(p, 1). The resulting (characterizing) Stein identity is E p [ρ p g] = -E p [∇g] for all g ∈ F(A p ). (33)
Equation ( 33) is reminiscent of the classical Stein identity (1); one can easily check that ρ φ (x) = ∇ log φ(x) = -Σ -1 (x -ν) when X ∼ φ the multivariate normal density with mean ν and covariance Σ. We shall see that for the Student-t distribution the choice of a non-constant f in (30) leads to more tractable operators and identities. We will discuss operator [START_REF] Ernst | First order covariance inequalities via Stein's method[END_REF] and Stein identity (33) -and variations thereon -in the more general context of score functions and Fisher information in Section 7.2.

Divergence based first order operators and Stein kernels

In this section, we focus on the product rule ( 12) rewritten as

T div,p (Fg) = (T div,p F)g + F∇g (34)
for all properly chosen F : R d → R m × R d and all g : R d → R. By construction, T div,p (Fg) ∈ R m for all F, g. From this we deduce a family of Stein operators for p obtained by fixing F in F(p) and considering A p := A F,p defined by g → A p g = (T div,p F)g + F∇g with domain F(A p ) := dom(p, F). We also deduce the Stein identities (one for each F)

E p [(T div,p F)g] = -E p [F∇g] for all g ∈ F(A p ). (35)
One particularly important choice of F in ( 35) is the Stein kernel first defined in dimension d = 1 in [START_REF] Stein | Approximate computation of expectations[END_REF]. Here we consider general dimension d ≥ 1 and propose the following definition.

x → τ p,i (x) ∈ R d belonging to F(p) such that T div,p τ p,i (x) = ν i -x i (36)
for all x ∈ Ω p . A Stein kernel is any square matrix τ τ τ p = (τ i,j ) 1≤i,j≤d such that each line τ τ τ i = (τ i1 , . . . , τ id ) is a kernel in the direction e i . The τ τ τ p -kernel-Stein operator is the

R d -valued operator A p g(x) = τ τ τ p (x)∇g(x) -(x -ν)g(x) (37)
with domain F(A p ) = dom(p, τ τ τ p ).

When the kernel is clear, we leave out the τ τ τ p -prefix in the kernel-Stein operator. To illustrate the concept we give two examples; Section 6 contains more constructions for some cases.

Example 3.9. A Stein kernel in the sense of Definition 3.8 is immediately obtained in the following way. For each i = 1, . . . , d let τ p,i = (0, . . . , τ i , . . . , 0) with

τ i (x) = 1 p(x) ∞ x i (ν i -u i )p(x \i , u i )du i where (x \i , u i ) = (x 1 , . . . , x i-1 , u i , x i+1 , . . . , x d ).
The diagonal matrix τ τ τ obtained by stacking up the above vectors is a Stein kernel. If X has finite first moment, this Stein kernel is defined for almost every x by Fubini's theorem. Such diagonal kernels are often not interesting (mainly because the Stein identity (40) below will not hold for a large enough class of functions).

Example 3.10. The score function and the divergence operator are linked through

T div,p F = Fρ p + div(F),
where F ∈ F(p). Hence for a density p with mean ν, a Stein kernel is any matrix τ τ τ satisfying τ τ τ (x)ρ p (x) + div(τ τ τ (x)) = -x for all x ∈ Ω p . A general construction of Stein kernels is thus given by finding an F ∈ F(p) such that for some constants α, β ∈ R, α + β = 0,

F(x)ρ p (x) = α(x -ν) + r(x) (38) div(F(x)) = β(x -ν) -r(x) (39)
(r(x) is any function) and setting

τ τ τ = 1 α + β F.
Remark 3.11. 1. With Definition 3.8, existence of the Stein kernel is readily checked, and the resulting Stein identity is

E [τ τ τ (X)∇g(X)] = E [(X -ν)g(X)] for all g ∈ F(A p ) = dom(p, τ τ τ ). (40
E [ τ p,i (X), ∇g(X) ] = E [(X i -ν i )g(X)] , (41)
for all g ∈ C ∞ c (Ω p ), the set of functions infinitely differentiable with compact support in Ω p . Thus, showing some matrix is a Stein kernel can be done equivalently by checking pointwisely [START_REF] Gallouët | Regularity of solutions of the Stein equation and rates in the multivariate central limit theorem[END_REF], or by checking the previous identity holds at least for infinitely differentiable functions with compact support included in Ω p (and that it is in F(p)). However this class of functions is narrow, and in most applications, one would want [START_REF] Goldstein | Stein's method and the zero bias transformation with application to simple random sampling[END_REF] to hold for more test functions g -e.g. for bounded g. In this work, we use the above definition to find Stein kernels, and we extend the class of functions for which the identity holds on a case-by-case basis. This will be done, for instance, for the Stein kernels we build for elliptical distributions. From the divergence theorem, when the boundary of Ω p is smooth enough, any g ∈ C 1 (Ω p ) such that τ gp → 0 on the boundary of Ω p will verify the identity (41).

Divergence based second order operators

The starting point here is the divergence product rule [START_REF] Barbour | Multivariate approximation in total variation, i: equilibrium distributions of markov jump processes[END_REF] and the corresponding Stein operator equation [START_REF] Chen | Normal approximation by Stein's method[END_REF], extended as follows. Choose A, B two d×d matrix valued functions in F(p). Then

T div,p (∇g) T AB T = ∇g, T div,p (B) A + ∇ (∇g) T A , B HS (42)
for all g ∈ C 1 (Ω p ). This leads to a family of second order scalar valued operators obtained by fixing A and/or B and considering A p g := A A,B,p g defined in [START_REF] Gómez-Sánchez-Manzano | A multivariate generalization of the power exponential family of distributions[END_REF] with domain the collection of g such that (∇g T A ∈ dom(p, B). We single out three particular choices for A, and B.

1. A = B = I d yields A p g = ∇ log p, ∇g + ∆g, (43)
(∆ being the Laplacian on R d ) with domain dom(A p );

2. A = I d and B = τ τ τ T p the (transpose of a) Stein kernel of X yields

B p g = ∇g, ν -• + ∇ 2 g, τ τ τ p HS (44)
(where ∇ 2 (f ) = ∇(∇f T ) is the Hessian of f and ν is the mean of X) with domain dom(B p ).

3.

A symmetric definite positive and B such that T div,p (B) = b yields

C p g = ∇g, b A + ∇ 2 A g, B HS (45)
(∇ 2 A g = ∇ (∇g) T A ) with domain dom(C p ). One recognizes in operators such as [START_REF] Gorham | Measuring sample quality with Stein's method[END_REF] and ( 45) the infinitesimal generators of multivariate diffusions, see e.g. [START_REF] Gorham | Measuring sample quality with diffusions[END_REF].

Stein equations and Stein factors

Let X ∼ p with Stein canonical class and operator (F(p), T D, p), where D denotes ∂ e , ∇, or div, and consider some standardization of the form

A p : F(A p ) → im(A p ) : g → A p g (46)
as detailed in Subsection 3.2.

Definition 3.12 (Stein's equations and (magic) factors). Instate all previous notations. Let H ⊂ L 1 (p) be a family of test functions, and suppose that h -E p h ∈ im(A p ). The (A p -H) Stein equation for X is the family of differential equations

A p g = h -E p h; h ∈ H. ( 47
)
For a given h, a solution to (47) is an absolutely continuous function g h such that there exists a version of the derivatives for which (47) is satisfied at all x ∈ R d . An (A p -H) Stein factor for X is any uniform bound on some moment of (derivatives of )

g h = A -1 p (h- E p h) over all h ∈ H of solutions to Stein equations (47).
Solving, and bounding the solution of, the Stein equation ( 47) is well studied in the univariate case under quite general conditions on p, see for example [START_REF] Stein | Approximate computation of expectations[END_REF][START_REF] Chen | Normal approximation by Stein's method[END_REF][START_REF] Nourdin | Normal approximations with Malliavin calculus : from Stein's method to universality[END_REF][START_REF] Döbler | Stein's method of exchangeable pairs for the beta distribution and generalizations[END_REF][START_REF] Döbler | An iterative technique for bounding derivatives of solutions of Stein equations[END_REF] and [START_REF] Ernst | Distances between distributions via Stein's method[END_REF]. Matters are much more complicated in the multivariate setting. When p is a multivariate Gaussian, then [START_REF] Barbour | Stein's method for diffusion approximations[END_REF] identified a solution of the Stein equation ( 2) to be given by the Mehler formula [START_REF] Adcock | Mean-variance-skewness efficient surfaces, Stein's lemma and the multivariate extended skew-Student distribution[END_REF]. Such explicit dependence of solution f h on the function h permits to study regularity properties of f h in terms of those of h, see for example [START_REF] Reinert | Multivariate normal approximation with Stein's method of exchangeable pairs under a general linearity condition[END_REF][START_REF] Raič | A multivariate CLT for decomposable random vectors with finite second moments[END_REF][START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF].

For solutions of the Stein equation here we focus on the general score equation based on [START_REF] Gorham | Measuring sample quality with diffusions[END_REF], namely

∆u + ∇ log p, ∇u = h -E p h. (48)
Bounds on the solution of this Stein equation are available (a) when p is strongly logconcave, or (b) when p admits a Poincaré constant. For the first condition, recall that a smooth density p is k-strongly log-concave for some k > 0 if

∀(x, y) ∈ R d × R d , y -x, ∇ log p(y) -∇ log p(x) < -k y -x 2 .
[59] proved that if p is k-strongly log-concave and h is 1-Lipschitz, then (48) admits a solution f = f h such that for p-almost every x,

|∇f h (x)| 1 k . ( 49
)
Similar bounds based on Malliavin calculus are obtained in [START_REF] Fang | Multivariate approximations in Wasserstein distance by Stein's method and Bismut's formula[END_REF]. [START_REF] Gorham | Measuring sample quality with diffusions[END_REF] computed uniform Stein factors for diffusions satisfying "distant dissipativity" and, more generally, for diffusions that couple rapidly. This class admits not only strong log-concave densities but also a large class of non-log concave, multimodal, and heavy-tailed distributions. We also point out that the authors manage to obtain bounds on higher-order derivatives of f h , under suitable assumptions. Related results can be found in [START_REF] Fang | Multivariate approximations in Wasserstein distance by Stein's method and Bismut's formula[END_REF].

For the second set of results, if

X ∼ p, recall that C p is a Poincaré constant associated to p if for every differentiable function ϕ ∈ L 2 (p) such that Eϕ(X) = 0, E ϕ 2 (X) C p E |∇ϕ(X)| 2
. Not all densities p admit a finite C p and probability distributions which possess a Poincaré constant are also referred to as having a spectral gap. Let

C ∞ c,0 (R d ) = {f ∈ C ∞ c (R d ) : f p = 0} ⊂ L 2 (p). Define on C ∞ c,0 (R d ) 2 the Dirichlet form E(f, g) = ∇f, ∇g p and the scalar product E 1 (f, g) = (f g + ∇f, ∇g )p. Assume that E is closable (i.e. for all u n ∈ C ∞ c,0 (R d ); n 0, if u n → 0 in L 2 (p) and (u n ) is E-Cauchy, then E(u n , u n ) → 0). This condition ensures that there exists a set of functions W 0 1,2 (p) ⊃ C ∞ c,0 (R d )
, which admit a gradient, and such that W 0 1,2 (p) is a Hilbert space for the scalar product E 1 . In particular, by passing to the limit, φ p = 0 for all φ ∈ W 0 1,2 (p), and the Poincaré inequality holds for such φ if it holds for any differentiable f ∈ L 2 (p) such that f p = 0. For background and more material on closable Dirichlet forms as well as sufficient conditions on p for a Dirichlet form to be closable, we refer to [START_REF] Ma | Introduction to the theory of (non-symmetric) Dirichlet forms[END_REF]. Proposition 3.13. Let h be a 1-Lipschitz function. Let X be a random vector with density p with Poincaré constant C p and satisfying the closability property. Let W 0 1,2 (p) be defined as above. Then there exists a weak solution u ∈ W 0 1,2 (p) to (48) such that

|∇u| 2 p C p .
The proof of the result follows exactly the lines of the proof in [START_REF] Courtade | Existence of Stein kernels under a spectral gap, and discrepancy bound[END_REF] and is hence omitted. We point out that here p does not need to satisfy Assumption A.

Remark 3.14. For strongly k-log concave p, [START_REF] Landsman | Stein's lemma for elliptical random vectors[END_REF] gives that there exists a (strong) solution u such that |∇u(x)| 1/k, for all x in the domain. Now it is known that when X has k-log-concave density, then the law of X possesses a Poincaré constant C p = 1/k see [START_REF] Brascamp | 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[END_REF] (it is not necessary to assume that the density is k-strongly log-concave). Hence we can also apply Proposition 3.13, which gives gives |∇u| 2 p 1/k. Thus, the constants in (49) and in Proposition (3.13) are the same; the bound in Proposition (3.13) is weaker only because the norm is weaker. Remark 3.15. In [START_REF] Courtade | Existence of Stein kernels under a spectral gap, and discrepancy bound[END_REF] it is shown that for any mean zero multivariate distribution p which is absolutely continuous with respect to the Lebesgue measure with finite second moment, and which satisfies a Poincaré inequality with constant C p , there exists a unique function g ∈ W 1,2 p such that τ τ τ p = ∇g is a Stein kernel for p. Moreover,

τ τ τ p 2 HS p ≤ C p |x| 2 p. ( 50
)
In the one-dimensional case, [START_REF] Ley | Distances between nested densities and a measure of the impact of the prior in Bayesian statistics[END_REF] (inspired by a similar result from [START_REF] Döbler | Stein's method of exchangeable pairs for the beta distribution and generalizations[END_REF]) show that there exists a solution u such that

|u (x)| τ p (x), (51)
τ p being the univariate Stein kernel associated to X. Under log concavity of the density, in the univariate case a stronger bound is available in [START_REF] Saumard | Weighted Poincaré inequalities, concentration inequalities and tail bounds related to the behavior of the Stein kernel in dimension one[END_REF]: Assume X is centered and has a smooth, non-vanishing k-log concave density p on R. Then τ p (x) 1/k, for all x ∈ R. Thus, in the univariate case where p is k-log concave, the bound (51) implies [START_REF] Landsman | Stein's lemma for elliptical random vectors[END_REF], which in turn implies the L 2 bound from Proposition 3.13.

Wasserstein distance between nested distributions

In this section we illustrate Stein's method in the multivariate setting in three examples which are based on a novel bound on the Wasserstein distance between nested distributions. The (1-)Wasserstein distance between two distributions F and G is Let P 1 and P 2 be two probability measures on R d , with respective pdfs p 1 and p 2 = π 0 p 1 , so that the two densities are nested. Assume p 1 is k-log concave. Consider as in (43) a Stein operator for P i , i = 1, 2 defined by

d W (F, G) = sup h∈W |F (h) -G(h)|
A i u = ∇ log p i , ∇u + ∆u.
Then A 2 u = A 1 u + ∇ log π 0 , ∇u . This relationship between generators makes it straightforward to bound the Wasserstein distance between the distributions, as follows.

Proposition 4.1. Assume p 2 = π 0 p 1 and that E[|∇π 0 (X 1 )|] < ∞. 1. Assume that p 1 is k-strongly log concave. Then with X 1 ∼ p 1 and X 2 ∼ p 2 , d W (X 1 , X 2 ) 1 k E[|∇π 0 (X 1 )|]. (52)
2. Assume that p 1 admits a Poincaré constant C p and π 0 ∈ W 1,2 (p 1 ). Then with

X 1 ∼ p 1 and X 2 ∼ p 2 , d W (X 1 , X 2 ) C p E[|∇π 0 (X 1 )| 2 ]. (53)
Proof. Let h : R d → R be a 1-Lipschitz function. To show [START_REF] Ledoux | Stein's method, logarithmic Sobolev and transport inequalities[END_REF], by [START_REF] Landsman | Stein's lemma for elliptical random vectors[END_REF], there exists a solution

u h to A 1 u h = h -hp 1 such that ∇u h ∞ 1/k. Let X 1 (X 2 ) have distribution P 1 (P 2 ). Then E[h(X 2 )] -E[h(X 1 )] = E[A 1 u h (X 2 )] = E [A 2 u h (X 2 ) -∇ log π 0 (X 2 ), ∇u h (X 2 ) ] = -E [ ∇ log π 0 (X 2 ), ∇u h (X 2 ) ] .
By [START_REF] Landsman | Stein's lemma for elliptical random vectors[END_REF] we have | ∇ log π 0 (X 2 ), ∇u h (X 2 ) | |∇ log π 0 (X 2 )|/k, and thus

|E[h(X 2 )] -E[h(X 1 )]| 1 k E[|∇ log π 0 (X 2 )|], (54)
and the same bound holds for the Wasserstein distance.

To show [START_REF] Ley | Stein's density approach and information inequalities[END_REF], by Proposition 3.13, there exists a weak solution

u h ∈ W 0 1,2 (p) to (48) such that |∇u| 2 p C p . It is straightforward to see that E [ ∇u h (X), ∇v(X) ] = E h(X)v(X) (55)
for any v ∈ W 1,2 (p i ). Indeed, by definition, Af = T ∇,p (∇f ) so that, by [START_REF] Ernst | Stein based goodness-of-fit tests[END_REF],

E [Af (X)g(X)] = E [T ∇,p (∇f (X))g(X)] = -E [ ∇f (X), ∇g(X) ] ; and (55)follows. Let h(x) = h(x) -E[h(X 1 )]. By 55, we have for any v ∈ W 1,2 (p), ∇u h , ∇v p 1 = -h v p 1 .

Applying this equation to

v = -π 0 ∈ W 1,2 (p 1 ), E[h(X 2 )-h(X 1 )] = ∇u h , ∇π 0 p 1 |∇u h | 2 p 1 |∇π 0 | 2 p 1 1/2 C p E[|∇π 0 (X 1 )| 2 ].
Remark 4.2. In dimension 1, in [START_REF] Ley | Distances between nested densities and a measure of the impact of the prior in Bayesian statistics[END_REF], Equation (4.2), gives

d W (X 1 , X 2 ) E[τ 1 (X 1 )|π 0 (X 1 )|], (56)
where τ 1 is the Stein kernel associated to X 1 . If p 1 admits a Poincaré constant C p , then from Proposition 3.13 (applied to h(x) = -x), we have that (τ (x)) 2 p 1 C 2 p . Thus by the Cauchy-Schwarz inequality, (56) is a stronger bound than (53).

Example 1: Copulas

Let (V 1 , V 2 ) be a 2-dimensional random vector, such that the marginals V 1 and V 2 have a uniform distribution on [0, 1]. We want to bound the Wasserstein distance between (V 1 , V 2 ) and its independent version (U 1 , U 2 ) (U 1 and U 2 are uniform and independent), in terms of the copula of (V 1 , V 2 ) defined as

C(x 1 , x 2 ) = P[V 1 x 1 , V 2 x 2 ], (x 1 , x 2 ) ∈ [0, 1] 2 .
(Note that the copula for

(U 1 , U 2 ) is (x 1 , x 2 ) → x 1 x 2 .) Assume that V 1 , V 2 has a pdf c = ∂ 2
x 1 x 2 C. An optimal Poincaré constant for the uniform distribution on [0, 1] 2 is C p = 2/π 2 , see [START_REF] Payne | An optimal Poincaré inequality for convex domains[END_REF]. Then, a simple application of (52) yields Corollary 4.3. Let (V 1 , V 2 ) have uniform marginals on [0, 1], and pdf c. Let (U 1 , U 2 ) also have uniform marginals, and let U 1 be independent of U 2 . Then

d W [(V 1 , V 2 ) , (U 1 , U 2 )] 2 π 2 [0,1] 2 |∇c(x 1 , x 2 )| 2 dx 1 dx 2 .
In some cases, one can compute the gradient of c in a closed form.

Example 4.4. The Ali-Mikhail-Haq copula [START_REF] Ali | A class of bivariate distributions including the bivariate logistic[END_REF] has pdf

c(x 1 , x 2 ) = (1 -θ){1 -θ(1 -x 1 )(1 -x 2 )} + 2θx 1 x 2 {1 -θ(1 -x 1 )(1 -x 2 )} 3 .
Here θ ∈ [-1, 1] is a measure of association between the two components V 1 and V 2 of the vector (V 1 , V 2 ) with uniform marginals each. If θ = 0 then the uniform copula (x 1 , x 2 ) → x 1 x 2 is recovered. Using Corollary 4.3 we can assess the Wasserstein distance between the Ali-Mikhail-Haq copula and the uniform copula in terms of θ. for -1 < θ < 1

[0,1] 2 |∇c(x 1 , x 2 )| 2 dx 1 dx 2 ≤ 128θ 2 1 {1 -|θ|} 8 . Hence d W ([(V 1 , V 2 ) , (U 1 , U 2 )]) ≤ 2.3 |θ|{1 -|θ|} -4 .
This bound shows the expected behaviour -it tends to 0 for θ → 0, whereas it diverges for |θ| → 1.

Example 2: Normal model with normal prior

Consider a normal N (θ, Σ 2 ) model with mean θ ∈ R d and positive definite covariance matrix Σ. The likelihood of a sample (x 1 , . . . , x n ) (where x i ∈ R d for all i) is given by L(x 1 , . . . , x n |θ) = i f (x i |θ) where

f (x|θ) = (2π) -d/2 det(Σ) -1/2 exp - 1 2 (x -θ) T Σ -1 (x -θ) .
To compare the posterior distribution P 1 of θ with uniform prior with the posterior P 2 with normal prior with parameters (µ, Σ 2 ) we employ the operator norm |||A||| = sup x =1 Ax . Corollary 4.5. Let P 1 denote the posterior distribution of θ with uniform prior and P 2 the posterior of θ with prior N (µ, Σ 2 ); Σ 2 is assumed positive definite. Then

d W (P 1 , P 2 ) |||Σ||| |||(Σ + nΣ 2 ) -1 ||| x -µ + √ 2Γ(d/2 + 1/2) Γ(d/2) |||Σ||| n |||(Σ 2 + nΣ 2 Σ -1 Σ 2 ) -1/2 |||.
Proof. It is a standard calculation that P 1 ∼ N (x, n -1 Σ), with x = 1 n n i=1 x i Thus (see Example 2.10 in [START_REF] Saumard | Log-concavity and strong log-concavity: a review[END_REF]) p 1 is strongly 1/λ-log concave,where λ is the greatest eigenvalue of Σ; λ is also the operator norm of Σ,. and ( 52) can be applied. From [START_REF] Ley | A general parametric Stein characterization[END_REF], we deduce

d W (P 1 , P 2 ) |||Σ||| n E[ ∇ log π 0 (X 2 ) ]
where X 2 is a r.v. with law P 2 . It remains to calculate this last expectation. It is again a standard calculation that P 2 ∼ N (μ, Σn ) with

μ = (Σ -1 2 + (n -1 Σ) -1 ) -1 ((n -1 Σ) -1 x + Σ -1 2 µ) = µ + n Σn Σ -1 (x -µ) (57) Σn = (Σ -1 2 + nΣ -1 ) -1 . Since p 2 (θ) ∝ p 1 (θ) exp[(θ -µ) T Σ -1 2 (θ -µ)], it holds that ∇ log π 0 (θ) = -Σ -1 2 (θ -µ).
Together with [START_REF] Liu | A kernelized Stein discrepancy for goodness-of-fit tests[END_REF], it follows that ∇ log π 0 (X 2 ) has a normal distribution with mean -nΣ -1 2 Σn Σ -1 (x -µ) and covariance matrix Σ -1 2 Σn Σ -1 2 . Note that

Σ -1 2 Σn Σ -1 = Σ -1 2 (Σ -1 2 + nΣ -1 ) -1 Σ -1 = (Σ + nΣ 2 ) -1
and, in the same way,

Σ -1 2 Σn Σ -1 2 = (Σ 2 + nΣ 2 Σ -1 Σ 2 ) -1 .
Let N stand for a d-dimensional standard normal vector. From the above, we deduce that ∇ log π 0 (X

2 ) ∼ -n(Σ+nΣ 2 ) -1 (x- µ) + (Σ 2 + nΣ 2 Σ -1 Σ 2 ) -1/2 N . Thus, |||Σ||| n E[ ∇ log π 0 (X 2 ) ] |||Σ||| |||(Σ + nΣ 2 ) -1 ||| x -µ + |||Σ||| n |||(Σ 2 + nΣ 2 Σ -1 Σ 2 ) -1/2 |||E[ N ].
Since N 2 ∼ χ 2 (d), N follows a chi distribution with d degrees of freedom and

E[ N ] = √ 2Γ(d/2+1/2) Γ(d/2)
. The assertion follows.

Remark 4.6. The bound in Corollary (4.5) has a contribution of order n -1 from the different covariance matrices, but also a contribution which depends on x -µ . By the law of large numbers, this term is order n -1 in probability, but crucially depends on the actual sample which was drawn.

Remark 4.7. When d = 1, we retrieve the bound of [START_REF] Ley | Distances between nested densities and a measure of the impact of the prior in Bayesian statistics[END_REF].

As in [START_REF] Ley | Distances between nested densities and a measure of the impact of the prior in Bayesian statistics[END_REF], our Stein framework lends itself naturally for assessing distributional distances between nested densities more generally than in the previous example. Proposition 4.8. Assume p 2 = p 1 π 0 , where p 1 and p 2 are densities which satisfy Assumption A. Let X 1 ∼ p 1 and X 2 ∼ p 2 denote two random variables with distributions having densities p 1 and p 2 , respectively. Then

E[τ τ τ 1 (X 1 )∇π 0 (X 1 )] d W (X 1 , X 2 ) sup f h |E[ ∇ log π 0 (X 2 ), ∇f h (X 2 ) )]| (58)
where the supremum is taken over all f h solving the p 1 -Stein equation (48) for h a 1-Lipschitz function.

Proof. Consider the lower bound first. Let e be a unit vector. Since x → x, e is 1-Lipschitz, using the nested structure we have

d W (X 1 , X 2 ) E[ X 1 , e -X 2 , e ] = E[ X 1 , e (1 -π 0 (X 1 ))] = -E[τ τ τ 1 (X 1 )∇π 0 (X 1 )], e . Taking e = -E[τ τ τ 1 (X 1 )∇(p 2 /p 1 )(X 1 )]/ E[τ τ τ 1 (X 1 )∇π 0 (X 1 )] gives d W (X 1 , X 2 ) E[τ τ τ 1 (X 1 )∇π 0 (X 1 )] .
For the upper bound, if h is 1-Lipschitz, using the score Stein operator for p 1 ,

E[h(X 2 ) -h(X 1 )] = E[(∆f h )(X 2 ) -∇ log p 1 (X 2 ), ∇f h (X 2 ) )] = E[ ∇ log π 0 (X 2 ), ∇f h (X 2 ) )].
The last equality follows from the nested structure p 2 = p 1 θ 0 .

Remark 4.9. The following argument shows that the gradient of the likelihood ratio between two densities arises naturally in the Stein framework. Suppose that q is another density on R d and that K p , the support of p, is a subset of K q , the support of q. Then the likelihood ratio r = p/q is well defined over R d (with the convention that r = 0 outside of K q ) and, for every f ∈ F(p) ∩ F(q), the product rule [START_REF] Barbour | Stein's method for diffusion approximations[END_REF] gives

T div,p f = div(f p) p = div(f q r) q 1 = div(f q) q + f , ∇r r = T div,q (f ) + f , ∇r r . ( 59
)
In combination with the ideas of Proposition 4.8, this leads to a general bound on the difference between densities in terms of the likelihood ratio.

Remark 4.10. Based on Remark 4.9, the comparisons of the posteriors arising from different priors can be carried out more generally. Here is a sketch. In the Bayesian framework, we aim to compare Θ 1 and Θ 2 obtained through different priors, where for θ (x) the likelihood function,

Θ 1 ∼ π 1 (θ) = κ 1 (x) θ (x) Θ 2 ∼ π 2 (θ) = κ 2 (x)π 0 (θ) θ (x)
for π 0 a nonnegative function with support a subset of R d , and κ i (x), i = 1, 2 the normalizing constants -here x = (x 1 , . . . , x J ) ∈ R J is a fixed sample of size J and θ ∈ R d is the variable. Then, in the notations of Remark 4.9, we have r(θ) = κ 2 (x)/κ 1 (x)π 0 (θ) so that, introducing ρ 0 = ∇π 0 /π 0 -which does not depend on the normalizing constants -we obtain, from (59),

T div,π 2 f = T div,π 1 f + f , ρ 0 ( 60
)
for any function f ∈ R n × R d ∈ F(π 1 ) ∩ F(π 2 )
, where n ≥ 1 is arbitrary. In particular,

E [T div,π 1 f (Θ 2 )] = E [ f (Θ 2 ), ρ(Θ 2 ) ] (61)
for all sufficiently regular vector-valued functions f . Next suppose that there exists some well chosen matrix valued function B for which π 1 is characterized by a second-order divergence Stein operator acting on real valued functions u : R d → R via

u → A 1 u = T div,π 1 (B∇u) = b, ∇u + B, ∇ 2 u HS (62)
with b = T div,π 1 B (recall (45)); then in light of (60), we know that

A 2 u = T div,π 2 (B∇u) = A 1 u + B∇u, ρ 0 (63)
satisfies E A 2 u(Θ 2 ) = 0 for all admissible functions u. In particular, if u h is a solution to the (second order) Stein equation

A 1 u h (θ) = h(θ) -E [h(Θ 1 )]
we get, with the help of (63) and under suitable conditions,

E [h(Θ 2 )] -E [h(Θ 1 )] = E [A 1 u h (Θ 2 )] = E [A 2 u h (Θ 2 )] -E [ B(Θ 2 )∇u h (Θ 2 ), ρ 0 (Θ 2 ) ] = -E [ B(Θ 2 )∇u h (Θ 2 ), ρ 0 (Θ 2 ) ] .
Investigating this approach in more detail will be part of future research.

Example 3: Comparing Azzalini-Dalla Valle skew-normal distributions to multivariate normal

The density of the Azzalini-Dalla Valle type r.v. X is given by

p α (x) = 2ω d (x -µ; Σ)Φ(α T x),
where ω d (x -µ : Σ) is the density of N (µ, Σ), Φ the c.d.f. of the standard normal on R, and α ∈ R d is a skew parameter, see [START_REF] Azzalini | The multivariate skew-normal distribution[END_REF]. Now we assume for simplicity Σ = I d and µ = 0.

In [START_REF] Ley | Stein's method for comparison of univariate distributions[END_REF], an exact expression for the Wasserstein distance is given for d = 1; here we extend this result to general d.

d W (X, Z) = 2 π α 1 + α 2 . ( 64
)
Proof. First we show the upper bound on the Wasserstein distance. As the multivariate standard normal distribution is 1-strongly log concave,(52) can be applied with π 0 (x) = 2Φ(α T x). As ∇π 0 (x) = 2αφ(α T x) with φ the one-dimensional standard normal pdf, it suffices to calculate that

E[φ(α T Z)] = (2π) -(d+1)/2 R d exp - 1 2 (α T x) 2 + x 2 dx = (2π) -(d+1)/2 (2π) d/2 det(M ) = 1 2π(1 + α 2 ) .
Thus,

d W (X, Z) 2 π α 1 + α 2 . ( 65
)
Now we prove we have actually equality. Consider the 1-Lipschitz test function h(x) = α α , x ; note that E[h(Z)] = 0. Then u h (x) = -h(x) is a solution to the Stein equation ∆u -x, ∇u(x) = h(x) and ∆u = 0, and

E[h(X) -h(Z)] = 2E[∆u h (X) -X, ∇u h (X) ] = 2E[ ∇u h (Z), αφ(α T Z) ] = 2E[ ∇h(Z), αφ(α T Z) ] = 2E 1 α α, αφ(α T Z) = 2 π α 1 + α 2 , so that d W (X, Z) 2 π α √ 1+ α 2 .
Thus we obtained (64).

Stein operators for elliptical distributions

In this section we detail the constructions explicitly for the entire family of elliptical distributions, defined as follows, see [START_REF] Landsman | Stein's lemma for elliptical random vectors[END_REF].

Definition 5.1. An absolutely continuous d-random vector has multivariate elliptical distribution E d (ν, Σ, φ) if its density is of the form

p(x) = κ|Σ| -1/2 φ 1 2 (x -ν) T Σ -1 (x -ν) , x ∈ R d , (66)
for φ : R + → R + a measurable function called density generator, ν ∈ R d the location parameter, κ the normalising constant and Σ = (σ ij ) a symmetric positive definite d × d dispersion matrix.

A particular important case is E d (0, I d , φ) called spherical distribution. Note that the matrix Σ in definition [START_REF] Oates | Control functionals for Monte Carlo integration[END_REF] is not necessarily the covariance matrix; also not all choices of φ lead to well-defined densities, see [START_REF] Landsman | Stein's lemma for elliptical random vectors[END_REF] or [START_REF] Paindaveine | Elliptical symmetry[END_REF] for a discussion and references. Prominent members of the elliptical family are 1. Gaussian distribution N d (ν, Σ), with φ(t) = e -t and κ = (2π) -d/2 .

Power exponential distribution

φ(t) = exp(-b p,ζ t ζ ) with ζ > 0, b p,ζ a scale factor
and κ defined accordingly, see [START_REF] Gómez-Sánchez-Manzano | A multivariate generalization of the power exponential family of distributions[END_REF][START_REF] Aerts | Robust asymptotic tests for the equality of multivariate coefficients of variation[END_REF] for details.

Multivariate Student-t distribution, with

φ(t) = (1 + 2t/k) -(k+d)/2 and κ = c k,d,Σ .
4. Symmetric generalized hyperbolic distribution with density

p(x) = √ χψ -λ ψ d/2 (2π) d/2 |Σ| K λ √ χψ K λ-d/2 (χ + (x -ν) T Σ -1 (x -ν)) ψ (χ + (x -ν) T Σ -1 (x -ν)) ψ (d/2)-λ (67)
where K λ denotes a modified Bessel function of the third kind and λ, χ, ψ are real parameters, see [START_REF] Mcneil | Quantitative risk management: Concepts, techniques and tools[END_REF]Example 3.8] for details. To put in the parameterization [START_REF] Oates | Control functionals for Monte Carlo integration[END_REF] we take φ(t

) = K λ-d/2 ( √ 2t)/( √ 2t) (d/2)-λ and κ = ( √ χψ) -λ ψ d/2 (2π) d/2 ( √ χψ) .
Throughout this section we let Ω φ be the image of Ω p through x → 1 2 (x-ν) T Σ -1 (x-ν). We make the following assumption.

Assumption B: There is a non-empty open set Ω φ such that φ > 0 and φ is C 1 on this open set, and such that the Lebesgue measure of K φ \Ω φ is zero.

Under Assumption B, it is readily checked that p satisfies Assumption A with Ω p defined by

Ω p = x ∈ R d : 1 2 (x -ν) T Σ -1 (x -ν) ∈ Ω φ = x ∈ R d : φ 1 2 (x -ν) T Σ -1 (x -ν) > 0 . ( 68
)
It is straightforward to calculate that if p is of the form (66) and satisfies Assumption B then

ρ p (x) = Σ -1 (x -ν) φ ((x -ν) T Σ -1 (x -ν)/2) φ((x -ν) T Σ -1 (x -ν)/2) ( 69
)
is the score function of p (defined on Ω p ). Hence the score-Stein operator ( 32) is easily obtained for this family of distributions.

Stein kernels for elliptical distributions

The following proposition shows that in order to find Stein kernels for elliptical distributions, it suffices to consider the case Σ = I d and ν = 0.

Proposition 5.2. The application

τ τ τ → (x → Σ 1/2 τ τ τ (Σ -1/2 (x -ν))Σ 1/2 ),
is a bijection between the set of Stein kernels of E d (0, I d , φ) and that of E d (ν, Σ, φ).

Proof. Let p (resp. q) be the density of

E d (0, I d , φ) (resp. E d (ν, Σ, φ)). If τ τ τ X is a Stein ker- nel for X ∼ E d (0, Id, φ), then for all f ∈ C ∞ c (Ω p ), E[τ τ τ X (X)∇f (X)] = E[Xf (X)]. Setting f (x) = g(Σ 1/2 x+ν), we deduce that for all g ∈ C ∞ c (Ω q ), we have E[τ τ τ X (X)Σ 1/2 ∇g(Σ 1/2 X+ ν)] = E[Xg(Σ 1/2 X + ν)]. Thus E[Σ 1/2 τ τ τ X (Σ -1/2 ((Σ 1/2 X +ν)-ν))Σ 1/2 ∇g(Σ 1/2 X +ν)] = E[(Σ 1/2 X +ν -ν)g(Σ 1/2 X +ν)]. It follows x → Σ 1/2 τ τ τ X (Σ -1/2 (x -ν))Σ 1/2 is a Stein kernel for Σ 1/2 X + ν.
The converse is shown in the same way. Now we state the following result, due to [START_REF] Landsman | Stein's lemma for elliptical random vectors[END_REF][START_REF] Landsman | Some Stein-type inequalities for multivariate elliptical distributions and applications[END_REF] but for which we give a new proof.

Lemma 5.3 (Proposition 2, [51]). If X ∼ E d (ν, Σ, φ) then the matrix τ τ τ 1 (x) = 1 φ((x -ν) T Σ -1 (x -ν)/2) +∞ (x-ν) T Σ -1 (x-ν)/2 φ(u)du Σ, ( 70
) is a Stein kernel for X if τ τ τ 1 ∈ F(p).
Proof. For transparency here we consider general ν and Σ. We will show the Stein identity holds for functions in C ∞ c (Ω p ). Consider any functions f : R → Ω φ ∈ F(p) and g : Ω p → R ∈ dom(p, f ). We start by inverting (69) to get

x -ν = Σρ p (x) φ((x -ν) T Σ -1 (x -ν)/2) φ ((x -ν) T Σ -1 (x -ν)/2) .
Fixing ν = 0 (see Proposition 5.2) and introducing the temporary notations ψ(t) = φ(t)/φ (t) and t = x T Σ -1 x/2, T = X T Σ -1 X/2, by ( 33),

E p [Xf (T )g(X)] = E p [Σρ p (X)ψ(T )f (T )g(X)] = -ΣE p [∇ {ψ(T )f (T )g(X)}] = -ΣE p [∇ {ψ(T )f (T )} g(X)] -ΣE p [ψ(T )f (T )∇g(X)] = -ΣE p ψ (T )f (T ) + ψ(T )f (T ) Σ -1 Xg(X) -ΣE p [ψ(T )f (T )∇g(X)]
and thus

E p f (T )(1 + ψ (T )) + ψ(T )f (T ) Xg(X) = -E p [ψ(T )f (T )Σ∇g(X)] . (71)
In order to obtain a Stein kernel, it suffices to choose f solution to the ODE

f (t)(1 + ψ (t)) + ψ(t)f (t) = -1
to ensure that the function x → ψ(t)f (t)Σ satisfies [START_REF] Goldstein | Zero biasing in one and higher dimensions, and applications[END_REF], and is a Stein kernel in the sense of Definition 3.8. Now note that the function u(t

) := 1 φ(t) +∞ t φ(u)du satisfies u = -1 ψ u -1; hence the choice f = u/ψ satisfies (f ψ) = (u) = -1 ψ u -1
τ τ τ p (x) = 1 φ((x -ν) 2 /2) +∞ (x-ν) 2 /2 φ(u)du = 1 p(x) +∞ (x-ν) 2 /2 p( √ 2u)du = 1 p(x) ∞ x (y -ν)p(y)dy.
The case x < 0 is treated similarly. Hence (70) indeed recovers the Stein kernel.

The identity (70) resulting from Lemma 5.3 has found many applications, [START_REF] Adcock | Extensions of Stein's lemma for the skew-normal distribution[END_REF][START_REF] Adcock | Asset pricing and portfolio selection based on the multivariate extended skew-Student-t distribution[END_REF][START_REF] Landsman | A note on Stein's lemma for multivariate elliptical distributions[END_REF][START_REF] Adcock | Mean-variance-skewness efficient surfaces, Stein's lemma and the multivariate extended skew-Student distribution[END_REF][START_REF] Landsman | Some Stein-type inequalities for multivariate elliptical distributions and applications[END_REF][START_REF] Vanduffel | A Stein type lemma for the multivariate generalized hyperbolic distribution[END_REF] and the references therein.

The following proposition gives a way of finding Stein kernels of a particular form which generalizes Lemma 5.3.

Proposition 5.5. Let a, b : Ω φ → R two C 1 functions such that, for all all u ∈ Ω φ , (a(t)φ(t)) φ(t) + 2t (b(t)φ(t)) φ(t) + (d + 1)b(t) + 1 = 0 ( 72
)
and set τ τ τ a,b (t) = a(t)Σ + b(t) (x -ν)(x -ν) T with t = 1 2 (x -ν) T Σ -1 (x -ν). If τ τ τ a,b ∈ F(p) then this function is a Stein kernel for X ∼ E d (ν, Σ, φ).
Moreover, if φ is continuous and positive on [0, +∞), then the Stein identity (41) holds for every test function g ∈ C 1 (R d ) such that g(x)a(t)φ(t) and g(x)b(t)t go to zero when |x| → +∞.

Proof. From Proposition 5.2, we can assume ν = 0 and Σ = I d . It is readily checked that div(I d ) = 0, div(xx T ) = (d + 1)x, xx T x = 2tx. Thus, by the chain rule and noting that ∇t = x,

div(φ(t)(a(t)I d + b(t)xx T )) = (a(t)φ(t)) I d x + a(t)φ(t)div(I d ) + (b(t)φ(t)) xx T x + b(t)φ(t)div(xx T ) = (a(t)φ(t)) x + 2t(b(t)φ(t)) x + (d + 1)b(t)φ(t)x.
Hence a(t)I d + b(t)xx T is a Stein kernel if the last quantity is equal to -φ(t)x, so that the result follows.

To see that a test function satisfying the stated conditions verifies [START_REF] Goldstein | Stein's method and the zero bias transformation with application to simple random sampling[END_REF], simply note that since (x -ν)(x -ν) T is of order t when |x| is large, the conditions imply that τ τ τ (x)p(x)f (x) → 0 when |x| is large, and the result follows from the divergence theorem (see discussion below Definition 3.8).

Remark 5.6. 1. A particular instance of Proposition 5.5 is given by the following expression:

τ τ τ 2,β (x) = (β + 2) -2 φ (t)/φ (t) φ (t)/φ(t) (β -2)(d -1)   2   d -1 (β + 2) φ (t) φ(t) -2 φ (t) φ (t) + t   Σ -(x -ν)(x -ν) T   (73
) is a Stein kernel for X for all β = 2 as long as τ τ τ 2,β ∈ F(p).

2. It is straightforward to generalize the previous proposition in the following way. Here without loss of generality we take Σ = I d and ν = 0. Assume we are given matrices U 1 , . . . , U m such that for every i = 1, . . . , m and some functions α i , β i : R → R,

divU i = α i (t)x, U i x = β i (t)x. If 1 + m i=1 a i α i + (a i φ) φ β i = 0, then a 1 (t)U 1 (t) + . . . + a m (t)U m (t) is a Stein kernel for X ∼ E d (0, I d , φ) if this function is in the class F(p).
By setting b ≡ 0, we obtain a(t) = 1 φ(t) +∞ t φ(u)du, and we retrieve Lemma 5.3. Setting a ≡ 0 leads to the following

Corollary 5.7. Set t = 1 2 (x -ν) T Σ -1 (x -ν). If +∞ 0 u d-1 2 φ(u)du < ∞, then t -d+1 2 2φ(t) +∞ t u d-1 2 φ(u)du (x -ν)(x -ν) T , ( 74
) is a Stein kernel for X ∼ E d (ν, Σ, φ). Moreover, if Ω p = R d , the Stein identity (41) holds for every function f ∈ C 1 (R d ) such that f (x)t d-1 2 is bounded.
Proof. With a ≡ 0, (72) becomes

(bφ) + d + 1 2t bφ = - φ 2t , which integrates to b(t) = t -d+1 2 2φ(t) +∞ t u d-1 2 φ(u)du. Now if f is as stated in the corollary, since (x -ν)(x -ν) T is of order t for large |x|, then τ (x)f (x)p(x) = O( +∞ t u d-1
2 φ(u)du) for large |x|, so that τ τ τ (x)f (x)p(x) goes to zero at infinity, and the Stein identity follows again from the divergence theorem.

Remark 5.8. For d = 1, (74) leads to the classical Stein kernel (4). Indeed, assuming ν = 0, for x > 0,

t -1 2φ(t) +∞ t φ(u)du = 1 x 2 p(x) +∞ x 2 /2 p(2 √ u)du = 1 x 2 p(x) +∞ x s p(s)ds,
and multiplying by x 2 yields the claim. The case x < 0 is treated similarly.

In the next three subsections we develop Stein kernels for three distributional families: the multivariate Gaussian, the power exponential, and the multivariate Student t-distribution. Similar computations are possible for the symmetric generalized hyperbolic distribution but are not pursued here. We refer to [START_REF] Landsman | Some Stein-type inequalities for multivariate elliptical distributions and applications[END_REF][START_REF] Vanduffel | A Stein type lemma for the multivariate generalized hyperbolic distribution[END_REF] and the references therein.

The multivariate Gaussian distribution

Consider a Gaussian d-dimensional random vector Z ∼ N d (ν, Σ) with pdf ϕ on R d and let µ(dx) = ϕ(x)dx be the corresponding probability measure. As Z ∼ E d (ν, Σ, φ) with φ(t) = e -t and φ (t)/φ(t) = -1, from (69), we recover that ρ ϕ (x) = -Σ -1 (x -ν) is the score function of ϕ. Since 1 φ(t) ∞ t φ(u)du = 1 for all t, Lemma 5.3 shows that τ τ τ 1 = Σ is, as is well-known, a Stein kernel for ϕ. Moreover, [START_REF] Ross | Fundamentals of Stein's method[END_REF] gives, after some simplifications, the following family of Stein kernels which are indexed by β = 2 (we set ν = 0 to save space):

τ τ τ 2,β (x) = β (β -2)(d -1) x T Σ -1 xΣ -xx T - 2 β -2 Σ.
It is easy to check that these functions are in the class F(p). Several interesting cases stand out. Sending β to 0, on the one hand, and to +∞ on the other hand, we obtain

τ τ τ 2,0 (x) = Σ and τ τ τ 2,∞ (x) = 1 d -1 x T Σ -1 xΣ -xx T . In dimension d ≥ 3, setting β = 2(d -1) we get τ τ τ 2,2(d-1) (x) = 1 d -2 x T Σ -1 xΣ -xx T -Σ
Moreover, we find that the Gaussian multivariate normal satisfies

β (β -2)(d -1) E X T Σ -1 XΣ -XX T ∇g(X) = E 2 β -2 Σ∇g(X) + Xg(X) (75)
for all g ∈ dom(ϕ, τ τ τ 2,β ) and all β = 2. In particular, every g ∈ C 1 (R d ) with at most polynomial growth at infinity lies in this domain.

Power exponential distribution

Consider a d random vector Z ∼ PE d,ζ (ν, Σ) distributed according to the multivariate power exponential distribution with power ζ > 0, location µ, scale b, shape

Σ ∈ R d × R d and pdf ϕ ζ (x) = a d,ζ |Σ| -1/2 exp -b (x -ν) T Σ -1 (x -ν) ζ (76)
on R d (a d,ζ is the normalizing constant), ζ ∈ (0, ∞), and let µ(dx) = ϕ ζ (x)dx be the corresponding probability measure.

Clearly Z ∼ E d (ν, Σ, φ) with φ(t) = e -bt ζ so that φ (t)/φ(t) = -bζt ζ-1 and φ (t)/φ (t) = -bζt ζ-1 + ζ-1 t . From (69), the score function of ϕ ζ is ρ ζ (x) = -2bζ((x -ν) T Σ -1 (x -ν)) ζ-1 Σ -1 (x -ν). (77)
While, except when ζ = 1,the kernel from (5.3) does not lead to palatable expressions, applying [START_REF] Ross | Fundamentals of Stein's method[END_REF] we obtain (for ζ = 1)

τ τ τ 2,ζ (x) = β + 2 ζ-1 bζt ζ (β -2)(d -1) 1 - d -1 βbζt ζ + 2(ζ -1) (x -ν) T Σ -1 (x -ν)Σ -(x -ν)(x -ν) T .
(78) These functions are Stein kernels : they are in the class F(p), since τ τ τ 2,ζ ϕ ζ → 0 when x → +∞. %eqrefeq:26. We do not provide details here.

Note again that the Stein identity (41) holds for every g ∈ C 1 (R d ) with at most polynomial growth at infinity, since in this case we have τ τ τ 2,ζ g ϕ ζ → 0 at infinity.

The multivariate Student t-distribution with k > 1

Consider a d random vector X ∼ t k (ν, Σ) distributed according to the multivariate Student-t distribution with k > 1 degrees of freedom, location ν ∈ R d , shape Σ ∈ R d × R d and pdf

t k (x) = c k,d |Σ| -1/2 1 + (x -ν) T Σ -1 (x -ν) k -(k+d)/2 (79)
with normalizing constant c k,d,Σ = Γ((k + d)/2)/(Γ(k/2)k d/2 π d/2
). Let µ(dx) = t k (x)dx be the corresponding probability measure. The assumption that k > 2 ensures that this distribution has finite mean and finite variance. This distribution is an elliptical distribution with φ(t) = (1 + 2t/k) -(k+d)/2 and hence φ (t)/φ(t) = -(k + d)/(k + 2t) leading to the score function

ρ t k (x) = -(k + d) Σ -1 (x -ν) k + (x -ν) T Σ -1 (x -ν) . (80)
Moreover,

1 φ(t) +∞ t φ(u)du = k + 2t d + k -2 (from k > 1 it follows that d + k > 2) and hence Lemma 5.3 gives that τ τ τ 1 (x) = (x -ν) T Σ -1 (x -ν) + k d + k -2 Σ (81)
is a Stein kernel for the multivariate Student distribution for k > 1, as then τ τ τ 1 ∈ F(t k ).

Similarly, using that φ (t)/φ (t) = -(d + k + 2)/(k + 2t), Lemma 73 gives a family of Stein kernels which are indexed by β ∈ R:

τ τ τ 2,β (x) = β(d + k) -4 (d + k)(β -2)(d -1) 2 (d -1)(k + 2t) 4 -β(d + k) + t Σ -(x -ν)(x -ν) T = β(d + k) -4 (d + k)(β -2)(d -1) 2 (d -1)k + t (2(d + 1) -β(d + k))) 4 -β(d + k) Σ -(x -ν)(x -ν) T .
It is easy to verify that τ τ τ 2,β ∈ F(t k ). If we choose β so that (d + k)β = 2(d + 1), i.e. β = 2(d + 1)/(d + k) then, after simplifications, we obtain for k > 2

τ τ τ 2 (x) = 1 k -1 (x -ν)(x -ν) T + kΣ . (82)
Since τ τ τ 1 and τ τ τ 2 are of order t for large |x|, and since φ(t) ∼ t→+∞ t -(k+d)/2 , the Stein identity (41) holds for every g ∈ C 1 (R d ) such that t -(k+d-2)/2 g(x) → 0 at infinity. In particular constant functions verify [START_REF] Goldstein | Stein's method and the zero bias transformation with application to simple random sampling[END_REF] and the Stein kernels are in F(p). Note that both τ τ τ 1 and τ τ τ 2 simplify to τ (x) = (x 2 + kσ 2 )/(k -1) when d = 1; this last quantity is well-known to be the univariate kernel for the Student-t distribution with k degrees of freedom and centrality parameter ν, see e.g. [56, page 30].

There are several types of operators and identities that can be obtained; below are some examples.

1. Vector valued operators. Applying the product rule [START_REF] Döbler | An iterative technique for bounding derivatives of solutions of Stein equations[END_REF]

directly with f (x) = k +(x-ν) T Σ -1 (x-ν) we obtain for g : R d → R ∈ dom(f, t k ) the operator A t k g(x) = (k + (x -ν) T Σ -1 (x -ν))∇g(x) + (2 -k -d)Σ -1 (x -ν)g(x).
Taking expectations for X ∼ t k (ν, Σ) we obtain the vector-Stein identity

E (k + (X -ν) T Σ -1 (X -ν))∇g(X) = (k + d -2)E Σ -1 (X -ν)g(X) (83)
By definition of the Stein kernel we also get new Stein operators and identities. Using τ τ τ 1 recovers (83), whereas using τ τ τ 2 we obtain

E ((X -ν)(X -ν) T + kΣ)∇g(X) = (k -1)E [(X -ν)g(X)] (84)
(still with X ∼ t k (ν, Σ)).

2. Scalar valued operators. Suppose for simplicity that Σ = I d and ν = 0. Taking B successively equal to τ τ τ 1 then τ τ τ 2 in (44) leads to

B 1 g(x) = -∇g(x), x + 1 d + k -2 x T x + 2k, ∇ 2 g(x) HS B 2 g(x) = -∇g(x), x + 1 k -1 xx T + kI d , ∇ 2 g(x) HS
acting on functions g such that ∇g ∈ F 1 (t k ).

3.

A covariance identity. Starting from (84) with

X = (X 1 , X 2 ) T ∈ R d 1 × R d 2
multivariate student with location (ν 1 , ν 2 ) and shape

Σ 11 Σ 12 Σ 21 Σ 22 , taking g(x) = g(x 2
) and considering only the first d 1 components of the resulting identity we obtain

E ((X 1 -ν 1 )(X 2 -ν 2 ) T + kΣ 12 )∇g(X 2 ) = (k -1)E [(X 1 -ν 1 )g(X 2 ).] (85)
Many more such covariance identities can be obtained by this approach, thus complementing those obtained in [START_REF] Adcock | Extensions of Stein's lemma for the skew-normal distribution[END_REF].

Generalities on Stein kernels

Let T div,p be the canonical Stein operator (34) acting on F(p), the corresponding Stein class. In this section we explore properties of the Stein kernels from Definition 3.8. Proposition 6.1 (Properties). Let τ p,i be a Stein kernel for p in the direction e i and τ τ τ p the matrix with i th raw being τ p,i . Then

1. For all j = 1, . . . , d we have

∂ ∂x j R d τ p,i (x)p(x)dx = R d ∂ ∂x j (τ p,i (x)p(x))dx = 0.
2. If p admits a second moment, and if x -ν ∈ dom(p, τ τ τ p ), then

E [τ τ τ p (X)] = Var(X).
Proof. The first statement follows by the requirement that the kernel belongs to F(p), which in particular imposes that all components of τ p,i belong to F 1 (p). To see the second claim, taking expectations in [START_REF] Gan | Dirichlet approximation of equilibrium distributions in Cannings models with mutation[END_REF] yields that the Stein kernel necessarily satisfies

E [τ τ τ p (X)∇g(X)] = E [(X -ν)g(X)] (86)
for all g : R d → R belonging to dom(p, τ τ τ p ). By assumption, g i (x) = x i -ν i belongs to dom(p, τ τ τ p ) for all i = 1, . . . , d. Plugging these functions in (86) leads to

E [(τ τ τ p,i (X)) j ] = E [(X i -ν i )(X j -ν j )]
for all i, j = 1, . . . , d. The claim follows.

Proposition 6.2. Given k ≤ d and {i 1 , . . . , i k } ⊂ {1, . . . , n} denote by V = e i 1 , . . . , e i k the space generated by e i 1 , . . . , e i k . Also, write any x ∈ R d as x = (x V , x V ⊥ ) and let p V = V ⊥ p be the marginal of p on V. Suppose that p admits a p-integrable Stein kernel τ i j = τ i j ,1 , . . . , τ i j ,d in each direction e i 1 , . . . , e i k . Then the vector τ

V i j = τ V i j ,1 , . . . , τ V i j ,k
with components

τ V i j , (x V ) = E τ i j ,i (X) | X V = x V , = 1, . . . , k (87)
is a Stein kernel for p V in the direction e i j .

Proof. Without loss of generality we suppose that p is centered. Fix

x V ∈ V. Then k =1 ∂ ∂ i τ V i j , (x V )p V (x V ) = k =1 ∂ ∂ i V ⊥ τ i j ,i (x V , x V ⊥ )p(x V , x V ⊥ )dx V ⊥ = k =1 V ⊥ ∂ ∂ i τ i j ,i (x V , x V ⊥ )p(x V , x V ⊥ ) dx V ⊥ = -x i j V ⊥ p(x V , x V ⊥ )dx V ⊥ + m∈{1,...,n}\{i 1 ,...,i k } V ⊥ ∂ ∂ m τ i j ,m (x V , x V ⊥ )p(x V , x V ⊥ ) dx V ⊥
where the second-last line is allowed thanks to integrability of the Stein kernel and the last follows from [START_REF] Gallouët | Regularity of solutions of the Stein equation and rates in the multivariate central limit theorem[END_REF], giving the identity

d m=1 ∂ ∂ m τ i j ,m (x V , x V ⊥ )p(x V , x V ⊥ ) = -x i j p(x V , x V ⊥ )
which is valid for any x V ⊥ such that (x V , x V ⊥ ) lies in the support of Ω p . Now using the

requirement τ i j ∈ F(p) ensures that V ⊥ ∂ ∂m τ i j ,m (x V , x V ⊥ )p(x V , x V ⊥ ) dx V ⊥ = 0 for all m / ∈ {i 1 , . . . , i k } so that k =1 ∂ ∂ i τ V i j , (x V )p V (x V ) = -x i j p V (x),
as required. It remains to check that τ V i j , (x V ) ∈ F(p V ), but this is a direct consequence of the definitions.

Next, we provide formulas for computing Stein kernels explicitly based on univariate Stein kernels. Proposition 6.3 (Bivariate Stein kernels). Let X = (X 1 , X 2 ) T ∼ p with p a continuous pdf on R 2 satisfying Assumption A. Let p 1 be the marginal of p in direction e 1 , ρ 1 the corresponding univariate score and τ 1 the corresponding univariate kernel (which we suppose to exist). Set τ 11 (x 1 , x 2 ) = τ 1 (x 1 ), and

τ 12 (x 1 , x 2 ) = τ 1 (x 1 ) p 1 (x 1 ) p(x 1 , x 2 ) ∂ 1 ∞ x 2 p(x 1 , v)dv/p 1 (x 1 ) . (88)
Then the vector (x 1 , x 2 ) → (τ 1 (x 1 , x 2 ), τ 12 (x 1 , x 2 )) 1≤i,j≤2 is a Stein kernel for p in the direction e 1 .

Note that

τ 1 (x 1 ) p 1 (x 1 ) p(x 1 , x 2 ) ∂ 1 ∞ x 2 p(x 1 , v)dv/p 1 (x 1 ) = τ 1 (x 1 ) p 1 (x 1 ) p(x 1 , x 2 ) - p 1 (x 1 ) (p 1 (x 1 )) 2 ∞ x 2 p(x 1 , v)dv + 1 p 1 (x 1 ) ∞ x 2 ∂ 1 p(x 1 , v)dv = 1 p(x 1 , x 2 ) τ 1 (x 1 ) ∞ x 2 -ρ 1 (x 1 )p(x 1 , v) + ∂ 1 p(x 1 , v) dv. (89)
This alternative form of (88) is often more convenient than (88).

Proof. We need to prove that

2 j=1 ∂ j (τ 1j (x)p(x)) = -(x 1 -ν 1 )p(x) (90)
for all x = (x 1 , x 2 ) T ∈ R 2 . Applying (89) we have

τ 12 (x)p(x 1 , x 2 ) = p 1 (x 1 )τ 1 (x 1 )∂ 1 ∞ x 2 p(x 1 , v)dv/p 1 (x 1 ) = -τ 1 (x 1 ) p 1 (x 1 ) p 1 (x 1 ) ∞ x 2 p(x 1 , v)dv + τ 1 (x 1 ) ∞ x 2 ∂ 1 p(x 1 , v)dv (91) so that ∂ 2 (τ 12 (x)p(x 1 , x 2 ))) = τ 1 (x 1 ) p 1 (x 1 ) p 1 (x 1 ) p(x 1 , x 2 ) -τ 1 (x 1 )∂ 1 p(x 1 , x 2 ). (92)
Also, as ∂ 1 (p 1 (x 1 )τ 1 (x 1 )) = (ν 1 -x 1 )p 1 (x 1 ),

∂ 1 (τ 1 (x)p(x 1 , x 2 )) = ∂ 1 p 1 (x 1 )τ 1 (x 1 ) p(x 1 , x 2 ) p 1 (x 1 ) = ∂ 1 (p 1 (x 1 )τ 1 (x 1 )) p(x 1 , x 2 ) p 1 (x 1 ) + p 1 (x 1 )τ 1 (x 1 )∂ 1 p(x 1 , x 2 ) p 1 (x 1 ) = -(x 1 -ν 1 )p(x 1 , x 2 ) + p 1 (x 1 )τ 1 (x 1 ) - p 1 (x 1 ) (p 1 (x 1 )) 2 p(x 1 , x 2 ) + ∂ 1 p(x 1 , x 2 ) p 1 (x 1 ) = -(x 1 -ν 1 )p(x 1 , x 2 ) -τ 1 (x 1 ) p 1 (x 1 ) p 1 (x 1 ) p(x 1 , x 2 ) + τ 1 (x 1 )∂ 1 p(x 1 , x 2 ). (93)
Adding up ( 92) and (93) we get (90).

Remark 6.4. The proof of Proposition 6.3 is of purely computational nature. The inspiration for formula (88) is [11, equation (9)], where a similar quantity is introduced via a transport argument. To see the connection, note that (89) gives

τ 12 (x 1 , x 2 )p(x 1 , x 2 ) = τ 1 (x 1 ) x 2 -∞ (ρ 1 (x 1 )p(x 1 , v) -∂ 1 p(x 1 , v)) dv. (94)
We introduce p X | X i =x i (x 1 , x 2 ) = p(x 1 , x 2 )/p i (x i ) the conditional density of X at X i = x i . Fix i = 1 and, for each t, t , x 2 let x 2 → T t,t (x 2 ) be the mapping transporting the conditional density at x 1 = t to that at x 1 = t , implicitly defined via

p X | X 1 =t (x 2 ) = p X | X 1 =t T t,t (x 2 ) ∂ x 2 T t,t (x 2 ). ( 95
)
Taking derivatives in (95) with respect to t and setting t = t = x 1 we deduce (using the fact that T t,t (x 2 ) = x 2 ) that

0 = ∂ 1 p(x 1 , x 2 ) p 1 (x 1 ) + ∂ 2 p(x 1 , x 2 ) p 1 (x 1 ) ∂ t T t,t (x 2 ) t =t=x 1 - p(x 1 , x 2 ) p 1 (x 1 ) p 1 (x 1 ) p 1 (x 1 ) + p(x 1 , x 2 ) p 1 (x 1 ) ∂ 2 (∂ t T t,t (x 2 )) t =t=x 1 ,
that is,

∂ 2 p(x 1 , x 2 )∂ t T t,t (x 2 ) t =t=x 1 = p(x 1 , x 2 ) p 1 (x 1 ) p(x 1 ) -∂ 1 p(x 1 , x 2 )
and we recognize from (94) that, up to a function which depends only on x 1 ,

τ 12 (x 1 , x 2 ) τ 1 (x 1 ) = ∂ t T t,t (x 2 ) t =t=x 1 . ( 96
)
This is not the only Stein kernel in connection with transport maps, see [START_REF] Fathi | Stein kernels and moment maps[END_REF] where yet another construction is introduced.

Further, inspired by [START_REF] Artstein | Solution of Shannon's problem on the monotonicity of entropy[END_REF][START_REF] Artstein | On the rate of convergence in the entropic central limit theorem[END_REF], we can directly postulate our next result which guarantees existence of Stein kernels under smoothness conditions. Theorem 6.5. Let p : R d → (0, ∞) be a continuously twice differentiable density on R

d with ∇p 2 p , ∇ 2 (p) < ∞. Let τ (1)
i , i = 1, . . . , d be the marginal Stein kernels. Then, for any direction e i , i = 1, . . . , d there exists a Stein kernel for p in direction e i τ (d)

p,i (x) = τ (1) i (x i ) τ (d) i,1 (x | x i ) • • • τ (d) i,i-1 (x | x i ) 1 τ (d) i,i+1 (x | x i ) • • • τ (d) i,d (x | x i ) with coefficients τ i• , i = 1, . . . , d solving the equations T div,p (τ i• (x | x i )) = ρ i (x i ). ( 97
)
Here ρ i (x i ) = p i (x i )/p i (x i )
is the score function of the ith marginal and x = (x 1 , . . . , x d ).

Proof. The result is almost immediate from [START_REF] Artstein | Solution of Shannon's problem on the monotonicity of entropy[END_REF]Theorem 4] where it is proved (see middle of page 978) that, under the stated conditions, there exist continuously differentiable vector fields R h such that

div(R h (x)p(x)) p(x) = h (u) h(u)
for any marginal u → h(u) of p, in any direction. Collecting these into a single vector and adapting the notations leads to (97). To see the connection with Stein kernels, write

τ ij (x) = τ i (x i )τ ij (x | x i ). (98)
Then

d j=1 ∂ j (τ ij (x)p(x)) = d j=1 ∂ j (τ ij (x | x i )p(x)τ i (x i )) = d j=1 ∂ j (τ ij (x | x i )p(x))τ i (x i ) + d j=1 τ ij (x | x i )p(x)∂ j (τ i (x i )) = ρ i (x i )p(x)τ i (x i ) + p(x)∂ i (τ i (x i ))
where in the last line we use (97) in the first sum and ∂ j (τ i (x i )) = 0 for all j = i in the second sum. Clearly by the definition of the univariate Stein kernel

∂ i (τ i (x i )) = -ρ i (x i )τ i (x i ) + E[X i ] -x i
and the claim follows.

Remark 6.6. Theorem 6.5 gives a mechanism allowing to generalize the bivariate construction from Proposition 6.3. Under the conditions of Theorem 6.5, for d = 3, we

set p ij (x i , x j ) = ∞ -∞ p(x)dx k and P i (x) = x -∞ p i (v)dv, i = 1, 2, 3. Then we can choose τ i,i = τ i and τ (3) 1,2 (x | x 1 )p(x 1 , x 2 , x 3 ) = x 2 -∞ (ρ 1 (x 1 )p(x 1 , v, x 3 ) -∂ 1 p(x 1 , v, x 3 )) dv -P 2 (x 2 ) ρ 1 (x 1 )p 13 (x 1 , x 3 ) -∂ 1 p 13 (x 1 , x 3 ) τ (3) 1,3 (x | x 1 )p(x 1 , x 2 , x 3 ) = p 2 (x 2 ) x 3 -∞ (ρ 1 (x 1 )p 13 (x 1 , w) -∂ 1 p 13 (x 1 , w)) dw
and similarly for τ

(3) i,j (x | x i ) for all i, j. Direct computations suffice for this claim. Also, setting τ

(2)

12 (x 1 , x 2 ) = τ (1) 1 (x i )E[τ (3) 12 (x | x 1 ) | X 1 = x 1 , X 2 = x 2 ], (99)
the vector (τ

1 (x 1 ), τ

12 (x 1 , x 2 )) forms a bivariate Stein kernel for (X 1 , X 2 ). Moreover,

+∞ -∞ τ (3) 1,2 (x | x 1 ) p(x 1 , x 2 , x 3 ) p(x 1 , x 2 ) dx 3 = 1 p(x 1 , x 2 ) x 2 -∞ (ρ 1 (x 1 )p(x 1 , v) -∂ 1 p(x 1 , v)) dv -P 2 (x 2 ) ρ 1 (x 1 )p 1 (x 1 ) -∂ 1 p 1 (x 1 ) = 1 p(x 1 , x 2 ) x 2 -∞ (ρ 1 (x 1 )p(x 1 , v) -∂ 1 p(x 1 , v)) dv
which is equivalent to the expression (89).

The k-variate extension can also be constructed, as follows. For all k ≥ 1, and under the same conditions, for all 1 ≤ j ≤ d -1 we can define

τ (k) 1,j (x | x 1 )p(x) = x j -∞ (ρ 1 (x 1 )p(x 1 , x j = v, . . . , x d ) -∂ 1 p(x 1 , x j = v, . . . , x d )) dv -P j (x j ) (ρ 1 (x 1 )p(x 1 , x j+1 . . . , x d ) -∂ 1 p(x 1 , x j+1 , . . . , x d )) (100) and for j = d τ (k) 1,d (x | x 1 )p(x) = P d-1 (x d 1 ) x d -∞ (ρ 1 (x 1 )p(x 1 , x d = v) -∂ 1 p(x 1 , x d = v)) dv (101) Example 6.7 (Multivariate Gaussian). If X ∼ N 2 (ν, Σ) is multivariate d-dimensional
Gaussian then direct computations of the kernel as provided by Proposition 6.3 leads to τ τ τ 2 given in (82). The expression is more complicated in dimension d ≥ 3, and so far we have not been able to give a probabilistic interpretation of it.

Example 6.8. If X = (X 1 , X 2 ) T ∼ t k (ν, Σ) follows the bivariate Student distribution then direct computations of the kernel as provided by Proposition 6.3 leads to τ τ τ 2 given in (82). Again, we have not been able to give a probabilistic interpretation of the expression in dimension d ≥ 3.

Stein discrepancies

Instead of using the Wasserstein metric which uses Lipschitz functions, more general classes functions G in sup g∈G |E[A p g(Y )]| can be useful to assess distributional distances, leading to the notion of Stein discrepancies.

Integral probability metrics and Stein discrepancies

Differences between distributions can be measured using probability metrics. For applying Stein's method, so-called integral probability metrics are well suited.

d H (F, G) := sup h∈H F (h) -G(h) (103)
defines a distance between probability distributions. This intuition leads to the following general definition.

Definition 7.2 (Stein discrepancy).

Let p be a density on R d and A p a Stein operator acting on some class dom(A p ). Then for any random Y ∼ q, any G ⊂ dom(A p ) and any norm • , the quantity

S • (q, A p , G) = sup g∈G E [A p g(Y )] (104) is the ( • -G -A p ) Stein discrepancy from Y to X.
Definition 7.2 is motivated by the reference [START_REF] Gorham | Measuring sample quality with Stein's method[END_REF] where, to the best of our knowledge, such a unified notation for general Stein-based discrepancies (with freedom of choice both in the operator and the class of functions) is first introduced. The choice of norm • is generally fixed by context such as dimensionality, basic properties of the operator and the random variables X, Y under study. In the sequel we will generally drop the indexation in the norm and simply write S(Y, A p , G) instead. The next subsection links Stein discrepancies to information metrics.

Our set-up allows for a combination of the operators from (105) and (106) which is sometimes more suitable for the problem at hand. All classes of functions are designed to ensure that, for all g ∈ G = F(A p ) ∩ F(A q ), we have

E [A p g(Y )] = E [A p g(Y )] -E [A q g(Y )] = E [(A p (Y ) -A q (Y )) ∇g(Y )] + E [(a p (Y ) -a q (Y ))g(Y )]
=: E A p/q (Y )∇g(Y ) + E a p/q (Y )g(Y ) .

(109)

Thus in particular we can take A p/q and a p/q in (105) (110)

-A p/q anda p/q in (106).

(111)

Two particular choices of input matrices A p and A q stand out:

• A p (x) = τ τ τ p (x) and A q (x) = τ τ τ q (x) for which a p/q = 0 and (109) becomes E [A p g(Y )] = E [(τ τ τ p (Y ) -τ τ τ q (Y ))∇g(Y )] =: E τ τ τ p/q (Y )∇g(Y ) (112)

• A p (x) = A q (x) = I d for which A p/q = 0 and a p/q = ρ p -ρ q and (109) becomes

E [A p g(Y )] = E [(ρ p (Y ) -ρ q (Y ))g(Y )] =: E ρ p/q (Y )g(Y ) . (113)
The Stein operators (105) and (106) can also be applied jointly to functions k(x, x ) with marginals in G.

A T p,1 A p,2 k(x, x ) = A T p (x)A q (x )∇ T x ∇ x k(x, x ) + A T p (x)a q (x )∇ T x k(x, x ) + a p (x) T ∇ x A q (x )k(x, x ) + a p (x) T a q (x )k(x, x ).

In [START_REF] Chwialkowski | A kernel test of goodness of fit[END_REF], such kernelized expressions are studied for the score function choice (113). Inspired by [START_REF] Chwialkowski | A kernel test of goodness of fit[END_REF][START_REF] Liu | A kernelized Stein discrepancy for goodness-of-fit tests[END_REF] we give the following result which follows immediately from (109) with the choices (110) and (111). This result shows that kernelized methods have a larger range of applicability than usually assumed. At the same time, it illustrates the power of our general set-up. Theorem 7.4. Let Y, Y be independently drawn from q on the same space and consider functions k : R d × R d → R such that both y → k(y, •) and y → k(•, y ) belong to F(A p ) ∩ F(A q ). Then E A T p,1 A p,2 k(Y, Y ) = E ∇ T y A p/q (Y )A p/q (Y ) ∇ y k(y, y )| y=Y,y =Y + E ∇ T y A p/q (Y )a p/q (Y )k(y, Y ) | y=Y + E a p/q (Y ) T ∇ y A p/q (Y )k(y, Y ) | y=Y + E a p/q (Y ) T k(Y, Y )a p/q (Y ) (114)

with L i , i = 1, 2 the operator L applied with respect to the ith variable of the function k(•, •).

Theorem 7.4, with the particular choice that A p (x) = A q (x) = I d , so that A p/q = 0 and a p/q = ρ p -ρ q , was proposed in [START_REF] Liu | A kernelized Stein discrepancy for goodness-of-fit tests[END_REF]. If both p and q are known up to a normalising constant, then (114) does not depend on the normalising constant and hence, again, is particularly attractive in high-dimensional Bayesian inference.

From ( 112) and (113) we introduce the following two special cases of particular interest.

1. If A p (x) = τ τ τ p (x) and A q (x) = τ τ τ q (x) for which a p/q = 0 then E ∇ T y τ τ τ p/q (Y )τ τ τ p/q (Y ) ∇ y k(t, t )| y=Y,y =Y = E A T p,1 A p,2 k(Y, Y ) .

(115)

2. If A p (x) = A q (x) = I d for which A p/q = 0 and a p/q = ρ p -ρ q then E ρ p/q (Y ) T k(Y, Y )ρ p/q (Y ) = E A T p,1 A p,2 k(Y, Y ) with J(Y /X) the classical Fisher Information Distance between X and Y , see [START_REF] Johnson | Information theory and the central limit theorem[END_REF]. Taking n = d, α i = 1 and e i (y) = y i and supposing that all marginals have unit variance leads to the natural measure of discrepancy

S n (Y, X) = 2 i<j (σ X ij -σ Y ij ) 2
with σ X ij (resp., σ Y ij ) the covariance between the marginals i and j of X (resp., of Y ).

In terms of potential applications, one of the most interesting aspects of identity (114) is the fact that the right-hand side justifies the use of the left-hand side as a discrepancy metric. Applications of (116) have begun to be explored [START_REF] Chwialkowski | A kernel test of goodness of fit[END_REF][START_REF] Liu | A kernelized Stein discrepancy for goodness-of-fit tests[END_REF], and more general versions have been touched upon in [START_REF] Gorham | Measuring sample quality with kernels[END_REF]. The freedom of choice in the input matrices A p , A q encourages us to be hopeful that these quantities will have numerous applications.

Inspired by [START_REF] Chwialkowski | A kernel test of goodness of fit[END_REF][START_REF] Liu | A kernelized Stein discrepancy for goodness-of-fit tests[END_REF], we conclude the section with an illustration of how to apply Stein discrepancies to obtain a goodness-of-fit test for a Student-t distribution.

Example 7.9 (Dimension 2). Fix d = 2 and let p be the centered Student-t distribution with degrees of freedom and Σ = Id the identity matrix. We only consider the operator A 1 f (y) = τ τ τ (y)∇f (y) -yf (y) with τ τ τ the Stein kernel matrix given in (82), with entries (τ ij ) 1≤i,j≤2 ; the resulting kernelized discrepancy is S ν (t , q, ) = 1 n 1 n 2 n 1 i=1 n 2 j=1 u q (y i , y j ). Our naive implementation of the bivariate test appears to have difficulties in distinguishing the bivariate Student from the bivariate Gaussian (obtained at = 1000). Such an observation is perhaps not so surprising, see e.g. [START_REF] Mcassey | An empirical goodness-of-fit test for multivariate distributions[END_REF] where a similar problem is tested (by different means) with low power for the case of Gaussian vs Student, see page 1126. The problem of devising tractable powerful goodness-of-fit tests for multivariate distributions seems to be difficult; we will concentrate on this in future publications.

Definition 3 . 1 (

 31 Multivariate Stein class and operator). Let X be a d-dimensional random vector with pdf p : R d → R satisfying Assumption A. 1. The canonical directional Stein class for p in direction e is the vector field F 1,e (p) of all functions φ : R d → R such that the restriction of φ to Ω p is C 1 , p φ has integrable gradient and satisfies R d ∂ e (p φ) = 0. 2. The canonical Stein derivative for p in the direction e is the operator T e,p on F 1,e (p)

Definition 3 . 3 (

 33 Stein adjoint class).To every (scalar/vector/matrix valued function) F ∈ F(p) we associate dom(p, F) the collection of all matrix-valued functions G with compatible dimensions such that G is C 1 on Ω p , FG ∈ F(p) and F(∂ e G) ∈ L 1 (p) for all e ∈ S d-1 .

Remark 3 . 7 .

 37 One can easily extend Definition 3.6 to introduce a score (and score-Stein operator) in any direction e ∈ S d-1 , by considering the gradient ∇ = (∂ e , ∂ e ⊥ ) along e.

Definition 3 . 8 (

 38 Stein kernel). Consider a density p : R d → R + satisfying Assumption A, and suppose p has finite mean ν ∈ R d . For each canonical direction e i ∈ R d , i = 1, . . . , d, a Stein kernel for p in the direction e i is any vector field

with W = Lip( 1 )

 1 the collection of Lipschitz functions h : R d → R with Lipschitz constant smaller than 1. Here F (h) = h dF is the expectation of h under F . Abusing notation, we also write d W (X, Y ) = d W (F, G) when X ∼ F and Y ∼ G.

Proposition 4 . 11 .

 411 Let X ∈ R d have pdf p α (x) = 2ω d (x; Id)Φ(α T x),, and let Z be a d-dimensional standard normal element. Then

Remark 5 . 4 .

 54 and thus is exactly what we need. Setting t = x T Σ -1 x/2 the claim follows. In dimension d = 1, the Stein kernel (4) of p is the function τ p (x) = 1 p(x) +∞ x (u -ν)p(u)du. Changing variables in (70) for x ≥ 0,

Example 7 . 5 (

 75 Fisher information distance). Pick k(y, y ) = δ y=y the Dirac delta on the diagonal. Then (108) with the choice (116) becomesS(Y, X, δ) = E ρ X/Y (Y ) T ρ X/Y (Y ) =: J(Y /X)(117)

Example 7 . 6 (

 76 Independent kernels). Let (e i ) i=1,...,n be a sequence of functions in F(X) and k(x, y) = d i=1 α i e i ⊗ e i (which belongs to G for any (α i ) i=1,...,n ). ThenS(Y, X, (α, e) n ) = n i=1 α i (E [A p e i (Y )]) T (E [A p e i (Y )]).(118) Example 7.7. (Kernelized Stein discrepancies for comparing Gaussian random vectors) Let X and Y be independent centered multivariate normal random variables in R d with variances Σ 1 and Σ 2 , respectively. Take A p g(x) = Σ∇g(x) -xg(x). Then, for any sufficiently regular function e : R d → R we have E [A p e(Y )] = E [(Σ 1 -Σ 2 )∇e(Y )] so that the kernelized discrepancy (118) becomes S(Y, X, (α, e) n ) = n i=1 α i E [∇e i (Y )] T (Σ 1 -Σ 2 ) 2 E [∇e i (Y )] .

  with discrepancy generating function u q (y, y ) =A T p,1 A p,2 k(y, y ) =[τ 11 (y)τ 11 (y ) + τ 12 (y)τ 12 (y )]∂y 1 ∂ y 1 k(y, y ) + [τ 11 (y)τ 12 (y ) + τ 12 (y)τ 22 (y )]∂ y 1 ∂ y 2 k(y, y ) + [τ 12 (y)τ 11 (y ) + τ 22 (y)τ 12 (y )]∂ y 2 ∂ y 1 k(y, y ) + [τ 12 (y)τ 12 (y ) + τ 22 (y)τ 22 (y )]∂ y 2 ∂ y 2 k(y, y ) -(τ 11 (y)y 1 + τ 12 (y)y 2 )∂ y 1 k(y, y ) -(τ 12 (y)y 1 + τ 22 (y)y 2 )∂ y 2 k(y, y ) -(y 1 τ 11 (y ) + y 2 τ 12 (y ))∂ y 1 k(y, y ) -(y 1 τ 12 (y ) + y 2 τ 22 (y ))∂ y 2 k(y, y ) + (y 1 y 1 + y 2 y 2 )k(y, y )As in the previous example, we present simulation results on a rather modest simulation study. We compare X ∼ p a bivariate (centered scaled) Student with = 5 degrees of freedom with Y ∼ q a bivariate (centered scaled) Student with degrees of freedom, via the kernelized discrepancies based on the RBF kernel k(x, y) = e -(x-y) 2 /2 with n 1 = n 2 = 100. The quantiles were estimated by simulation, with J = 10 3 experiments; we obtained 2.5% 97.5% -0.07256331 0.08441458 (which indicates some asymmetry in the sample distribution). The results for 10 3 simulations run for each value of degrees of freedom ∈ {0.1, 1, 5, 10, 100, 1000} (with = 5 corresponding to the null hypothesis) are reported below: 339 0.690 0.056 0.043 0.046 0.038 0.048 0.052

  . Many authors call this score-Stein operator the Stein operator for multivariate p; see e.g. [57, definition 2.1], [26, Section 2.1], [66, Section 2.2.2].

  Equation[START_REF] Ernst | First order covariance inequalities via Stein's method[END_REF] is often called the Stein operator for multivariate p; see e.g. [57, definition 2.1], [26, Section 2.1], [66, Section 2.2.2]. This operator is particularly interesting when constant functions 1 ∈ F 1 (p), an assumption which holds if p is a differentiable density such that

)

  Equation[START_REF] Goldstein | Zero biasing in one and higher dimensions, and applications[END_REF] is reminiscent of the classical Stein identity (1); one can easily check that τ τ τ p (x) = Σ is a Stein kernel for the multivariate standard Gaussian distribution.

2. The question of whether a density p admits a Stein kernel is non trivial, if one wants the Stein identity to hold for a larger class of functions than the regular ones with compact support in Ω p (see Section 6). Definition (3.8) agrees with the one of [27], except for the class of test functions we impose. 3. The Stein kernel τ p,i ∈ F(p) in the direction e i in Definition (3.8) could be defined equivalently by requiring

  Definition 7.1 (Integral Probability Metrics). Let F(R d ) be a collection of cumulative distribution functions on R d and denote L 1 (F(R d )) the class of Borel measurable functions h : R d → R (for some q ≥ 1) such thatF (|h|) = |h|dF < ∞ for all F ∈ F(R d ). A metric on F(R d) is an integral probability metric (IPM) if it can be written in the form for some class of real-valued bounded measurable test functionsH ⊂ L 1 (F(R d )) (|•| is the Euclidean norm).The expression on the right-hand side of (102) is called an IPMdiscrepancy.Many important probability metrics can be represented as integral probability metrics; classical references are[START_REF] Zolotarev | Probability metrics[END_REF][START_REF] Gibbs | On choosing and bounding probability metrics[END_REF]. The Wasserstein distance between X and Y , which we have already used in this paper, takes H = W the collection of Lipschitz functions h : R

	d H (F, G) := sup	|F (h) -G(h)|	(102)
	h∈H

d → R with Lipschitz constant smaller than 1. The Kolmogorov distance between two random vectors X ∼ F and

Y ∼ G is Kol(X, Y ) = d H Kol (F, G) with H Kol = I (-∞,z] , z ∈ R d

. The Wasserstein distance between X and Y takes H = W the collection of Lipschitz functions h : R d → R with Lipschitz constant smaller than 1. The total variation distance takes H TV the collection of Borel measurable functions h : R d → [0, 1]. For other examples and references see [62, Appendix E]. In Definition 7.1 the Euclidean norm | • | is used, but the definition generalises easily to other norms • as long as

Acknowledgements. GM and YS gratefully acknowledges support by the Fonds de la Recherche Scientifique -FNRS under Grant MIS F.4539. 16. GR acknowledges partial support from EPSRC grant EP/K032402/1 and the Alan Turing Institute. We also thank Christophe Ley and Guillaume Poly for interesting discussions, as well as Lester Mackey and Steven Vanduffel for suggesting some references which we had overlooked.

Information metrics and kernelized Stein discrepancies

In principle, the Stein discrepancy (7.2) can be used as a basis for goodness-of-fit tests, and could be estimated numerically. In high-dimensional problems, the class G is often too large to allow numerical integration. In high-dimensional goodness of fit tests, restricting the class of functions to a ball in a reproducing kernel Hilbert space associated with a positive definite kernel k(x, x ) has been shown in [START_REF] Chwialkowski | A kernel test of goodness of fit[END_REF][START_REF] Liu | A kernelized Stein discrepancy for goodness-of-fit tests[END_REF] to be an efficient way of estimating Stein discrepancies. In this context these discrepancies are called kernelized Stein discrepancies, with the kernel k(x, x ) in mind. The framework of the present paper shows how to generalise their approach, as follows.

Let X ∼ p and Y ∼ q be two random variables on R d with differentiable densities and respective Stein classes F(p) and F(q). Suppose, for simplicity, that both share the same mean and the same support S, satisfying Assumption A. Fix d × d matrix valued functions A p ∈ F(p) and A q ∈ F(q), set a p = T div,p (A p ), a q = T div,q (A q ) and introduce the divergence based vector valued standardizations

A q g(x) = A q (x)∇g(x) + a q (x)g(x), g :

as in [START_REF] Fathi | Stein kernels and moment maps[END_REF].

The Stein heuristic that if p and q are close, then E q⊗q [A p⊗p k(Y, Y )] should be small still holds, where A p⊗p is the concatenated operator A p operating on functions g : R 2d → R by treating the first d and the last d components independently. It turns out that iterating the operator is a more elegant way to obtain a kernelized discrepancy. Write A T p for the transpose of the operator A p in the matrix transpose sense. Then for any positive definite symmetric kernel k with marginals in F(p), using (105),

By conditioning on X it is easy to see that

The operator (107) has been used in [START_REF] Chwialkowski | A kernel test of goodness of fit[END_REF] for the particular choice A p (x) = I d for which a p = ρ p , the score operator. In this case, evaluating (107) does not require knowledge of the normalising constant for the density p and is hence particularly attractive for applications in Bayesian inference. Equation (107) motivates our general definition of kernelized Stein discrepancies. We use the convention that L i , i = 1, 2 denotes the operator L applied with respect to the ith variable of the function k(•, •).

Definition 7.3 (Kernelized Stein discrepancies).

Let A p (resp., A p ) be a Stein operator for p (resp., for q) with class F(p) (resp., F(q)). Let k be some kernel k :

When X ∼ p and Y ∼ q, in abuse of notation we also write S(X, Y, k) = S(p, q, k).

Example 7.8 (Dimension 1). Let p be the centered Student-t distribution with degrees of freedom, with score function given by

The operator obtained from (119) is

The preceding developments lead to postulating the sample-based discrepancy

with y 1 , . . . , y n 1 and y 1 , . . . , y n 2 two i.i.d. samples independently drawn from q and

and k(y, y ) a well chosen kernel. Particularizing (121) we get

Let Y, Y ∼ q be two independent copies. Since under natural conditions, E q [u ρ q (Y, Y )] = 0 if and only if q = p, a natural goodness-of-fit test in this context is to reject the null assumption H 0 : q = p whenever S ν (t , q, k) is too large.

For the sake of proof of concept rather than anything else, here are the result of simulations comparing X ∼ p a Student with = 5 degrees of freedom with Y ∼ q a Student with degrees of freedom, via the kernelized discrepancies based on the RBF kernel k(x, y) = e -(x-y) 2 /2 . The quantiles for S ρ (t , q, k) under the null hypothesis were estimated by simulation, with J = 10