Year: 2019

Abstract:
We propose a novel supervised learning method to optimize the kernel in maximum mean discrepancy generative adversarial networks (MMD GANs). Specifically, we characterize a distributionally robust optimization problem to compute a good distribution for the random feature model of Rahimi and Recht to approximate a good kernel function. Due to the fact that the distributional optimization is infinite dimensional, we consider a Monte-Carlo sample average approximation (SAA) to obtain a more tractable finite dimensional optimization problem. We subsequently leverage a particle stochastic gradient descent (SGD) method to solve the derived finite dimensional optimization problem. Based on a mean-field analysis, we then prove that the empirical distribution of the interactive particles system at each iteration of the SGD follows the path of the gradient descent flow on the Wasserstein manifold. We also establish the non-asymptotic consistency of the finite sample estimator. Our empirical evaluation on synthetic data-set as well as MNIST and CIFAR-10 benchmark data-sets indicates that our proposed MMD GAN model with kernel learning indeed attains higher inception scores well as Frèchet inception distances and generates better images compared to the generative moment matching network (GMMN) and MMD GAN with untrained kernels.

heading: 1 INTRODUCTION
text: A fundamental and long-standing problem in unsupervised learning systems is to capture the underlying distribution of data. While deep generative models such as Boltzmann machines Salakhutdinov & Hinton (2009) and auto-encoding variational Bayes Kingma & Welling (2013) accomplish this task to some extent, they are inadequate for many intractable probabilistic computations that arise in maximum likelihood estimation. Moreover, in many machine learning tasks such as caption generation Xu et al. (2015), the main objective is to obtain new samples rather than to accurately estimate the underlying data distribution. Generative adverserial network (GAN) Goodfellow et al. (2014) provides a framework to directly draw new samples without estimating data distribution. It consists of a deep feedforward network to generate new samples from a base distribution (e.g. Gaussian distribution), and a discriminator network to accept or reject the generated samples. However, training GAN requires finding a Nash equilibrium of a non-convex minimax game with continuous, high-dimensional parameters. Consequently, it is highly unstable and prone to miss modes Salimans et al. (2016); Che et al. (2016). To obtain more stable models, the generative moment matching networks (GMMNs) Li et al. (2015) are proposed, wherein instead of training a discriminator network, a non-parametric statistical hypothesis test is performed to accept or reject the generated samples via the computation of the kernel maximum mean discrepancy Gretton et al. (2007). While leveraging a statistical test simplifies the loss function for training GMMN, in practice, the diversity of generated samples by GMMN is highly sensitive to the choice of the kernel. Thus, to improve the sampling performance, the kernel function also needs to be jointly optimized with the generator. Rather than
i
optimizing the kernel directly, the MMD GAN model Li et al. (2017) is proposed in which an embedding function is optimized in conjunction with a fixed user-defined kernel (e.g. RBF Gaussian kernel). However, there are no theoretical guarantees that the user-defined kernel is the ‘right’ kernel for embedded features.
Contributions. To address the kernel model selection problem in MMD GAN Li et al. (2017), in this paper we put forth a novel framework to learn a good kernel function from training data. Our kernel learning approach is based on a distributional optimization problem to learn a good distribution for the random feature model of Rahimi and Recht Rahimi & Recht (2008; 2009) to approximate the kernel. Since optimization with respect to the distribution of random features is infinite dimensional, we consider a Monte Carlo approximation to obtain a more tractable finite dimensional optimization problem with respect to the samples of the distribution. We then use a particle stochastic gradient descent (SGD) to solve the approximated finite dimensional optimization problem. We provide a theoretical guarantee for the consistency of the finite sample-average approximations. Based on a mean-field analysis, we also show the consistency of the proposed particle SGD. In particular, we show that when the number of particles tends to infinity, the empirical distribution of the particles in SGD follows the path of the gradient descent flow of the distributional optimization problem on the Wasserstein manifold.

heading: 2 PRELIMINARIES OF MMD GANS
text: Assume we are given data {vi}ni=1 that are sampled from an unknown distribution PV with the support V . In many unsupervised tasks, we wish to attain new samples from the distribution PV without directly estimating it. Generative Adversarial Network (GAN) Goodfellow et al. (2014) provides such a framework. In vanilla GAN, a deep network G(·;W ) parameterized by W ∈ W is trained as a generator to transform the samples Z ∼ PZ ,Z ∈ Z from a user-defined distribution PZ (e.g. Gaussian) into a new sample G(Z;W ) ∼ PW , such that the distributions PW and PV are close. In addition, a discriminator network D(·; δ) parameterized by δ ∈ ∆ is also trained to reject or accept the generated samples as a realization of the data distribution. The training of the generator and discriminator networks is then accomplished via solving a minimax optimization problem as below
min W∈W max δ∈∆ IEPV [D(X; δ)] + IEPZ [log(1−D(G(Z;W ); δ))]. (1)
In the high dimensional settings, the generator trained via the min-max program of equation 1 can potentially collapse to a single mode of distribution where it always emits the same point Che et al. (2016). To overcome this shortcoming, other adversarial generative models are proposed in the literature, which propose to modify or replace the discriminator network by a statistical two-sample test based on the notion of the maximum mean discrepancy which is defined below:
Definition 2.1. (MAXIMUM MEAN DISCREPANCY GRETTON ET AL. (2007)) Let (X , d) be a metric space, F be a class of functions f : X → IR, and P,Q ∈ B(X ) be two probability measures from the set of all Borel probability measures B(X ) onX . The maximum mean discrepancy (MMD) between the distributions P and Q with respect to the function class F is defined below
DF [P,Q] def = sup
f∈F ∫ X f(x)(P −Q)(dx). (2)
Different choices of the function class F in equation 2 yield different adversarial models such as Wasserstein GANs (WGAN) Arjovsky et al. (2017), f -GANs Nowozin et al. (2016), and GMMN and MMD GAN Li et al. (2017; 2015). In the latter two cases, the function class F def= {f : X → IR : ‖f‖HX ≤ 1}, where HX is a reproducing kernel Hilbert space (RKHS) of functions with a kernel K : X × X → IR, denoted by (HX ,K). Then, the squared MMD loss in equation 2 as a
ii
measure of the distance between the distributions PV and PW has the following expression
DK [PV , PW ] def = sup
f :X→IR:‖f‖HX≤1
∫ X f(x)(PV − PW )(dx) (3)
= IEV ,V ′∼PV [K(V ;V ′)] + IEW ,W ′∼PW [K(W ,W ′)]− 2IEV ∼PV ,W∼PW [K(V ;W )], (4)
where X = V ∪ W . Instead of training the generator via solving the minimax optimization in equation 1, the MMD GAN model of Li et al. (2015) proposes to optimize the discrepancy between two distributions via optimization of an embedding function ι : IRD 7→ IRp, p ≤ D, i.e.,
min W∈W max ι∈Q MMDk◦ι[PV , PW ], (5)
where k : IRp × IRp → IR is a user-defined fixed kernel. In Li et al. (2015), the proposal for the kernel k : IRp × IRp → IR is a mixture of the Gaussians,
k ◦ ι(x,y) = k(ι(x), ι(y)) = m∑ i=1 ( ‖ι(x)− ι(y)‖22 σ2i ) , (6)
where the bandwidth parameters σ1, · · · , σm > 0 are manually selected. Nevertheless, in practice there is no guarantee that the user-defined kernel k(ι(x), q(y)) can capture the structure of the embedded features ι(x).

heading: 3 PROPOSED APPROACH: KERNEL LEARNING WITH RANDOM FEATURES FOR
text: MMD GANS
In this section, we first expound our kernel learning approach. Then, we describe a novel MMD GAN model based on the proposed kernel learning approach.

heading: 3.1 ROBUST DISTRIBUTIONAL OPTIMIZATION FOR KERNEL LEARNING
text: To address the kernel model selection issue in MMD GAN Li et al. (2017), we consider a kernel optimization scheme with random features Rahimi & Recht (2008; 2009). Let ϕ : IRd × IRD → [−1, 1] denotes the explicit feature maps and µ ∈ M(IRD) denotes a probability measure from the space of probability measuresM(IRD) on IRD. The kernel function is characterized via the explicit feature maps using the following integral equation
Kµ(x,y) = IEµ[ϕ(x; ξ)ϕ(y; ξ)] = ∫ Ξ ϕ(x; ξ)ϕ(y; ξ)µ(dξ). (7)
Let MMDµ[PV , PW ] def = MMDKµ [PV , PW ]. Then, the kernel optimization problem in can be formulated as a distribution optimization for random features, i.e,
min W∈W sup µ∈P MMDµ[PV , PW ]. (8)
Here, P is the set of probability distributions corresponding to a kernel class K. In the sequel, we consider P to be the distribution ball of radius R as below
P def= IBpR(µ0) def = {µ ∈M(IRD) : Wp(µ, µ0) ≤ R}, (9)
where µ0 is a user-defined base distribution, and To establish the proof, we consider the Wp(·, ·) is the p-Wasserstein distance defined as below
Wp(µ1, µ2) def = ( inf
π∈Π(µ1,µ2) ∫ IRD×IRD ‖ξ1 − ξ2‖p2dπ(ξ1, ξ2) ) 1 p , (10)
iii
where the infimum is taken with respect to all couplings π of the measures µ, µ0 ∈ M(IRD), and Π(µ, µ0) is the set of all such couplings with the marginals µ and µ0.
The kernel MMD loss function in equation 8 is defined with respect to the unknown distributions of the data-set PV and the model PW . Therefore, we construct an unbiased estimator for the MMD loss function in equation 8 based on the training samples. To describe the estimator, sample the labels from a uniform distribution y1, · · · , yn ∼i.i.d. Uniform{−1,+1}, where we assume that the number of positive and negative labels are balanced. In particular, consider the set of positive labels I = {i ∈ {1, 2, · · · , n} : yi = +1}, and negative labels J = {1, 2, · · · , n}/I, where their cardinality is |I| = |J | = n2 . We consider the following assignment of labels:
• Positive class labels: If yi = +1, sample the corresponding feature map from data-distribution xi = vi ∼ PV .
• Negative class labels: If yi = −1, sample from the corresponding feature map from the generated distribution xi = G(Zi,W ) ∼ PW ,Zi ∼ PZ .
By this construction, the joint distribution of features and labels PY,X has the marginals PX|Y=+1 = PV , and PX|Y=−1 = PW . Moreover, the following statistic, known as the kernel alignment in the literature (see, e.g., Sinha & Duchi (2016); Cortes et al. (2012)), is an unbiased estimator of the MMD loss in equation 8,
min W∈W sup µ∈P M̂MDµ
[ PV , PW ] def =
8 n(n− 1) ∑
1≤i<j≤n
yiyjKµ(xi,xj). (11)
See Appendix C.1 for the related proof. The kernel alignment in equation 11 can also be viewed through the lens of the risk minimization
min W∈W inf µ∈P
M̂MD α
µ
[ PV , PW ] def =
8 n(n− 1)α ∑
1≤i<j≤n
(αyiyj −Kµ(xi,xj))2 (12a)
= 8 n(n− 1)α ∑
1≤i<j≤n
(αyiyj − IEµ[ϕ(xi; ξ)ϕ(xj ; ξ)])2 . (12b)
Here, α > 0 is a scaling factor that determines the separation between feature vectors, and K∗ def = αyyT is the ideal kernel that provides the maximal separation between the feature vectors over the training data-set, i.e., K∗(xi,xj) = α when features have identical labels yi = yj , and K∗(xi,xj) = −α otherwise. Upon expansion of the risk function in equation 12, it can be easily shown that it reduces to the kernel alignment in equation 11 when α → +∞. Intuitively, the risk minimization in equation 12 gives a feature space in which pairwise distances are similar to those in the output space Y = {−1,+1}.

heading: 3.2 SAA FOR DISTRIBUTIONAL OPTIMIZATION
text: The distributional optimization problem in equation 8 is infinite dimensional, and thus cannot be solved directly. To obtain a tractable optimization problem, instead of optimizing with respect to the distribution µ of random features, we optimize the i.i.d. samples (particles) ξ1, · · · , ξN ∼i.i.d. µ generated from the distribution. The empirical distribution of these particles is accordingly defined as follows
µ̂N (ξ) def =
1
N N∑ k=1 δ(ξ − ξk), (13)
where δ(·) is the Dirac’s delta function concentrated at zero. In practice, the optimization problem in equation 12 is solved via the Monte-Carlo sample average approximation of the objective function,
min W∈W min µ̂N∈PN
M̂MD α
µ̂N [ PV , PW ] =
8 n(n− 1)α ∑
1≤i<j≤n
( αyiyj − 1
N N∑ k=1 ϕ(xi; ξ k)ϕ(xj ; ξ k) )2 , (14)
iv
where PN def = IBNR (µ̂ N 0 ) = { µ̂N ∈ M(IRD) : Wp(µ̂N , µ̂N0 ) ≤ R } , and µ̂N0 is the empirical
measure associated with the initial samples ξ10, · · · , ξN0 ∼i.i.d. µ0. The empirical objective function in equation 14 can be optimized with respect to the samples ξ1, · · · , ξN using the particle stochastic gradient descent. For the optimization problem in equation 14, the (projected) stochastic gradient descent (SGD) takes the following recursive form,1
ξkm+1 = ξ k m −
η
N
( ymỹm − 1
αN N∑ k=1 ϕ(xm; ξ k m)ϕ(x̃m; ξ k m)
) ∇ξ ( ϕ(xm; ξ k m)ϕ(x̃m; ξ k m) ) , (15a)
for k = 1, 2, · · · , N , where (ym,xm), (ỹm, x̃m) ∼i.i.d Px,y and η ∈ IR>0 denotes the learning rate of the algorithm, and the initial particles are ξ10, · · · , ξN0 ∼i.i.d. µ0. At each iteration of the SGD dynamic in equation 15, a feasible solution for the inner optimization of the empirical risk function in equation 14 is generated via the empirical measure
µ̂Nm(ξ) = 1
N N∑ k=1 δ(ξ − ξkm). (16)
Indeed, we prove in Section 4 that for an appropriate choice of the learning rate η > 0, the empirical measure in equation 16 remains inside the distribution ball µ̂Nm ∈ PN for all m ∈ [0, NT ] ∩ IN, and is thus a feasible solution for the empirical risk minimization equation 14 (see Corollary 4.2.1 in Section 4).

heading: 3.3 PROPOSED MMD GAN WITH KERNEL LEARNING
text: In Algorithm 1, we describe the proposed method MMD GAN model with the kernel learning approach described earlier. Algorithm 1 has an inner loop for the kernel training and an outer loop for training the generator, where we employ RMSprop Tieleman & Hinton (2012). Our proposed MMD GAN model is distinguished from MMD GAN of Li et al. (2017) in that we learn a good kernel function in equation 17 of the inner loop instead of optimizing the embedding function that is implemented by an auto-encoder. However, we mention that our kernel learning approach is compatible with the auto-encoder implementation of Li et al. (2017) for the dimensionality reduction of features (and particles). In the case of including an auto-encoder, the inner loop in Algorithm 1 must be modified to add an additional step for training the auto-encoder. However, to convey the main ideas more clearly, the training step of the auto-encoder is omitted from Algorithm 1.

heading: 3.4 COMPUTATIONAL COMPLEXITY ANALYSIS
text: Sampling the labels y, ỹ ∼i.i.d. Uniform{−1,+1} require O(1) complexity, while sampling x̃|ỹ = +1,x|y = +1 ∼ PV and x̃|ỹ = +1,x|y = −1 ∼ PW has a complexity of O(d). The computation of the stochastic gradient step in equation 17 requires computing the sum 1 N ∑N i=1 ϕ(x; ξ k)ϕ(x̃; ξk). This can be done as a separate preprocessing step prior to executing the SGD in equation 17, and requires preparing the vectors ϕ def= [ϕ(x; ξ1), · · · , ϕ(x; ξN )] and ϕ̃ def = [ϕ(x̃; ξ1), · · · , ϕ(x̃; ξN )]. Using the random feature model of Rahimi and Recht Rahimi & Recht (2008; 2009), where ϕ(x; ξ) = √
2 cos(xT ξ+ b). Here b ∼ Uniform[−π,+π], the complexity of computing the vectors ϕ and ϕ̃ is of the order O(Nd) on a single processor. However, this construction is trivially parallelizable. Furthermore, computation can be sped up even further for certain distributions µ0. For example, the Fastfood technique can approximateϕ and ϕ̃ inO(N log(d))
1To avoid clutter in our subsequent analysis, the normalization factor 16 n(n−1) of the gradient is omitted by
modifying the step-size η.
v
Algorithm 1 MMD GAN with a supervised kernel learning Method (Monte-Carlo Approach) Inputs: The learning rates η̃, η > 0 , the number of iterations of discriminator per generator update T ∈ N, the batch-size n, the number of random featuresN ∈ N. Regularization parameter α > 0. whileW has not converged do
for t = 1, 2, · · · , T do Sample the labels y, ỹ ∼i.i.d Uniform{−1, 1}. Sample the features x|y = +1 ∼ PV , and x|y = −1 ∼ PW . Similarly, x̃|ỹ = +1 ∼ PV , and x̃|ỹ = −1 ∼ PW . For all k = 1, 2, · · · , N , update the particles,
ξ k ← ξk −
η
N
( αyỹ − 1
N N∑ k=1 ϕ(x; ξ k )ϕ(x̃; ξ k )
) ∇ξ ( ϕ(x; ξ k )ϕ(x̃; ξ k ) ) , (17)
end for Sample a balanced minibatch of labels {yi}ni=1 ∼i.i.d. Uniform{−1,+1}. Sample the minibatch {x}ni=1 such that xi|yi = +1 ∼ PV , and xi|yi = −1 ∼ PW for all i = 1, 2, · · · , n. Update the generator
gW ← ∇W D̂αµ̂N [ PV , PW ] , µ̂ N = 1
N N∑ k=1 δ(ξ − ξk). (18a)
W ←W − η̃RMSprop(gW ,W ). (18b)
end while
time for the Gaussian kernel Le et al. (2013). Updating each particle in equation 17, involves the computation of the gradient ∇ξ(ϕ(ξ; ξk)ϕ(ξ̃; ξk))) which is O(d). Thus, the complexity of one iteration of SGD for all the particles is O(Nd). Overall, one step of the kernel learning has a complexity of O(Nd). On the other hand, to attain ε-suboptimal solution maxk=1,2,··· ,N ‖ξk − ξk∗‖ ≤ ε, the SGD requires has the sample complexity O ( log( 1ε ) ) . Consequently, the computational complexity of the kernel learning is of the order of
O(Nd log( 1ε )). To compare this complexity with that of MMD GAN is of the orderO ( B2` log( 1ε ) ) , where B is the batch size for approximation of the population MMD, and ` is the number of kernel mixtures.

heading: 4 CONSISTENCY AND A MEAN-FIELD ANALYSIS
text: In this section, we provide theoretical guarantees for the consistency of various approximations we made to optimize the population MMD loss function in equation 8. We defer the proofs of the following theoretical results to AppendixC. The main assumptions ((A.1),(A.2), and (A.3)) underlying our theoretical results are also stated in the same section.
Consistency of finite-sample estimate: In this part, we prove that the solution to finite sample optimization problem in equation 14 approaches its population optimum in equation 8 as the number of data points as well as the number of random feature samples tends to infinity.
Theorem 4.1. (NON-ASYMPTOTIC CONSISTENCY OF FINITE-SAMPLE ESTIMATOR) Suppose conditions (A.1)-(A.3) of Appendix C are satisfied. Consider the distribution balls P and PN that are defined with respect to the 2-Wasserstein distance (p = 2). Furthermore, consider the optimal MMD values of the population optimization and its finite sample estimate
(W∗, µ∗) def = arg min
W∈W arg sup µ∈P MMDµ[PV , PW ]. (19a)
(ŴN∗ , µ̂ N ∗ )
def = arg min
W∈W arg inf µ̂N∈PN M̂MD
α µ̂N [PV , PW ], (19b)
vi
respectively. Then, with the probability of (at least) 1 − 3% over the training data samples {(xi, yi)}ni=1 and the random feature samples {ξk0}Nk=1, the following non-asymptotic bound holds∣∣∣MMDµ∗ [PV , PW∗ ]−MMDµ̂N∗ [PV , PŴN∗ ]∣∣∣ (20) ≤ √ L2(d+ 2)
N ln
1 2
( 28Ndiam2(X )
%
) + 2 max { c1L 2
n ln
1 2
( 4
%
) , c2RL 4
n2 ln
( 4e L4 9
%
)} + 8L2
α ,
where c1 = 3 1 4 × 24, and c2 = 9× 211.
The proof of Theorem 4.1 is presented in Appendix C.1.
Notice that there are three key parameters involved in the upper bound of Theorem 4.1. Namely, the number of training samples n, the number of random feature samples N , and the regularization parameter α. The upper bound in equation 20 thus shows that when n,N, α → +∞, the solution obtained from solving the empirical risk minimization in equation 12 yields a MMD population value tending to the optimal value of the distributional optimization in equation 8.
Consistency of particle SGD for solving distributional optimization. The consistency result of Theorem 4.1 is concerned with the MMD value of the optimal empirical measure µ̂N∗ (ξ) = 1 N ∑N k=1 δ(ξ − ξk∗) of the empirical risk minimization equation 14. In practice, the particle SGD is executed for a few iterations and its values are used as an estimate for (ξ1∗, · · · , ξN∗ ). Consequently, it is desirable to establish a similar consistency type result for the particle SGD estimates (ξ1m, · · · , ξNm) at the m-th iteration. To reach this objective, we define the scaled empirical measure as follows
µNt = µ̂ N bNtc =
1
N N∑ k=1 δ(ξ − ξbNtc), 0 ≤ t ≤ T. (21)
At any time t, the scaled empirical measure µNt is a random element, and thus (µ N t )0≤t≤T is a measured-valued stochastic process. Therefore, we characterize the evolution of its Lebesgue density pNt (ξ) def = µNt (dξ)/dξ in the following theorem: Theorem 4.2. (MCKEAN-VLASOV MEAN-FIELD PDE) Suppose conditions (A.1)-(A.3) of Appendix C are satisfied. Further, suppose that the Radon-Nikodyme derivative q0(ξ) = µ0(dξ)/dξ exists. Then, there exists a unique solution (p∗t (ξ))0≤t≤T to the following non-linear partial differential equation∂pt(ξ)∂t = − ηα ∫∫ X×Y (∫ IRp ϕ(x, ξ̃)ϕ(x̃, ξ̃)pt(ξ̃)dξ̃ − αyỹ ) ∇ξ(pt(ξ)∇ξ(ϕ(x; ξ)ϕ(x̃; ξ))dP⊗2X,Y ,
p0(ξ) = q0(ξ).
(22)
Moreover, the measure-valued process {(µNt )0≤t≤T }N∈IN defined in equation 21 converges (weakly) to the unique solution µ∗t (ξ) = p ∗ t (ξ)dξ as the number of particles tend to infinityN →∞. 2.
Due to the mean-field analysis of Theorem 4.2, we can prove that the empirical measure µ̂Nm of the particles in SGD dynamic equation 15 remains inside the feasible distribution ball PN :
Corollary 4.2.1. Consider the learning rate η = O (
Rp
T √ NT log(2/δ)
) for the SGD in equation 15.
Then, the empirical measure µ̂Nm of the particles remains inside the distributional ball µ̂ N m ∈ PN = {µ̂N ∈ M(IRD) : Wp(µ̂N , µ̂N0 ) ≤ R} for all m ∈ [0, NT ] ∩ IN, with the probability of (at least) 1− δ.
2The notion of the weak convergence of a sequence of empirical measures is formally defined in Supplementary.
vii
Let us make two remarks about the PDE in equation 22. First, the seminal works of Otto Otto (2001), and Jordan, et al. Jordan et al. (1998) establishes a deep connection between the McKeanVlasov type PDEs specified in equation 22 and the gradient flow on the Wasserstein manifolds. More specifically, the PDE equation in equation 22 can be thought of as the minimization of the energy functional
inf µ∈M(IRD)
Eα(pt(ξ)) def =
1
α ∫ IRD Rα(ξ, pt(ξ))pt(ξ)dξ (23a)
Rα(ξ, pt(ξ)) def = −α(IEPX,Y [yϕ(x; ξ)])2 + IEξ̃∼pt [( IEPX [ϕ(x; ξ)ϕ(x; ξ̃)] )2] , (23b)
using the following gradient flow dynamics dpt(ξ)
dt = −η · gradptEα(pt(ξ)), p0(ξ) = q0(ξ), (24)
where gradptE(pt(ξ)) = ∇ξ · (pt(ξ)∇ξRα(pt(ξ))) is the Riemannian gradient of Rα(µt(ξ)) with respect to the metric of the Wasserstein manifold . This shows that when the number of particles in particle SGD equation 15 tends to infinity (N → +∞), their empirical distribution follows a gradient descent path for minimization of the population version (with respect to data samples) of the distributional risk optimization in equation 12. In this sense, the particle SGD is the ‘consistent’ approximation algorithm for solving the distributional optimization.

heading: 4.1 RELATED WORKS
text: The mean-field description of SGD dynamics has been studied in several prior works for different information processing tasks. Wang et al. Wang et al. (2017) consider the problem of online learning for the principal component analysis (PCA), and analyze the scaling limits of different online learning algorithms based on the notion of finite exchangeability. In their seminal papers, Montanari and co-authors Mei et al. (2018); Javanmard et al. (2019); Mei et al. (2019) consider the scaling limits of SGD for training a two-layer neural network, and characterize the related Mckean-Vlasov PDE for the limiting distribution of the empirical measure associated with the weights of the input layer. Similar mean-field type results for two-layer neural networks are also studied recently in Rotskoff & Vanden-Eijnden (2018); Sirignano & Spiliopoulos (2018). Our work is also related to the unpublished work of Wang, et al. Wang et al., which proposes a solvable model of GAN and analyzes the scaling limits. However, our GAN model is different from Wang et al. and is based on the notion of the kernel MMD. Our work is also closely related to the recent work of Li, et al Li et al. (2019) which proposes an implicit kernel learning method based on the following kernel definition
Kh(ι(x), ι(y)) = IEξ∼µ0
[ e(ih(ξ)(ι(x)−ι(y))) ] , (25)
where µ0 is a user defined base distribution, and h ∈ H is a function that transforms the base distribution µ0 into a distribution µ that provides a better kernel. Therefore, the work of Li, et al Li et al. (2019) implicitly optimizes the distribution of random features via transforming a random variable with a function. In contrast, the proposed distributional optimization framework in this paper optimizes the distribution of random feature explicitly, via optimizing their empirical measures. Perhaps more importantly from a practical perspective is the fact that our kernel learning approach does not require the user-defined function class H. Moreover, our particle SGD method in equation 15 obviates tuning hyper-parameters related to the implicit kernel learning method such as the gradient penalty factor and the variance constraint factor (denoted by λGP and λh, respectively, in Algorithm 1 of Li et al. (2019)).

heading: 5 EMPIRICAL EVALUATION

heading: 5.1 SYNTHETIC DATA-SET
text: Due to the space limitation, the experiments on the synthetic data are deferred to Appendix A.
viii

heading: 5.2 PERFORMANCE ON BENCHMARK DATASETS
text: We evaluate our kernel learning approach on large-scale benchmark data-sets. We train our MMD GAN model on two distinct types of datasets, namely on MNIST LeCun et al. (1998) and CIFAR-10 LeCun et al. (1998), where the size of training instances are 60 × 103 and 50 × 103, respectively. All the generated samples are from a fixed noise random vectors and are not singled out. Implementation and hyper-parameters. We implement Algorithm 1 as well as MMD GAN Li et al. (2017) in Pytorch using NVIDIA Titan V100 32GB graphics processing units (GPUs). The source code of Algorithm 1 is built upon the code of Li et al. (2017), and retains the auto-encoder implementation. In particular, we use a sequential training of the auto-encoder and kernel as explained in the Synthetic data in Section A of Supplementary. For a fair comparison, our hyper-parameters are adjusted as in Li et al. (2017), i.e., the learning rate of 0.00005 is considered for RMSProp Tieleman & Hinton (2012). Moreover, the batch-size for training the generator and auto-encoder is n = 64. The learning rate of particle SGD is tuned to η = 10. Random Feature Maps. To approximate the kernel, we use the the random feature model of Rahimi and Recht Rahimi & Recht (2008; 2009), where ϕ(x; ξ) = √ 2 cos(xT ξ + b). Here b ∼ Uniform{−π,+π} is a random bias term. Practical considerations. When data-samples {Vi} ∈ IRd are high dimensional (as in CIFAR10), the particles ξ1, · · · , ξN ∈ IRD, D = d in SGD equation 15 are also high-dimensional. To reduce the dimensionality of the particles, we apply an auto-encoder architecture similar to Li et al. (2017), and train our kernel on top of learned embedded features. More specifically, in our simulations, we train an auto-encoder where the dimensionality of the latent space is h = 10 for MNIST, and h = 128 (thus D = d = 128) for CIFAR-10. Therefore, the particles ξ1, · · · , ξN in subsequent kernel training phase have the dimension of D = 10, and D = 128, respectively. Choice of the scaling parameter α. There is a trade-off in the choice of α. While for large values of α, the kernel is better able to separate data-samples from generated samples, in practice, it slows down the convergence of particle SGD. This is due to the fact that the coupling between the particle dynamics in equation 15 decrease as α increase. The scaling factor is set to be α = 1 in all the following experiments. Qualitative comparison. We now show that without the bandwidth tuning for Gaussian kernels and using the particle SGD to learn the kernel, we can attain better visual results on benchmark data-sets. In Figure 1, we show the generated samples on CIFAR-10 and MNIST data-sets, using our Algorithm 1, MMD GAN Li et al. (2017) with a mixed and homogeneous Gaussian RBF kernels, and GMMN Li et al. (2015).
Figure 1(a) shows the samples from Algorithm 1, Figure 1(b) shows the samples from MMD GAN Li et al. (2017) with a mixture RBF Gaussian kernel κ(x,y) = ∑5 k=1 κσk(x,y), where σk ∈ {1, 2, 4, 8, 16} are the bandwidths of the Gaussian kernels that are fine tuned and optimized. We observe that our MMD GAN with automatic kernel learning visually attains similar results to MMD GAN Li et al. (2017) which requires manual tuning of the hyper-parameters. In Figure 1(c), we show the MMD GAN result with a single kernel RBF Gaussian kernel whose bandwidth is manually tuned at σ = 16. Lastly, in Figure 1(d), we show the samples from GMMN Li et al. (2015) which does not exploit an auto-encoder or kernel training. Clearly, GMMN yield a poor results compared to other methods due to high dimensionality of features, as well as the lack of an efficient method to train the kernel.
On MNIST data-set in Figure 1,(e)-(h), the difference between our method and MMD GAN Li et al. (2017) is visually more pronounced. We observe that without a manual tuning of the kernel bandwidth and by using the particle SGD equation 15 to optimize the kernel, we attain better generated images in Figure 1(e), compared to MMD GAN with mixed RBF Gaussian kernel and manual bandwidth tuning in Figure 1(f). Moreover, using a single RBF Gaussian kernel yields a poor result regardless of the choice of its bandwidth. The generated images from GMMN is also shown in Figure 1(h).
ix
Quantitivative comparison. To quantitatively measure the quality and diversity of generated samples, we compute the inception score (IS) Salimans et al. (2016) as well as Frèchet Inception Distance (FID) Heusel et al. (2017) on CIFAR-10 images. Intuitively, the inception score is used for GANs to measure samples quality and diversity. Accordingly, for generative models that are collapsed into a single mode of distribution, the inception score is relatively low. The FID improves on IS by actually comparing the statistics of generated samples to real samples, instead of evaluating generated samples independently.
In Table 1, we report the quantitative measures for different MMD GAN model using different scoring metric. Note that in Table 1 lower FID scores and higher IS scores indicate a better performance. We observe from Table 1 that our approach attain lower FID score, and higher IS score compared to MMD GAN with single Gaussian kernel (bandwidth σ = 16), and a mixture Gaussian kernel (bandwidths {1, 2, 4, 8, 16}).
x

heading: APPENDIX
text: We provide additional material to support the content presented in the paper. This Appendix is organized as follows:
Appendix A. We provide the results of our experiments on the synthetic data-set described in Section 5.1 of the main text.
Appendix B. We provide the results of our experiments on the LSUN and CelebA data-sets.
Appendix C. We present the proofs of the main theoretical results of Section 4 in the main text.
Appendix D. We present the proofs of auxiliary lemmas used to support the proof of main results.
Appendix E. We prove additional theoretical results regarding the so-called chaoticity of the particle SGD in equation 15.

heading: A EXPERIMENTAL RESULTS ON THE SYNTHETIC DATA-SET
text: The synthetic data-set we consider is as follows:
• The distribution of training data is PV = N(0, (1 + λ)Id×d), • The distribution of generated data is PW = N(0, (1− λ)Id×d).
xiii
To reduce the dimensionality of data, we consider the embedding ι : IRd 7→ IRp,x 7→ ι(x) = Σx, where Σ ∈ IRp×d and p < d. In this case, the distribution of the embedded features are PX|Y=+1 = N(0, (1 + λ)ΣΣT ), and PX|Y=−1 = N(0, (1− λ)ΣΣT ).
Note that λ ∈ [0, 1] is a parameter that determines the separation of distributions. In particular, the Kullback-Leibler divergece of the two multi-variate Gaussian distributions is controlled by λ ∈ [0, 1],
DKL(PX|Y=−1, PX|Y=+1) = 1
2
[ log ( 1− λ 1 + λ ) − p+ p(1− λ2) ] . (26)
In Figure 2, we show the distributions of i.i.d. samples from the distributions PV and PW for different choices of variance parameter of λ = 0.1, λ = 0.5, and λ = 0.9. Notice that for larger λ the divergence is reduced and thus performing the two-sample test is more difficult. From Figure 2, we clearly observe that for large values of λ, the data-points from the two distributions PV and PW have a large overlap and conducting a statistical test to distinguish between these two distributions is more challenging.

heading: A.0.1 KERNEL LEARNING APPROACH
text: Figure 4 depicts our two-phase kernel learning procedure which we also employed in our implementations of Algorithm 1 on benchmark data-sets of Section 5.2 in the main text. The kernel learning approach consists of training the auto-encoder and the kernel optimization sequentially, i.e.,
sup µ̂N∈PN sup ι∈Q
M̂MD α
Kµ̂N ◦ι [PV , PW ]. (27)
where the function class is defined Q def= {ι(z) = Σz,Σ ∈ IRD×d}, and (Kµ̂N ◦ ι)(x1,x2) = Kµ̂N (ι(x1), ι(x2)). Now, we consider a two-phase optimization procedure:
• Phase (I): we fix the kernel function, and optimize the auto-encoder to compute a covariance matrix Σ for dimensionality reduction • Phase (II): we optimize the kernel based on the learned embedded features.
This two-phase procedure significantly improves the computational complexity of SGD as it reduces the dimensionality of random feature samples ξ ∈ IRD, D d. When the kernel function K is fixed, optimizing the auto-encoder is equivalent to the kernel learning step of Li et al. (2017).
xiv

heading: A.0.2 STATISTICAL HYPOTHESIS TESTING WITH THE KERNEL MMD
text: Let V1, · · · ,Vm ∼i.i.d. PV = N(0, (1 + λ)Id×d), and W1, · · · ,Wn ∼i.i.d. PW = N(0, (1 − λ)Id×d). Given these i.i.d. samples, the statistical test T ({Vi}mi=1, {Wi}nj=1) : Vm×Wn → {0, 1} is used to distinguish between these hypotheses:
• Null hypothesis H0 : PV = PW (thus λ = 0), • Alternative hypothesis H1 : PV 6= PW (thus λ > 0).
To perform hypothesis testing via the kernel MMD, we require thatHX is a universal RKHS, defined on a compact metric space X . Universality requires that the kernel K(·, ·) be continuous and, HX be dense in C(X ). Under these conditions, the following theorem establishes that the kernel MMD is indeed a metric: Theorem A.1. (METRIZABLITY OF THE RKHS) Let F denotes a unit ball in a universal RKHS HX defined on a compact metric space X with the associated continuous kernel K(·, ·). Then, the kernel MMD is a metric in the sense that MMDK [PV , PW ] = 0 if and only if PV = PW .
To design a test, let µ̂Nm(ξ) = 1 N ∑N k=1 δ(ξ − ξkm) denotes the solution of SGD in equation 15 for solving the optimization problem. Consider the following MMD estimator consisting of two
xv
U -statistics and an empirical function
M̂MDKµ̂Nm◦ι [ {Vi}mi=1, {Wi}ni=1 ] =
1
m(m− 1) N∑ k=1 ∑ i 6=j ϕ(ι(Vi), ξ k m)ϕ(ι(Vj), ξ k m)
+ 1
n(n− 1) N∑ k=1 ∑ i 6=j ϕ(ι(Wi), ξ k m)ϕ(ι(Wj), ξ k m)
− 2 nm N∑ k=1 m∑ i=1 n∑ j=1 ϕ(ι(Wi), ξ k m)ϕ(ι(Vj), ξ k m). (28)
Given the samples {Vi}mi=1 and {Wi}ni=1, we design a test statistic as below
T ({Vi}mi=1, {Wi}ni=1) def =
{ H0 if M̂MDKµ̂Nm◦ι [ {Vi}mi=1, {Wi}ni=1 ] ≤ τ
H1 if M̂MDKµ̂Nm◦ι [ {Vi}mi=1, {Wi}ni=1 ] > τ, . (29)
where τ ∈ IR is a threshold. Notice that the unbiased MMD estimator of equation 28 can be negative despite the fact that the population MMD is non-negative. Consequently, negative values for the statistical threshold τ equation 29 are admissible. In the following simulations, we only consider non-negative values for the threshold τ .
A Type I error is made when H0 is rejected based on the observed samples, despite the null hypothesis having generated the data. Conversely, a Type II error occurs when H0 is accepted despite the alternative hypothesis H1 being true. The significance level of a test is an upper bound on the probability of a Type I error: this is a design parameter of the test which must be set in advance, and is used to determine the threshold to which we compare the test statistic. The power of a test is the probability of rejecting the null hypothesis H0 when it is indeed incorrect. In particular,
Power def = IP(reject H0|H1 is true). (30)
In this sense, the statistical power controls the probability of making Type II errors.

heading: A.0.3 EMPIRICAL RESULTS
text: In Figure 3, we show the evolution of the empirical measure µNm(ξ) of SGD particles by plotting the 2D histogram of the particles ξ1m, · · · , ξNm ∈ IR
D at different iterations of SGD (D = d). Clearly, starting with a uniform distribution in 3(a), the empirical measure seemingly evolves into a Gaussian measure in Figure 3(d). The evolution to a Gaussian distribution demonstrates that the RBF Gaussian kernel corresponding to a Gaussian distribution for the random features indeed provides a good kernel function for the underlying hypothesis test with Gaussian distributions.
In Figure 4, we evaluate the power of the test for 100 trials of hypothesis test using the test statistics of equation 29. To obtain the result, we used an autoencoder to reduce the dimension from d = 100 to p = 50. Clearly, for the trained kernel in Panel (a) of Figure 4, the threshold τ for which Power = 1 increases after learning the kernel via the two phase procedure described earlier. In comparison, in Panel (b), we observe that training an auto-encoder only with a fixed standard Gaussian kernel K(x,y) = e−‖x−y‖ 2 2 attains lower thresholds compared to our two-phase procedure. In Panel (c), we demonstrate the case of a fixed Gaussian kernel without an auto-encoder. In this case, the threshold is significantly lower due to the large dimensionality of the data.
From Figure 4, we also observe that interestingly, the phase transition in the statistical threshold τ is not sensitive to the parameter λ. This phenomenon can be justified by the fact that the kernel learning indeed improves the MMD value more significantly for smaller values of λ (i.e., more difficult hypothesis testing problems) than larger values of λ. See Figure 5.
xvi

heading: B EXPERIMENTAL RESULTS ON LSUN BEDROOM AND CELEBA DATA-SETS
text: In this section, we present additional simulations using CelebA Liu et al. (2015), and LSUN bedroom Yu et al. (2015) data-sets. The LSUN dataset of bedroom pictures resized to 64× 64, and the CelebA dataset of celebrity face images resized and cropped to 160 × 160. The sample generated images are shown in Figure 6.
The inception score for LSUN bedroom data set is 3.860 ± 0.0423 using a single Gaussian kernel with the bandwidth of σ = 16, and 4.064± 0.061 using our proposed method in Algorithm 1.

heading: C PROOFS OF MAIN THEORETICAL RESULTS
text: Before we delve into proofs, we state the main assumptions underlying our theoretical results:
Assumptions:
(A.1) The feature space X = V ∪ W ⊂ IRd is compact with a finite diameter diam(X ) < ∞, where V = support(PV ) and W = support(PW ) are the supports of the distributions PV and PW respectively.
xvii
(A.2) The feature maps are bounded and Lipchitz almost everywhere (a.e.) ξ ∈ IRD. In particular, supx∈X |ϕ(x; ξ)| < L0, supx∈X ‖∇ξϕ(x; ξ)‖2 ≤ L1, and supξ∈IRD ‖∇xϕ(x; ξ)‖ < L2. Let L def = max{L0, L1, L2} < +∞.
(A.3) Let µ̂N0 (ξ) def = 1N ∑N k=1 δ(ξ − ξk0 ) denotes the empirical measure for the initial particles
ξ10, · · · , ξN0 . We assume that µ̂N0 (ξ) converges (weakly) to a deterministic measure µ0 ∈ M(IRD). Furthermore, we assume the limiting measure µ0 is absolutely continuous w.r.t. Lebesgue measure and has a compact support support(µ0) = Ξ ⊂ IRD.
Notation: We denote vectors by lower case bold letters, e.g. x = (x1, · · · , xn) ∈ IRn, and matrices by the upper case bold letters, e.g., M = [Mij ] ∈ IRn×m. The Frobenius norm of a matrix is denoted by ‖M‖F = ∑n i=1 ∑m j=1 |Mij |2. Let IBr(x) def = {y ∈ IRd : ‖y − x‖2 ≤ r} denote the Euclidean ball of radius r centered at x. For a given metric space X , Let Cb(IRd) denote the space of bounded and continuous functions on X equipped with the usual supremum norm
‖f‖∞ def = sup x∈X |f(x)|. (31)
Further, Ckb (X ) the space of all functions in Cb(X ) whose partial derivatives up to order k are bounded and continuous, and Ckc (X ) the space of functions whose partial derivatives up to order k are continuous with compact support.
xviii
We denote the class of the integrable functions f with f(t) ≥ 0 a.e., on 0 ≤ t ≤ T by L1+[0, T ]. Similarly, L∞+ [0, T ] will denote the essentially bounded functions with f(t) ≥ 0 almost everywhere. For a given metric space X , we denote the Borel σ-algebra by B(X ). For a Borel set B ∈ B(X ), the measure value of the set B with respect to the measure is given by µ(B). The space of finite non-negative measures defined on X is denoted byM(X ). The Dirac measure with the unit mass at x ∈ X is denoted by δ(x). For any measure µ ∈ M(X ) and any bounded function f ∈ Cb(X ), we define
〈µ, f〉 def= ∫ X f(x)µ(dx). (32)
The space M(X ) is equipped with the weak topology, i.e., a (random) sequence {µN}N∈IN converges weakly to a deterministic measure µ ∈ M(X ) if and only if 〈µN , f〉 → 〈µ, f〉 for all f ∈ Cb(X ). We denote the weak convergence by µNt
weakly→ µ. Notice that when X is Polish, thenM(X ) equipped with the weak topology is also Polish.3 For a Polish space X , let DX ([0, T ]) denotes the Skorokhod space of the cádlág functions that take values in X defined on [0, T ]. We assume that DX ([0, T ]) is equipped with the Skorokhod’s J1-topology Billingsley (2013), which in that case DX ([0, T ]) is also a Polish space. We use asymptotic notations throughout the paper. We use the standard asymptotic notation for sequences. If an and bn are positive sequences, then an = O(bn) means that lim supn→∞ an/bn < ∞, whereas an = Ω(bn) means that lim infn→∞ an/bn > 0. Furthermore, an = Õ(bn) implies an = O(bnpoly log(bn)). Moreover an = o(bn) means that limn→∞ an/bn = 0 and an = ω(bn) means that limn→∞ an/bn = ∞. Lastly, we have an = Θ(bn) if an = O(bn) and an = Ω(bn). Finally, for positive a, b > 0, denote a . b if a/b is at most some universal constant. Definition C.1. (ORLICZ NORM) The Young-Orlicz modulus is a convex non-decreasing function ψ : IR+ → IR+ such that ψ(0) = 0 and ψ(x)→∞ when x→∞. Accordingly, the Orlicz norm of an integrable random variable X with respect to the modulus ψ is defined as
‖X‖ψ def = inf{β > 0 : IE[ψ(||X| − IE[|X|]|/β)] ≤ 1}. (33)
In the sequel, we consider the Orlicz modulus ψν(x) def = (xν)− 1 . Accordingly, the cases of ‖ · ‖ψ2 and ‖ ·‖ψ1 norms are called the sub-Gaussian and the sub-exponential norms and have the following alternative definitions: Definition C.2. (SUB-GAUSSIAN NORM) The sub-Gaussian norm of a random variable Z, denoted by ‖Z‖ψ2 , is defined as
‖Z‖ψ2 = sup q≥1
q−1/2(IE|Z|q)1/q. (34)
For a random vector Z ∈ IRn, its sub-Gaussian norm is defined as follows
‖Z‖ψ2 = sup x∈Sn−1 ‖〈x,Z〉‖ψ2 . (35)
Definition C.3. (SUB-EXPONENTIAL NORM) The sub-exponential norm of a random variable Z, denoted by ‖Z‖ψ1 , is defined as follows
‖Z‖ψ1 = sup q≥1
q−1(IE[|Z|q])1/q. (36)
For a random vector Z ∈ IRn, its sub-exponential norm is defined below
‖Z‖ψ1 = sup x∈Sn−1 ‖〈Z,x〉‖ψ1 . (37)
3A topological space is Polish if it is homeomorphic to a complete, separable metric space.
xix

heading: C.1 PROOF OF THEOREM 4.1
text: By the triangle inequality, we have that∣∣∣MMDµ∗ [PV , PW∗ ]−MMDµ̂N∗ [PV , PŴN∗ ]∣∣∣ ≤ A1 + A2 + A3 + A4, (38) where the terms Ai, i = 1, 2, 3, 4 are defined as follows
A1 def = ∣∣∣MMDµ∗ [PV , PW∗ ]− min
W∈W sup µ∈P
M̂MDµ[PV , PW ] ∣∣∣
A2 def = ∣∣∣ min W∈W sup µ∈P M̂MDµ[PV , PW ]− min W∈W sup µ̂N∈PN M̂MDµ̂N [ PV , PW ]∣∣∣ A3 def = ∣∣∣ min W∈W
sup µ̂N∈PN
M̂MDµ̂N [ PV , PW ] − M̂MDµ̂N∗ [ PV , PW ]∣∣∣ A4 def = ∣∣∣M̂MDµ̂N∗ [PV , PŴN∗ ]−MMDµ̂N∗ [PV , PŴN∗ ]∣∣∣.
In the sequel, we compute an upper bound for each term on the right hand side of equation 38:
Upper bound on A1:
First, notice that the squared kernel MMD loss in equation 4 can be characterized in terms of class labels and features defined in Section 3.1 as follows
MMDµ[PV , PW ] = 4IEP⊗2x,y [yŷKµ(x, x̂)] . (39)
To see this equivalence, we first rewrite the right hand side of equation 39 as follows IEP⊗2y,x [yŷKµ(x, x̂)] = IP{y = +1}IP{ŷ = +1}IEx,x̂∼P⊗2x|y=+1 [Kµ(x, x̂)]
+ IP{y = −1}IP{ŷ = −1}IEx,x̂∼P⊗2 x|y=−1 [Kµ(x, x̂)]
− IP{y = −1}IP{ŷ = +1}IEx∼Px|y=−1,x̂∼Px|y=+1 [Kµ(x, x̂)] − IP{y = +1}IP{ŷ = −1}IEx∼Px|y=+1,x̂∼Px|y=−1 [Kµ(x, x̂)]. (40)
Now, recall from Section 3.1 that Px|y=+1 = PV , and Px|y=−1 = PW by construction of the labels and random features. Moreover, y, ŷ ∼i.i.d. Uniform{−1,+1}, and thus IP{y = −1} = IP{y = +1} = 12 . Therefore, from equation 40, we derive
IEP⊗2y,x [yŷKµ(x, x̂)] = 1
4 IEP⊗2V
[Kµ(x; x̂)] + 1
4 IEP⊗2W
[Kµ(x; x̂)]− 1
2 IEPV ,PW [Kµ(x; x̂)]
= 1
4 MMDµ[PV , PW ].
For any givenW ∈ W , we have that∣∣∣ sup µ∈P M̂MDµ[PV , PW ]− sup µ∈P MMDµ[PV , PW ] ∣∣∣
≤ sup µ∈P ∣∣M̂MDµ[PV , PW ]−MMDµ[PV , PW ]∣∣ = 4 sup
µ∈P ∣∣∣∣∣ 1n(n− 1) ∑ i 6=j yiyjKµ(xi,xj)− IEP⊗2y,x [yŷKµ(x, x̂)] ∣∣∣∣∣ = 4 sup
µ∈P ∣∣IEµ[En(ξ)]∣∣ ≤ 4
∣∣∣∣sup µ∈P IEµ[En(ξ)] ∣∣∣∣ , xx
where the error term is defined using the random features
En(ξ) def =
1 n(n− 1) ∑ i 6=j yiyjϕ(xi; ξ)ϕ(xj ; ξ)− IEP⊗2x,y [yŷϕ(x; ξ)ϕ(x̂, ξ)]. (41)
Now, we invoke the following strong duality theorem Gao & Kleywegt (2016): Theorem C.4. (STRONG DUALITY FOR ROBUST OPTIMIZATION, (GAO & KLEYWEGT, 2016, THEOREM 1)) Consider the general metric space (Ξ, d), and any normal distribution ν ∈ M(Ξ), whereM(Ξ) is the set of Borel probability measures on Ξ. Then,
sup µ∈M(Ξ)
{ IEµ[Ψ(ξ)] : Wp(µ, ν) ≤ R } = min
λ≥0
{ λRp − ∫ Ξ inf ξ∈Ξ [λdp(ξ, ζ)−Ψ(ξ)]ν(dζ) } , (42)
provided that Ψ is upper semi-continuous in ξ.
Under the strong duality of Theorem C.4, we obtain that∣∣∣ sup µ∈P M̂MDµ[PV , PW ]− sup µ∈P MMDµ[PV , PW ] ∣∣∣
≤ 4 ∣∣∣∣minλ≥0 { λRp − ∫ IRD inf ζ∈IRD [ λ‖ξ − ζ‖p2 − En(ζ) ] µ0(dξ) }∣∣∣∣ . (43) In the sequel, let p = 2. The Moreau’s envelope Parikh & Boyd (2014) of a function f : X → IR is defined as follows
Mβf (y) def = inf x∈X
{ 1
2β ‖x− y‖22 + f(x)
} , ∀y ∈ X , (44)
where β > 0 is the regularization parameter. When the function f is differentiable, the following lemma can be established: Lemma C.5. (MOREAU’S ENVELOPE OF DIFFERENTIABLE FUNCTIONS) Suppose the function f : X → IR is differentiable. Then, the Moreau’s envelope defined in equation 44 has the following upper bound and lower bounds
f(y)− β 2 ∫ 1 0 sup x∈X ‖∇f(y + s(x− y))‖22ds ≤M β f (y) ≤ f(y). (45)
In particular, when f is Lf -Lipschitz, we have
f(y)− βL2f
2 ≤Mβf (y) ≤ f(y). (46)
The proof is presented in Appendix D.1.
Now, we return to Equation equation 43. We leverage the lower bound on Moreau’s envelope in equation 45 of Lemma C.5 as follows∣∣∣ sup
µ∈P M̂MDµ[PV , PW ]− sup µ∈P MMDµ[PV , PW ] ∣∣∣ ≤ 4 ∣∣∣∣minλ≥0 { λR2 − ∫ IRD M 1 2λ −En(ξ)µ0(dξ)
}∣∣∣∣ ≤ 4 ∣∣∣∣∣minλ≥0 { λR2 + IEµ0 [En(ξ)] + 1 4λ IEµ0 [∫ 1 0 sup ζ∈IRD ‖∇En((1− s)ξ + sζ)‖22ds
]}∣∣∣∣∣ ≤ 4|IEµ0 [En(ξ)]|+ 4RIEµ0 [∫ 1 0 sup ζ∈IRD ‖∇En((1− s)ξ + sζ)‖22ds ] . (47)
xxi
Let ζ∗ = ζ∗(ξ, s) = arg supζ∈IRD ‖∇En(1 − s)ξ + sζ‖2. Then, applying the union bound in conjunction with Inequality equation 47 yields
IP (∣∣∣ sup µ∈P M̂MDµ[PV , PW ]− sup µ∈P MMDµ[PV , PW ] ∣∣∣ ≥ δ)
≤ IP (∣∣∣∣∣ ∫ IRD En(ξ)µ0(dξ) ∣∣∣∣∣ ≥ δ8 ) + IP (∫ IRD ∫ 1 0 ‖∇En((1− s)ξ + sζ∗)‖22dsµ0(dξ) ≥ δ 8R ) .
(48)
Now, we state the following lemma: Lemma C.6. (TAIL BOUNDS FOR THE FINITE SAMPLE ESTIMATION ERROR) Consider the estimation error En defined in equation 41. Then, the following statements hold:
• Z = ‖∇En(ξ)‖22 is a sub-exponential random variable with the Orlicz norm of ‖Z‖ψ1 ≤ 9×29×L4
n2 for every ξ ∈ IR D. Moreover,
IP (∫ IRD ∫ 1 0 ‖∇En((1− s)ξ + sζ∗)‖22dsµ0(dξ) ≥ δ ) ≤ 2e− n2δ 9×29×L4 +L 4 9 ,
(49)
• En(ξ) is zero-mean sub-Gaussian random variable with the Orlicz norm of ‖En(ξ)‖ψ2 ≤ 16 √ 3L4
n for every ξ ∈ IR D. Moreover,
IP (∣∣∣∣∫ IRD En(ξ)µ0(dξ) ∣∣∣∣ ≥ δ) ≥ 2e− n2δ216√3L4 . (50) The proof of Lemma C.6 is presented in Appendix C.2.
Now, we leverage the probability bounds equation 49 and equation 50 of Lemma C.6 to upper bound the terms on the right hand side of equation 48 as below
IP (∣∣∣ sup µ∈P M̂MDµ[PV , PW ]− sup µ∈P MMDµ[PV , PW ] ∣∣∣ ≥ δ) ≤ 2e− n2δ2√3×211×L4 + 2e− n2δ9×212×RL4 +L49
≤ 4max { e − n 2δ2√ 3×211×L4 , e − n 2δ 9×212×RL4 +L 4 9 } ,
(51)
where the last inequality comes from the basic inequality a+ b ≤ 2 max{a, b}. Therefore, with the probability of at least 1− %, we have that∣∣∣ sup
µ∈P M̂MDµ[PV , PW ]− sup µ∈P MMDµ[PV , PW ] ∣∣∣ ≤ max 3 1 4 × 2 11 2 × L2 n ln 1 2 ( 4 % ) , 9× 212 ×RL4 n2 ln 4eL49 %
 , (52) for allW ∈ W . Lemma C.7. (DISTANCE BETWEEN MINIMA OF ADJACENT FUNCTIONS) Let Ψ(W ) : W → IR and Φ(W ) :W → IR. Further, suppose ‖Ψ(W )− Φ(W )‖∞ ≤ δ for some δ > 0. Then,∣∣∣∣ minW∈WΨ(W )− minW∈W Φ(W )
∣∣∣∣ ≤ δ. (53) xxii
See Appendix D.4 for the proof.
Let Ψ(W ) def= supµ∈P M̂MDµ[PV , PW ], and Φ(W ) def = supµ∈P MMDµ[PV , PW ]. Then, from Inequality equation 53, we have the following upper bound on A1
A1 = ∣∣∣MMDµ∗ [PV , PW∗ ]− min
W∈W sup µ∈P
M̂MDµ[PV , PW ] ∣∣∣
= ∣∣∣ min W∈W sup µ∈P MMDµ[PV , PW ]− min W∈W sup µ∈P M̂MDµ[PV , PW ] ∣∣∣
≤ max
{ 3 1 4 × 24 × L2
n ln
1 2
( 4
%
) , 9× 211 ×RL4
n2 ln
( 4e L4 9
%
)} . (54)
with the probability of (at least) 1− %.
Upper bound on A2:
To establish the upper bound on A2, we recall that
M̂MDµ̂N [PV , PW ] = 1 n(n− 1) 1 N ∑ i 6=j N∑ k=1 yiyjϕ(xi; ξ k)ϕ(xj ; ξ k) (55a)
M̂MDµ[PV , PW ] = 1 n(n− 1) ∑ i 6=j yiyjIEµ[ϕ(xi; ξ)ϕ(xj ; ξ)]. (55b)
Therefore,∣∣∣ sup µ∈P M̂MDµ[PV , PW ]− sup µ̂N∈PN M̂MDµ̂N [PV , PW ] ∣∣∣
≤ ∣∣∣ sup µ∈P IEµ[ϕ(xi; ξ)ϕ(xj ; ξ)]− sup µ̂N∈PN IEµ̂N [ϕ(xi; ξ)ϕ(xj ; ξ)] ∣∣∣. (56)
Here, the last inequality is due to Theorem C.4 and the following duality results hold
sup µ∈P IEµ[ϕ(x; ξ)ϕ(x̂; ξ)] = inf λ≥0
{ λR2 − ∫ IRD inf ζ∈IRD {λ‖ξ − ζ‖22 − ϕ(xi; ζ)ϕ(xj ; ζ)}µ0(dξ) }
sup µ̂N∈PN IEµ̂N [ϕ(xi; ξ)ϕ(xj ; ξ)] = inf λ≥0
{ λR2 − 1
N N∑ k=1 inf ζ∈IRD {λ‖ξk0 − ζ‖22 − ϕ(xi; ζ)ϕ(xj ; ζ)}
} .
Now, in the sequel, we establish a uniform concentration result for the following function
Tλ(x,x̂) : IR N×D 7→ IR
(ξ10, · · · , ξN0 ) 7→ Tλ(x,x̂)(ξ10, · · · , ξN0 ) = 1
N N∑ k=1 M 1 2λ −ϕ(x,·)ϕ(x̂,·)(ξ k 0 )− ∫ IRD M 1 2λ −ϕ(x,·)ϕ(x̂,·)(ξ)µ0(dξ).
Then, from equation 56 we have∣∣∣ sup µ∈P M̂MDµ[PV , PW ]− sup µ̂N∈PN M̂MDµ̂N [PV , PW ] ∣∣∣ ≤ sup λ≥0 sup x,x̂∈X |Tλ(x,x̂)(ξ 1 0, · · · , ξN0 )|. (57)
We now closely follow the argument of Rahimi & Recht (2008) to establish a uniform concentration result with respect to the data points x, x̂ ∈ X . In particular, consider an -net cover of X ⊂ IRd.
xxiii
Then, we require N = ( 4diam(X ) )d balls of the radius > 0, e.g., see (Pollard, 1990, Lemma
4.1, Section 4). Let Z = {z1, · · · , zN } ⊂ X denotes the center of the covering net. Now, let (ξ10, · · · , ξk0 , · · · , ξN0 ) ∈ IR N×D and (ξ10, · · · , ξ̃k0 , · · · , ξN0 ) ∈ IR N×D be two sequences that differs in the k-th coordinate for 1 ≤ k ≤ N . Then,∣∣T(zi,zj)(ξ10, · · · , ξk0 , · · · , ξN0 )− T(zi,zj)(ξ10, · · · , ξ̃k0 , · · · , ξN0 )∣∣ = 1
N ∣∣∣M 12λ−ϕ(zi;·)ϕ(zj ;·)(ξk0 )−M 12λ−ϕ(zi,·)ϕ(zj ,·)(ξ̃k0 )∣∣∣. (58) Without loss of generality suppose M 1 2λ
−ϕ(zi;·)ϕ(zj ;·)(ξ k 0 ) ≥M
1 2λ
−ϕ(zi;·)ϕ(zj ;·)(ξ̃ k 0 ). Then,
M 1 2λ
−ϕ(zi;·)ϕ(zj ;·)(ξ k 0 )−M
1 2λ
−ϕ(zi,·)ϕ(zj ,·)(ξ̃ k 0 )
= inf ζ∈IRD
{ λ‖ζ − ξk0‖22 − ϕ(zi; ζ)ϕ(zj ; ζ) } − inf ζ∈IRD { λ‖ζ − ξ̃k0‖22 − ϕ(zi; ζ)ϕ(zj ; ζ) } (a) ≤ −ϕ(zi; ξk0 )ϕ(zj ; ξk0 )− inf ζ∈IRD { λ‖ζ − ξ̃k0‖22 − ϕ(zi; ζ)ϕ(zj ; ζ)
} (b) ≤ −ϕ(zi; ξk0 )ϕ(zj ; ξk0 ) + sup ζ∈IRD {ϕ(zi; ζ)ϕ(zj ; ζ)}
(c) ≤ 2L2, (59)
where (a) follows by letting ζ = ξk0 in the first optimization problem, (b) follows by using the fact that −λ‖ζ − ξ̃k0‖2 is non-positive for any ζ ∈ IR
D and can be dropped, and (c) follows from Assumption (A.2).
Now, plugging the upper bound in equation 59 into equation 58 yields∣∣Tλ(zi,zj)(ξ10, · · · , ξk0 , · · · , ξN0 )− Tλ(zi,zj)(ξ10, · · · , ξ̃k0 , · · · , ξN0 )∣∣ ≤ 2L2N . From McDiarmid’s Martingale inequality McDiarmid (1989) and the union bound, we obtain that
IP ( ∪zi,zj∈Z |Tλ(zi,zj)(ξ 1 0, · · · , ξN0 )| ≥ δ ) ≤ ( 4diam(X ) )d · ( −Nδ 2
L2
) , (60)
for all λ ≥ 0. Now, consider arbitrary points (x, x̂) ∈ X ×X . Let the center of the balls containing those points be zi, zj ∈ Z , i.e., x ∈ IBε(zi) and x̂ ∈ IBε(zj) for some zi, zj ∈ Z . Then, by the triangle inequality, we have that
|Tλ(x,x̂)(ξ10, · · · , ξN0 )− T(zi,zj)(ξ 1 0, · · · , ξN0 )|
≤ |Tλ(x,x̂)(ξ10, · · · , ξN0 )− Tλ(zi,x̂)(ξ 1 0, · · · , ξN0 )|
+ |Tλ(zi,x̂)(ξ 1 0, · · · , ξN0 )− Tλ(zi,zj)(ξ 1 0, · · · , ξN0 )|
≤ ‖∇xTλ(x,x̂)(ξ10, · · · , ξN0 )‖2‖x− zi‖2 + ‖∇x̂Tλ(x,x̂)(ξ10, · · · , ξN0 )‖2‖x̂− zj‖2 ≤ 2LT , (61)
where LT = LT (ξ10, · · · , ξN0 ) def = supx,x̂∈X ‖∇xTλ(x,x̂)(ξ 1 0, · · · , ξN0 )‖2 is the Lipschitz constant of the mapping T . Note that the Lipschitz constant LT is a random variable with respect to the random feature samples ξ0, · · · , ξN . Let (x∗, x̂∗) def = arg supx,x̂∈X ‖∇xTλ(x,x̂)(ξ 1 0, · · · , ξN0 )‖2.
xxiv
We compute an upper bound on the second moment of the random variable LT as follows
IEµ0 [ L2T ] = IEµ0 [ ‖∇xTλ(x∗,x̂∗)(ξ 1 0, · · · , ξN0 )‖22 ] = IEµ0 ∥∥∥∥∥ 1N N∑ k=1 ∇xM 1 2λ −ϕ(x∗;·)ϕ(x̂∗;·)(ξ k 0 )− ∫ IRD ∇xM 1 2λ −ϕ(x∗;·)ϕ(x̂∗;·)(ξ)µ0(dξ) ∥∥∥∥∥ 2
2  = IEµ0 ∥∥∥∥∥ 1N N∑ k=1 ∇xM 1 2λ −ϕ(x∗;·)ϕ(x̂∗;·)(ξ k 0 ) ∥∥∥∥∥ 2
2
− IEµ0 [∥∥∥∥∫
IRD ∇xM
1 2λ −ϕ(x∗;·)ϕ(x̂∗;·)(ξ)µ0(dξ) ∥∥∥∥2 2 ]
≤ 1 N2 IEµ0 ∥∥∥∥∥ N∑ k=1 ∇xM 1 2λ −ϕ(x∗;·)ϕ(x̂∗;·)(ξ k 0 ) ∥∥∥∥∥ 2
2  . We further proceed using the triangle inequality as well as the basic inequality (a1 + a2 + · · · + aN ) 2 ≤ N(a21 + a22 + · · ·+ a2N ),
IEµ0 [ L2T ] = 1
N2 IEµ0 ∥∥∥∥∥ N∑ k=1 ∇xM 1 2λ −ϕ(x∗;·)ϕ(x̂∗;·)(ξ k 0 ) ∥∥∥∥∥ 2
2  ≤ 1 N2 IEµ0 ( N∑ k=1 ∥∥∥∇xM 12λ−ϕ(x∗;·)ϕ(x̂∗;·)(ξk0 )∥∥∥2 )2
≤ 1 N N∑ k=1 IEµ0 [∥∥∥∇xM 12λ−ϕ(x∗;·)ϕ(x̂∗;·)(ξk0 )∥∥∥22 ] . (62)
To proceed from equation 62, we leverage the following lemma: Lemma C.8. (MOREAU’S ENVELOP OF PARAMETRIC FUNCTIONS) Consider the parametric function f : X ×Θ→ IR and the associated Moreau’s envelope for a given θ ∈ Θ ⊂ IRd:
Mβf(·;θ)(x) = infy∈X
{ 1
2β ‖x− y‖22 + f(y;θ)
} . (63)
Furthermore, define the proximal operator as follows
Proxβf(·;θ)(x) = arg infy∈X
{ 1
2β ‖x− y‖22 + f(y;θ)
} . (64)
Then, Moreau’s envelope has the following upper bound∥∥∥∇θMβf(·;θ)(x)∥∥∥ 2 ≤ ∥∥∥∇θf (Proxβf(·;θ)(x);θ)∥∥∥ 2 . (65)
The proof is presented in Appendix D.2.
Equipped with Inequality equation 65 of Lemma C.8, we now compute an upper bound on the right hand side of equation 62 as follows
IEµ0 [ L2T ] ≤ 1 N N∑ k=1 IEµ0 [|ϕ(x̂∗; ξ)|2 · ‖∇xϕ(x∗; ξk0 )‖22] ≤ L4, (66)
where the last inequality is due to (A.2).
xxv
Invoking Markov’s inequality now yields IP ( |Tλ(x,x̂)(ξ 1 0, · · · , ξN0 )− Tλ(zi,zj)(ξ 1 0, · · · , ξN0 )| ≥ δ ) = IP ( LT ≥ δ
2 ) ≤ ( 2
δ
)2 IEµ0 [L 2 T ]
≤ ( 2
δ
)2 L4. (67)
Now, using the union bound, for every arbitrary pair of data points (x, x̂) ∈ X × X the following inequality holds
IP ( |Tλ(x,x̂)(ξ 1 0, · · · , ξN0 )| ≥ δ ) ≤ IP ( |Tλ(zi,zj)(ξ 1 0, · · · , ξN0 )| ≥ δ/2 ) + IP ( |Tλ(x,x̂)(ξ 1 0, · · · , ξN0 )− Tλ(zi,zj)(ξ 1 0, · · · , ξN0 )| ≥ δ/2
) (68)
≤ ( 2
δ
)2 L4 + ( 4diam(X ) )d · ( −Nδ 2
L2
) .
Following the proposal of Rahimi & Recht (2008), we select = (κ1/κ2) 1 d+2 , where κ1 def = (4diam(X ))d · e− 2Nλδ2 2L2λ+L4 and κ2 def = (2/δ)2L4. Then,
IP ( sup x,x̂∈X |Tλ(x,x̂)(ξ 1 0, · · · , ξN0 )| ≥ δ ) ≤ 28 ( L2diam(X ) δ )2 · ( − Nδ 2 L2(d+ 2) ) .
Thus, with the probability of at least 1− %, the following inequality holds
sup λ≥0 sup x,x̂∈X
|Tλ(x,x̂)(ξ 1 0, · · · , ξN0 )| ≤
( L2(d+ 2)
N W
( 28Ndiam2(X )
%
)) 1 2
, (69)
where W(·) is the Lambert W -function.4 Since W(x) ≤ ln(x) for x > e, we can rewrite the upper bound in terms of elementary functions
sup λ≥0 sup x,x̂∈X
|Tλ(x,x̂)(ξ 1 0, · · · , ξN0 )| ≤
√ L2(d+ 2)
N ln
1 2
( 28Ndiam2(X )
%
) , (70)
provided thatN is sufficiently large and/or % is sufficiently small so that 2 8Ndiam2(X )
% ≥ e. Plugging Inequality equation 70 in equation 57 now results in the following inequality∣∣∣ sup
µ∈P M̂MDµ[PV , PW ]− sup µ̂N∈PN M̂MDµ̂N [PV , PW ]
∣∣∣ ≤√L2(d+ 2) N ln 1 2 ( 28Ndiam2(X ) % ) , (71)
for allW ∈ W . Employingequation 53 from Lemma C.7 now yields the following upper bound A2 = ∣∣∣ min W∈W sup µ∈P M̂MDµ[PV , PW ]− min W∈W sup µ̂N∈PN M̂MDµ̂N [ PV , PW ]∣∣∣ ≤ √ L2(d+ 2)
N ln
1 2
( 28Ndiam2(X )
%
) . (72)
4Recall that the lambert W -function is the inverse of the function f(W ) =WeW .
xxvi
Upper bound on A3:
Recall that the solution of the empirical risk function of equation 14 is denoted by
(ŴN∗ , µ̂ N ∗ )
def = arg min
W∈W arg inf µ̂N∈PN M̂MD
α µ̂N [PV , PW ]
= arg min W∈W arg sup µ̂N∈PN
8 n(n− 1) ∑
1≤i<j≤n
yiyjIEµ̂N∈PN [ϕ(xi; ξ)ϕ(xj ; ξ)]
− 8 n(n− 1)α ∑ 1≤i<j≤n (IEµ̂N [ϕ(xi; ξ)ϕ(xj ; ξ)]) 2. (73)
We also define the solution of the empirical kernel alignment as follows
(ŴN , µ̂ N ) def = arg min
W∈W arg sup µ̂N∈PN M̂MDµ̂N [PV , PW ]
= 8 n(n− 1) ∑
1≤i<j≤n
yiyjIEµ̂N [ϕ(xi; ξ)ϕ(xj ; ξ)]. (74)
Due to the optimality of the empirical measure µ̂N∗ for the inner optimization in equation 73, the following inequality holds
M̂MD α
µ̂N [ PV , PŴN∗ ] ≤ M̂MD α µ̂N∗ [PV , PŴN∗ ]
≤ 8 n(n− 1) ∑ 1≤i<j≤n yiyjIEµ̂N∗ [ϕ(xi; ξ)ϕ(xj ; ξ)]. (75)
Upon expansion of M̂MD α
µ̂N [ PV , PŴN∗ ] , and after rearranging the terms in equation 75, we arrive
at M̂MDµ̂N [ PV , PŴN∗ ] − M̂MDµ̂N∗ [ PV , PŴN∗ ] = 8 n(n− 1) ∑
1≤i<j≤n
yiyj(IEµ̂N [ϕ(xi; ξ)ϕ(xj ; ξ)]− IEµ̂N∗ [ϕ(xi; ξ)ϕ(xj ; ξ)])
≤ 8 n(n− 1)α ∑ 1≤i<j≤n ( IEµ̂N [ϕ(xi; ξ)ϕ(xj ; ξ)] )2 ≤ 8L 4
α , (76)
where the last inequality is due to the fact that ‖ϕ‖∞ < L by (A.1). Now, due to optimality of ŴN for the outer optimization problem in equation 74, we have
M̂MDµ̂N [ PV , PŴN ] ≤ M̂MDµ̂N [ PV , PŴN∗ ] . (77)
Putting together Inequalities equation 76 and equation 77 yields M̂MDµ̂N [ PV , PŴN ] − M̂MDµ̂N∗ [ PV , PŴN∗ ] ≤ 8L 4 α . (78)
Similarly, due to the optimality of the empirical measure µ̂N for the optimization in equation 74 we have that
M̂MDµ̂N∗ [ PV , PŴN∗ ] ≤ M̂MDµ̂N∗ [ PV , PŴN ] ≤ M̂MDµ̂N [ PV , PŴN ] . (79)
xxvii
Combining equation 78 and equation 79 then yields A3 = ∣∣∣M̂MDµ̂N∗ [PV , PŴN∗ ]− M̂MDµ̂N [PV , PŴN ]∣∣∣ ≤ 8L4α . (80)
Upper bound on A4:
The upper bound on A4 can be obtained exactly the same way as A1. Specifically, from equation 52 it follows directly that
A4 = ∣∣∣M̂MDµ̂N∗ [PV , PŴN∗ ]−MMDµ̂N∗ [PV , PŴN∗ ]∣∣∣
≤ sup µ̂N∈PN ∣∣∣M̂MDµ̂N [PV , PŴN∗ ]−MMDµ̂N [PV , PŴN∗ ]∣∣∣ (81) ≤ max { 3 1 4 × 2 112 × L2
n ln
1 2
( 4
%
) , 9× 212 ×RL4
n2 ln
( 4e L4 9
%
)} . (82)
Now, plugging the derived upper bounds in A1-A4 in equation 38 and employing the union bound completes the proof.

heading: C.2 PROOF OF THEOREM 4.2
text: The proof has three main ingredients and follows the standard procedure in the literature, see, e.g., Wang et al. (2017); Luo & Mattingly (2017). In the first step, we identify the mean-field limit of the particle SGD in equation 15. In the second step, we prove the convergence of the measured-valued process {(µNt )0≤t≤T } to the mean-field solution by establishing the pre-compactness of Sokhorhod space. Lastly, we prove the uniqueness of the mean-field solution of the particle SGD.
Step 1-Identification of the scaling limit: First, we identify the weak limit of converging subsequences via the action of the empirical measure µ̂Nm(ξ) = 1 N ∑N k=1 δ(ξ − ξkm) on a test function f ∈ C3b (IR D). In particular, we use the standard techniques of computing the scaling limits from Luo & Mattingly (2017).
Recall that the action of an empirical measure on a bounded function is defined as follows
〈f, µ̂Nm〉 def =
1
N N∑ k=1 f(ξkm). (83)
We analyze the evolution of the empirical measure µ̂Nm via its action on a test function f ∈ C3b (IR D). Using Taylor’s expansion, we obtain 〈f, µ̂Nm+1〉 − 〈f, µ̂Nm〉 = 〈f, µ̂Nm+1〉 − 〈f, ν̂Nm+1〉
= 1
N N∑ k=1 f(ξkm+1)− f(ξkm)
= 1
N N∑ k=1 ∇f(ξkm)(ξkm+1 − ξkm)T +RNm.
where RNm is a remainder term defined as follows
RNm def =
1
N N∑ k=1 (ξkm+1 − ξkm)T∇2f(ξ̃k)(ξkm+1 − ξkm), (84)
xxviii
where ξ̃k def= (ξ̃k(1), · · · , ξ̃k(p)), and ξ̃k(i) ∈ [ξkm(i), ξkm+1(i)], for i = 1, 2, · · · , p.
Plugging the difference term (ξkm+1 − ξkm) from the SGD equation in equation 15 results in
〈f, µ̂Nm+1〉 − 〈f, µ̂Nm〉 (85)
= η
N2α N∑ k=1 ∇f(ξkm) ·
(( 1
N N∑ `=1 ϕ(xm; ξ ` m)ϕ(x̃m; ξ ` m)− αymỹm
) ∇ξ ( ϕ(xm; ξ k m)ϕ(x̃m; ξ k m) )) +RNm.
Now, we define the drift and Martingale terms as follows
DNm def =
η
Nα ∫∫ X×Y ( 〈ϕ(x, ξ)ϕ(x̃, ξ), µ̂Nm〉 − αyỹ ) (86a)
× 〈∇f(ξ)(ϕ(x̃; ξ)∇ξϕ(x; ξ) + ϕ(x; ξ)∇ξϕ(x̃; ξ)), µ̂Nm〉dP⊗2x,y((x, y), (x̃, ỹ))
MNm def =
η
Nα
( 〈ϕ(xm, ξ)ϕ(x̃m, ξ), µ̂Nm〉 − αymỹm ) (86b)
× 〈∇f(ξ)(ϕ(x̃m; ξ)∇ξϕ(xm; ξ) + ϕ(x̃m; ξ)∇ξϕ(xm; ξ)), µ̂Nm〉 − DNm.
respectively. Using the definitions ofDNm andM N m in equation 86a-equation 86b, we recast Equation equation 85 as follows
〈f, µ̂Nm+1〉 − 〈f, µ̂Nm〉 = DNm +MNm +RNm. (87) Summation over ` = 0, 1, 2 · · · ,m− 1 and using the telescopic sum yields
〈f, µ̂Nm〉 − 〈f, µ̂N0 〉 = m−1∑ `=0 DN` + m−1∑ `=0 MN` + m−1∑ `=0 RN` . (88)
We also define the following continuous embedding of the drift, martingale, and the remainder terms as follows
DNt def = bNtc∑ `=0 DN` (89a)
MNt def = bNtc∑ `=0 MN` (89b)
RNt def = bNtc∑ `=0 RN` , t ∈ (0, T ]. (89c)
The scaled empirical measure µNt def = µ̂NbNtc then can be written as follows
〈f, µNt 〉 − 〈f, µN0 〉 = DNt +MNt +RNt . (90)
Since the drift process (DNt )0≤t≤T is a piecewise cádlág process, we have
DN` =
∫ `+1 N
` N
R[µs]ds, (91)
where the functional R[µs] is defined as follows
R[µs] def = η
α ∫∫ X×Y (〈ϕ(x, ξ)ϕ(x̃, ξ), µs〉 − αyỹ) (92)
× 〈∇f(ξ)(ϕ(x̃; ξ)∇ξϕ(x; ξ) + ϕ(x; ξ)∇ξϕ(x̃; ξ))T , µs〉P⊗2x,y((dx,dx̃), (dy,dỹ)).
xxix
Therefore, the expression in equation 90 can be rewritten as follows 〈f, µNt 〉 − 〈f, µN0 〉 = ∫ t
0
R[µs]ds+MNt +RNt . (93)
In the following lemma, we prove that the remainder term sup0≤t≤T |RNt | vanishes in probabilistic sense as the number of particles tends to infinity N →∞: Lemma C.9. (LARGE N -LIMIT OF THE REMAINDER PROCESS) Consider the remainder process (RNt )0≤t≤T defined via scaling in equation 84-equation 89c. Then, there exists a constant C0 > 0 such that
sup 0≤t≤T
|RNt | ≤ C0T
N
( ηL2 + 2ηL4
α
) . (94)
and thus lim supN→∞ sup0≤t≤T |RNt | = 0.
Proof. The proof is relegated to Appendix D.5.
We can also prove a similar result for the process defined by the remainder term:
Lemma C.10. (LARGE N -LIMIT OF THE MARTINGALE PROCESS) Consider the Martingale process (MNt )0≤t≤T defined via scaling in equation 86b-equation 89b. Then, for some constant C1 > 0, the following inequality holds
IP ( sup
0≤t≤T |MNt | ≥ ε ) ≤ 1 Nαε 4 √ 2L2 √ bNT cηC1(L2 + α)2. (95)
In particular, with the probability of at least 1− ρ, we have
sup 0≤t≤T
|MNt | ≤ 1 Nαρ 4 √
2L2 √ bNT cηC1(L2 + α)2. (96)
and thus lim supN→∞ sup0≤t≤T |MNt | = 0 almost surely.
Proof. The proof is deferred to Appendix D.6.
Now, using the results of Lemmata C.10 and C.9 in conjunction with equation 93 yields the following mean-field equation as N →∞,
〈µt, f〉 = 〈µ0, f〉+ η
α ∫ t 0 (∫∫ X×Y (〈ϕ(x, ξ)ϕ(x̃, ξ), µs〉 − αyỹ) (97)
× 〈∇f(ξ)(ϕ(x̃; ξ)∇ξϕ(x; ξ) + ϕ(x; ξ)∇ξϕ(x̃; ξ)), µs〉P⊗2x,y((dx,dx̃), (dy,dỹ))
) ds.
Notice that he mean-field equation in equation 97 is in the weak form. When the Lebesgue density pt(ξ) = dµt/dξ exists, the McKean-Vlasov PDE in equation 22 can be readily obtained from equation 97.
Step 2: Pre-compactness of the Skorkhod space: To establish our results in this part of the proof, we need a definition and a theorem:
xxx
Definition C.11. (TIGHTNESS) A setA of probability measures on a metric space S is tight if there exists a compact subset S0 ⊂ S such that
ν(S0) ≥ 1− ε, for all ν ∈ A, (98)
for all ε > 0. A sequence {XN}N∈IN of random elements of the metric space S is tight if there exists a compact subset S0 ⊂ S such that
ν(XN ∈ S0) > 1− ε, (99)
for all ε > 0, and all N ∈ IN.
Now, to show the tightness of the measured valued process (µNt )0≤t≤T , we must verify Jakubowski’s criterion (Jakubowski, 1986, Thm. 1):
Theorem C.12. (JAKUBOWSKI’S CRITERION (JAKUBOWSKI, 1986, THM. 1)) A sequence of measured-valued process {(ζNt )0≤t≤T }N∈IN is tight in DM(IRD)([0, T ]) if and only if the following two conditions are satisfied:
(J.1) For each T > 0 and γ > 0, there exists a compact set UT,γ such that
lim N→∞
inf IP ( ζNt ∈ UT,γ ,∀t ∈ (0, T ] ) > 1− γ. (100)
This condition is referred to as the compact-containment condition.
(J.2) There exists a family H of real-valued functions H : M(IRD) 7→ IR that separates points inM(IRD) and is closed under addition such that for every H ∈ H, the sequence {(H(ξNt ))0≤t≤T }N∈IN is tight in DIR([0, T ]).
To establish (J1), we closely follow the proof of (Giesecke et al., 2013, Lemma 6.1.). In particular, for each L > 0, we define SL = [0, B]p. Then, SB ⊂ IRp is compact, and for each t ≥ 0, and N ∈ IN, we have
IE[µNt (IR p/SB)] =
1
N N∑ k=1 IP ( ‖ξkbNtc‖2 ≥ B ) (101)
(a) ≤ 1 N N∑ k=1 IE[‖ξkbNtc‖2] B
(102)
(b) ≤ c0 + ηαL 2T + 2ηL4T
B , (103)
where (a) follows from Markov’s inequality, and (b) follows from the upper bound on the norm of the particles in equation 177 of Appendix D. We now define the following set
UB = { µ ∈M(IRp) : µ(IRp/S(B+j)2) <
1√ B + j
for all j ∈ IN } . (104)
xxxi
We let UT,γ = UB , where UB is the completion of the set UB . By definition, UT,γ is a compact subset ofM(IRD). Now, we have
IP ( µNt 6∈ UT,γ ) ≤ ∞∑ j=1 IP ( µNt (IR p/S(B+j)2) > 1√ B + j )
≤ ∞∑ j=1 IE[µNt (IR p/S(B+j)2)] 1/ √ B + j
≤ ∞∑ j=1 c0 + ηL 2T + 2(η/α)L4T (B + j)2/ √ B + j
= ∞∑ j=1 c0 + ηL 2T + 2(η/α)L4T (B + j)3/2 . (105)
Now, since
lim B→∞ ∞∑ j=1 c0 + ηL 2T + 2(η/α)L4T (B + j)3/2 = 0, (106)
this implies that for any γ > 0, there exists a B > 0, such that
lim N→∞
inf IP ( µNt ∈ UB ,∀t ∈ (0, T ] ) > 1− γ. (107)
This completes the proof of (J.1). To verify (J.2), we consider the following class of functions
H def= {H : ∃f ∈ C3b (IR D) such that H(µ) = 〈µ, f〉,∀µ ∈M(IRD)}. (108)
By definition, every function H ∈ H is continuous with respect to the weak topology ofM(IRD) and further the class of functions H separate points inM(IRD) and is closed under addition. Now, we state the following sufficient conditions to establish (J.2). The statement of the theorem is due to (Robert, 2013, Thm. C.9): Theorem C.13. (TIGHTNESS IN DIR([0, T ]), (ROBERT, 2013, THM. C.9)) A sequence {(ZNt )0≤t≤T }N∈IN is tight in DIR([0, T ]) iff for any δ > 0, we have
(T.1) There exists > 0, such that
IP(|ZN0 | > ) ≤ δ, (109) for all N ∈ IN.
(T.2) For any ρ > 0, there exists σ > 0 such that
IP ( sup
t1,t2≤T,|t1−t2|≤ρ |ZNt1 − Z N t2 | > σ
) ≤ δ, (110)
This completes the tightness proof of the of the laws of the measured-valued process {(µNt )0≤t≤T }N∈IN. Now, we verify the condition (J.2) by showing that the sufficient conditions (T.1) and (T.2) hold for function values {(H(µNt ))0≤t≤T }N∈IN, where H ∈ H and H is defined in Eq. equation 108. Now, condition (T.1) is readily verified since
H(µN0 ) = 〈µN0 , f〉 = ∫
IRD f(ξ)µN0 (dξ) (111) ≤ ‖f‖∞ ∫
IRD µN0 (dξ) (112)
≤ b, (113)
xxxii
where in the last step, we used the fact that f ∈ C3b (IR D), and hence, ‖f‖∞ ≤ b. Thus, IP(H(µN0 ) ≥ b) = 0 for all N ∈ IN, and the condition (T.1) is satisfied. Now, consider the condition (T.2). From Equation equation 93, and with 0 ≤ t1 < t2 ≤ T we have
|H(µNt1)−H(µ N t1)| = |〈f, µ N t1〉 − 〈f, µ N t2〉| ≤ ∫ t2 t1 |R[µs]|ds+ |MNt1 −M N t2 |+ |R N t1 −R N t2 |. (114)
To bound the first term, recall the definition of R[µs] from equation 92. The following chain of inequalities holds,
|R[µs]| ≤ η
α IE P⊗2x,y
[|〈ϕ(x, ξ)ϕ(x̃, ξ), µs〉 − αyỹ||〈∇f(ξ)(ϕ(x̃; ξ)∇ξϕ(x; ξ) + ϕ(x; ξ)∇ξϕ(x̃; ξ))T , µs〉|]
≤ η α IE P⊗2x [(|〈ϕ(x, ξ)ϕ(x̃, ξ), µs〉|+ α)|〈∇f(ξ)(ϕ(x̃; ξ)∇ξϕ(x; ξ) + ϕ(x; ξ)∇ξϕ(x̃; ξ))T , µs〉|]. (115)
Let I : IRD → IR, I(ξ) = 1 denotes the identity function. Notice that 〈I, µs〉 = ∫
IRD µs(ds) = 1.
From equation 115, we proceed as follows
|R[µs]| ≤ η
α IE P⊗2 X
[(‖ϕ‖2∞ · |〈I, µs〉|+ α) · ‖∇f(ξ)(ϕ(x̃; ξ)∇ξϕ(x; ξ) + ϕ(x; ξ)∇ξϕ(x̃; ξ))T ‖∞ · |〈I, µs〉|]
≤ η α IE P⊗2 X [(‖ϕ‖2∞ + α) · ‖∇f(ξ)(ϕ(x̃; ξ)∇ξϕ(x; ξ) + ϕ(x; ξ)∇ξϕ(x̃; ξ))T ‖∞]
≤ 2η α (L2 + α)L2C1, (116)
where the last inequality is due to (A.1). Therefore,∫ t2 t1 |R[µs]|ds ≤ s0|t2 − t1|, (117)
where s0 def =
2η α (L2 + α)L2C1.
Consider the middle term of equation 114. Using the definition of the martingale term in equation 89b, we obtain that
|MNt1 −M N t2 | = ∣∣∣∣∣∣ bNt1c∑ `=0 MN` − bNt2c∑ `=0 MN` ∣∣∣∣∣∣ ≤
∣∣∣∣∣∣ bNt2c∑ `=bNt1c MN` ∣∣∣∣∣∣ . (118) In Equation of Section D, we have proved the following concentration bound
IP(|MNm | ≥ ε) ≤ 2 ( − N 2α2ε2
8mL4η2C21 (L 2 + α)2
) , ∀m ∈ [0, NT ] ∩ IN. (119)
Now, recall the alternative definition of the sub-Gaussian random variables: Definition C.14. (SUB-GAUSSIAN RANDOM VARIABLES BOUCHERON ET AL. (2013)) A random variable X is σ2-sub-Gaussian if
IE[(λ(X − IE[X]))] ≤ (λ2σ2
2
) . (120)
xxxiii
We enumerate a few standard consequences of sub-Gaussianity Boucheron et al. (2013). If Xi are independent and σ2i -sub-Gaussian, then ∑n i=1Xi is ∑n i=1 σ 2 i -sub-Gaussian. Moreover, X is σ2-sub-Gaussian if and only if
IP(|X − IE[X]| ≥ ε) ≤ ( − ε 2
2σ2
) . (121)
Now, it is clear from equation 119 andthat MNm is sub-Gaussian random variable with a zero mean, and with the parameter σ2m = 4mL4η2C21 (L 2 + α)2
N2α2 . Therefore, ∑bNt2c `=bNt1cM N ` is sub-Gaussian
with the parameter σ2(t1, t2) def =
2L4η2C21 (L 2 + α)2
N2α2 (bNt1c− bNt2c+ 1)(bNt1c+ bNt2c). Con-
sequently, from Inequality equation 118 and the concentration inequality in equation 121, we have
IP ( sup
t1,t2≤T,|t1−t2|≤ρ |MNt1 −M N t2 | ≥ ε
) ≤ IP  sup t1,t2≤T,|t1−t2|≤ρ ∣∣∣∣∣∣ bNt2c∑ `=bNt1c MN` ∣∣∣∣∣∣ ≥ ε  (122)
= IP ∣∣∣∣∣∣ bNt∗2c∑ `=bNt∗1c MN` ∣∣∣∣∣∣ ≥ ε  (123)
≤ 2 ( − ε 2
σ2(t∗1, t ∗ 2)
) (124)
≤ 2 ( − α 2ε2
4L4η2C21 (L 2 + α)2(ρ+ 1)T
) , (125)
where (t∗1, t ∗ 2) def = arg supt1,t2≤T,|t1−t2|≤ρ ∣∣∣∑bNt2c`=bNt1cMN` ∣∣∣. We first compute a bound for the last term of equation 114 using the definition of the scaled term RNt from equation 89c. We have
|RNt1 −R N t2 | = ∣∣∣∣∣∣ bNt1c∑ `=0 RN` − bNt2c∑ `=0 RN` ∣∣∣∣∣∣ =
∣∣∣∣∣∣ bNt2c∑ `=bNt1c R` ∣∣∣∣∣∣ ≤ bNt2c∑ `=bNt1c |R`|
(a) ≤ |bNt2c − bNt1c| C0 N2 (ηL2 + (L4/α))
(b) ≤ s1|t2 − t1|, (126)
where (a) follows from the upper bound in equation 178 of Section D, and in (b) we define s1 def = C0 N (ηL2 + (L4/α)).
xxxiv
Putting together equation 117, equation 122, and equation 126, we conclude from Inequality equation 114 that
IP (
sup t1,t2≤T,|t1−t2|≤ρ
|H(µNt1)−H(µ N t1)| ≥ σ
) ≤ IP ( sup
t1,t2≤T,|t1−t2|≤ρ |MNt1 −M N t2 |+ (s0 + s1)ρ ≥ σ ) ≤ 2 ( − α 2(σ − (s0 + s1)ρ)2
4L4η2C21 (L 2 + α)2(ρ+ 1)T
) .
Therefore, condition (T.2) is also satisfied. Since the sufficient conditions (T.1) and (T.2) are satisfied, the condition (J.2) is satisfied. This completes the tightness proof of the measured-valued sequence {µNt }N∈IN. Now, we prove its convergence to a mean-field solution (µ∗t )0≤t≤T . Theorem C.15. (PROKHOROV’S THEOREM PROKHOROV (1956)) A subset of probability measures on a complete separable metric space is tight if and only if it is pre-compact.
According to Theorem C.15, the tightness of the Skorkhod Space DM(IRD)([0, T ]) implies its precompactness which in turn implies the existence of a converging sub-sequence {(µNt )0≤t≤T }Nk of {µNt }N∈IN . Notice that {(µNt )0≤t≤T }Nk is a stochastic process defined on the Skorkhod space. Therefore, let πNk denotes the law of the converging sub-sequence {(µNt )0≤t≤T }Nk . By definition, πNk is an element of the measure spaceM(D[0,T ](M(IRD))). In the sequel, we closely follow the argument of (Wang et al., 2017, Proposition 4) to show that the limiting measure π∞ is a Dirac’s delta function concentrated at a mean-field solution µ∗t ∈ D[0,T ](M(IR
D)). We define the following functional
Ft : D[0,T ](M(IRD))→ IR, µt 7→ Ft[µt] = ∣∣∣∣〈µt, f〉 − 〈µ0, f〉 − ∫ t
0
R[µs]ds ∣∣∣∣ . (127) We compute the expectation of the functional Ft with respect to πNk . We then have
IEπNk [Ft(µ)] = IE[Ft[µ N t ]]
= IE [∣∣∣∣〈µNkt , f〉 − 〈µN0 , f〉 − ∫ t 0 R[µNks ]ds ∣∣∣∣ .] . (128) Now, from Equation equation 93, we have that
〈µNkt , f〉 − 〈µ Nk 0 , f〉 − ∫ t 0 R[µNks ]ds =M Nk t +R Nk t . (129)
Plugging equation 129 in equation 128 gives
IEπNk [Ft(µ)] = IE[Ft[µ Nk t ]] = IE [∣∣∣MNkt +RNkt ∣∣∣]
≤ IE [
sup 0≤t≤T
|MNkt | ] + IE [
sup 0≤t≤T
|RNkt | ]
= 1 Nαρ 4 √
2L2 √ bNT cηC1(L2 + α)2 + C0T
N (ηαL2T + 2ηL4T ), (130)
where the last equality is due to the bounds in equation 94 and equation 95 of Lemmata C.9 and C.10, respectively. Taking the limit of N →∞ from equation 130 yields
lim Nk→∞ IEπNk [|Ft[µ]|] = 0. (131)
xxxv
It can be shown that the functional Ft[·] is continuous and bounded. Therefore, due the weak convergence of the sequence {πNk}Nk∈IN to π∞, equation 131 implies that
IEπ∞ [|Ft(µ)|] = 0. (132)
Since the identity equation 132 holds for all bounded test functions f ∈ C3b (IR D) and for all t ∈ (0, T ], it follows that π∞ is a Dirac’s delta function concentrated at a solution (µ∗t )0≤t≤T of the mean-field equation.
Step 3: Uniqueness of a mean-field solution: Before we establish the uniqueness result we make two remarks:
First, we make it clear that from the compact-containment condition (J.1) of Jakubowski’s criterion in Theorem C.12, the support of the measured-valued process (µNt )0≤t≤T = (µ̂ N bNtc)0≤t≤T is compact for all 0 ≤ t ≤ T . Moreover, in Step 2 of the proof, we established that the measure valued process (µNt )0≤t≤T converges weakly to a mean-field solution as the number of particles tends to infinity (i.e., N → ∞). Thus, all the possible solutions of the mean-field equation also have compact supports. Let Ξ̂ ⊂ IRD denotes a compact set containing the supports of all such solutions at 0 ≤ t ≤ T . In the sequel, it suffices to establish the uniqueness of the mean-field solution for the test functions with a compact domain, i.e., let f ∈ C3b (Ξ̂).
Second, for all bounded continuous test functions f ∈ C3b (Ξ̂), the operator f → 〈µt, f〉 is a linear operator with µt(IRD) = 1. Hence, from Riesz-Markov-Kakutani representation theorem Rudin (1987); Varadarajan (1958) by assuming µt ∈ M(IRD) , existence of unique operator implies f 7→ 〈f, µt〉 implies the existence of the unique probability measure µt. Now, we equip the measure spaceM(IRD) with the following norm
‖µ‖ def= sup f∈C3b (Ξ̂) ‖f‖∞ 6=0 |〈f, µ〉| ‖f‖∞ . (133)
Given an initial measure µ0, we next prove that there exists at most one mean-field model solution by showing that there exists at most one real valued process 〈µt, f〉 corresponding to the mean-field model. Suppose (µ∗,1t )0≤t≤T , (µ ∗,2 t )0≤t≤T are two solutions satisfying the mean-field equations equation 97 with the initial distributions µ10, µ 2 0 ∈ M(IR
D), respectively. For any test function f ∈ C3b (Ξ̂) we have that
〈µ∗,1t − µ∗.2t , f〉 = 〈µ10 − µ20, f〉+ η
α ∫ t 0 (∫∫ X×Y ( 〈ϕ(x, ξ)ϕ(x̃, ξ), µ∗,1s − µ∗,2s 〉 − αyỹ
) (134)
× 〈∇f(ξ)(∇ξ(ϕ(x̃; ξ)ϕ(x; ξ)))T , µ∗,1s − µ∗,2s 〉P⊗2x,y((dz,dz̃)
) ds.
We bound the first term on the right side of Equation equation 134 as follows
〈µ10 − µ20, f〉 ≤ ‖µ10 − µ20‖ · ‖f‖∞ (135) ≤ b‖µ10 − µ20‖, (136)
where used the definition of the norm ‖ · ‖ on the measure spaceM(IRD) from equation 133.
xxxvi
Furthermore, let∫∫ X×Y αyỹ〈∇f(ξ)(∇ξ(ϕ(x̃; ξ)ϕ(x; ξ)))T , µ∗,1s − µ∗,2s 〉P⊗2x,y(d(x, y),d(x̃, ỹ))
≤ ∫∫ X×Y α|yỹ| · |〈∇f(ξ)(∇ξ(ϕ(x̃; ξ)ϕ(x; ξ)))T , µ∗,1s − µ∗,2s 〉|P⊗2x (d(x, y),d(x̃, ỹ))
≤ α‖µ∗,1s − µ∗,2s ‖ ∫ X ‖∇f(ξ)(∇ξ(ϕ(x̃; ξ)ϕ(x; ξ)))T ‖P⊗2x (dx,dx̃)
≤ α‖µ∗,1s − µ∗,2s ‖ ∫ X ‖∇f(ξ)‖∞ · ‖∇ξϕ(x̃; ξ)ϕ(x; ξ)‖∞P⊗2x (dx,dx̃) ≤ αL2C1‖µ∗,1s − µ∗,2s ‖, (137)
where in the last inequality, we used the fact that ‖∇f(ξ)‖ ≤ C1 since the test function is threetimes continuously differentiable f ∈ C3b (Ξ̂) on a compact support. Similarly, we have∫
X 〈ϕ(x, ξ)ϕ(x̃, ξ), µ∗,1s − µ∗,2s 〉〈∇f(ξ)(∇ξϕ(x, ξ)ϕ(x, ξ)), µ∗,1s − µ∗,2s 〉P⊗2x (dx,dx̃) ≤ ‖µ∗,1s − µ∗,2s ‖2 ∫ X ‖ϕ(x, ξ)ϕ(x̃, ξ)‖∞‖∇f(ξ)(∇ξϕ(x, ξ)ϕ(x, ξ))T ‖∞P⊗2x (dx,dx̃)
≤ L4C1‖µ∗,1s − µ∗,2s ‖2. (138) Putting together the inequalities in equation 136,equation 137, and equation 138 yield
〈µ∗,1t − µ ∗,2 t , f〉 ≤ b‖µ 1 0 − µ20‖+ L2C1η ∫ t 0 ‖µ∗,1s − µ∗,2s ‖ds+ ηL4C1 α ∫ t 0 ‖µ∗,1s − µ∗,2s ‖2ds. (139)
The above inequality holds for all bounded functions f ∈ C3b (Ξ̂). Thus, by taking the supremum with respect to f we obtain
‖µ∗,1t − µ ∗,2 t ‖ = sup f∈C3b (Ξ̂) 〈µ∗,1t − µ ∗,2 t , f〉 (140)
≤ b‖µ10 − µ20‖+ L2C1η ∫ t
0
‖µ∗,1s − µ∗,2s ‖ds+ L4C1η
α ∫ t 0 ‖µ∗,1s − µ∗,2s ‖2ds.
(141)
Now, we employ the following result which generalizes Gronewall’s inequality when higher order terms are involved: Lemma C.16. (EXTENDED GRONEWALL’S INEQUALITY, (WEBB, 2018, THM 2.1.)) Let p ∈ IN and suppose that for a.e. t ∈ [0, T ], u ∈ L∞+ [0, T ] satisfies
ut ≤ c0(t) + ∫ t
0
(c1(s)us + c2(s)u 2 s + · · ·+ cp+1(s)up+1s )ds, (142)
where c0 ∈ L∞[0, T ] is non-decreasing, and cj ∈ L1+[0, T ] for j ∈ {1, · · · , p+ 1}. Then, if∫ T 0 cj+1(s)u j sds ≤Mj , j ∈ {1, 2, · · · , p}. (143)
It follows that for a.e. t ∈ [0, T ] ut ≤ c0(t) (∫ t
0
c1(s)ds ) (M1 + · · ·+Mp). (144)
xxxvii
We now apply the extended Gronewall’s Inequality equation 142 with p = 1, c0(t) = b‖µ10 − µ20‖, 0 ≤ t ≤ T , c1(t) = ηL2C1, 0 ≤ t ≤ T , c2(t) = ηL4C1 α , 0 ≤ t ≤ T , and us = ‖µ ∗,1 s −µ∗,2s ‖. In this case, it is easy to see that M1 = 2bηTL 4C1
α . Consequently, from equation 140 and equation 144, we obtain that
‖µ∗,1t − µ ∗,2 t ‖ ≤ b‖µ10 − µ20‖ · ( ηL2C1t+
2bTηL4C1 α
) , 0 ≤ t ≤ T. (145)
Thus, starting from an initial measure µ10 = µ 2 0 = µ0, there exists at most one solution for the mean-field model equations equation 97.

heading: C.3 PROOF OF COROLLARY 4.2.1
text: To establish the proof, we recall from equation 88 that
〈f, µ̂Nm〉 − 〈f, µ̂N0 〉 = m−1∑ `=0 DN` + m−1∑ `=0 MN` + m−1∑ `=0 RN` , (146)
for all f ∈ Cb(IRD). Recall the definition of the total variation distance TV(·, ·) on a metric space (X , d).
TV(µ, ν) = 1
2 sup A⊂X
|µ(A)− ν(A)|. (147)
The total variation distance admits the following variational form TV(µ, ν) = sup
f :‖f‖∞≤ 12 〈f, µ〉 − 〈f, ν〉. (148)
Now, using the variation form and by taking the supremum of equation 146 with respect to the functions from the function class Fc def = {f ∈ C1/2(IRD)}, we obtain the following upper bound on the TV distance
TV(µ̂Nm, µ̂ N 0 ) ≤
1
2 m−1∑ `=0 |DN` |+ 1 2 m−1∑ `=0 |MN` |+ 1 2 m−1∑ `=0 |RN` |. (149)
Based on the upper bound equation 178 on the remainder term, we have
|RN` | ≤ C0 N2
( ηL2 + 2 η α L4 ) , ` ∈ [0,m− 1], (150)
for some constant C0 > 0. Moreover, from the concentration inequality equation 189, we also have that with the probability of at least 1− δ, the following inequality holds
|MN` | ≤ 8 √ `ηL2C1(L 2 + α)
Nα log
( 2
δ
) . (151)
Lastly, recall the definition of the drift term in equation 86a. By carrying out a similar bounding method leading to equation 116, it can be shown that
|DN` | ≤ 2η
Nα (L2 + α)L2C1. (152)
By plugging equation 150, equation 151, and equation 152 into equation 149, we derive that TV ( µ̂Nm, µ̂ N 0 ) ≤ mηC0
2N2
( L2 + 2 L4
α
) + 8m √ mηL2C1(L 2 + α)
Nα log
( 2
δ
) + mη
Nα (L2 + α)L2C1,
(153)
with the probability of 1− δ. We now leverage the following lemma:
xxxviii
Lemma C.17. (BOUNDED EQUIVALENCE OF THE WASSERSTEIN AND TOTAL VARIATION DISTANCES, SINGH & PÓCZOS (2018)) Suppose (X , d) is a metric space, and suppose µ and ν are Borel probability measures on X with countable support; i.e., there exists a countable set X ′ ⊆ X such that µ(X ′) = ν(X ′) = 1. Then, for any p ≥ 1, we have
Sep(X ′)(2TV(µ, ν)) 1 p ≤Wp(µ, ν) ≤ Diam(X ′)(2TV(µ, ν)) 1 p , (154)
where Diam(X ′) def= supx,y∈X ′ d(x, y), and Sep(X ′) def = infx 6=y∈X ′ d(x, y).
Consider the metric space (IRD, ‖ · ‖2). Note that the empirical measures µ̂Nm, µ̂N0 have a countable support X ′ = {ξkm}Nk=1 ∪ {ξk0}Nk=1 ⊂ IR
D. Therefore, using the upper bounds in equation 154 of Lemma C.17 and 153, we conclude that when the step-size is of the order
η = O (
Rp
T √ NT log(2/δ)
) , (155)
then Wp(µ̂Nm, µ̂ N 0 ) ≤ R for all m ∈ [0, NT ] ∩ IN.

heading: D PROOFS OF AUXILIARY RESULTS

heading: D.1 PROOF OF LEMMA C.5
text: The upper bound follows trivially by letting x = y in the optimization problem equation 44.
Now, consider the lower bound. Define the function g : [0, 1] → IR, t 7→ g(t) = f(y + t(x − y)). Then, when f is differentiable, we have g′(t) = 〈x − y,∇f(y + t(x − y))〉. In addition, g(0) = f(y), and g(1) = f(x). Based on the basic identity g(1) = g(0) + ∫ 1 0 g′(s)ds, we derive
f(x) = f(y) + ∫ 1 0 〈x− y,∇f(y + s(x− y))〉ds
≥ f(y)− ‖x− y‖2 ∫ 1
0
‖∇f(y + s(x− y))‖2ds, (156)
where the last step is due to the Cauchy-Schwarz inequality. Using Inequality equation 156 yields the following lower bound on Moreau’s envelope
Mβf (y) ≥ f(y) + infx∈X
{ 1
2β ‖x− y‖22 − ‖x− y‖2 ∫ 1 0 ‖∇f(y + s(x− y))‖2ds }
= f(y) + inf x∈X {( 1√ 2β ‖x− y‖2 − √ β 2 ∫ 1 0 ‖∇f(y + s(x− y))‖2ds )2
− β 2 (∫ 1 0 ‖∇f(y + s(x− y))‖2ds )2}
≥ f(y)− β 2 sup x∈X (∫ 1 0 ‖∇f(y + s(x− y))‖2ds )2 (a) ≥ f(y)− β 2 sup x∈X ∫ 1 0 ‖∇f(y + s(x− y))‖22ds
≥ f(y)− β 2 ∫ 1 0 sup x∈X ‖∇f(y + s(x− y))‖22ds, (157)
where (a) is due to Jensen’s inequality.
xxxix

heading: D.2 PROOF OF LEMMA C.8
text: Let u ∈ IRd denotes an arbitrary unit vector ‖u‖2 = 1. From the definition of the gradient of a function, we have that 〈
∇θMβf(·;θ)(x),u 〉 = Du[Mf(·;θ)(x)], (158)
where Du[Mf(·;θ)(x)] is the directional derivative
Du[Mf(·;θ)(x)] def = lim
δ→0 Mf(·;θ+δu)(x)−Mf(·;θ)(x) δ . (159)
We now have
Mf(·;θ+δu)(x) = inf x∈X
{ 1
2β ‖x− y‖22 + f(y;θ + δu) } = inf x∈X { 1 2β ‖x− y‖22 + f(y;θ) + δ〈∇θf(y;θ),u〉 } +O(δ2)
(a) ≤ 1 2β ‖x− Proxf(·;θ)(x)‖22 + f(Proxf(·;θ)(x);θ)
+ δ〈∇θf(Proxf(·;θ)(x);θ),u〉+O(δ2), (160) where the inequality in (a) follows by letting y = Proxf(·;θ)(x) in the optimization problem. Now, recall that
Mf(·;θ)(x) = inf y∈X
{ 1
2β ‖x− y‖22 + f(y;θ)
} ,
Proxf(·;θ)(x) = arg min y∈X
{ 1
2β ‖x− y‖22 + f(y;θ)
} .
Therefore,
Mf(·;θ)(x) = 1
2β ‖x− Proxf(·;θ)(x)‖22 + f(Proxf(·;θ)(x);θ). (161)
Substitution of equation 161 in equation 160 yields Mf(·;θ+δu)(x) ≤Mf(·;θ)(x) + δ〈∇θf(Proxf(·;θ)(x);θ),u〉+O(δ2). (162)
Hence, Du[Mf(·;θ)(x)] ≤ 〈∇θf(Proxf(·;θ)(x);θ),u〉. From equation 158 and by using CauchySchwarz inequality, we compute the following bound on the inner product of the gradient with the unit vectors u ∈ IRd, ‖u‖2 = 1,〈 ∇θMβf(·;θ)(x),u 〉 ≤ ‖∇θf(Proxf(·;θ)(x);θ)‖2 · ‖u‖2 = ‖∇θf(Proxf(·;θ)(x);θ)‖2. (163)
Since the preceding upper bound holds for all the unit vectors u ∈ IRd, we let u = ∇θMβf(·;θ)(x) ‖∇θMβf(·;θ)(x)‖2 to get Inequality equation 65.

heading: D.3 PROOF OF LEMMA C.6
text: Let z ∈ Sd−1 denotes an arbitrary vector on the sphere. Define the random variable Qz ( (y1,x1), · · · , (yn,xn) ) def = 〈z,∇En(ξ)〉
= 1 n(n− 1) ∑ i 6=j yiyj ( ϕ(xi; ξ)〈z,∇ξϕ(xj ; ξ)〉+ ϕ(xj ; ξ)〈z,∇ξϕ(xi; ξ)〉 ) − IE
P⊗2x,y
[ yŷ ( ϕ(xi; ξ)〈z,∇ξϕ(xj ; ξ)〉+ ϕ(xj ; ξ)〈z,∇ξϕ(xi; ξ)〉 )] .
(164)
xl
Clearly, IEPx,y [Qz] = 0. Now, let (ŷm, x̂m) ∈ Y × X , 1 ≤ m ≤ n. By repeated application of the triangle inequality, we obtain that∣∣∣Qz((y1,x1), · · · , (ym,xm), · · · , (yn,xn))−Qz((y1,x1), · · · , (ŷm, x̂m), · · · , (yn,xn))∣∣∣
≤ 1 n(n− 1) ∣∣∣∑ i 6=m yiϕ(xi; ξ)〈z, ym∇ξϕ(xm; ξ)− ŷm∇ξϕ(x̂m; ξ)〉 ∣∣∣
+ 1
n(n− 1) ∣∣∣∑ i 6=m yi〈z,∇ξϕ(xi; ξ)〉(ymϕ(xm; ξ)− ŷmϕ(x̂m; ξ)) ∣∣∣
≤ 1 n(n− 1) ∑ i 6=m |ϕ(xi; ξ)| · ‖z‖2 · ‖ym∇ξϕ(xm; ξ)− ŷm∇ξϕ(x̂m; ξ)‖2
+ 1 n(n− 1) ∑ i 6=m ‖z‖2 · ‖∇ξϕ(xi; ξ)‖2 · |ymϕ(xm; ξ)− ŷmϕ(x̂m; ξ)|
≤ 4L 2
n , (165)
where the last inequality is due to assumption (A.2) and the fact that ‖z‖2 = 1 for z ∈ Sd−1. In particular, to derive Inequality equation 165, we employed the following upper bounds
|ϕ(xi; ξ)| ≤ L, |ymϕ(xm; ξ)− ŷmϕ(x̂m; ξ)| ≤ |ϕ(xm; ξ)|+ |ϕ(x̂m; ξ)| ≤ 2L,
‖∇ξϕ(xi; ξ)‖2 ≤ L, ‖ym∇ξϕ(xm; ξ)− ŷm∇ξϕ(x̂m; ξ)‖2 ≤ ‖∇ξϕ(xm; ξ)‖2 + ‖∇ξϕ(x̂m; ξ)‖2 ≤ 2L.
Using McDiarmid Martingale’s inequality McDiarmid (1989) then gives us IP (∣∣∣Qz((y1,x1), · · · , (yn,xn))∣∣∣ ≥ u) ≤ 2(− nu2
16L4
) , (166)
for x ≥ 0. Now, for every p ∈ IN, the 2p-th moment of the random variable Qz is given by IE [ Q2pz ((y1,x1), · · · , (yn,xn)) ] = ∫ IR+ 2pu2p−1IP(Qz((y1,x1), · · · , (yn,xn)) ≥ u)du
(a) ≤ ∫
IR+
4pu2p−1 ( − u 2
16nL4
) du
= 2(16L4/n)2pp!, (167)
where (a) is due to the concentration bound in equation 166. Now Therefore,
IE [( Q2z((y1,x1), · · · , (yn,xn))/γ2 )] = ∞∑ p=0 1 p!γ2p IE [ φ2pz ((y1,x1), · · · , (yn,xn)) ] = 1 + 2
∑ p∈IN ( 16L4 nγ )2p = 2
1− (16L4/nγ)2 − 1.
For γ = 16 √ 3L4/n, we obtain IE [( Q2z((y1,x1), · · · , (yn,xn))/γ2 )] ≤ 2. Therefore, ‖Qz‖ψ2 = ‖〈z,∇En(ξ)〉‖ψ2 ≤ 16 √ 3L4/n for all z ∈ Sn−1 and ξ ∈ IRD. Consequently, by the definition
xli
of the sub-Gaussian random vector in equation 35 of Definition C.2, we have ‖∇En(ξ)‖ψ2 ≤ 16 √
3L4/n for every ξ ∈ IRD. We invoke the following lemma proved by the first author in (Khuzani & Li, 2017, Lemma 16):
Lemma D.1. (THE ORLICZ NORM OF THE SQUARED VECTOR NORMS, (KHUZANI & LI, 2017, LEMMA 16)) Consider the zero-mean random vector Z satisfying ‖Z‖ψν ≤ β for every ν ≥ 0. Then, ‖‖Z‖22‖ψ ν
2 ≤ 2 · 3 2ν · β2.
Using Lemma D.1, we now have that ‖‖∇En(ξ)‖22‖ψ1 ≤ 4608L4/n2 for every ξ ∈ IR D. Applying the exponential Chebyshev’s inequality with β = 4608L4/n2 yields
IP (∫ IRD ∫ 1 0 ∣∣∣‖∇En((1− s)ξ + sζ∗)‖22 − IEx,y[‖∇En((1− s)ξ + sζ∗)‖22]∣∣∣µ0(dξ) ≥ δ )
≤ e− n2δ
4608L4 IEx,y
[ e ( n2 4608L4 ∫ IRD ∫ 1 0 ∣∣‖∇En((1−s)ξ+sζ∗)‖22−IEx,y [‖∇En((1−s)ξ+sζ∗)‖22]∣∣dsµ0(dξ))] (a) ≤ e− n2δ 4608L4
∫ IRD ∫ 1 0 IEx,y [ e n2 4608L4 (|‖∇En((1−s)ξ+sζ∗)‖22−IEx,y [‖∇En((1−s)ξ+sζ∗)‖ 2 2])| ] dsµ0(dξ)
(b) ≤ 2e− n2δ 4608L4 ,
where (a) follows by Jensen’s inequality, and (b) follows from the fact that
IEx,y
[ e n2 4608L4 (|‖∇En((1−s)ξ+sζ∗)‖22−IEx,y [‖∇En((1−s)ξ+sζ∗)‖ 2 2])| ] ≤ 2, (168)
due to Definition C.1. Therefore,
IP (∫ IRD ∫ 1 0 ‖∇En((1− s)ξ + sζ∗)‖22dsµ(dξ) ≥ δ )
≤ 2 ( − n2(δ − ∫ 1 0 ∫ IRD
IEx,y[‖∇En((1− s)ξ + sζ∗)‖22]dsµ0(dξ)) 4608L4
) . (169)
It now remains to compute an upper bound on the expectation IEx,y[‖∇En((1− s)ξ+ sζ∗)‖22]. But this readily follows from equation 167 by letting p = 1 and z = ∇En((1−s)ξ+sζ∗)‖∇En((1−s)ξ+sζ∗)‖2 as follows
IEx,y[‖∇En((1− s)ξ + sζ∗)‖22] = IEx,y [〈 ∇En((1− s)ξ + sζ∗) ‖∇En((1− s)ξ + sζ∗)‖2 ,∇En((1− s)ξ + sζ∗) 〉2]
= IEx,y [ Q2 ∇En((1−s)ξ+sζ∗) ‖∇En((1−s)ξ+sζ∗)‖2 ] ≤ 29L 8
n2 . (170)
Plugging the expectation upper bound of equation 170 into equation 169 completes the proof of the first part of Lemma C.6.
The second part of Lemma C.6 follows by a similar approach and we thus omit its proof.
xlii

heading: D.4 PROOF OF LEMMA C.7
text: Let W∗ def = arg minW∈W Ψ(W ) and W def = arg minW∈W Φ(W ). Then, since |Ψ(W ) − Φ(W )| ≤ δ for allW ∈ W , we have that
|Ψ(W∗)− Φ(W∗)| = ∣∣∣∣ minW∈WΨ(W )− Φ(W∗) ∣∣∣∣ ≤ δ. (171) Therefore,
min W∈W Φ(W ) ≤ Φ(W∗) ≤ min W∈W Ψ(W ) + δ. (172)
Similarly, it can be shown that
min W∈W Ψ(W ) ≤ Ψ(W ) ≤ min W∈W Φ(W ) + δ. (173)
Combining equation 172 and equation 173 yields the desired inequality.

heading: D.5 PROOF OF LEMMA C.10
text: We recall the expression of the remainder term {RNm}0≤m≤NT from equation 84. For each 0 ≤ m ≤ NT , N ∈ IN, we can bound the absolute value of the remainder term as follows
∣∣RNm∣∣ = 1N ∣∣∣∣∣ N∑ k=1 (ξkm+1 − ξkm)∇2f(ξ̃k)(ξkm+1 − ξkm)T ∣∣∣∣∣
≤ 1 N N∑ k=1 ∣∣∣(ξkm+1 − ξkm)∇2f(ξ̃k)(ξkm+1 − ξkm)T ∣∣∣ ≤ 1 N N∑ k=1 ‖ξkm+1 − ξkm‖22 · ∥∥∇2f(ξ̃k)∥∥ F . (174)
Next, we characterize a bound on the difference term ‖ξkm+1− ξkm‖2. To attain this goal, we use the iterations of the particle SGD in Equation equation 15. We have that
‖ξkm+1 − ξkm‖2
≤ η N ∥∥∥∥∥ ( ymỹm − 1 Nα N∑ k=1 ϕ(xm; ξ k m)ϕ(x̃m; ξ k m) ) ∇ξ ( ϕ(xm; ξ k m)ϕ(x̃m; ξ k m) )∥∥∥∥∥
2
≤ η N |ymỹm|
( |ϕ(xm; ξkm)| · ‖∇ξϕ(x̃m; ξkm)‖2 + |ϕ(x̃m; ξkm)| · ‖∇ξϕ(xm; ξkm)‖2 ) + η
N
( 1
Nα N∑ k=1
∣∣∣ϕ(xm; ξkm)ϕ(x̃m; ξkm)∣∣∣ )( |ϕ(xm; ξkm)|‖∇ξϕ(x̃m; ξkm)‖2 + |ϕ(x̃m; ξkm)|‖∇ξϕ(xm; ξkm)‖2 ) (a)
≤ ηL 2
N +
2ηL4
Nα , (175)
where in (a), we used the fact that ‖ϕ‖∞ < L and ‖∇ξϕ(x, ξ)‖2 < L due to (A.1), and ym, ỹm ∈ {−1, 1}. Plugging the last inequality in equation 175 yields
|RNm| ≤ 1
N3
( ηL2 + 2ηL4
α ) N∑ k=1 ‖∇2f(ξ̃k)‖F . (176)
xliii
We next compute an upper bound on the Frobenious norm of the Hessian matrix ∇2f(ξ̃k). To this end, we first show that there exists a compact set C ⊂ IRD such that ξkm ∈ C for all k = 1, 2, · · · , N and all m ∈ [0, NT ] ∩ IN. For each k = 1, 2, · · · , N , from Inequality equation 175 we obtain that
‖ξkm‖2 ≤ ‖ξkm−1‖2 + ηL2
N +
2ηL4
Nα
= ‖ξk0‖2 + mηL2
N +
2mηL4
Nα
≤ ‖ξk0‖2 + ηL2T + 2(η/α)L4T. (177)
Now, ‖ξk0‖2 < c0 for some constant c0 > 0 since the initial samples ξ10, · · · , ξN0 are drawn from the measure µ0 whose support support(µ0) = Ξ is assumed to be compact by (A.3). From upper bound in equation 177, it thus follows that ‖ξkm‖2 < C for some constant C > 0, for all m ∈ [0, NT ]∩ IN. Now, recall that ξ̃k = (ξ̃k(1), · · · , ξ̃k(p)), where ξ̃k(i) ∈ [ξkm(i), ξkm+1(i)], i = 1, 2, · · · ,m + 1, for i = 1, 2, · · · , p. Therefore, ξ̃k ∈ C. Since all the test function f ∈ C3b (IR 3) are three-times continuously differentiable, it follows that there exists a constant C0 def = C0(T ) > 0 such that supξ̃∈C ‖∇ 2f(ξ̃)‖F < C0. From Inequality equation 176, it follows that
|RNm| ≤ C0 N2
( ηL2 + 2ηL4
α
) , m ∈ [0, NT ] ∩ IN. (178)
Now, recall the definition of the scaled term RNt from equation 89c. Using the Inequality equation 178 as well as the definition ofRNt , we obtain
sup 0≤t≤T
|RNt | ≤ C0T
N
( ηL2 + 2ηL4
α
) . (179)

heading: D.6 PROOF OF LEMMA C.10
text: Let Fm−1 = σ((xk, yk)0≤k≤m−1, (x̃k, ỹk)0≤k≤m−1) denotes the σ-algebra generated by the samples up to time m− 1. We define F−1 def = ∅. Further, define the following random variable
∆Nm def = ( 〈ϕ(xm, ξ)ϕ(x̃m, ξ), µ̂Nm〉 − αymỹm ) × 〈∇f(ξ)∇ξ(ϕ(x̃m; ξ)ϕ(xm; ξ)), µ̂Nm〉. (180)
Notice that 1N IE[∆ N m|Fm−1] = DNm . We now rewrite the martingale term in equation 86b in term of ∆Nm,
MNm def =
η
Nα m∑ `=0 (∆N` − IE[∆N` |F`−1]), (181)
with MN0 = 0.
By construction of MNm in equation 181, it is a Martingale IE[M N m |Fm−1] = MNm−1. We now prove that MNm has also bounded difference. To do so, we define the shorthand notations
aNm def = 〈ϕ(xm, ξ)ϕ(x̃m, ξ), µ̂Nm〉 − αymỹm, (182)
bNm def = 〈∇f(ξ)(∇ξ(ϕ(x̃m; ξ)ϕ(xm; ξ)))T , µ̂Nm〉. (183)
xliv
Then, we compute
|MNm −MNm−1| = η
Nα ∣∣∆Nm − IE[∆Nm|Fm−1]∣∣ ≤ η Nα |∆Nm|+ η Nα IE[|∆Nm||Fm−1] ≤ η Nα |aNm| · |bNm|+ η Nα IE [ |aNm| · |bNm||Fm−1 ] . (184)
For the difference terms, we derive that |aNm| = ∣∣∣〈ϕ(xm, ξ)ϕ(x̃m, ξ), µ̂Nm〉 − αymỹm∣∣∣
≤ 1 N N∑ k=1 ∣∣ϕ(xm, ξkm)ϕ(x̃m, ξkm)∣∣+ α|ymỹm| ≤ L2 + α, (185)
where the last step follows from the fact that ‖ϕ‖∞ ≤ L due to (A.1). Similarly, we obtain that |bNm| = |〈∇f(ξ)(∇ξ(ϕ(x̃m; ξ)ϕ(xm; ξ)))T , µ̂Nm〉|
≤ 1 N N∑ k=1 |ϕ(x̃m; ξkm)| · ∣∣∇f(ξkm)(∇ξϕ(x̃m; ξkm))T ∣∣
+ 1
N N∑ k=1 |ϕ(xm; ξkm)| · |∇f(ξkm)(∇ξϕ(x̃m; ξkm))T |
(a) ≤ L N N∑ k=1 ∣∣∇f(ξkm)(∇ξϕ(x̃m; ξkm))T ∣∣+ LN N∑ k=1 |∇f(ξkm)(∇ξϕ(x̃m; ξkm))T |
(b) ≤ L N N∑ k=1 ‖∇f(ξkm)‖2 · ‖∇ξϕ(x̃m; ξkm)‖2 + L N N∑ k=1 ‖∇f(ξkm)‖2 · ‖∇ξϕ(x̃m; ξkm)‖2
(c) ≤ 2L 2
N N∑ k=1 ‖∇f(ξkm)‖2, (186)
where (a) and (c) follows from (A.1), and (b) follows from the Cauchy-Schwarz inequality. From Inequality equation 177 and the ensuing disucssion in Appendix D.5, we recall that ‖ξkm‖2 < C for some constant and for all m ∈ [0, NT ] ∩ IN, and k = 1, 2, · · · , N . For the two times continuously test function f ∈ C3b (IR
D), it then follows that |∇f(ξkm)‖2 ≤ C1 for some constant C1 > 0. The following bound can now be computed from equation 186,
|bNm| ≤ 2C1L
2
N . (187)
Plugging the upper bounds on |aNm| and |bNm| from equation 185-equation 187 into equation 184 we obtain that
|MNm −MNm−1| ≤ 4ηC1 Nα L2(L2 + α). (188)
Thus, (MNm )m∈[0,NT ]∩IN is a Martingale process with bounded difference. From the AzumaHoeffding inequality it follows that
IP(|MNm | ≥ ε) = IP(|MNm −MN0 | ≥ ε) ≤ 2 ( − N 2α2ε2
8mL4η2C21 (L 2 + α)2
) , ∀m ∈ [0, NT ] ∩ IN.
(189)
xlv
Therefore, sinceMNt = MNbNtc, we have IP(|MNT | ≥ ε) ≤ 2 ( − N 2α2ε2
8L4bNT cη2C21 (L2 + α)2
) . (190)
Then,
IE [ |MNT | ] = ∫ ∞ 0 IP(|MNT | ≥ ε)dε
≤ 2 ∫ ∞
0
( − N 2α2ε2
8L4bNT cη2C21 (L2 + α)2 ) = 1
Nα 4 √
2L2 √ bNT cηC1(L2 + α)2. (191)
where the inequality follows from equation 190.
By Doob’s Martingale inequality Doob (1953), the following inequality holds
IP ( sup
0≤t≤T |MNt | ≥ ε
) ≤ IE[|M N T |]
ε (192)
≤ 1 Nαε
4 √ 2L2 √ bNT cηC1(L2 + α)2. (193)
In particular, with the probability of at least 1− ρ, we have
sup 0≤t≤T
|MNt | ≤ 1 Nαρ 4 √
2L2 √ bNT cηC1(L2 + α)2. (194)

heading: E CHAOTICITY AND PROPAGATION OF CHAOS IN PARTICLE SGD
text: In this appendix, we establish the so called ‘propagation of chaos’ property of particle SGD. We now establish the so called ‘propagation of chaos’ property of particle SGD. At a high level, the propagation of chaos means that when the number of samples {ξk}Nk=1 tends to infinity (N → +∞), their dynamics are decoupled. Definition E.1. (EXCHANGABLITY) Let ν be a probability measure on a Polish space S and. For N ∈ IN, we say that ν⊗N is an exchangeable probability measure on the product space Sn if it is invariant under the permutation π def= (π(1), · · · , π(N)) of indices. In particular,
ν⊗N (π ·B) = ν⊗N (B), (195) for all Borel subsets B ∈ B(Sn).
An interpretation of the exchangablity condition equation 195 can be provided via De Finetti’s representation theorem which states that the joint distribution of an infinitely exchangeable sequence of random variables is as if a random parameter were drawn from some distribution and then the random variables in question were independent and identically distributed, conditioned on that parameter.
Next, we review the mathematical definition of chaoticity, as well as the propagation of chaos in the product measure spaces: Definition E.2. (CHAOTICITY) Suppose ν⊗N is exchangeable. Then, the sequence {ν⊗N}N∈IN is ν-chaotic if, for any natural number ` ∈ IN and any test function f1, f2, · · · , fk ∈ C2b (S), we have
lim N→∞ 〈∏̀ k=1 fk(s k), ν⊗N (ds1, · · · ,dsN ) 〉 = ∏̀ k=1 〈fk, ν〉 (196)
xlvi
According to equation 196 of Definition E.2, a sequence of probability measures on the product spaces S is ν-chaotic if, for fixed k the joint probability measures for the first k coordinates tend to the product measure ν(ds1)ν(ds2) · · · ν(dsk) = ν⊗k on Sk. If the measures ν⊗N are thought of as giving the joint distribution of N particles residing in the space S, then {ν⊗N} is ν-chaotic if k particles out of N become more and more independent as N tends to infinity, and each particles distribution tends to ν. A sequence of symmetric probability measures on SN is chaotic if it is ν-chaotic for some probability measure ν on S.
If a Markov process on SN begins in a random state with the distribution ν⊗N , the distribution of the state after t seconds of Markovian random motion can be expressed in terms of the transition function KN for the Markov process. The distribution at time t > 0 is the probability measure UNt ν ⊗N is defined by the kernel
UNt ν ⊗N (B) def = ∫ SN KN (s,B, t)ν⊗N (ds). (197)
Definition E.3. (PROPOGATION OF CHAOS) A sequence functions{ KN (s,B, t) } N∈IN
(198)
whose N -th term is a Markov transition function on SN that satisfies the permutation condition
KN (s,B, t) = KN (π · s,π ·B, t), (199)
propagates chaos if whenever {ν⊗N}N∈IN is chaotic, so is {UNt } for any t ≥ 0, where UNt is defined in equation 197.
We note that for finite systems size N ,the states of the particles are not independent of each other. However, as we prove in the following result, in the limiting system N → +∞, the particles are mutually independent. This phenomena is known as the propagation of chaos (a.k.a. asymptotic independence): Theorem E.4. (CHAOTICITY IN PARTICLE SGD) Consider Assumptions (A.1)− (A.3). Furthermore, suppose that {ξk0}1≤k≤N ∼i.i.d. µ0 is exchangable in the sense that the joint law is invariant under the permutation of indices. Then, at each time instant t ∈ (0, T ], the scaled empirical measure µNt ∈M(IR D) defined via scaling
µNt (dξ 1, · · · ,dξN ) def= µ̂NbNtc(dξ 1, · · · ,dξN ) = IP{ξ1bNtc ∈ dξ 1, · · · , ξNbNtc ∈ dξ N}, (200)
is µ∗t -chaotic, where µ ∗ t is mean-field solution of equation 97.
Proof. To establish the proof, it suffices to show that for every integer ` ∈ IN, and for all the test functions f1, · · · , fk ∈ C3b (IR D), we have
lim sup N→∞ ∣∣∣∣∣IE [∏̀ k=1 fk(ξ k bNtc) ] − ∏̀ k=1 〈µ∗t , fk〉 ∣∣∣∣∣ = 0. (201) Using the triangle inequality, we now have that∣∣∣∣∣IE [∏̀ k=1 fk(ξ k bNtc) ] − ∏̀ k=1 〈µ∗t , fk〉
∣∣∣∣∣ ≤
∣∣∣∣∣IE [∏̀ k=1 〈µ̂Nt , fk〉 ] − ∏̀ k=1 〈µ∗t , fk〉 ∣∣∣∣∣+ ∣∣∣∣∣IE [∏̀ k=1 〈µ̂Nt , fk〉 ] − IE [∏̀ k=1 fk(ξ k bNtc) ]∣∣∣∣∣ . (202) xlvii
For the first term on the right side of equation 202 we have
lim sup N→∞ ∣∣∣∣∣IE [∏̀ k=1 〈µ̂Nt , fk〉 ] − ∏̀ k=1 〈µ∗t , fk〉 ∣∣∣∣∣ (a)≤ lim supN→∞ IE [∣∣∣∣∣∏̀ k=1 〈µ̂Nt , fk〉 − ∏̀ k=1 〈µ∗t , fk〉 ∣∣∣∣∣ ]
(b) ≤ IE [ lim sup
N→∞ ∣∣∣∣∣∏̀ k=1 〈µ̂Nt , fk〉 − ∏̀ k=1 〈µ∗t , fk〉 ∣∣∣∣∣ ]
(c) ≤ b`−1IE [∑̀ k=1 lim sup N→∞ ∣∣〈µ̂Nt , fk〉 − 〈µ∗t , fk〉∣∣ ]
(d) = 0, (203)
where (a) is by Jensen’s inequality, (b) is by Fatou’s lemma, (c) follows from the basic inequality∣∣∣∏Ni=1 ai −∏Ni=1 bi∣∣∣ ≤ ∑Ni=1 |ai − bi| for |ai|, |bi| ≤ 1, i = 1, 2, · · · , N , as well as the fact that 〈µ∗t , fk〉 ≤ b and 〈µ̂Nt , fk〉 ≤ b for all k = 1, 2, · · · , N due to the boundedness of the test functions f1, · · · , f` ∈ C3b (IR D), and (d) follows from the weak convergence µ̂Nt weakly→ µ∗t to the mean-field solution equation 97.
Now, consider the second term on the right hand side of equation 202. Due to the exchangeability of the initial states (ξk0 )1≤k≤N , the law of the random variables (ξ k 0 )1≤k≤N is also exchangeable. Therefore, we obtain that
IE [∏̀ k=1 fk(ξ k bNtc) ] = `! N ! IE  ∑ π∈Π(`,N) ∏̀ k=1 fk(ξ π(k) bNtc)  , (204) where Π(`,N) is the set of all permutations of ` numbers selected from {1, 2, · · · , N}. Notice that the right hand side of equation 204 is the symmetrized version of the left hand side equation 205.
Further, by the definition of the empirical measure µ̂Nt we obtain that
IE [∏̀ k=1 〈µ̂Nt , fk〉 ] = 1 N ` IE [∏̀ k=1 ( N∑ m=1 fk(ξ m bNtc) )] (205)
= 1
N ` IE  ∑ π∈Π̃(`,N) (∏̀ k=1 fk(ξ π(k) bNtc) ) . (206) Therefore, subtracting equation 204 and equation 205 yields∣∣∣∣∣IE [∏̀ k=1 〈µ̂Nt , fk〉 ] − IE [∏̀ k=1 fk(ξ k bNtc) ]∣∣∣∣∣ ≤ b` ( 1− N ! `!N ` ) . (207)
Hence,
lim sup N→∞ ∣∣∣∣∣IE [∏̀ k=1 〈µ̂Nt , fk〉 ] − IE [∏̀ k=1 fk(ξ k bNtc) ]∣∣∣∣∣ = 0. (208) Combining equation 203-equation 208 yields the desired result.
xlviii

Review summaries:
This paper proposes to learn a kernel for training MMD-GAN by optimizing over the probability distribution that defines the kernel by means of random features. This is unlike the usual setting of MMD-GAN where the kernel is parametrized by composing a fixed top-kernel with a discriminator network that is optimized during training. The main motivation for this approach is to avoid having to 'manually' fix some parameters of the top-level kernel like the bandwidth of the RBF kernel. They provide an algorithm to achieve such optimization in probability space along with some consistency results and perform experiments on MNIST and Cifar10 to demonstrate empirically the advantage of such an approach over those that fix the top-level kernel.

This paper aims to improve the kernel selection issue of the MMD-based generative models. The author formulates the kernels via inverse Fourier transform and the goal is to learn the optimal N finite random Fourier features (RFF). The RFF samples are optimized by the proposed kernel alignment loss where the positive and negative labels are defined as samples coming from real and negative data distributions, respectively. Some theoretical analysis regarding the consistency of the learned kernel is provided. Experiment results on the IS score and FID on CIFAR-10 show improvement of the proposed methods over MMD-GAN baselines, while the results are not comparable to the original MMD-GAN due to unknown results.