Stein Points
Wilson Ye Chen Lester Mackey Jackson Gorham François-Xavier Briol Chris J. Oates


[
Stein Points

















Wilson Ye Chenuts
Lester Mackeymic
Jackson Gorhamod
François-Xavier Briolwar,imp,ati
Chris. J. Oatesncl,ati


utsSchool of Mathematical and Physical Sciences, University of Technology Sydney, Australia
micMicrosoft Research New England, USA
odOpendoor Labs, Inc., USA
warDepartment of Statistics, University of Warwick, UK
impDepartment of Mathematics, Imperial College London, UK
nclSchool of Mathematics, Statistics and Physics, Newcastle University, UK
atiAlan Turing Institute, UK

Wilson Ye Chenye.chen@uts.edu.au
Lester Mackeylmackey@microsoft.com




Discrepancy, Kernel, Stein's Method

0.3in
]









  



An important task in computational statistics and machine learning is to approximate a posterior distribution p(x) with an empirical measure supported on a set of representative points {x_i}_i=1^n.
This paper focuses on methods where the selection of points is essentially deterministic, with an emphasis on achieving accurate approximation when n is small.
To this end, we present Stein Points.
The idea is to exploit either a greedy or a conditional gradient method to iteratively minimise a kernel Stein discrepancy between the empirical measure and p(x).
Our empirical results demonstrate that Stein Points enable accurate approximation of the posterior at modest computational cost.
In addition, theoretical results are provided to establish convergence of the method.





§ INTRODUCTION


This paper is motivated by approximation of a Borel distribution P, defined on a topological space X, with deterministic point sets or sequences {x_i}_i=1^n ⊂ X for n ∈ℕ, such that

    1/n∑_i=1^n h(x_i)    →   ∫ h  dP

as n →∞ for all functions h : X →ℝ in a specified set ℋ.
Throughout it will be assumed that P admits a density p, with respect to a reference measure, available in a form that is un-normalised (i.e., we know q(x) in closed form where p(x) = q(x)/C for some C>0).
Such problems occur in Bayesian statistics where P represents a posterior distribution, and the integral represents a posterior expectation of interest.
Markov chain Monte Carlo (MCMC) methods are extensively used for this task but suffer (in terms of accuracy) from `clustering' of the points {x_i}_i=1^n when n is small.
This observation motivates us to instead consider a range of goal-oriented discrete approximation methods that are designed with un-normalised densities in mind.
   
The problem of discrete approximation of a distribution, given its normalised density, has been considered in detail and relevant methods include quasi-Monte Carlo (QMC) <cit.>, kernel herding <cit.>, support points <cit.>, transport maps <cit.>, and minimum energy methods <cit.>.
On the other hand, the question of how to proceed with un-normalised densities has been primarily answered with increasingly sophisticated MCMC.


At the same time, recent work had led to theoretically-justified measures of sample quality in the case of an un-normalised target.
In <cit.> it was shown that Stein's method can be used to construct discrepancy measures that control weak convergence of an empirical measure to a target.
This was later extended in <cit.> to encompass a family of discrepancy measures indexed by a reproducing kernel.
In the latter case, the discrepancy measure can be recognised as a maximum mean discrepancy <cit.>.
As such, one can consider discrete approximation as an optimisation problem in a Hilbert space and attempt to optimise this objective with either a greedy or a conditional gradient method.
The resulting method – Stein Points – and its variants are proposed and studied in this work.




  
Our Contribution

This paper makes the following contributions:


  * Two algorithms are proposed for minimisation of the kernel Stein discrepancy <cit.>; a greedy algorithm and a conditional gradient method. 
In each case, a convergence result of the form in Eqn. <ref> is established.

  * Novel kernels are proposed for the KSD, and we prove that, with these kernels, the KSD controls weak convergence of the empirical measure to the target.
In other words, the test functions h for which our results hold constitute a rich set ℋ.




  
Outline

The paper proceeds as follows.
In Section <ref> we provide background, and in Section <ref> we present the approximation methods that will be studied.
Section <ref> applies these methods to both simulated and real approximation problems and provides a extensive empirical comparison.
All technical material is contained in Section <ref>, where we derive novel theoretical results for the methods we proposed.
Finally we summarise our findings in Section <ref>.



§ BACKGROUND
 

Throughout this section it will be assumed that X is a metric space, and we let 𝒫(X) denote the collection of Borel distributions on X.
In this context, weak convergence of the empirical measure to P corresponds to taking the set ℋ in Eqn. <ref> to be the set ℋ_CB of functions which are continuous and bounded.
In this work we also consider sets ℋ that correspond to stronger modes of convergence in 𝒫(X). 

First, in <ref>, we recall how discrepancy measures are constructed.
Then we recall the use of Stein's method in this context in <ref>.
Formulae for KSD are presented in <ref>.








 §.§ Discrepancy Measures
 

A discrepancy is a quantification of how well the points {x_i}_i=1^n cover the domain X with respect to the distribution P.
This framework will be developed below in reproducing kernel Hilbert spaces <cit.>, but the general theory of discrepancy can be found in <cit.>.
Note that we focus on unweighted point sets for ease of presentation, but our discussions and results generalise straightforwardly to point sets that are weighted.

Let k : X × X →ℝ be the reproducing kernel of a RKHS 𝒦 of functions 𝒳→ℝ.
That is, 𝒦 is a Hilbert space of functions with inner product ⟨·,·⟩_𝒦 and induced norm ·_𝒦 such that, for all x ∈ X, k(x,·) ∈𝒦 and f(x) = ⟨ f , k(x,·) ⟩_𝒦 whenever f ∈𝒦. The Cauchy-Schwarz inequality in 𝒦 gives that

    | 1/n∑_i=1^n f(x_i) - ∫ f dP | 
       ≤   f_𝒦  D_𝒦,P({x_i}_i=1^n)

where the final term

    D_𝒦,P({x_i}_i=1^n )   :=  1/n∑_i=1^n k(x_i,·) - ∫ k(x,·) dP(x) _𝒦

is the canonical discrepancy measure for the RKHS. 
The Bochner integral k_P := ∫ k(x,·) dP(x) ∈𝒦 is known as the mean embedding of P into 𝒦 <cit.>.
Thus, if ℋ = B(𝒦) := {f ∈𝒦 : f_𝒦≤ 1} is the unit ball in 𝒦, then D_𝒦,P({x_i}_i=1^n) → 0 implies the convergence result in Eqn. <ref>.

The RKHS framework is now standard for QMC analysis <cit.>.
Its popularity derives from the fact that, when both k_P and k_P,P :=
∫ k_P dP are explicit, the canonical discrepancy measure is
also explicit:

    D_𝒦,P({x_i}_i=1^n)     =    
          √( k_P,P - 2/n∑_i=1^n k_P(x_i) + 1/n^2∑_i,j=1^n k(x_i,x_j) )

Table 1 in <cit.> collates pairs (k,P) for which k_P and k_P,P are explicit.

If P is a posterior distribution, so that p has unknown normalisation
constant, it is unclear how the terms k_P and k_P,P can be computed in closed form,
and so similarly for the discrepancy D_𝒦,P. 
This has so far prevented QMC and related methods such as kernel herding <cit.> from being used to compute posterior integrals. 
A solution to this problem can be found in Stein's method, presented next.




 §.§ Kernel Stein Discrepancy
 

The method of <cit.> was introduced as an analytical tool for establishing convergence in distribution of random variables, but its potential for generating and analyzing computable discrepancies was developed in <cit.>.
In what follows, we recall the kernelised version of the Stein discrepancy, first presented for an optimally-weighted point set in 2.3.3 of <cit.> and later generalised to an arbitrarily-weighted point set in <cit.>.

Suppose that X carries the structure of a smooth manifold, and consider a linear differential operator 𝒯_P on X, together with a set ℱ of sufficiently differentiable functions, with the following property:

    ∫𝒯_P[f]  dP     =     0   ∀ f ∈ℱ.

Then 𝒯_P is called a Stein operator and ℱ a Stein set. In the
kernelised version of Stein's method, the set ℱ is either an
RKHS 𝒦 with reproducing kernel k : X × X →ℝ, or the product 𝒦^d, which contains vector-valued
functions f = (f_1,…,f_d) with f_j ∈𝒦 and is equipped
with a norm[For what follows, any vector norm can be used to combine the component norms f_j_𝒦 <cit.>.] f_𝒦^d = (∑_j=1^d f_j_𝒦^2
)^1/2.
For the case ℱ = 𝒦, the image of 𝒦 under a Stein operator 𝒯_P is denoted 𝒦_0 = 𝒯_P 𝒦.
The notation can be justified since, under appropriate regularity assumptions, the set 𝒯_P 𝒦 admits structure from the reproducing kernel k_0(x,x') = 𝒯_P 𝒯_P k(x,x') <cit.>.
Here 𝒯_P is the adjoint of the operator 𝒯_P and acts on the second argument x' of the kernel.
If instead ℱ = 𝒦^d, then we suppose that 𝒯_P f = ∑_j=1^d 𝒯_P,j f_j so that the set 𝒦_0 = 𝒯_P 𝒦^d admits structure from the reproducing kernel k_0(x,x') = ∑_j=1^d 𝒯_P,j𝒯_P,j k(x,x').
In either case, we will call the reproducing kernel k_0 of 𝒦_0 a Stein reproducing kernel.

Stein reproducing kernels possess the useful property that k_0,P = ∫ k_0(x,·) dP = 0 and k_0,P,P = ∫ k_0,PdP = 0, so in particular both are explicit.
Thus, if k_0 is a Stein reproducing kernel, then Eqn. <ref> can be simplified:

    D_𝒦_0,P({x_i}_i=1^n )     =    √(1/n^2∑_i,j = 1^n k_0(x_i,x_j)) .

We call this quantity a kernel Stein discrepancy (KSD).
Next, we exhibit some differential operators for which Eqn. <ref> is satisfied and Eqn. <ref> can be computed.




 §.§ Stein Operators and Their Reproducing Kernels
 

The divergence theorem can be used to construct Stein operators on a manifold.
For P supported on X = ℝ^d, <cit.> considered the Langevin Stein operator

    𝒯_P f    :=    ∇· (p f) /p

where ∇· is the usual divergence operator and f ∈𝒦^d.
Thus, for the Langevin Stein operator, we obtain a Stein reproducing kernel

    k_0(x,x')     =    ∇_x ·∇_x' k(x,x') 
           + ∇_x k(x,x') ·∇_x'log p(x') 
           + ∇_x' k(x,x') ·∇_x log p(x) 
           + k(x,x') ∇_x log p(x) ·∇_x'log p(x').

To evaluate this kernel, the normalisation constant for p is not required.
Other Stein operators for the Euclidean case were developed in <cit.>.
For P supported on a closed Riemannian manifold X, <cit.> proposed the second order Stein operator 𝒯_P f := 1/p∇· (p∇ f) where ∇ and ∇· are, respectively, the gradient and divergence operators on the manifold and f ∈𝒦.
Other Stein operators for the general case are proposed in the supplement of <cit.>.

The theoretical results in <cit.> established that certain combinations of Stein operator 𝒯_P and base kernel k ensure that KSD controls weak convergence; that is, D_𝒦_0 , P({x_i}_i=1^n) → 0 implies that Eqn. <ref> holds with ℋ = ℋ_CB.
This important result motivates our next contribution, where numerical optimisation methods are used to select points {x_i}_i=1^n to approximately minimise KSD.
Theoretical analysis of the proposed methods is reserved for Section <ref>.







§ METHODS
 

In this paper, two algorithms to select points {x_i}_i=1^n are studied in detail.
The first of these is a greedy algorithm, which at each iteration attempts to minimise the KSD, whilst the second is a conditional gradient algorithm, known as herding, which also targets the KSD.
In <ref> and <ref> the two algorithms are described, whilst in <ref> some alternative approaches are briefly discussed.



 §.§ Greedy Algorithm
 

The simplest algorithm that we consider follows a greedy strategy, whereby
the first point x_1 is taken to be a global maximum of p (an operation
which does not require the normalisation constant) and each subsequent point
x_n is taken to be a global maximum of D_𝒦_0 ,
  P({x_i}_i=1^n), with the KSD being viewed as a function of x_n holding {x_i}_i=1^n-1 fixed.
Equivalently, at iteration n > 1 of the greedy algorithm, we select

    x_n    ∈   _x ∈ X  k_0(x,x)/2 + ∑_i=1^n-1 k_0(x_i,x) .

Note that each iteration of the algorithm requires the solution of a global optimisation problem over X; in practice we employed a numerical optimisation method, and this choice is discussed in detail in connection with the empirical results in Section <ref> and the theoretical results in Section <ref>.

If a user has a budget of at most n points, the greedy algorithm can be run for n iterations and thereafter improved using (block) coordinate descent on the KSD objective to update an existing point x_i instead of introducing a new point. 
The cost of each update is equal to the cost of adding the n-th greedy Stein Point.
This budget-constrained variant of the method will be called Stein Greedy-n in the sequel (see Section <ref> for more details).




 §.§ Herding Algorithm
 

The definition of discrepancy in Section <ref> suggests that selection of {x_i}_i=1^n can be elegantly formulated as a single global optimisation problem over 𝒦_0.
Let M(𝒦_0) be the marginal polytope of 𝒦_0; i.e. the convex hull of the set {k_0(x,·)}_x ∈ X <cit.>.
The mean embedding Q ↦ k_Q, as a map 𝒫(X) → M(𝒦), is injective whenever the kernel k is universal and X is compact <cit.>, so that in this case k_Q fully characterises Q. 
Results in a similar direction for Stein reproducing kernels were established in <cit.> and <cit.>.
Thus, as P is mapped to 0 under the embedding, we are motivated to consider non-trivial solutions to

    _f ∈ M(𝒦_0) J(f),      J(f)   :=  1/2 f _𝒦_0^2  .

As might be expected, the objective function is closely related to KSD; for f(·) = 1/n∑_i=1^n k_0(x_i,·) we have J(f) = 1/2 D_𝒦_0,P({x_i}_i=1^n)^2.
An iterative algorithm, called kernel herding, was proposed in <cit.> to solve problems in the form of Eqn. <ref>.
This was later shown to be equivalent to a conditional gradient algorithm, the Frank-Wolfe algorithm, in <cit.>.
The canonical Frank-Wolfe algorithm, which results in an unweighted point set <cit.>, is presented next.

The first point x_1 is again taken to be a global maximum of p; this corresponds to an element f_1 = k_0(x_1,·) ∈ M(𝒦_0).
Then, at iteration n > 1, the convex combination f_n = n-1/n f_n-1 + 1/n f_n^* ∈ M(𝒦_0) is constructed where the element f_n^* encodes a direction of steepest descent:

    f_n    ∈   _f ∈ M(𝒦_0)⟨ f , DJ(f_n-1) ⟩_𝒦_0 ,

where DJ(f) is the representer of the Fréchet derivative of J at f.
Given that minimisation of a linear objective over a convex set can be restricted to the boundary of that set, it follows that f_n^* = k(x_n,·) for some x_n ∈ X.
Thus, at iteration n > 1 of the proposed algorithm, we select

    x_n    ∈   _x ∈ X  ∑_i=1^n-1 k_0(x_i,x)

to obtain f_n(·) = 1/n∑_i=1^n k_0(x_i,·), the embedding of the empirical distribution of {x_i}_i=1^n.
As in the standard kernel herding algorithm of <cit.>, each iteration in practice requires the solution of a global optimisation problem over X.

Compared to Eqn. <ref>, the greedy algorithm is seen to be a regularised version of herding with regulariser 1/2 k_0(x,x).
The two algorithms coincide if k_0(x,x) is independent of x; however, this is typically not true for a Stein reproducing kernel.
The computational cost of either method is O(n^2); thus we anticipate applications in which evaluation of p(x) (and its gradient) constitute the principal computational bottleneck. 
The performance of both algorithms is studied empirically in Section <ref> and theoretically in Section <ref>.
In a similar manner to Stein Greedy-n, a budget-constrained variant of the above method can be considered, which we call Stein Herding-n in the sequel.



 §.§ Other Algorithms
 

The output of either of our algorithms will be called Stein Points.
These are extensible point sequence S_n = (x_i)_i=1^n, meaning that S_n can be incrementally extended S_n = (S_n-1,x_n) as required.
Another recently proposed extensible method is the (sequential) minimum energy design (MED) of <cit.>, here used as a benchmark.

For some problems the number of points n will be fixed in advance and the aim will instead be to select a single optimal point set {x_i}_i=1^n.
This alternative problem demands different methodologies, and a promising method in this direction is Stein variational gradient descent <cit.>.
A natural point set analogue of our approach would be to optimise KSD for n fixed.
This approach was considered for other discrepancy measures in <cit.>, where the Newton method was used.
We instead employ our budget-constrained algorithms Stein Greedy-n and Stein Herding-n for this use case.














§ RESULTS
 

In this section, the proposed greedy and herding algorithms are empirically assessed and compared.
In <ref> a Gaussian mixture problem is studied in detail, whilst in <ref> and <ref>, respectively, the methods are applied to approximate the parameter posterior in a non-parametric regression model and an IGARCH model.
First, in <ref> we provide details on the experimental protocol.




 §.§ Experimental Protocol
 

Here we describe the parameters and settings that were varied in the experiments that are presented.



  
Stein Operator

To limit scope, we focus on the case X = ℝ^d and always take 𝒯_P to be the Langevin Stein operator in Eqn. <ref>.





  
Choice of Kernel

For the kernel k in Eqn. <ref> we considered one standard choice – the inverse multi-quadratic (IMQ) kernel – together with two novel alternatives:


  *  (IMQ) k_1(x,x')   =  (α + x - x'_2^2 )^β

  * (inverse log) k_2(x,x')   =  (α + log(1 + x - x'_2^2))^-1

  * (IMQ score) 
 k_3(x,x')   =  (α + ∇log p(x) - ∇log p(x')_2^2)^β.

In all cases α > 0 and β∈ (-1,0).
To limit scope, in what follows we considered a finite number of judiciously selected configurations for α,β, though in principle these could be optimised as in <cit.>.
The best set of parameter values was selected for each algorithm and each target distribution, where the possible values were α∈{0.1η, 0.5η, η, 2η, 4η, 8η} and β∈{0.1, 0.3, 0.5, 0.7, 0.9}, with η > 0 problem-dependent (see the Supplement).
The IMQ kernel, together with the Langevin Stein operator, was proven in <cit.> to provide a KSD that controls weak convergence.
Similar results for novel kernels k_2 and k_3 are established in Section <ref>.



  
Numerical Optimisation Method
Any optimisation procedure could be used to (approximately) solve the global optimisation problem embedded in each iteration of the proposed algorithms.
In our experiments, we considered the following numerical methods, for which full details appear in the Section <ref>.


  * Nelder-Mead (): At iteration n, parallel runs of Nelder-Mead were employed, initialised at draws from a Gaussian mixture proposal centred on the current point set Π = 1/n-1∑_i=1^n-1𝒩(x_i,λ I) with problem-specific λ > 0.

  * Monte Carlo (): The optimisation problem at iteration n was solved over a sample of points drawn from the same Gaussian mixture proposal Π.

  * Grid search (): Through brute force, the optimisation problem at iteration n was solved over a regular grid of width 1/√(n). This required O(n^-d/2) points; if required, the domain was first truncated with a large bounding box.





  
Performance Assessment
To obtain a reasonably objective assessment, we focused on the 1-Wasserstein distance between the empirical measure and P:

    W_P({x_i}_i=1^n)     =    sup_h ∈ℋ_Lip| 1/n∑_i=1^n h(x_i) - ∫ h dP | ,

where ℋ_Lip is the set of all function h : X →ℝ with Lipschitz constant Lip(h) ≤ 1.
By replacing P with the empirical measure P_N = 1/N∑_i=1^N
δ_y_i for y_i iid∼ P, the expected error
from using W_P_N({x_i}_i=1^n) in lieu of W_P({x_i}_i=1^n)
converges at a N^-1/2log N rate for d=2 and N^-1/d rate for
d>2 <cit.>.
By employing L_1-spanners, the approximation W_P_N({x_i}_i=1^n) can be
computed in O((n+N)^2 log^2d-1(n+N)) time <cit.>.

For all reported results, the {y_i}_i=1^N were obtained by brute-force Monte Carlo methods applied to P, with N sufficiently large that approximation error can be neglected.

The computational cost associated to any given method was quantified as the total number n_eval of times either the log-density log p or its gradient ∇log p were evaluated.
This can be justified since in most applications the `parameter to data' map dominates the computational cost associated with the likelihood.





  
Benchmarks

Two existing methods were used as a benchmark:


  * The MED method of <cit.> relies on numerical optimisation methods to minimise an energy measure ℰ_δ,P({x_i}_i=1^n), adapted to P.
This measure has one tuning parameter δ∈ [1,∞).
See Section <ref> of the Supplement for full detail.

  * The SVGD method of <cit.> performs a version of gradient descent on the Kullback-Leibler divergence, described in Section <ref> of the Supplement.


To avoid confounding of the empirical results by incomparable algorithm parameters, (1) the collection of numerical optimisation methods used for KSD were also used for MED, and (2) the same collection of kernels k_1, …, k_3 was considered for SVGD as was used for KSD.
Note that, apart from standard Monte Carlo, none of the methods considered in these experiments are re-parametrisation invariant.







 §.§ Gaussian Mixture Test
 

For our first test, we considered a Gaussian mixture model

    P     =    1/2𝒩(μ_1 , Σ_1) + 1/2𝒩(μ_2 , Σ_2)

defined on X = ℝ^2.
Full settings for each of the methods considered are detailed in Section <ref> in the Supplement.
Typical point sets are displayed  over the contours of P for μ_1 = (-1.5,0), μ_2 = (1.5,0), Σ_1 = Σ_2 = I in Figure <ref>.
Additionally, point sets for the n point budget-constrained algorithms Stein Greedy-n and Stein Herding-n are presented in Figure <ref> in the Supplement.
For each of the methods shown in Figures <ref> and <ref>, tuning parameters were varied and the overall performance was captured in Figure <ref>.
It was observed that for (a-c) the choice of numerical optimisation method was the most influential tuning parameter, with the simpler Monte Carlo-based method being most successful.
The kernels k_1,k_2 were seen to perform well, but in (a,b,d) the kernel k_3 was sometimes seen to fail.

A subjectively-selected exemplar was extracted for each method, and these `best' results for each method are overlaid in Figure <ref>.
The total number of points was limited to n = 100.
In terms of our proposed methods, two qualitative regimes were observed:
(i) For low computational budget log n_eval≤ 7, the standard Monte Carlo method performed best.
(ii) For a larger computational budget 7 < log n_eval, greedy Stein points were not out-performed.

Note that KSD and SVGD are based on the log target and its gradient, whilst for MED the target p(x) itself is required.
As a result, numerical instabilities were sometimes encountered with MED.

Next, we turned our attention to two important posterior approximation problems that occur in the real world.
















 §.§ Gaussian Process Regression Model
 

The Gaussian process (GP) model is a popular choice for uncertainty quantification in the non-parametric regression context <cit.>.
The data 𝒟 = {(x_i,y_i)}_i=1^n that we considered are from a light detection and ranging (LIDAR) experiment <cit.>. They consist of 221 realisations of an independent scalar variable x_i (distances travelled before the light is reflected back to its source) and a dependent scalar variable y_i (log-ratios of received light from two laser sources); these were modelled as
y_i  =  g(x_i) + ϵ_i,  for ϵ_i i.i.d.∼𝒩(0,σ^2)
and a known value of σ.
The unknown regression function g is modelled as a centred GP with covariance function cov(x,x') = θ_1 exp( - θ_2 (x - x')^2).
The hyper-parameters θ_1,θ_2 > 0 determine the suitability of the GP model, but appropriate values will be unknown in general.
In this experiment we re-parametrised ϕ_i = logθ_i and placed a standard multivariate Cauchy prior on ϕ = (ϕ_1, ϕ_2), defined on X = ℝ^2.
The task is thus to approximate the conditional distribution p(ϕ | 𝒟).
This problem is motivated by the computation of posterior predictive marginal distributions p(y^* | x^* , 𝒟) for a new input x^*, which is defined as the integral ∫ p(y^* | x^* , ϕ, 𝒟) p (ϕ | 𝒟) dϕ.
Note that the density p(ϕ | 𝒟) can be differentiated, and an explicit formula is provided in <cit.>.

For each class of method, `best' tuning parameters were selected and these are presented on the same plot in Figure <ref>.
In addition, typical point sets provided by each method are presented in Figures <ref> and <ref> in the Supplement.
MED was not included because the method exhibited severe numerical instability on this task, as earlier discussed.
Results indicated three qualitative regimes where, respectively, Monte Carlo, greedy Stein points and SVGD provided the best performance for fixed cost.













 §.§ IGARCH Model
 

The integrated generalised autoregressive conditional heteroskedasticity (IGARCH) model is widely-used to describe financial time series (y_t) with time-varying volatility (σ_t) <cit.>.
The model is as follows:

    y_t     =    σ_t ϵ_t, ϵ_t i.i.d.∼𝒩(0,1) 
    σ_t^2     =    θ_1 + θ_2 y_t-1^2 + (1-θ_2) σ_t-1^2

with parameters θ = (θ_1,θ_2), θ_1 > 0 and 0 < θ_2 < 1.
The data y = (y_t) that we considered were 2,000 daily percentage returns of the S&P 500 stock index (from December 6, 2005 to November 14, 2013), and an improper uniform prior was placed on θ.
Thus the task was to approximate the posterior p(θ | y).
Note that, whilst the domain X = ℝ_+ × (0,1) is bounded, for these data the posterior density is negligible on the boundary ∂ X.
This ensures that Eqn. <ref> holds essentially to machine precision; see also the discussion in <cit.>.
For the IGARCH model, gradients ∇log p(θ | y) can be obtained as the solution of a recursive system of equations for ∂σ_t / ∂θ_2.

As before, the `best' performing of each class of method was selected and these are presented on the same plot in Figure <ref>.
In addition, typical point sets provided by each method are presented in Figures <ref> and <ref> in the Supplement.
(Numerical instability again prevented results for MED from being obtained.)
Results were consistent with the Gaussian mixture experiment, favouring either Monte Carlo or greedy Stein points depending on the computational budget.











#1#1 
#1#1
#1#1
#1#1
#1#1
#1#1
#1#1

#1I[#1] 
E 
#1[#1]
#1#2_#1[#2]
#1𝒪(#1) 
P 
#1(#1)
#1#2_#1[#2]

#1#1 
 
 
 
 
 
#1#1 
 
 
 

 
≜ 

 
 
 
 
 

#1supp(#1)







Grönwall
Hölder
Itô
Nyström
Schätten
Matérn




ℝ 
ℝ
ℤ 
ℤ
ℚ 
ℚ
ℕ 
ℕ
ℂ 

lemma[theorem]Lemma

Δ   
Re
Im




























































#1#2KL(#1‖#2)

#1#2#1#24mu#1#2










§ THEORETICAL RESULTS
 
In this section we establish two important forms of theoretical guarantees: 
(1) discrepancy control, i.e., D_𝒦_0,P({x_i}_i=1^n) → 0 as n→∞ for our extensible Stein Point sequences
and
(2) distributional convergence control, i.e., for our kernel choices and appropriate choices of target, D_𝒦_0,P({x_i}_i=1^n) → 0 
implies that the empirical distribution 1/n∑_i=1^nδ_x_i converges in distribution to P.




 §.§ Discrepancy Control

Earlier work has shown that, when a kernel is uniformly bounded (i.e., sup_x∈ X k_0(x,x) ≤ R^2), 
the greedy and kernel herding algorithms decrease the associated discrepancy
D_𝒦_0,P at an O(n^-1/2) rate <cit.>.

We extend these results to cover all growing, P-sub-exponential kernels. 


We say a reproducing kernel k_0 is P-sub-exponential if

    Z∼ P k_0(Z,Z) ≥ t ≤ c_1 e^-c_2t

for some constants c_1, c_2 > 0 and all t ≥ 0.

Notably, any uniformly bounded reproducing kernel is P-sub-exponential, and, when P is a sub-Gaussian distribution, 
any kernel with at most quadratic growth (i.e., k_0(x,x) = O(x^2)) is also P-sub-exponential.
Our first result, proved in herding-convergence-proof, shows that if we truncate the search domain suitably in each step, Stein Herding decreases the discrepancy at an O(√(log(n)/n)) rate.


This result holds even if each point x_i is selected suboptimally with error δ/2.
This extra degree of freedom allows a user to conduct a grid search or a search over appropriately generated random points on each step <cit.> and still obtain a rate of convergence.



Suppose k_0 with k_0,P = 0 is a P-sub-exponential reproducing kernel.
Then there exist constants c_1, c_2 > 0 depending only on k_0 and P 
such that any point sequence {x_i}_i=1^n satisfying 

    ∑_i=1^j-1 k_0(x_i, x_j) ≤δ/2 
    	+ min_x∈ X: k_0(x,x)≤ R_j^2∑_i=1^j-1 k_0(x_i, x)

with k_0(x_j,x_j)≤ R_j^2∈[2log(j)/c_2,2log(n)/c_2] for each 1 ≤ j ≤ n also satisfies 

    D_K_0,P({x_i}_i=1^n) 
    	≤ e^π/2√(2log(n)/c_2n +c_1/n + δ/n).


Our next result, proved in greedy-convergence-proof, shows that Stein Greedy decreases the discrepancy at an O(√(log(n)/n)) rate whether we choose to truncate (R_j < ∞) or not (R_j = ∞).  This highlights an advantage of the Stein Greedy algorithm over Stein Herding: the extra k_0(x,x)/2 term acts as a regularizer ensuring that no truncation is necessary. The result also accommodates points x_i selected suboptimally with error δ/2.



Suppose k_0 with k_0,P = 0 is a P-sub-exponential reproducing kernel.
Then there exist constants c_1, c_2 > 0 depending only on k_0 and P 
such that any point sequence {x_i}_i=1^n satisfying 

    k_0(x_j,x_j)/2   +∑_i=1^j-1 k_0(x_i, x_j) 
       ≤δ/2 + 
    	min_x∈ X: k_0(x,x)≤ R_j^2k_0(x,x)/2+∑_i=1^j-1 k_0(x_i, x)


with √(2log(j)/c_2)≤ R_j ≤∞ for each 1 ≤ j ≤ n also satisfies 

    D_K_0,P({x_i}_i=1^n) 
    	≤ e^π/2√(2log(n)/c_2n +c_1/n + δ/n).







 §.§ Distributional Convergence Control

To present our final results, we overload notation to define the KSD associated with any probability measure μ:

    D_𝒦_0,P(μ)     =    √((Z,Z')∼μ×μk_0(Z,Z')) .

Our original D_𝒦_0,P definition (Eq. <ref>) for a point set {x_i}_i=1^n is recovered when μ is the empirical measure 1/n∑_i=1^n δ_x_i.
We also write μ_m ⇒ P to indicate that a sequence of probability measures (μ_m)_m=1^∞ converges in distribution to P.

<cit.> showed that KSDs with IMQ base kernel <ref> and Langevin Stein operator T_P control distributional convergence 
 whenever P belongs to the set  of distantly dissipative distributions (i.e., log p(x) - log p(y)x-y≤ -κx-y^2 + C for some C≥ 0, κ >0) with Lipschitz log p.
Surprisingly, Gaussian, , and other kernels with light tails do not satisfy this property <cit.>.

Our next theorem establishes distributional convergence control for our newly introduced log inverse kernel <ref>.


Suppose P ∈.
Consider a Stein reproducing kernel k_0 = T_P𝒯_Pk_2
with Langevin operator T_P and base kernel
k_2(x,x') = (α + log(1 + x-x'^2))^β for
α > 0 and β < 0. If D_𝒦_0,P(μ_m ) →
0, then μ_m ⇒ P.


Our final theorem, proved in imq-score-implies-weak-convergence-proof, guarantees distributional convergence control for the new IMQ score kernel <ref> under the additional assumption that log p is strictly concave.


Suppose P ∈ has strictly concave log density. 
Consider a Stein reproducing kernel k_0 = T_P𝒯_Pk_3
with Langevin operator T_P and base kernel
k_3(x,x') = (c^2 + log p(x)-log p(x')^2)^β for c>0 and
β∈ (-1,0). If D_𝒦_0,P(μ_m ) → 0, then μ_m ⇒ P.





§ CONCLUSION
 

This paper proposed and studied Stein Points, extensible point sequences rooted in minimisation of a KSD, building on the recent theoretical work of <cit.>.
Although we focused on KSD to limit scope, our methods could in fact be applied to any computable Stein discrepancy, even those not based on reproducing kernels <cit.>.
Stein Points provide an interesting counterpoint to other recent work focussing on point sequences <cit.> and point sets <cit.>.
Moreover, when X is a finite set {y_i}_i=1^N (e.g., an inexpensive initial point set generated by MCMC), Stein Points provide a compact and convergent approximation to the optimally weighted probability measure ∑_i=1^N w_iδ_y_i with minimum KSD (see Section <ref> for more details).

Theoretical results were provided which guarantee the asymptotic correctness of our methods.
However, we were only able to establish an O(√(log(n)/n)) rate, which leaves a theoretical gap between the faster convergence that was sometimes empirically observed.
Relatedly, the O(n^2) computational cost could be reduced to O(n) by using finite-dimensional kernels <cit.>, but the associated distributional convergence control results must first be developed.

Our experiments were relatively comprehensive, but we did not consider other Stein operators, nor higher-dimensional or non-Euclidean manifolds X.
Related methods not considered in this work include those based on optimal transport <cit.> and self-avoiding particle-based samplers <cit.>.
The comparison against these methods is left for future work.







§ ACKNOWLEDGEMENTS

WYC was supported by the ARC Centre of Excellence for Mathematical and Statistical Frontiers.
FXB was supported by EPSRC [EP/L016710/1, EP/R018413/1].
CJO was supported by the Lloyd's Register Foundation programme on data-centric engineering at the Alan Turing Institute, UK.
This material was based upon work partially supported by the National Science Foundation under Grant DMS-1127914 to the Statistical and Applied Mathematical Sciences Institute. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.







langley00


icml2018





Supplement

This electronic supplement is organised as follows:
In Section <ref> proofs for the theoretical results in the main text are provided.
In Section <ref> we provide details for the two existing methods (MED, SVGD) that formed our experimental benchmark.
Then, in Section <ref>, we provide additional numerical results that elaborate on those reported in the main text.



  
Code

Code to reproduce these experiments is available from:


<github.com/wilson-ye-chen/stein_points>





§ PROOF OF THEORETICAL RESULTS IN THE MAIN TEXT
 








 §.§ Proofs of herding-convergence and greedy-convergence: Stein Herding and Stein Greedy Convergence

We will show that both herding-convergence and greedy-convergence follow from the following unified Stein Point convergence result, proved in stein-point-convergence-proof.



Suppose k_0 with k_0,P = 0 is a P-sub-exponential reproducing kernel.
Then there exist constants c_1, c_2 > 0 depending only on k and P 
such that any point sequence {x_i}_i=1^n satisfying 

    k_0(x_j,x_j)/2   +∑_i=1^j-1 k_0(x_i, x_j) 
    	≤δ/2 + S_j^2/2
    	+ min_x∈ X: k_0(x,x)≤ S_j^2∑_i=1^j-1 k_0(x_i, x)

with S_j ∈ [√(2log(j)/c_2), √(2log(n)/c_2)] for each 1 ≤ j ≤ n and δ≥ 0 also satisfies 

    D_K_0,P({x_i}_i=1^n) 
    	≤ e^π/2√(2log(n)/c_2n +c_1/n + δ/n).






  §.§.§ Proof of herding-convergence: Stein Herding Convergence


Instantiate the constants c_1,c_2>0 from stein-point-convergence, and 
consider any point sequence {x_i}_i=1^n satisfying 

    ∑_i=1^j-1 k_0(x_i, x_j) ≤δ/2 
    	+ min_x∈ X: k_0(x,x)≤ R_j^2∑_i=1^j-1 k_0(x_i, x)

with k_0(x_j,x_j)≤ R_j^2∈[2log(j)/c_2,2log(n)/c_2].
We immediately have

    k_0(x_j,x_j)/2   +∑_i=1^j-1 k_0(x_i, x_j) 
    	≤δ/2 + R_j^2/2
    	+ min_x∈ X: k_0(x,x)≤ R_j^2∑_i=1^j-1 k_0(x_i, x)

so the desired conclusion follows from stein-point-convergence.




  §.§.§ Proof of greedy-convergence: Stein Greedy Convergence


Instantiate the constants c_1,c_2>0 from stein-point-convergence, and 
consider any point sequence {x_i}_i=1^n satisfying 

    k_0(x_j,x_j)/2   +∑_i=1^j-1 k_0(x_i, x_j) 
    	≤δ/2 + 
    	min_x∈ X: k_0(x,x)≤ R_j^2k_0(x,x)/2+∑_i=1^j-1 k_0(x_i, x)


with S_j = √(2log(j)/c_2)≤ R_j ≤∞ for each 1 ≤ j ≤ n.
We immediately have

    k_0(x_j,x_j)/2   +∑_i=1^j-1 k_0(x_i, x_j) 
    	≤δ/2 + 
    	min_x∈ X: k_0(x,x)≤ S_j^2k_0(x,x)/2+∑_i=1^j-1 k_0(x_i, x) 
    	≤δ/2 + S_j^2/2
    	+ min_x∈ X: k_0(x,x)≤ S_j^2∑_i=1^j-1 k_0(x_i, x),

so the desired conclusion follows from stein-point-convergence.




  §.§.§ Proof of stein-point-convergence: Stein Point Convergence


Our high-level strategy is to show that, when k_0 is P-sub-exponential, 
optimizing over a suitably truncated search space on each step is sufficient to optimize the discrepancy globally.
To obtain an explicit rate of convergence, we adapt the greedy approximation error analysis of <cit.>, which applies to uniformly bounded kernels.
We begin by fixing any sequence of truncation levels (S_j)_j=1^∞ with each S_j ∈ [0,∞), defining the truncation sets B_j = { x∈ X: k_0(x,x) ≤ S_j^2},
and letting ℳ_j denote the convex hull
of { k_0(x,·) }_x∈ B_j.
Next we identify a truncation-optimal h_j ∈_f ∈ℳ_j J(f).
Now, fix any point sequence {x_i}_i=1^n satisfying 

    k_0(x_j,x_j)/2   +∑_i=1^j-1 k_0(x_i, x_j) 
    	≤δ/2 + S_j^2/2
    	+ min_x∈ X: k_0(x,x)≤ S_j^2∑_i=1^j-1 k_0(x_i, x)

for some approximation level δ≥ 0 and each 1 ≤ j ≤ n.
In the remainder, we will recursively bound the discrepancy of this point sequence in terms of each S_j and h_j__0, bound each h_j__0 in terms of S_j using the P-sub-exponential tails of k_0, and show that an appropriate setting of each S_j delivers the advertised claim.



  
Bounding discrepancy
For each j, let f_j = 1/j∑_i=1^j k_0(x_i,·) and _j = D_K_0,P({x_i}_i=1^j) =  f_j_K_0.
By Cauchy-Schwarz and the arithmetic-geometric mean inequality, we have the estimates

    n^2 _n^2 - δ   = k_0(x_n,x_n) + (n-1)^2 _n-1^2 + 2(n-1)f_n-1(x_n) - δ
       ≤ S_n^2 + (n-1)^2 _n-1^2 + 2(n-1)min_x ∈ B_n f_n-1(x)
       = S_n^2 + (n-1)^2 _n-1^2 + 2(n-1)inf_f ∈M_nff_n-1_K_0
       ≤ S_n^2 + (n-1)^2 _n-1^2 + 2(n-1)h_nf_n-1_K_0
       ≤ S_n^2 + (n-1)^2 _n-1^2 + n^2h_n_K_0^2 + (n-1)^2/n^2_n-1^2.

Unrolling the recursion, we obtain

    n^2 _n^2 
    	   ≤∑_i=0^n-1 (S_n-i^2 + δ + h_n-i_K_0^2 (n-i)^2) ∏_j=1^i (1+1/(n-j+1)^2).

Moreover, the products in this expression are uniformly bounded in i as

    log(∏_j=1^i (1+1/(n-j+1)^2))
    	= ∑_j=1^ilog( (1+1/(n-j+1)^2))
    	≤∫_0^∞log(1+1/x^2) = π.

Therefore, 

    n^2 _n^2 
    	   ≤ e^π∑_i=1^nS_i^2 + δ + i^2h_i_K_0^2.




  
Bounding h_i__0
To bound each h_i__0, we consider the truncated mean embeddings

    k_i^-    := ∫ k_0(x,·) k_0(x,x) ≤ S_i^2dP(x) and
    
    k_i^+    := ∫ k_0(x,·) k_0(x,x) > S_i^2dP(x) = k_P - k_i^-.

Since k_P = 0, we have k_i^+_K_0 = k_i^-_K_0.
Moreover, since k_i^- ∈M_i, we deduce that

    h_i_K_0^2
    	   ≤k_i^-_K_0^2 = k_i^+_K_0^2
       = ∬ k_0(x,y) k_0(x,x) > S_i^2dP(x) k_0(y,y) > S_i^2dP(y) 
       ≤(∫√(k_0(x,x))k_0(x,x) > S_i^2dP(x) )^2 
       ≤∫ k_0(x,x)k_0(x,x) > S_i^2dP(x)

where the final two inequalities follow by Cauchy-Schwarz and Jensen's inequality.

Let Y = k_0(Z,Z) for Z ∼ P.
We will bound the tail expectation in the final display by considering the biased random variable Y^* = k_0(Z^*,Z^*)
for Z^* with density ρ(z^*) = k_0(z^*,z^*) p(z^*)/[Y].
By <cit.>, since Y is sub-exponential, there exists c_0> 0 such that
[e^λ Y] < ∞ for all |λ| ≤ c_0.
For any λ≠ 0 with |λ| ≤ c_0/2, we have, by the relation x ≤ e^x,

    [e^λ Y^*] 
    	= [e^λ k_0(Z^*,Z^*)] 
    	= [k_0(Z,Z)e^λ k_0(Z,Z)]/[Y]
    	= [λ Y e^λ Y] /λ[Y]≤[e^2λ Y] /λ[Y]
    	< ∞.

Hence, by <cit.>, Y^* is also sub-exponential and satisfies, for some c̃_1, c_2 > 0,
Y^* ≥ t ≤c̃_1 e^-c_2t for all t > 0.

Applying this finding to the bounding of h_i, we obtain

    h_i_K_0^2
    	≤∫ k_0(x,x)k_0(x,x) > S_i^2dP(x)
    	= 𝔼[Y]∫k_0(x,x) > S_i^2ρ(x)
    	= 𝔼[Y]Y^* ≥ S_i^2≤ c_1 e^-c_2S_i^2

where c_1 = c̃_1 𝔼[Y].
Hence

    D_K_0,P({x_i}_i=1^n) 
    	   ≤ e^π/2√(1/n^2∑_i=1^nS_i^2 + δ + i^2c_1 e^-c_2S_i^2).




  
Setting each S_i
By choosing S_i ∈ [√(2log(i)/c_2),√(2log(n)/c_2)] for each i we obtain 

    D_K_0,P({x_i}_i=1^n) 
    	   ≤ e^π/2√(1/n^2∑_i=1^n2log(n)/c_2 + δ + c_1)≤ e^π/2√(2log(n)/c_2 n + δ/n + c_1/n).







 §.§ Proof of log-inverse-implies-weak-convergence: Log Inverse KSD Controls Convergence

Fix any α > 0 and β < 0.
Our proof will leverage <cit.>. This requires demonstrating two
separate properties for the log inverse kernel: first, the log inverse function Φ(z)  (α + log(1 +
z^2))^β has a nonvanishing
generalized Fourier transform, and second, whenever
D_𝒦_0,P(μ_m )→ 0, the measures μ_m are
uniformly tight. We will repeatedly use the notation γ(r)  (α +
log(1 + r))^β and ϕ(r) γ(r^2) throughout the
proof. Moveover, we will use f̂ to denote the
(generalized) Fourier transform of a function f, and V_d will represent
the volume of the unit Euclidean ball in d dimensions. Finally, we write f^(m) for the m-th derivative of any sufficiently differentiable function f:→.

To demonstrate the first property, we begin with the following lemma.



Fix any α > 0 and β < 0.
The function γ(r)  (α +log(1 + r))^β is completely monotone, i.e.,
γ∈ C^∞ and (-1)^m γ^(m)(r) ≥ 0 for all
m∈N_0 and all r ≥ 0,
and hence
the function k_2:R^d×R^d→R given by
k_2(x,x') γ(x-x'^2) is a
kernel function for all dimensions d∈N.




By <cit.> we know that Φ is positive semidefinite for
all dimensions d∈N if and only if γ is completely monotone. Thus it remains to show that γ is
completely monotone.

Since α > 0, γ(r) > 0 for all r ≥ 0. To verify (-1)^mγ^(m)(r)≥ 0 for all m≥
1, we will proceed by induction. Let us suppose that for some m≥ 1,

    γ^(m)(r) = (-1)^m∑_l=1^m c_l,m (α + log(1 + r))^β-l
    (1 + r)^-m

where each c_l,m∈R is positive. Taking another derivative
yields

    γ^(m+1)(r)
      
      
      
      
      
      
         = (-1)^m+1∑_l=1^m+1 c_l,m+1 (α + log(1 + r))^β-l (1 + r)^-m-1,

where c_1, m+1 m  c_1,m, c_l,m+1 m  c_l, m +
(l-β-1)  c_l-1,m for l > 1 and c_l,m 0 for all l > m,
completing the induction step.

As for the base case, notice γ'(r) = β(α + log(1 + r))^β-1 (1 +
r)^-1, which establishes the identity for l = 1 by setting
c_1,1 -β. The conclusion of this
proof by induction implies (-1)^m γ^(m)(r)≥ 0 for all m and all
r ≥ 0. By inverse-log-radial-derivs, γ∈ C^∞,
establishing the lemma.


Knowing that γ is a completely monotone function, we can now demonstrate
Φ̂ has a nonvanishing generalized Fourier transform.



Consider the function Φ:^d→ given by Φ(z) = (α + log(1 +
z^2))^β for some α > 0 and β < 0. Its
generalized Fourier transform Φ̂(w) is radial, nonvanishing, and
continuous for w≠ 0. Moreover, Φ̂(w)→ 0 as w→∞.





We will first use induction to prove an intermediate result
that states for any m∈N_0,

    ^m Φ(z) = ∑_(u,v)∈ S_mτ_u,vz^2vγ^(u)(z^2)

where τ_u,v > 0 are positive reals, S_m = {(u, v)∈N_0^2
 |  v≤ u - m, u≤ 2m} and γ(r) (α + log(1 + r))^β.

Note for the base case m=0, the claim above for ^0 Φ = Φ
clearly holds. Now suppose it holds from some m∈N_0.
If A:R^d→R is a function that can decomposed as
A(z)  f(z^2) g(z^2)
where f,g ∈ C^∞([0,∞)), then we have

    A(z)
      
      
      
      
      
         =  [2d g'(z^2) + 4z^2 g”(z^2) ]
      f(z^2) +
       [2d g(z^2) + 4z^2 g'(z^2) ]
      f'(z^2) +
      4z^2 g(z^2) f”(z^2).

Consider each term in the decomposition of ^m Φ(z) from the induction
hypothesis. If we let g(r) = r^v and f(r) = ϕ^(u)(r), we see that
each term from laplacian-product-rule is of the form τ_u,v'
z^2v'ϕ^(u')(z^2) where the values for (u',
v') are (u, v-1), (u, v-1), (u+1, v), (u+1, v), (u+2, v+1)
respectively. Notice that when v=0 or v=1, the first or second
derivative of g will be zero and these terms may disappear altogether. Thus all
these tuples will lie in S_m+1 for any (u,v)∈ S_m, and so we must have
^m+1Φ(z) satisfies the induction hypothesis as well, completing
the proof by induction.








Now we can prove the lemma. Suppose 2m ≥ d. Then by the triangle inequality and a
radial substitution <cit.>,

    ∫_R^d |^m Φ(z)|≤∑_(u,v)∈ S_m∫_R^dτ_u,vz^2v
        |ϕ^(u)(z^2)| 
        = d V_d ∑_(u,v)∈ S_m∫_0^∞τ_u,v r^2v + d - 1
        |ϕ^(u)(r^2)| .

Because |ϕ^(u)(r)| = O(r^-ulog^β-1(r)) as r→∞ for u∈N by
inverse-log-radial-derivs, we see that each integrand above is O(r^2(v-u)
  + d - 1log^β-1(r)). But since v ≤ u - m, this will imply that
each integrand is O(r^-2m + d - 1log^β-1(r)), which is integrable for
large r yielding ^m Φ∈ L^1(R^d).

By <cit.> and the fact that positive definiteness is preserved by summation,
we have ^m Φ is a positive definite
function. This along with the fact that ^m Φ∈ L^1(^d)
allows us to invoke <cit.> and <cit.> to obtain ^m Φ is continuous,
radial and nonvanishing. Moreover, ^m Φ belonging to
L^1(^d) implies its Fourier transform belongs to L^∞(^d).
The lemma follows by noticing ^m
  Φ(w) = w^2mΦ̂(w), i.e.,
Φ̂(w) = w^-2m^m Φ(w)
for all w≠ 0.



We now need to demonstrate the second property to complete the proof of
log-inverse-implies-weak-convergence, but in order to do so, we
first will establish the lemma below.  By
log-inverse-nonvanishing-gft, we know Φ̂ is radial and
thus can write Φ̂(w) = ϕ_∧(w) for some
continuous function ϕ_∧:(0,∞)→ (0, ∞). Our first
priority will be to lower bound ϕ_∧ near the origin.



If Φ is the log inverse function on ^d from
log-inverse-nonvanishing-gft,
then lim inf_r→ 0^+
r^d(α + log(1+1/r^2))^-β+1ϕ_∧(r) > 0 where
Φ̂(w) = ϕ_∧(w) for all w≠ 0.



First we will show that ϕ_∧ is strictly decreasing. Since
r ↦ (α + log(1 + r))^β was shown to be completely
monotone in log-inverse-is-completely-monotone, by <cit.> we must have
Φ(z) = ∫_0^∞ e^-t z^2∂ v(t)
for some finite, non-negative Borel measure v on [0,∞) that is not
  concentrated at zero. Let (φ_m)_m=1^∞ be a sequence of Schwartz functions <cit.> defined on
  ^d. Then, for each m, both φ̂_m and Φφ̂_m are also Schwartz
functions, and thus

    ∫_^d [∫_0^∞ |e^-t x^2φ̂_m(x)|
    ∂ v(t)  ] 
      = ∫_^d [∫_0^∞ e^-t x^2 |φ̂_m(x)|
        ∂ v(t)  ] 
      = ∫_^dΦ(x) |φ̂_m(x)|  < ∞,

as all Schwartz functions are integrable.  This allows us to use Fubini's
theorem in conjunction with Plancherel's Theorem to argue

    ∫_^dΦ̂(w)φ_m(w)
      = ∫_^dΦ(x)φ̂_m(x)   = ∫_^d [∫_0^∞ e^-t x^2∂
        v(t) ] φ̂_m(x)
       = ∫_^d∫_0^∞ e^-t x^2φ̂_m(x) ∂
        v(t) 
       = ∫_0^∞∫_^d  e^-t x^2φ̂_m(x)∂ v(t) 
       = ∫_0^∞∫_^d (2t)^d/2 e^-1/4tw^2φ_m(w)∂ v_+(t) + v_0 φ_m(0),

where we have used the decomposition v  v_+ + v_0 δ_0 for  v_0
≥ 0 and v_+ non-zero and absolutely continuous with respect to
Lebesgue measure on [0,∞). Let ℬ:^d→ be a bump
function, e.g., ℬ(x) 𝒵^-1exp{-1/(1-x^2)}x<1 where 𝒵 is the
normalization constant chosen such that ∫_^dℬ(x)
= 1. Then let us define φ_m:^d→
via the mapping φ_m(w)  m^d ℬ(m(w-w_0e_1)) - m^d ℬ(m(w-w_1e_1)),
where 0 < w_0 < w_1 and e_1∈^d is the first standard basis vector. Then ∫_^dΦ̂(w)φ_m(w)→Φ̂(w_0e_1) - Φ̂(w_1e_1) =
ϕ_∧(w_0) - ϕ_∧(w_1)
since Φ̂ is a continuous in neighborhoods of w_0e_1 and w_1e_1
<cit.>.

Because v_+ cannot be the zero measure, there must be some finite interval
[a_0,b_0]⊂ (0,∞) such that v_+([a_0,b_0])>0.
For each t > 0 and m > max(1/w_0,2/w_1 - w_0), we have

    A_m(t)∫_^d (2t)^d/2 e^-1/4tw^2φ_m(w)
    	=  ∫_^d (2t)^d/2 (e^-1/4tw-w_0e_1^2 - e^-1/4tw-w_1e_1^2)m^d ℬ(mw) 
    	> 0,

since w-w_0e_1 < w-w_1e_1 when w < min(w_0 - w_1/2, w_0).
Using <cit.> again,
we have A_m(t)→ (2t)^d/2(e^-1/4tw_0^2 -
e^-1/4tw_1^2) as m→∞ for any t > 0.
Moreover, for all t ∈ [a_0,b_0] and m≥ 1, we have

    |A_m(t)| ≤∫_^d (2t)^d/2 e^-1/4tw-w_0e_1^2m^d ℬ(mw)    ≤ (2b_0)^d/2sup_w < 1 e^-1/4b_0w-w_0e_1^2 < ∞.

Hence, the dominated convergence theorem allows us to exchange the limit over m and integral
over t below to conclude

    ϕ_∧(w_0) - ϕ_∧(w_1)
      = lim_m→∞∫_0^∞ A_m(t)  ∂ v_+(t) + v_0 φ_m(0)
      ≥lim_m→∞∫_a_0^b_0 A_m(t)  ∂ v_+(t) 
      = ∫_a_0^b_0lim_m→∞ A_m(t)  ∂ v_+(t) 
       = ∫_a_0^b_0 (2t)^d/2(e^-1/4tw_0^2 -
    e^-1/4tw_1^2) ∂ v_+(t)
      ≥ v_+([a_0,b_0]) min_t∈ [a_0,b_0]{(2t)^d/2
      (e^-1/4tw_0^2 -
      e^-1/4tw_1^2)} > 0,

showing ϕ_∧ is strictly decreasing as claimed.

Suppose ψ:[0,∞)→ is a C^∞ function
with support [a,b] for 0 < a < b such that ψ(r) > 0 for all r ∈ (a,b) and
∫_0^∞ψ(r)  = 1. Then because ϕ_∧ is strictly decreasing,
by the mean value theorem we have

    ϕ_∧(b/λ) ≤∫_0^∞λϕ_∧(r)ψ(λ r)≤ϕ_∧(a/λ)

for all λ > 0. If we assign Ψ(w)ψ(w) to be
the radial continuation of ψ, by <cit.> the quantity
sandwiched above becomes

    ∫_0^∞λϕ_∧(r)ψ(λ r)
      = ∫_0^∞ϕ_∧(s/λ)ψ(s)
      = 1/d V_d∫_^dΦ̂(w/λ) Ψ(w)/w^d-1.

Next suppose that ξ:[0, ∞)→ is a Schwartz function
satisfying ξ^(k)(0) = 0 for all integral k≥ 0, and
let Ξ:^d→ given by Ξ(x) ξ(x) be
the radial continuation of ξ. Then by Plancherel's Theorem,
scaling the input of a Fourier transform as in <cit.>,
and the change to spherical coordinates in <cit.>, for any
λ > 0, we have

    ∫_^dΦ̂(w/λ) Ξ(w)
      = ∫_^dΦ(w) Ξ̂(w/λ)
      = d V_d ∫_0^∞ r^d-1ϕ(r) ξ_∧(r/λ) 
      = d V_d λ^d ∫_0^∞ s^d-1ϕ(λ s) ξ_∧(s)  d s,

where s = r/λ and ξ_∧ is the radial function associated with Ξ̂,
i.e., Ξ̂(w) = ξ_∧(w) for all w. 
















Let us define ω :[0,∞)→ by the mapping ω(t) 
  (α + t)^β. Then by the mean value theorem and the fact that
  ω' is increasing, we have for all s > 1

    -ω'(log(1 + λ^2s^2)) ≤ -ω(log(1 + λ^2s^2)) -
      ω(log(1 + λ^2))/log(1 + λ^2s^2) - log(1 + λ^2)≤ -ω'(log(1 + λ^2)).

By rearranging terms, this implies for all λ > 0

    (-β)(α + log(1 + λ^2s^2)/α + log(1 +
        λ^2) )^β-1log (1 + λ^2s^2/1 + λ^2 )
      ≤ -ω(log(1 + λ^2s^2)) -
        ω(log(1 + λ^2))/ω(log(1 + λ^2))(α + log(1 + λ^2))^-1≤ (-β)log (1 + λ^2s^2/1 + λ^2 ).

Since log(1+λ^2s^2/1+λ^2)→ 2log s as
λ→∞, and the sandwiched term above is -(α +
log(1+λ^2))(ϕ(λ s)/ϕ(λ) - 1), we have
(α + log(1+λ^2))(ϕ(λ s)/ϕ(λ) - 1)→ 2βlog s
as λ→∞ for all s > 1. The case for s ∈ (0,1] is
analogous and yields the same asymptotic limit.

With this new asymptotic expansion in hand, we will revisit phi-xi.
We have

    λ^-dϕ(λ)^-1(α + log(1+λ^2))∫_^dΦ̂(w/λ) Ξ(w)    = d  V_d ϕ(λ)^-1(α + log(1+λ^2))∫_0^∞ϕ(λ s) s^d-1ξ_∧(s)
       = d  V_d  (α + log(1+λ^2))∫_0^∞ϕ(λ s)/ϕ(λ) s^d-1ξ_∧(s)
       = d  V_d ∫_0^∞
     (α + log(1+λ^2))  [ϕ(λ
         s)/ϕ(λ) - 1  ] s^d-1ξ_∧(s).

Notice that final integrand converges to 2β s^d-1 (log s) ξ_∧(s)
pointwise for all s ≥ 0 as λ→∞. Since ξ_∧ is
a Schwartz function on [0,∞), we can utilize the fact that s↦log s is integrable near the origin to reason that s^d-1 (log s)
ξ_∧(s) is a Schwartz function as well, and thus integrable. Hence 
by the dominated convergence theorem, we have the integral above converges to
2 β  d  V_d ∫_0^∞  s^d-1 (log s) ξ_∧(s)
as λ→∞.

Now suppose we choose Ξ(x)x^1-dΨ(x). By
sandwich-gft-near-origin we have

    lim_λ→∞λ^-dϕ(λ)^-1(α + log(1+λ^2)) ϕ_∧(b/λ) ≤
    2β∫_0^∞ s^d-1(log s)ξ_∧(s)≤lim_λ→∞λ^-dϕ(λ)^-1(α + log(1+λ^2)) ϕ_∧(a/λ).

By log-inverse-nonvanishing-gft, we know ϕ_∧(r) > 0 for
all r > 0, and thus the left-hand side above must be non-negative. Hence if we can
show for some choice of ψ that the sandwiched term is non-zero, then the proof of
the lemma will follow from choosing r = a/λ.

Let us define L(x) = logx with generalized Fourier transform L̂. 
As usual, let
l:[0,∞)→ and l_∧:[0,∞)→ be the radial
    functions associated with L and L̂.  Notice that again by
    Plancherel's Theorem

    ∫_0^∞ s^d-1(log s)ξ_∧(s)
      = 1/d V_d∫_^dΞ̂(w)L(w) 
      = 1/d V_d∫_^dΞ(x)L̂(x)    = 1/d V_d∫_^dΨ(x)/x^d-1L̂(x) 
       = ∫_0^∞ψ(r) l_∧(r).

Since we are free to choose ψ to be any Schwartz function with support
[a,b], if we could not find a function ψ such that
the quantity in multivariate-log-gft-limit is non-zero, this would
imply the support of l_∧ is a subset of {0}. But this would
mean l_∧ is some multiple of a point mass at zero, which would
imply l is a constant function, a contradiction. Thus we must be able to
find some ψ such that the integral above is non-zero, completing the
lemma.



Fix any a_0 > 0
and α_0∈ (0,1/2).
Our strategy for showing the KSD controls tightness will mimic
<cit.>: we will
show that a bandlimited approximation of the function g_j(x) = 2α_0 x_j(a_0^2 +
x^2)^α_0-1 belongs to the inverse log RKHS and thus enforces tightness.

First note that in the proof of <cit.>, it was shown
h = T_P g was a coercive, Lipschitz, and bounded-below function for P∈. Moreover,
in the proof of <cit.>, a random vector Y with
density ρ(y) is constructed such that
the support of ρ̂ belongs to [-4,4]^d and also Y is
integrable. Consider the new function g^∘(x) g(x +
  Y) for all x∈^d. By the convolution theorem, ĝ^∘_j
= ĝ_j ρ̂ and so g^∘_j is bandlimited for all j.
In the proof of <cit.>, ĝ_j was shown to grow
asymptotically at the rate (iw_j) w^-d-2α_0 as
w→ 0. Thus

    ∑_j=1^d ∫_^dĝ^∘_j(w)ĝ^∘_j(w)/Φ̂(w)
        = ∑_j=1^d ∫_[-4,4]^dĝ_j(w)ĝ_j(w)ρ̂(w)^2/Φ̂(w)   ≤κ_0 ∫_[-4,4]^dw^-2d-4α_0+2/Φ̂(w)
       ≤κ_1 ∫_0^4√(d) r^-4α_0+1log^-β+1(1 + r^-2),

for some constants κ_0,κ_1 > 0 where we used
log-inverse-gft-lower-bound in the final inequality.
This integral is finite for
all α_0 ∈ (0,1/2) and any β < 0, which implies
g^∘ is in the log inverse RKHS by <cit.>.

Finally, notice that via the argument proving <cit.>,

    sup_x∈^d |T_P g^∘(x) - h(x)|
      ≤3dlog 2/π(sup_x∈^d h(x) +
      sup_x∈^d^2log p(x)·sup_x∈^dg(x)) < ∞.

Since h is bounded below and coercive, these properties are inherited by
T_Pg^∘. This allows us to apply <cit.> to argue D_𝒦_0,P(μ_m )→ 0
implies the measures μ_m are uniformly tight. Combining this with 
log-inverse-nonvanishing-gft allows us to utilize <cit.> for the log inverse kernel, thereby concluding the proof.






 §.§ Proof of imq-score-implies-weak-convergence: IMQ Score KSD Convergence Control


For b = log p, introduce the alias k_b = k_3, let _b denote the RKHS of k_b, and let C_c represent the set of continuous compactly supported functions on X.
Since P∈, the proof of Thm. 13 in <cit.> shows that 
if, for each h∈ C^1∩ C_c and > 0, there exists h_∈_b
such that sup_x ∈ X |(T_Ph)(x) - (T_Ph_)(x)| ≤, then μ_m ⇒ P
whenever D_𝒦_0,P(μ_m ) → 0 and (μ_m)_m=1^∞ is uniformly tight.
Hence, to establish our result, it suffices to show 
(1) that, for each h∈ C^1∩ C_c and > 0, there exists h_∈_b
such that sup_x ∈ Xmax( (h-h_)(x), b(x) (h-h_)(x)) ≤ϵ
and 
(2) that D_𝒦_0,P(μ_m ) → 0 implies (μ_m)_m=1^∞ is uniformly tight.



  §.§.§ Approximating C^1 ∩ C_c with _b

Fix any f∈ C^1∩ C_c and > 0, and let  denote the RKHS of k(x,y) = (c^2 + x-y^2)^β.
Since p is strictly log-concave, b is invertible with ( b(x)) never zero.
Since P ∈, b is Lipschitz.
By the following theorem, proved in composition-kernel-properties-proof, it therefore suffices to show that there exists f_∈ such that sup_x ∈ Xmax( (f-f_)(x),  x (f-f_)(x))≤ϵ.


For b : X → X invertible and k a reproducing kernel on X with induced RKHS , define the composition kernel k_b(x,y) = k(b(x), b(y))
with induced RKHS _b.




Suppose that, for each f∈ C^1∩ C_c and > 0, there exists f_∈ such that 
    sup_x ∈ Xmax( (f-f_)(x),  x (f-f_)(x))≤ϵ.

If b is Lipschitz and ( b(x)) is never zero, then, for each h∈ C^1∩ C_c and > 0, there exists h_∈_b
such that 
    sup_x ∈ Xmax( (h-h_)(x), b(x) (h-h_)(x)) ≤ϵ.




Since the identity map x↦ x is Lipschitz and f ∈ L^2 because it is continuous and compactly supported, <cit.> provides an explicit construction of f_∈ satisfying our desired property whenever k(x,y) = Φ(x-y) for Φ∈ C^2 with non-vanishing Fourier transform.  Our choice of IMQ k satisfies these properties by <cit.>.




  §.§.§ Controlling Tightness


Since P is distantly dissipative,

    -b(x)x≤b(x)x≤ -κx^2+C + b(0)x≤ -κx^2+C + b(0)x

by Cauchy-Schwarz.
Hence, b is norm-coercive, i.e., b(x)→∞ whenever x→∞.
Since b is bounded, our desired result follows from the following lemma
which guarantees tightness control on b under weaker conditions.


If b : X → X is norm coercive and differentiable, and _j b_j(x) = o(b(x)^2) as x→∞, then lim sup_m D_𝒦_0,P(μ_m ) < ∞ implies (μ_m)_m=1^∞ is tight.



Fix any a >c/2 and α∈(0,1/2(β + 1)). 
The proof of <cit.> showed that the function g_j(x) = 2α x_j (a^2 + x^2)^α - 1∈ for each j ∈{1,…, d}.
Hence g_b,j(x)  g_j(b(x)) ∈_b for each j ∈{1,…, d} by composition-rkhs.
By our assumptions on b, we have

    (T_P g_b)(x) 
    	   = 2α (b(x)^2 (a^2 + b(x)^2)^α - 1 
    	+ ∑_j=1^d _j b_j(x)  (a^2 + b(x)^2)^α - 1
    	+ b_j(x)^2 2(α - 1)(a^2 + b(x)^2)^α - 2_j b_j(x))
       = 2αb(x)^2 (a^2 + b(x)^2)^α - 1  + o(b(x)^2α),

so T_P g_b is coercive, and the proof of <cit.> therefore gives the result (μ_m)_m=1^∞ is uniformly tight whenever lim sup_m D_𝒦_0,P(μ_m ) finite. 






 §.§ Proof of composition-kernel-properties: Composition Kernel Approximation






Let c=b^-1 represent the inverse of b, and for any function f on X, let f_c(y) = f(c(y)) denote the composition of f and c so that f_c(b(x)) = f(x).
The following lemma shows that f_c inherits many of the properties of f under suitable restrictions on b.

 
For any function f on X and invertible function b on X, define f_c(y) = f(c(y)) for c = b^-1.
The following properties hold.

	
  * If f has compact support and b is continuous, then f_c has compact support.
	
  * If f∈ C^1, b ∈ C^1, and ( b(x)) is never zero, then f_c ∈ C^1. 







 



We prove each claim in turn.



  * If f is compactly supported and b is continuous, then f_c = b(f) is also compact, since continuous functions are compact-preserving <cit.>.


  * If f∈ C^1, b ∈ C^1, and ( b(x)) is never zero, then c is continuous by the inverse function theorem <cit.>, x↦ ( b(x))^-1 is continuous,
and hence f_c(y) = ( c(y)) ( f)(c(y)) = (( b)(c(y)))^-1 ( f)(c(y)) is continuous.

   

Our next lemma exposes an important relationship between the RKHSes  and _b.


Suppose f is in the RKHS  of a reproducing kernel k on X and b : X → X is invertible.
Then f_b is in the RKHS _b of k_b for f_b(x) = f(b(x)) and k_b(x,y) = k(b(x), b(y)).




Since f∈, there exist f_m = ∑_j=1^J_m a_m,j k(x_m,j,·) for m∈, a_m,j∈, and x_m,j∈ X such that lim_m→∞f_m - f_ = 0 and lim_m→∞ f_m(x) = f(x) for all x∈ X.
Now let c=b^-1, and define 

    f_m,b(x) = f_m(b(x)) = ∑_j=1^J_m a_m,j k(x_m,j,b(x)) = ∑_j=1^J_m a_m,j k_b(c(x_m,j),x).

Since _b = {∑_j=1^J a_j k_b(y_j,·) : J∈, a_j∈, y_j∈ X},
each f_m,b∈_b. 
Since (f_m)_m=1^∞ is a Cauchy sequence, and f_m,bf_m',b__b^2 =  ∑_j=1^J_m a_m,j∑_j'=1^J_m' a_m',j'k_b(c(x_m,j),c(x_m',j'))
= f_mf_m'_ so that f_m,b-f_m',b__b =  f_m-f_m'_ for all m,m',
the sequence (f_m,b)_m=1^∞ is also Cauchy and converges in ·__b to its pointwise limit f_b.
Since an RKHS is complete, f_b ∈_b.


With our lemmata in hand, we now prove the advertised claim. 






Suppose b is Lipschitz, ( b(x)) is never zero, and for each f∈ C^1∩ C_c and > 0 there exists f_∈ such that sup_x ∈ Xmax( (f-f_)(x),  x (f-f_)(x))≤ϵ.
Select any h ∈ C^1∩ C_c and any > 0.
By composition-properties, h_c ∈ C^1∩ C_c,
and hence there exists h_c,∈ such that 
sup_y ∈ Xmax( (h_c-h_c,)(y),  y (h_c-h_c,)(y))≤/max(1,M_1(b)).
Now define h_(x) = h_c,(b(x)) so that h_∈_b by composition-rkhs.
We have sup_x∈ Xb(x)(h_ -h)(x)≤sup_y∈ Xy(h_c,-h_c)(y)≤, and  

    sup_x∈ X h_(x) -  h(x) 
    	= sup_x∈ X( b(x)) (( h_c,)(b(x)) - ( h_c)(b(x)))≤ M_1(b) /max(1,M_1(b)) ≤.






§ IMPLEMENTATIONAL DETAIL
 



 §.§ Benchmark Methods


In this section we briefly describe the MED and SVGD methods used as our empirical benchmark, as well as the (block) coordinate descent method that was used in conjunction with Stein Points.



  §.§.§ Minimum Energy Designs
 

The first class of method that we consider is due to <cit.>.
That work restricted attention to X = [0,1]^d and constructed an energy functional:

    ℰ_δ,P({x_i}_i=1^n)     :=    ∑_i ≠ j[  p(x_i)^-1/2d p(x_j)^-1/2d/x_i - x_j_2 ]^δ

for some tuning parameter δ∈ [1,∞) to be specified.
In <cit.> the rule-of-thumb δ = 4d was recommended.
A heuristic argument in <cit.> suggests that the points {x_i}_i=1^n that minimise ℰ_δ,P({x_i}_i=1^n) form an empirical approximation that converges weakly to P.
The argument was recently made rigorous in <cit.>.

Minimisation of ℰ_δ,P does not require knowledge of how p is normalised.
However, the actual minimisation of ℰ_δ,P can be difficult.
In <cit.> an extensible (greedy) method was considered, wherein the first point is selected as

    x_1 ∈_x ∈ X   p(x)

and subsequent points are selected as

    x_n    ∈   _x ∈ X   p(x)^-δ/2d∑_i=1^n-1p(x_i)^-δ/2d/x_i - x_2^δ .

However, alternative approaches could easily be envisioned.
For instance, if n were fixed then one could consider e.g. applying the Newton method for optimisation over the points {x_i}_i=1^n.

Remark: There is a connection between certain minimum energy methods and discrepancy measures in RKHS; see <cit.>.

Remark: Several potential modifications to ℰ_δ,P were suggested in <cit.>, but that report appeared after this work was completed.
These could be explored in future work.

Remark: The MED objective function is typically numerically unstable due to the fact that the values of the density p(·) can be very small.
In contrast, our proposed methods operate on log p(·) and its gradient, which is more numerically robust.






  §.§.§ Stein Variational Gradient Descent
 

The second method that we considered was due to <cit.> and recently generalised in <cit.>.
The idea starts by formulating a continuous version of gradient descent on 𝒫(X) with the Kullback-Leibler divergence KL(· || P) as a target.
To this end, restrict attention to X = ℝ^d and consider the dynamics

    S_f(x)     =     x + ϵ f(x)

parametrised by a function f ∈𝒦^d.
For infinitesimal values of ϵ we can lift S_f to a pushforward map on 𝒫(X); i.e. Q ↦ S_f Q.
It was then shown in <cit.> that

    - . d/dϵKL(S_f Q || P) |_ϵ = 0    =    ∫𝒯_P f  dQ

where 𝒯_P is the Langevin Stein operator in Eqn. <ref>.
Recall that this operator can be decomposed as 𝒯_P f = ∑_j=1^d 𝒯_P,j f_j with 𝒯_P,j = ∇_j + ∇_j log p, where ∇_j denotes differentiation with respect to the jth coordinate in X.
Then the direction of fastest descent

    f^*(·)     :=    _f ∈ B(𝒦^d)  - . d/dϵKL(S_f Q || P) |_ϵ = 0

has a closed-form, with jth coordinate

    f_j^*(· ; Q)     =    ∫𝒯_P,j k(x,·)  dQ(x) .

The algorithm proposed in <cit.> discretises this dynamics in both space X, through the use of n points, and in time, through the use of a positive step size ϵ > 0, leading to a sequence of empirical measures based on point sets {x_i^m}_i=1^n for m ∈ℕ.
Thus, given an initialisation {x_i^0}_i=1^n of the points, at iteration m ≥ 1 of the algorithm we update 

    x_i^m     =     x_i^m-1 + ϵ f^*(x_i^m-1 ; Q_n^m)

in parallel, where

    Q_n^m     =    1/n∑_i=1^n δ_x_i^m-1

is the empirical measure, at a computational cost of O(n).
The output is the empirical measure Q_n^m.

Remark: The step size ϵ is a tuning parameter of the method.

Remark: At present there are not theoretical guarantees for this method.
Initial steps toward this goal are presented in <cit.>.








  §.§.§ Block Coordinate Descent
 

The Stein Point methods developed in the main text can be adapted to return a fixed number n of points for a given finite computational budget by first iteratively generating a size n point set, as described in the main text, and then performing (block) coordinate descent on this point set.
The (block) coordinate descent procedure is now described:

Fix an initial configuration {x_i^0}_i=1^n.
Then at iteration m ≥ 1 of the algorithm, perform the following sequence of operations:

    ∀ i        x_i^m ← x_i^m-1  then:
    for  i = 1,…,n        x_i^m ←_x ∈ X D_𝒦_0,P({x_j^m}_j ≠ i∪{x})

The output is the point set {x_i^m}_i=1^n.

Remark: The block coordinate descent method can equally be applied to MED; this was not considered in our empirical work.

Remark: Any numerical optimisation method can be used to solve the global optimisation problem in the inner loop.
In this work we considered the same three candidates in the main text; Monte Carlo, Nelder-Mead and grid search.
These are described next.








 §.§ Numerical Optimisation Methods
  

Computation of the nth term in the proposed Stein Point sequences, given the previous n-1 terms, requires that a global optimisation is performed over x_n ∈ X.
The same is true for both MED and KSD in the coordinate descent context.
For all experiments reported in the main text, three different numerical methods were considered for this task, denoted , ,  in the main text.
In this section we provide full details for how these methods were implemented.



  §.§.§ Nelder-Mead
  

The Nelder-Mead () method <cit.> proceeds as in Algorithm <ref>.
The function NM takes the following inputs: f is the objective function; t is the iteration count; n_init is the number of initial points to be drawn from a proposal distribution; n_delay is the number of iterations after which the proposal distribution becomes adaptive; μ_0 and Σ_0 are the mean vector and the covariance matrix of the initial proposal distribution; {x^curr_j}_j=1^n_curr is the set of existing points; λ is the variance of each mixture component of the adaptive proposal distribution; l and u are the lower- and upper-bounds of the search space. The non-adaptive initial proposal distribution is a truncated multivariate Gaussian 𝒩(μ_0, Σ_0) whose support is bounded by the hypercube [l,u]. The adaptive proposal distribution is a truncated Gaussian mixture Π({x^curr_j}_j=1^n_curr, λ) := 1/n_curr-1∑_j=1^n_curr-1𝒩(x^curr_j,λ I) with λ > 0 and support [l,u]. The expression NelderMead_x[f(x), x^init_i, l, u] denotes the standard Nelder-Mead procedure for objective function f, initial point x^init_i, and bound constraint x ∈ [l,u]. We use the symbol  to denote the assignment of a realised independent draw. The operator trunc_l^u[·] bounds the support of a distribution by the hypercube [l,u].







  §.§.§ Monte Carlo
  

The Monte Carlo () optimisation method proceeds as in Algorithm <ref>.
The function MC takes the following inputs: f is the objective function; t is the iteration count; n_test is the number of test points to be drawn from a proposal distribution; n_delay is the number of iterations after which the proposal distribution becomes adaptive; μ_0 and Σ_0 are the mean vector and the covariance matrix of the initial proposal distribution; {x^curr_j}_j=1^n_curr is the set of existing points; λ is the variance of each mixture component of the adaptive proposal distribution; l and u are the lower- and upper-bounds of the search space. The non-adaptive initial proposal distribution is a truncated multivariate Gaussian 𝒩(μ_0, Σ_0) whose support is bounded by the hypercube [l,u]. The adaptive proposal distribution is a truncated Gaussian mixture Π({x^curr_j}_j=1^n_curr, λ) := 1/n_curr-1∑_j=1^n_curr-1𝒩(x^curr_j,λ I) with λ > 0 and support [l,u].







  §.§.§ Grid Search
  

The grid search () optimisation method proceeds as in Algorithm <ref>.
The function GS takes the following inputs: f is the objective function; t is the iteration count; l and u are the lower- and upper-bounds of the grid; n_0 is the initial grid size.







 §.§ Remark on Application to a Reference Point Set


It is interesting to comment on the behaviour of our proposed methods in the case where X is a finite set or the global optimisation over X is replaced by a discrete optimisation over a pre-determined fixed set Y = {y_i}_i=1^N ⊆ X.
In this case it can be shown that:


  * The algorithm after n iterations will have selected n points {y_π(i)}_i=1^n with replacement from Y.
(Here π(i) indexes the point that was selected at iteration i of the algorithm.)

  * The empirical measure 1/n∑_i=1^n δ_y_π(i) can be expressed as ∑_i=1^N w_i y_i for some weights w_i.

  * The weights w_i converge to 

    (*)   =  _w ≥ 0 
     w_1 + … + w_N = 1 √(1/N^2∑_i,j = 1^N w_i w_j k_0(y_i,y_j) ) .


  * At iteration n, it holds that D_𝒦_0,P({y_π(i)}_i=1^n) = (*) + O(√(log(n) / n)).

Thus in this scenario the algorithms that we have proposed act to ensure that these points are optimally weighted in the sense just described.







§ EXPERIMENTAL PROTOCOL AND ADDITIONAL NUMERICAL RESULTS
  

This section contains additional numerical results that elaborate on the three experiments reported in the main text.



 §.§ Gaussian Mixture Test
 

Recall from the main text that the kernels k_1, k_2 and k_3 contain either one or two hyper-parameters that must be selected.
For each of the methods (a)-(f) reported in Figure <ref> in the main text we optimised these parameters over a discrete set, with respect to an objective function of W_P based on a point set of size n = 100 and the Nelder-Mead optimisation method. The set of possible values for α was {0.1η, 0.5η, η, 2η, 4η, 8η}, where η is a problem dependent “base scale" and chosen to be 1 for the Gaussian mixture test. The set of possible values for β was {-0.1, -0.3, -0.5, -0.7, -0.9}.
The sensitivity of the reported results to the variation in hyper-parameters is shown, for the Gaussian mixture test, in Figure <ref>.
Point sets obtained under representatives of each method class are shown in Figure <ref>.

For all the global optimisation methods we imposed a bounding box (-5, 5) × (-5,5); for the Nelder-Mead method, we set n_init = 3, n_delay = 20, μ_0 = (0, 0), Σ_0 = 25I, and λ = 1; for the Monte Carlo method, we set n_test = 20, n_delay = 20, μ_0 = (0, 0), Σ_0 = 25I, and λ = 1; for the grid search, we set n_0 = 100.

For MED the tuning parameter δ was considered for δ = 4, δ = 8 or δ = 16, with δ = 4d = 8 being the recommendation in <cit.>.

For SVGD we set the initial point-set to be an equally spaced rectangular grid over the bounding box. Following <cit.>, the step-size ϵ for SVGD was determined by AdaGrad with a master step-size of 0.1 and a momentum factor of 0.9.











 §.§ Gaussian Process Test
 

For the Gaussian process test, the base scale η is also set to 1. The sensitivity of results to the selection of kernel parameters was reported in Figure <ref>.
Point sets obtained under representatives of each method class are shown in Figures <ref> and <ref>.
Detailed results for each method considered are contained in Figure <ref>.

For all the global optimisation methods we imposed a bounding box of (-5, 5) × (-13,-7); for the Nelder-Mead method, we set n_init = 3, n_delay = 20, μ_0 = (0, -10), Σ_0 = 25I, and λ = 1; for the Monte Carlo method, we set n_test = 20, n_delay = 20, μ_0 = (0, -10), Σ_0 = 25I, and λ = 1; for the grid search, we set n_0 = 100.

For SVGD we set the initial point-set to be an equally spaced rectangular grid over the bounding box. Following <cit.>, the step-size ϵ for SVGD was determined by AdaGrad with a master step-size of 0.1 and a momentum factor of 0.9.
















 §.§ IGARCH Test
 

For the IGARCH test, we choose the base scale η to be 1e-5. The sensitivity of results to the selection of kernel parameters was reported in Figure <ref>.
Point sets obtained under representatives of each method class are shown in Figures <ref> and <ref>.
Detailed results for each method considered are contained in Figure <ref>.

For all the global optimisation methods we impose a bounding box of (0.002, 0.04) × (0.05, 0.2); for the Nelder-Mead method, we set n_init = 3, n_delay = 20, μ_0 = (0.021, 0.125), Σ_0 = diag[(1e-4, 1e-3)], and λ = 1e-5; for the Monte Carlo method, we set n_test = 20, n_delay = 20, μ_0 = (0.021, 0.125), Σ_0 = diag[(1e-4, 1e-3)], and λ = 1e-5; for the grid search, we set n_0 = 100.

For SVGD we set the initial point-set to be an equally spaced rectangular grid over the bounding box. Following <cit.>, the step-size ϵ for SVGD was determined by AdaGrad with a master step-size of 1e-3 and a momentum factor of 0.9.


















