{'fraction_non_alphanumeric': 0.0548598011544953, 'fraction_numerical': 0.01748762534216461, 'mean_word_length': 4.083489545782264, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 35, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Mean-Field Theory for Kernel Alignment with Random Features in Generative and Discriminative Models

We propose a novel supervised learning method to optimize the kernel in the maximum mean discrepancy generative adversarial networks (MMD GANs), and the kernel support vector machines (SVMs). Specifically, we characterize a distributionally robust optimization problem to compute a good distribution for the random feature model of Rahimi and Recht. 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. We evaluate our kernel learning method for the hypothesis testing problem by evaluating the kernel MMD statistics, and show that our learning method indeed attains better power of the test for larger threshold values compared to an untrained kernel. Moreover, our empirical evaluation on benchmark data-sets shows the advantage of our kernel learning approach compared to alternative kernel learning methods.


Introduction
A fundamental step in devising efficient learning algorithms is the design of good feature maps. Such feature maps embed data from an input space to some potentially highdimensional non-linear feature space, where it is often easier to discriminate data.
In the kernel methods of machine learning, such feature maps and thus the underlying kernel function, must be specified by the user. However, in many practical machine learning applications, it is unclear how to select an appropriate kernel function. In particular, choosing a good kernel often entails some domain knowledge about the underlying learning applications. In practice, it can be difficult for the practitioner to find prior justification for the use of one kernel instead of another to generate random feature maps. Moreover, a bad kernel model yields a poor classification accuracy. a good distribution for the random feature model of Rahimi and Recht [40,41]. The main contributions of this paper are summarized below: • Based on the notion of the kernel-target alignment, we characterize a distributional optimization problem to learn a good distribution for Rahimi and Recht's random feature model of a kernel function [40,41]. • To obtain a tractable optimization problem, we leverage the sample average approximation (SAA) method to transform the infinite dimensional distributional optimization problem to a finite dimensional optimization problem. • We establish the consistency of the finite sample average approximation. In particular, we show that as the number of samples in SAA tends to infinity, the optimal value of the finite sample estimates tend to their population values. • We propose a particle stochastic gradient descent method to solve the finite dimensional optimization problem associated with the Monte-Carlo sampling method. Using a mean-field analysis, we then show that the interactive particle system with SGD dynamics follows the gradient decent path on the Wasserstein manifold to minimize the distributional optimization problem. In this sense, we establish the consistency of the proposed particle SGD method for solving the distributional optimization problem. • We validate our kernel learning approach for generative modeling and discriminative analysis on benchmark data-sets.
1.1. Related works. The mean-field description of SGD dynamics has been studied in several prior works for different information processing tasks. Wang et al. [55] 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 [34,22,33] 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. They also establish the uniqueness and existence of the solution for the PDE using the connection between Mckean-Vlasov type PDEs and the gradient flows on the Wasserstein manifolds established by Otto [36], and Jordan, Kinderlehrer, and Otto [23]. Similar mean-field type results for two-layer neural networks are also studied recently in [43,49]. The present application of the mean-field theory to the kernel optimization is partly motivated by the work of Sinha and Duchi [48], which studied a distributional optimization for kernel learning with random features in the context of classification problems. Therein, the authors have proposed a robust optimization framework for the importance sampling of the random features. In contrast to the work of [48] that assign a weight (importance) to each sample and optimizes the weights, in this paper we directly optimize the samples.
For generative modeling using MMD GANs, kernel learning has been studied in several prior works [29,57,19]. Nevertheless, those methods are either based on heuristics, or are difficult to characterize theoretically. Our work is also related to the unpublished work of Wang, et al. [54], which proposes a solvable model of GAN and analyzes the scaling limits. However, our GAN model is significantly different from [54] and is based on the notion of the kernel MMD. Our work for kernel learning in generative modeling is also closely related to the recent work of Li, et al [29] which proposes an implicit kernel learning method based ii on the following kernel definition K h (ι(x), ι(y)) = IE ξ∼µ0 e (ih(ξ)(ι(x)−ι(y))) , where µ 0 is a user defined base distribution, and h ∈ H is a functions that transforms the base distribution µ 0 into a distribution µ that provides a better kernel. Therefore, the work of Li, et al [29] implicitly optimizes the distribution of random features through 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 (16) 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 [29]).
For classification problems using kernel SVMs, Cortes, et al. studied a kernel learning procedure from the class of mixture of base kernels. They have also studied the generalization bounds of the proposed methods. The same authors have also studied a two-stage kernel learning in [8] based on a notion of alignment. The first stage of this technique consists of learning a kernel that is a convex combination of a set of base kernels. The second stage consists of using the learned kernel with a standard kernel-based learning algorithm such as SVMs to select a prediction hypothesis. In [26], the authors have proposed a semi-definite programming for the kernel-target alignment problem. An alternative approach to kernel learning is to sample random features from an arbitrary distribution (without tuning its hyper-parameters) and then apply a supervised feature screening method to distill random features with high discriminative power from redundant random features. This approach has been adopted by Shahrampour et al. [46], where an energy-based exploration of random features is proposed. The feature selection algorithm in [46] employs a score function based on the kernel polarization techniques of [2] to explore the domain of random features and retains samples with the highest score.
The rest of this paper is organized as follows. In Section 2, we propose the kernel learning approach for the generative modeling. In Section 3, we study the kernel learning approach for the discriminative modeling. In Section 4, we present the main theoretical results, while deferring their proofs to Appendix. In Section 5, we evaluate our kernel learning approach for generative modeling using synthetic data-set and CIFAR-10 and MNIST benchmark data-sets. In Section 6, we provide empirical evaluations of our proposed kernel learning approach for classification via the kernel SVMs, where we evaluate our method on some benchmark data-sets. Lastly, in Section 8, we conclude this paper.

Kernel learning in generative models
Assume we are given data {v i } n i=1 that are sampled from an unknown distribution P V with the support V. In many unsupervised tasks, we wish to attain new samples from the distribution P V without directly estimating it. Generative Adversarial Network (GAN) [17] 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 ∼ P Z , Z ∈ Z from a userdefined distribution P Z (e.g. Gaussian distribution) into a new sample G(Z; W ) ∼ P W , such that the distributions P W and P V are close under some specified metric. 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 iii discriminator networks is then accomplished via solving a minimax optimization problem as below min W ∈W max δ∈∆ IE P V [D(X; δ)] + IE P Z [log(1 − D(G(Z; W ); δ))]. (2) In the high dimensional regimes, the generator trained via the minimax program of Eq. (2) can potentially collapse to a single mode of distribution where it always emits the same point [6]. To overcome this shortcoming, other adverserial generative models are developed 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 (MMD). Below, we formalize the notion of MMD: Definition 2.1. (Maximum Mean Discrepancy, [18]) 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 ) on X . The maximum mean discrepancy (MMD) between the distributions P and Q with respect to the function class F is defined below Adapting different function classes F in Eq. (3) of Definition (2.1) yield different adversarial models such as Wasserstein GANs (WGAN) [1], f -GANs [35], GMMN [28], and MMD GAN [30]. In the latter two cases, the function class F corresponds to a reproducing kernel Hilbert space (RKHS) of functions with a kernel K : X × X → IR, denoted by (H X , K). For RKHS function class, the squared MMD in Eq. (3) between the distributions P = P V and Q = P W has the following expression 1 where X = V ∪ W.
Instead of training the generator via solving the Min-Max optimization in Eq. (2), the MMD GAN model of [30] proposes to optimize the discrepancy between two distributions via optimization of an embedding function ι : IR d → IR d0 , d 0 d, i.e., where k : IR d0 × IR d0 → IR is a user-defined fixed kernel. In [30], the proposal for the kernel k : IR d0 × IR d0 → IR is a mixture of the Gaussians, where the bandwidth parameters σ 1 , · · · , σ m > 0. However, there are two challenges associated with the optimization of the embedding map with a per-defined Gaussian mixture kernel: 1 We note that the notation MMD K [P V , P W ] is conventionally used to denote the square root of the expression on the right hand side of Eq. (4). However, to simplify the arguments in subsequent results and their proofs, we adapt the notation of Eq. (4). iv (i ) The Gaussian mixture kernel is not expected to be optimal. Even if the Gaussian mixture kernel yields a good performance on a given data-set, there are no theoretical guarantees to suggest that it performs well across all data-sets that are encountered in practice. (ii ) When the Gaussian mixture kernel model is admissible, the choice of the Gaussian bandwidths σ 1 , · · · , σ m poses a statistical model selection problem. Although such model selection issues can be addressed using the cross-validation or jackknife, in practice such methods slow down the training process due to repeatedly re-fitting the model. To address the aforementioned issues, in this paper we propose to optimize the MMD loss with respect to the underlying kernel over a suitable kernel class K. Notice that the optimization of the embedding function in Eq. (5) is a special case of the general kernel learning problem formulated in Eq. (7) when the kernel class is selected as K def = {K : K(x, y) = k(ι(x), ι(y)), ι ∈ Q}.
2.1. Robust distributional optimization for kernel learning. To address the kernel model selection issue in MMD GAN [28], we consider a kernel optimization scheme with random features [40,41]. Let ϕ : IR d × IR D → [−1, 1] denotes the explicit feature maps and µ ∈ M(IR D ) denotes a probability measure from the space of probability measures M(IR D ) on IR D . The kernel function is characterized via the explicit feature maps using the following integral equation Then, the kernel optimization problem in can be formulated as a distribution optimization for random features, i.e, 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 where µ 0 is a user-defined base distribution, and To establish the proof, we consider the W p (·, ·) is the p-Wasserstein distance defined as below where the infimum is taken with respect to all couplings π of the measures µ, µ 0 ∈ M(IR D ), and Π(µ, µ 0 ) is the set of all such couplings with the marginals µ and µ 0 .
The kernel MMD loss function in Eq. (9) is defined with respect to the unknown distributions of the data-set P V and the model P W . Therefore, we construct an unbiased estimator for the MMD loss function in Eq. (9) based on the training samples. To describe the estimator, sample the labels from a uniform distribution y 1 , · · · , y n ∼ 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} : y i = +1}, and negative labels v J = {1, 2, · · · , n}/I, where their cardinality is |I| = |J | = n 2 . We consider the following assignment of labels: • Positive class labels: If y i = +1, sample the corresponding feature map from datadistribution x i = v i ∼ P V . • Negative class labels: If y i = −1, sample from the corresponding feature map from the generated distribution x i = G(Z i , W ) ∼ P W , Z i ∼ P Z . By this construction, the joint distribution of features and labels P Y,X has the marginals P X|Y =+1 = P V , and P X|Y =−1 = P W . Moreover, the following statistic, known as the kernel alignment in the literature (see, e.g., [48,9]), is an unbiased estimator of the MMD loss in Eq. (9), See Appendix A.1 for the related proof. The kernel alignment in (12) can also be viewed through the lens of the risk minimization Here, α > 0 is a scaling factor that determines the separation between feature vectors, and K * def = αyy T is the ideal kernel that provides the maximal separation between the feature vectors over the training data-set, i.e., K * (x i , x j ) = α when features have identical labels y i = y j , and K * (x i , x j ) = −α otherwise. Upon expansion of the risk function in (13), it can be easily shown that it reduces to the kernel alignment in (12) when α → +∞. Intuitively, the risk minimization in (13) gives a feature space in which pairwise distances are similar to those in the output space Y = {−1, +1}.
2.2. SAA for distributional optimization. The distributional optimization problem in Eq. (9) 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 where δ(·) is the Dirac's delta function concentrated at zero. In practice, the optimization problem in Eq. (13) is solved via the Monte-Carlo sample average approximation of the objective function, ≤ R , and µ N 0 is the empirical measure associated with the samples ξ 1 0 , · · · , ξ N 0 ∼ i.i.d. µ 0 . The empirical objective function vi in Eq. (15) can be optimized with respect to the samples ξ 1 , · · · , ξ N using the particle stochastic gradient descent. For the optimization problem in Eq. (15), the (projected) stochastic gradient descent (SGD) takes the following recursive form, 2 for k = 1, 2, · · · , N , where (y m , x m ), ( y m , x m ) ∼ i.i.d P x,y and η ∈ IR >0 denotes the learning rate of the algorithm, and the initial particles are ξ 1 0 , · · · , ξ N 0 ∼ i.i.d. µ 0 . At each iteration of the SGD dynamic in Eq. (16), a feasible solution for the inner optimization of the empirical risk function in Eq. (15) is generated via the empirical measure Indeed, we prove in Section 4 that for an appropriate choice of the learning rate η > 0, the empirical measure in Eq. (17) [51]. Our proposed MMD GAN model is distinguished from MMD GAN of [28] in that we learn a good kernel function in Eq. (26) 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 [28] for 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.

Kernel learning in discriminative models
In this section, we now consider the kernel learning problem in discriminative models. We focus on the classical setting of supervised learning, whereby we are given are given n feature vectors and their corresponding uni variate class labels ( In this paper, we are concerned with the binary classification and so the target spaces is given by Given a loss function : Y × IR → IR + , where (·, z) is convex (e.g. hinge function for SVMs, or quadratic for linear regression), and a reproducing Hilbert kernel space (H X , K), a classifier f ∈ H X ⊕ {1} is learned by minimizing the empirical risk minimization with a quadratic regularization,
In the particular case of the soft margin SVMs classifier (y, z) = max{0, 1 − yz}, the primal and dual optimizations take the following forms Primal : min ω∈IR n min where is the Hadamard (element-wise) product. The last term in the dual objective function in Eq. (22b), suggest that for a fixed dual vector β ∈ IR n , the optimal choice of the kernel can be obtained by maximizing the kernel target alignment, i.e., where in the above optimization we used the random feature model of Eq. (8) to transform the optimization of the kernel over the kernel class K to a distributional optimization over the distributional ball P. In Equation (23)  Inputs: The learning ratesη, η > 0 , the number of iterations of discriminator per generator update T ∈ N, the batch-size n, the number of random features N ∈ N. Regularization parameter α > 0. for t = 1, 2, · · · , T do Sample the labels y, y ∼ i.i.d Uniform{−1, 1}. Sample the features x|y = +1 ∼ P+, and x|y = −1 ∼ P− from the training data-set {(y1, x1), · · · , (yn, xn)}. Similarly, x| y = +1 ∼ P+, and x| y = −1 ∼ P−.
For all k = 1, 2, · · · , N , update the particles, end for Solve the following dual optimization problem max β∈IR n : β,y =0,0 β C1 Using a Monte-Carlo sample average approximation similar to Eq. (13), yields the following optimization problem where we recall from Eq. (13) that α > 0 is a regularization parameter. Due to similarity of the kernel-target alignment in the generative model of Eq. (13) and the underlying optimization of (24), we apply the SGD method of Eq. (16) to solve Eq. (24).
3.1. Proposed SVMs classifier with SGD kernel learning. In Algorithm 2, we describe the proposed SVMs classifer with the kernel learning approach. Instead of a joint optimization over the dual vector β ∈ IR n and the kernel K ∈ K, we consider a two phase optimization problem. In the first phase (the for loop), we optimize the kernel using the SGD optimization of Eq. (16). In the second phase, we construct the Gram matrix ΦΦ T [n] , and suppose ξ 1 , · · · , ξ N ∼ i.i.d. µ. Then, the following concentration results are established by some of the authors of this paper (cf. [25,Lemma 5.3]), where · F and · 2 are Frobenius and spectral norms, respectively, Diam(X ) is the diameter of the feature space, d is the dimension of features x 1 , · · · , x n ∈ X ⊂ IR d , and σ p def = IE µ [ ξ 2 2 ] is the second moment of the random vector ξ ∼ µ. From Eq. (25a), we observe that as the number of particles N tends to infinity, the Gram matrix ΦΦ T concentrates in various norms around the (scaled) kernel matrix K. Due to the concentration phenomenon of Eq. (25a), the substitution of the kernel matrix K in the dual optimization of Eq. (22b) with the scaled Gram matrix ΦΦ T /N yields a good approximation. Equation (27) of Algoirthm 2 is the result of such an approximation.

Main results: consistency and mean-field analysis
In this section, we provide theoretical guarantees for the consistency of various approximations we made to optimize the population MMD loss function in (9). We defer the proofs of the following theoretical results to Section A of Appendix. Due to similarity of the kernel learning approach for generative and classification tasks, we state the underlying assumptions and theoretical results only for the generative model. Nevertheless, it is easy to see that our theoretical results generalizes to the classification problem.

4.1.
Assumptions. Before we delve into technical results, we state the main assumptions underlying them: where V = support(P V ) and W = support(P W ) are the supports of the distributions P V and P W respectively. (A.2) The feature maps are bounded and Lipchitz almost everywhere (a.e.) ξ ∈ IR D .

4.2.
Consistency of finite-sample estimator. n this part, we prove that the solution to finite sample optimization problem in (15) approaches its population optimum in (9) as the number of data points as well as the number of random feature samples tends to infinity.
In this part, we prove that the solution to finite sample optimization problem in (15) approaches its population optimum in (9) as the number of data points as well as the number of random feature samples tends to infinity. Consider the distribution balls P and P N 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 respectively. Then, with the probability of (at least) 1 − 3 over the training data samples {(x i , y i )} n i=1 and the random feature samples {ξ k 0 } N k=1 , the following non-asymptotic bound holds where c 1 = 3 1 4 × 2 4 , and c 2 = 9 × 2 11 .
x The proof of Theorem 4.1 is presented in Appendix A.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 (29) thus shows that when n, N, α → +∞, the solution obtained from solving the empirical risk minimization in (13) yields a MMD population value tending to the optimal value of the distributional optimization in (9). 4.3. Consistency of particle SGD algorithm. 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 (15). In practice, the particle SGD is executed for a few iterations and its values are returned as an estimate for (ξ 1 * , · · · , ξ N * ). Consequently, it is desirable to establish a consistency type result for the particle SGD estimates (ξ 1 m , · · · , ξ N m ) at the m-th iteration, where the notion of consistency will be made precise shortly. To prove such a consistency result, we define the scaled empirical measure as follows At any time t, the scaled empirical measure µ N t 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 p N t (ξ) def = µ N t (dξ)/dξ in the following theorem: in Section A are satisfied. Further, suppose that the Radon-Nikodyme derivative q 0 (ξ) = µ 0 (dξ)/dξ exists. Then, there exists a unique solution (p * t (ξ)) 0≤t≤T to the following nonlinear partial differential equation Moreover, the measure-valued process {(µ N t ) 0≤t≤T } N ∈IN defined in Eq. (30) converges (weakly) to the unique solution µ * t (ξ) = p * t (ξ)dξ as the number of particles tends to infinity N → ∞. 3 .
The proof of Theorem 4.2 is presented in Appendix A.1. Due to the mean-field analysis of Theorem 4.2, we can prove that the empirical measure µ N m of the particles in SGD dynamic (16) remains inside the feasible distribution ball P N : for the SGD in (16).
Let us make two remarks about the PDE in Eq. (31). First, the seminal works of [36], and [23] establishes an intriguing connection between the McKean-Vlasov type PDEs specified in (31) and the gradient flow on the Wasserstein manifolds. More specifically, the PDE equation in Eq. (31) can be thought of as the minimization of the energy functional using the following gradient flow dynamics with respect to the metric of the Wasserstein manifold . This shows that when the number of particles in particle SGD (16) 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 Eq. (13). In this sense, the particle SGD is a 'consistent' approximation algorithm for solving the distributional optimization.
Second, notice that as the scaling parameter tends to infinity α → ∞, the energy func- Interestingly, this limiting energy functional is precisely the kernel polarization of Baram [2], measuring the correlation between the class labels and random features 4.4. Kernel Learning by Solving A PDE. The kernel selection methods in the literature focuses on optimization methods based on kernel alignment optimization problem; see, e.g., [9,26,10].
The PDE in Eq. (31) puts forth an alternative method to finding good kernel functions. Namely, a good kernel function can be computed in two stages: first compute the stationary solution p * (ξ) corresponding to a solution of the PDE in Eq. (31) with ∂p * t (ξ)/∂t = 0, where the expectation with respect to the unknown data distribution is replaced by its finite-sample average. Then approximating the kernel via a Monte-Carlo sampling ξ 1 , · · · , ξ N ∼ i.i.d. p * (ξ), While such a method can be successful when the random features ξ ∈ IR D are low dimensional, in the high-dimensional settings (D 1), solving the PDE in Eq. (31) numerically and sampling from a high dimensional density function p * (ξ) appears to be computationally more expensive than solving the kernel-target alignment via SGD iterative optimization.

Propagation of Chaos.
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 } N k=1 tends to infinity (N → +∞), their dynamics are decoupled. Definition 4.3. (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 S n if it is invariant under the permutation π def = (π(1), · · · , π(N )) of indices. In particular, for all Borel subsets B ∈ B(S n ).
xii An interpretation of the exchangablity condition (35) 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: According to Eq. (36) of Definition 4.4, 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 ν(ds 1 )ν(ds 2 ) · · · ν(ds k ) = ν ⊗k on S k . 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 S N is chaotic if it is ν-chaotic for some probability measure ν on S.
If a Markov process on S N 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 K N for the Markov process. The distribution at time t > 0 is the probability measure U N t ν ⊗N is defined by the kernel  (38) whose N -th term is a Markov transition function on S N that satisfies the permutation condition 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): µ 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 µ N t ∈ M(IR D ) defined via scaling

Empirical Results for Generative Tasks
We now turn to empirical evaluations. We test the performance of Algorithm 1 on synthetic data-set, as well as on benchmark data-sets.

Experiments on Synthetic Data-Set.
The synthetic data-set we consider is as follows: To reduce the dimensionality of data, we consider the embedding ι : In this case, the distribution of the embedded features are P X|Y =+1 = N(0, (1 + λ)ΣΣ T ), and P X|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], In Figure 1, we show the distributions of i.i.d. samples from the distributions P V and P W 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 1, we clearly observe that for large values of λ, the data-points from the two distributions P V and P W have a large overlap and conducting a statistical test to distinguish between these two distributions is more challenging. 5.1.1. Kernel Learning Approach. Figure 3 depicts our two-phase kernel learning procedure which we also employed in our implementations of Algorithm 1 on benchmark data-sets of Section 5.3 in the main text. The kernel learning approach consists of training the autoencoder and the kernel optimization sequentially, i.e., where the function class is defined Q ). Here, σ(·) is the sigmoid non-linearity. Now, we consider a two-phase optimization procedure: • Phase (I): we fix the kernel function, and optimize the auto-encoder to compute a co-variance matrix Σ and the bias term b for the dimensionality reduction. • Phase (II): we optimize the kernel based from the learned embedded features ι(x). This two-phase procedure significantly improves the computational complexity of SGD as it reduces the dimensionality of random feature samples ξ ∈ IR D , D = d 0 d. When the kernel function K is fixed, optimizing the auto-encoder is equivalent to the kernel learning step of [28].

Statistical Hypothesis Testing with the Kernel
1} is used to distinguish between these hypotheses: • Null hypothesis H 0 :  To perform hypothesis testing via the kernel MMD, we require that H X is a universal RKHS, defined on a compact metric space X . Universality requires that the kernel K(·, ·) be continuous and, H X be dense in C(X ). Under these conditions, the following theorem establishes that the kernel MMD is indeed a metric: Theorem 5.1. (Metrizablity of the RKHS) Let F denotes a unit ball in a universal RKHS H X defined on a compact metric space X with the associated continuous kernel K(·, ·). Then, the kernel MMD is a metric in the sense that MMD K [P V , P W ] = 0 if and only if To design a test, let µ N m (ξ) = 1 N N k=1 δ(ξ − ξ k m ) denotes the solution of SGD in (16) for solving the optimization problem. Consider the following MMD estimator consisting of two xvi U -statistics and an empirical function we design a test statistic as below where τ ∈ IR is a threshold. Notice that the unbiased MMD estimator of (43) can be negative despite the fact that the population MMD is non-negative. Consequently, negative values for the statistical threshold τ (44) are admissible. In the following simulations, we only consider non-negative values for the threshold τ .
A Type I error is made when H 0 is rejected based on the observed samples, despite the null hypothesis having generated the data. Conversely, a Type II error occurs when H 0 is accepted despite the alternative hypothesis H 1 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 H 0 when it is indeed incorrect. In particular, In this sense, the statistical power controls the probability of making Type II errors.

5.2.
Empirical results on benchmark data-sets. In Figure 2, we show the evolution of the empirical measure µ N m (ξ) of SGD particles by plotting the 2D histogram of the particles ξ 1 m , · · · , ξ N m ∈ IR d at different iterations of SGD. Clearly, starting with a uniform distribution in 2(a), the empirical measure seemingly evolves into a Gaussian measure in Figure  2(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 3, we evaluate the power of the test for 100 trials of hypothesis test using the test statistics of (44). 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 3, 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 autoencoder only with a fixed standard Gaussian kernel K(x, y) = exp(− 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 3, we also observe that interestingly, the phase transition in the statistical threshold τ is less sensitive to the parameter λ. This phenomenon is also observed in Figure  [28] with an autoencoder for dimensionality reduction in conjunction with a mixed RBF Gaussian kernel whose bandwidths are optimized via cross-validation, Panels (c)-(g): MMD GAN in [28] with a single RBF Gaussian kernel with an auto-encoder for dimensionality reduction in conjunction with a single RBF Gaussian kernel whose bandwidth is optimized via cross-validation, and Panel (d)-(g): GMMN without an auto-encoder [30]. 4, where the improvements due to SGD for more difficult two sample tests (i.e., smaller values of λ) is larger than the case of distinguishable distributions.

5.3.
Performance on benchmark datasets. We evaluate our kernel learning approach on large-scale benchmark data-sets. We train our MMD GAN model on two distinct types of data-sets, namely on MNIST [27] and CIFAR-10 [27], where the size of training instances are 60 × 10 3 and 50 × 10 3 , respectively. All the generated samples are from a fixed noise random vectors and are not singled out.
5.3.1. Implementation and hyper-parameters. We implement Algorithm 1 as well as MMD GAN [28] in Pytorch using NVIDIA Titan V100 32GB graphics processing units (GPUs). The source code of Algorithm 1 is built upon the code of [28], 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-set. For a fair comparison, our hyper-parameters are adjusted as in [28], i.e., the learning rate of 0.00005 is considered for RMSProp [51]. Moreover, the batch-size for training the generator and auto-encoder is n = 64. The learning rate of particle SGD is tuned to η = 10. To approximate the kernel, we use the the random feature model of Rahimi and Recht [40,41], where ϕ(x; ξ) = 2/N cos(x T ξ + b). Here b ∼ Uniform[−π, π] is a random bias term. 5.3.3. Practical considerations. When data-samples {V i } ∈ IR d are high dimensional (e.g. CIFAR-10), the particles ξ 1 , · · · , ξ N ∈ IR p , p = d in SGD (16) are also high-dimensional. To reduce the dimensionality of the particles, we apply an auto-encoder architecture similar to [28], 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 p = 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.

5.3.4.
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, a large value of α slows down the convergence of particle SGD. This is due to the fact that the coupling strength between the particles in Eq. (16) decrease as α increase. The scaling factor is set to be α = 1 in all the following experiments.

5.3.5.
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 8, we show the generated samples on CIFAR-10 and MNIST data-sets, using our Algorithm 1, MMD GAN [28] with a mixed and homogeneous Gaussian RBF kernels, and GMMN [30]. Figure 8(a) shows the samples from Algorithm 1, Figure 8(b) shows the samples from MMD GAN [28] with the RBF Gaussian mixture kernel of Eq. , 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 [28] which requires manual tuning of the hyper-parameters. In Figure 8(c), we show the MMD GAN result with a single kernel RBF Gaussian kernel whose bandwidth are optimized via cross-validation and is σ = 16. Lastly, in Figure 8(d), we show the samples from GMMN [30] 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 8(e)-(h), the difference between our method and MMD GAN [28] is visually more pronounced. We observe that without a manual tuning of the kernel bandwidth and by using the particle SGD (16) to optimize the kernel, we attain better generated images in Figure 8(e), compared to MMD GAN with mixed RBF Gaussian kernel and manual bandwidth tuning in Figure 8(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 8(h). 5.3.6. Quantitivative comparison. To quantitatively measure the quality and diversity of generated samples, we compute the inception score (IS) [45] as well as Frèchet Inception Distance (FID) [20] on CIFAR-10 images.
Intuitively, the inception score is used for GANs to measure samples quality and diversity. This score is based on the Inception-v3 Network [50] which is a deep convolutional architecture designed for classification tasks on ImageNet [11], a dataset consisting of 1.2 million RGB images from 1000 classes. Given an image x, the task of the network is to output a class label y in the form of a vector of probabilities. The inception score uses an Inception-v3 MMD GAN (Gaussian) [28] 67.244 ± 0.134 5.608±0.051 MMD GAN (Mixture Gaussian) [28]  Network pre-trained on ImageNet and calculates a statistic of the networks outputs when applied to generated images. More precisely, the inception score of a generative model is where D KL (·||·) is the Kullback-Leibler divergence. The definition of the inception score is motivated by the following two observations: (i ) The images generated should contain meaningful objects, i.e., P Y |X should be low entropy. In other words, images with meaningful objects are supposed to have low label (output) entropy, that is, they belong to few object classes (ii ) The generative algorithm should generate diverse images from all the different classes in ImageNet, i.e., the distribution of labels P Y should have a high entropy.
The FID improves on IS by actually comparing the statistics of generated samples to real samples, instead of evaluating generated samples independently. In particular, Heusel, et al. [20] propose to use the Fréchet distance between two multivariate Gaussians N(µ 1 , Σ 1 ) and N(µ 2 , Σ 2 ) as follows 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}).

Empirical Results for Classification Tasks
We now turn to empirical evaluations for classification on synthetic and benchmark datasets. We compare our method with two alternative kernel learning techniques, namely, the importance sampling of Sinha and Duchi [48] , and the Gaussian bandwidth optimization via k nearest neighbor (kNN) [7]. In the sequel, we provide a brief description of each method: • Importance Sampling: the importance sampling of Sinha and Duchi [48] which proposes to assign a weight w m to each sample ξ m for m = 1, 2, · · · , N . The weights are then optimized via the following standard optimization procedure  where Q N def = {w ∈ IR N + : w, 1 = 1, χ(w||1/N ) ≤ R} is the distributional ball, and χ(·||·) is the χ 2 −divergence. Notice that the complexity of solving the optimization problem in Eq. (48) is insensitive to the dimension of random features x i ∈ IR d .
• Gaussian Kernel with kNN Bandwidth Selection Rule: In this method, we fix the Gaussian kernel K(x, y) = exp( x − y] 2 2 /σ 2 ). To select a good bandwidth for the kernel, we use the deterministic k nearest neighbor (kNN) approach of [7]. In particular, we choose the bandwidth according to the following rule where x kNN i is defined as k nearest neighbor of x i . In our experiments, we let k = 3.
We then generate random Fourier features ϕ(x; 6.1. Empirical Results on the Synthetic Data-Set. For experiments with the synthetic data, we follow the model of Sinha and Duchi [48] and compare our results with the importance sampling method proposed therein and the publicly available codes on Github. 4 In particular, we generate the features from the standard Gaussian distribution x ∼ N(0, I d×d ) with the class labels y = sign( The underpinning idea behind this data generation model is that the Gaussian kernel is ill-suited for this classification task, as the Euclidean distance used in this kernel does not capture the underlying structure of the classes. Nevertheless, we compare our results with the Gaussian kernel with kNN bandwidth selection procedure. 6.2. Quantitiative Comparison. In Table 2 and Table 3, we present the training and test errors using kernel SVMs in conjunction with the trained kernel of Algorithm 2, the trained kernel via auto-encoder embedding only, the trained kernel via the importance sampling optimization of Eq. (48), and the standard Gaussian kernel. 567 In Table 2, the number of random feature samples is fixed N = 1000, while the data dimension d changes from d = 5 to d = 20. We observe that already at dimension d = 10, the SGD optimization method outperforms the importance sampling method, albeit with a larger computational cost. In Table 3, we the test and training errors of the kernel learning methods for different number of random feature samples N . Clearly, already at N = 1000 random feature samples, the training and test errors for SGD is comparable to that of importance sampling using N = 10000 samples. We observe the improvement in the accuracy of both SGD and importance sampling as the number of random features increase.

Implementation on Benchmark Data-Sets
We now apply our kernel learning method for classification and regression tasks on realworld data-sets.    Table 4. Description of the benchmark data-sets used in this paper.

SGD
7.1. Data-Sets Description. We apply our kernel learning approach to classification and regression tasks of real-world data-sets. In Table 4, we provide the characteristics of each data-set. All of these datasets are publicly available at UCI repository. The 4097 data points are then divided and shuffled every into 23 segments, each segment contains 178 data points for 1 second, and each data point is the value of the EEG recording at a different point in time. Figure 9, we present the training and test results for regression and classification tasks on benchmark data-sets, using top d = 35 features from each data-set and for different number of random feature samples N . In all the experiments, the SGD method provides a better accuracy in both the training and test phases. Nevertheless, in the case of seizure detection, we observe that for a small number of random feature samples, the importance sampling and Gaussian kernel with kNN for bandwidth outperforms SGD.

Quantitative Comparison. In
In Figure 10, we illustrate the time consumed for training the kernel using SGD and importance sampling methods versus the number of random features N . For the linear regression on Buzz, and Online news popularity the difference between the run-times are negligible. However, for classification tasks using Adult and Seizure, the difference in run-times are more pronounced.
Those visual observations are also repeated in a tabular form in Table 5, where the training and test errors as well algorithmic efficiencies (run-times) are presented for N = 100, N = 1000, N = 2000, and N = 3000 number of features.

Conclusions
We have presented a novel supervised method to optimize the kernel in generative and discriminative machine learning models. Our kernel learning approach is based on a distributionally robust optimization for the distribution of random features over the Wasserstein distributional ball. To obtain a tractable finite dimensional optimization problem, we used a sample average approximation (SAA) for random features, and subsequently applied a particle stochastic gradient descent (SGD). We provided theoretical guarantees for the consistency the finite sample (Monte-Carlo) approximations. In addition, using methods of mean-field theory, we proved that SGD is a consistent algorithm in the sense that the empirical distribution of the particles is an admissible solution for the distributional optimization problem. Our empirical evaluations on synthetic and benchmark data-sets shows the efficacy of our kernel learning method for generative and discriminative machine learning tasks.   Table 5. Comparison between the SGD kernel learning method (Algorithm 2), the importance Sampling of [48], and a regular Gaussian kernel in conjunction with the kernel SVMs for different random feature samples N on benchmark data-sets. xxvii 9. Acknowledgement M.B.K. would like to thank his colleague Dr. Chuang Wang from Harvard University for his intriguing theoretical work on the scaling limits of the online algorithms [55] which motivated this work. M.B.K. would also like to thank Prof. Ravi Mazumdar from The University of Waterloo for his theoretical paper [53] which inspired the proof of tightness of Skhorhod spaces in this paper. We also thank Dr. Maxime Bassenne for useful discussion regarding the MacKean-Vlasov PDE.
The research of MBK, LS, and LX is supported by NIH under the grant 1R01 CA176553 and R01E0116777. The content of this article are solely the responsibility of the authors and do not necessarily reflect the official NIH views.

Appendix Appendix A. Proofs of Main Theoretical Results
Notation: We denote vectors by lower case bold letters, e.g. x = (x 1 , · · · , x n ) ∈ IR n , and matrices by the upper case bold letters, e.g., M = [M ij ] ∈ IR n×m . The Frobenius norm of a matrix is denoted by denote the Euclidean ball of radius r centered at x. For a given metric space X , Let C b (IR d ) denote the space of bounded and continuous functions on X equipped with the usual supremum norm Further, C k b (X ) the space of all functions in C b (X ) whose partial derivatives up to order k are bounded and continuous, and C k c (X ) the space of functions whose partial derivatives up to order k are continuous with compact support.
We denote the class of the integrable functions f with f (t) ≥ 0 a.e., on 0 ≤ t ≤ T by L 1 + [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 by M(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 ∈ C b (X ), we define 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 ∈ C b (X ). We denote the weak convergence by µ N t weakly → µ. Notice that when X is Polish, then M(X ) equipped with the weak topology is also Polish. 9 For a Polish space X , let D X ([0, T ]) denotes the Skorokhod space of the cádlág functions that take values in X defined on [0, T ]. We assume that D X ([0, T ]) is equipped with the Skorokhod's J 1 -topology [4], which in that case D X ([0, T ]) is also a Polish space.
To establish the concentration results in this paper, we require the following definitions: Definition A.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 In the sequel, we consider the Orlicz modulus ψ ν (x) def = exp(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: The sub-Gaussian norm of a random variable Z, denoted by Z ψ2 , is defined as For a random vector Z ∈ IR n , its sub-Gaussian norm is defined as follows Definition A.3. (Sub-exponential Norm) The sub-exponential norm of a random variable Z, denoted by Z ψ1 , is defined as follows For a random vector Z ∈ IR n , its sub-exponential norm is defined below We use asymptotic notations throughout the paper. We use the standard asymptotic notation for sequences. If a n and b n are positive sequences, then a n = O(b n ) means that lim sup n→∞ a n /b n < ∞, whereas a n = Ω(b n ) means that lim inf n→∞ a n /b n > 0. Furthermore, a n = O(b n ) implies a n = O(b n poly log(b n )). Moreover a n = o(b n ) means that lim n→∞ a n /b n = 0 and a n = ω(b n ) means that lim n→∞ a n /b n = ∞. Lastly, we have a n = Θ(b n ) if a n = O(b n ) and a n = Ω(b n ).
A.1. Proof of Theorem 4.1. By the triangle inequality, we have that where the terms A i , i = 1, 2, 3, 4 are defined as follows In the sequel, we compute an upper bound for each term on the right hand side of Eq. (57): Upper bound on A 1 : xxix First, notice that the squared kernel MMD loss in Eq. (4) can be characterized in terms of class labels and features defined in Section 2.1 as follows To see this equivalence, we first rewrite the right hand side of Eq. (58) as follows Now, recall from Section 2.1 that P x|y=+1 = P V , and P x|y=−1 = P W 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} = 1 2 . Therefore, from Eq. (59), we derive For any given W ∈ W, we have that where the error term is defined using the random features Now, we invoke the following strong duality theorem [14]: provided that Ψ is upper semi-continuous in ξ.
xxx Under the strong duality of Theorem A.4, we obtain that In the sequel, let p = 2. The Moreau's envelope [37] of a function f : X → IR is defined as follows where β > 0 is the regularization parameter. When the function f is differentiable, the following lemma can be established: In particular, when f is L f -Lipschitz, we have The proof is presented in Appendix B.1. Now, we return to Equation (62). We leverage the lower bound on Moreau's envelope in Eq. (64) of Lemma A.5 as follows Let ζ * = ζ * (ξ, s) = arg sup ζ∈IR D ∇E n (1 − s)ξ + sζ 2 . Then, applying the union bound in conjunction with Inequality (66) yields Now, we state the following lemma: xxxi Lemma A.6. (Tail Bounds for the Finite Sample Estimation Error) Consider the estimation error E n defined in Eq. (60). Then, the following statements hold: • Z = ∇E n (ξ) 2 2 is a sub-exponential random variable with the Orlicz norm of Z ψ1 ≤ 9×2 9 ×L 4 n 2 for every ξ ∈ IR D . Moreover, • E n (ξ) is zero-mean sub-Gaussian random variable with the Orlicz norm of E n (ξ) ψ2 ≤ 16 √ 3L 4 n for every ξ ∈ IR D . Moreover, The proof of Lemma A.6 is presented in Appendix A.2. Now, we leverage the probability bounds (68) and (69) of Lemma A.6 to upper bound the terms on the right hand side of Eq. (67) as below 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 with the probability of (at least) 1 − .

Upper bound on A 2 :
xxxii To establish the upper bound on A 2 , we recall that Here, the last inequality is due to Theorem A.4 and the following duality results hold Now, in the sequel, we establish a uniform concentration result for the following function Then, from Eq. (75) we have We now closely follow the argument of [40] to establish a uniform concentration result with respect to the data points x, x ∈ X . In particular, consider an -net cover of X ⊂ IR d .
We compute an upper bound on the second moment of the random variable L T as follows We further proceed using the triangle inequality as well as the basic inequality (a 1 + a 2 + · · · + a N ) 2 ≤ N (a 2 1 + a 2 2 + · · · + a 2 N ), To proceed from (81), we leverage the following lemma: Lemma A.8. (Moreau's Envelop of Parametric Functions) Consider the parametric function f : X × Θ → IR and the associated Moreau's envelope for a given θ ∈ Θ ⊂ IR d : Furthermore, define the proximal operator as follows Then, Moreau's envelope has the following upper bound The proof is presented in Appendix B.2. Equipped with Inequality (84) of Lemma A.8, we now compute an upper bound on the right hand side of Eq. (81) as follows where the last inequality is due to (A.2).
Invoking Markov's inequality now yields xxxv Now, using the union bound, for every arbitrary pair of data points (x, x) ∈ X × X the following inequality holds Following the proposal of [40], we choose = ( .
Thus, with the probability of at least 1 − , the following inequality holds where W(·) is the Lambert W -function. 10 Since W(x) ≤ ln(x) for x > e, we can rewrite the upper bound in terms of elementary functions provided that N is sufficiently large and/or is sufficiently small so that 2 8 N diam 2 (X ) ≥ e. Plugging Inequality (89) in Eq. (76) now results in the following inequality for all W ∈ W. Employing(72) from Lemma A.7 now yields the following upper bound Upper bound on A 3 : Recall that the solution of the empirical risk function of Eq. (15) is denoted by We also define the solution of the empirical kernel alignment as follows Due to the optimality of the empirical measure µ N * for the inner optimization in Eq. (92), the following inequality holds Upon expansion of MMD α µ N P V , P W N * , and after rearranging the terms in Eq. (94), we arrive at 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 Eq. (93), we have Putting together Inequalities (95) and (96) yields Similarly, due to the optimality of the empirical measure µ N for the optimization in Eq. (93) we have that Combining Eqs. (97) and (98) then yields Upper bound on A 4 : xxxvii The upper bound on A 4 can be obtained exactly the same way as A 1 . Ideed, from Eq. (71) it follows directly that Now, plugging the derived upper bounds in A 1 -A 4 in Eq. (57) and employing the union bound completes the proof.
A.2. Proof of Theorem 4.2. The proof has three main ingredients and follows the standard procedure in the literature, see, e.g., [55,31]. In the first step, we identify the mean-field limit of the particle SGD in Eq. (16). In the second step, we prove the convergence of the measured-valued process {(µ N t ) 0≤t≤T } to the mean-field solution by establishing the precompactness 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 sub-sequences via the action of the empirical measure µ N m (ξ) = 1 . In particular, we use the standard techniques of computing the scaling limits from [31].
Recall that the action of an empirical measure on a bounded function is defined as follows We analyze the evolution of the empirical measure µ N m via its action on a test function f ∈ C 3 b (IR p ). Using Taylor's expansion, we obtain where R N m is a remainder term defined as follows where ξ k def = ( ξ k (1), · · · , ξ k (p)), and ξ k (i) ∈ [ξ k m (i), ξ k m+1 (i)], for i = 1, 2, · · · , p.
Plugging the difference term (ξ k m+1 − ξ k m ) from the SGD equation in Eq. (16) results in Now, we define the drift and Martingale terms as follows respectively. Using the definitions of D N m and M N m in Eqs. (105a)-(105b), we recast Equation (104) as follows Summation over = 0, 1, 2 · · · , m − 1 and using the telescopic sum yields We also define the following continuous embedding of the drift, martingale, and the remainder terms as follows The scaled empirical measure µ N t def = µ N N t then can be written as follows Since the drift process (D N t ) 0≤t≤T is a piecewise cádlág process, we have xxxix where the functional R[µ s ] is defined as follows d x), (dy, d y)). Therefore, the expression in Eq. (109) can be rewritten as follows In the following lemma, we prove that the remainder term sup 0≤t≤T |R N t | vanishes in probabilistic sense as the number of particles tends to infinity N → ∞: In particular, with the probability of at least 1 − ρ, we have and thus lim sup N →∞ sup 0≤t≤T |M N t | = 0 almost surely. Proof. The proof is deferred to Appendix B.6. Now, using the results of Lemmata A.10 and A.9 in conjunction with Eq. (112) yields the following mean-field equation as N → ∞, × ∇f (ξ)(ϕ( x; ξ)∇ ξ ϕ(x; ξ) + ϕ(x; ξ)∇ ξ ϕ( x; ξ)), µ s P ⊗2 x,y ((dx, d x), (dy, d y)) ds.
Notice that he mean-field equation in Eq. (116) is in the weak form. When the Lebesgue density p t (ξ) = dµ t /dξ exists, the McKean-Vlasov PDE in Eq. (31) can be readily obtained from Eq. (116).

xl
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: Definition A.11. (Tightness) A set A of probability measures on a metric space S is tight if there exists a compact subset S 0 ⊂ S such that for all ε > 0. A sequence {X N } N ∈IN of random elements of the metric space S is tight if there exists a compact subset S 0 ⊂ S such that This condition is referred to as the compact-containment condition. To establish (J1), we closely follow the proof of [15,Lemma 6.1.]. In particular, for each L > 0, we define S L = [0, B] p . Then, S B ⊂ IR p is compact, and for each t ≥ 0, and N ∈ IN, we have where (a) follows from Markov's inequality, and (b) follows from the upper bound on the norm of the particles in Eq. (203) of Appendix B. We now define the following set xli We let U T,γ = U B , where U B is the completion of the set U B . By definition, U T,γ is a compact subset of M(IR D ). Now, we have this implies that for any γ > 0, there exists a B > 0, such that This completes the proof of (J.1). To verify (J.2), we consider the following class of functions This completes the tightness proof of the of the laws of the measured-valued process {(µ N H is defined in Eq. (127). Now, condition (T.1) is readily verified since where in the last step, we used the fact that f ∈ C 3 b (IR p ), and hence, f ∞ ≤ b. Thus, IP(H(µ N 0 ) ≥ b) = 0 for all N ∈ IN, and the condition (T.1) is satisfied. Now, consider the condition (T.2). From Equation (112), and with 0 ≤ t 1 < t 2 ≤ T we have To bound the first term, recall the definition of R[µ s ] from Eq. (111). The following chain of inequalities holds, Let I : IR p → IR, I(ξ) = 1 denotes the identity function. Notice that I, µ s = IR p µ s (ds) = 1. From (134), we proceed as follows where the last inequality is due to (A.1). Therefore, Consider the middle term of (133). Using the definition of the martingale term in Eq. (108b), we obtain that In Equation of Section B, we have proved the following concentration bound xliii Now, recall the alternative definition of the sub-Gaussian random variables: Definition A.14. (Sub-Gaussian Random Variables [5]) A random variable X is We enumerate a few standard consequences of sub-Gaussianity [5]. If X i are independent and σ 2 i -sub-Gaussian, then Now, it is clear from (138) andthat M N m is sub-Gaussian random variable with a zero mean, and with the parameter σ 2 m = where (t * 1 , t * 2 ) def = arg sup t1,t2≤T,|t1−t2|≤ρ N t2 = N t1 M N . We first compute a bound for the last term of (133) using the definition of the scaled term R N t from (108c). We have Putting together (136), (141), and (145), we conclude from Inequality (133) that 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 {µ N t } N ∈IN . Now, we prove its convergence to a mean-field solution (µ * t ) 0≤t≤T .
Theorem A.15. (Prokhorov's theorem [39]) A subset of probability measures on a complete separable metric space is tight if and only if it is pre-compact.
According to Theorem A.15, the tightness of the Skorkhod Space D M(IR D ) ([0, T ]) implies its pre-compactness which in turn implies the existence of a converging sub-sequence } N k is a stochastic process defined on the Skorkhod space. Therefore, let π N k denotes the law of the converging sub-sequence {(µ N t ) 0≤t≤T } N k . By definition, π N k is an element of the measure space M(D [0,T ] (M(IR D ))). In the sequel, we closely follow the argument of [55,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 We compute the expectation of the functional F t with respect to π N k . We then have Now, from Equation (112), we have that Plugging (148) in Eq. (147) gives xlv where the last equality is due to the bounds in Eqs. (113) and (114) of Lemmata A.9 and A.10, respectively. Taking the limit of N → ∞ from Eq. (149) yields It can be shown that the functional F t [·] is continuous and bounded. Therefore, due the weak convergence of the sequence {π N k } N k ∈IN to π ∞ , Eq. (150) implies that Since the identity (151) holds for all bounded test functions f ∈ C 3 b (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 A.12, the support of the measured-valued process (µ N t ) 0≤t≤T = ( µ N N t ) 0≤t≤T is compact for all 0 ≤ t ≤ T . Moreover, in Step 2 of the proof, we established that the measure valued process (µ N t ) 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 meanfield equation also have compact supports. Let Ξ ⊂ IR D 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 ∈ C 3 b ( Ξ). Second, for all bounded continuous test functions f ∈ C 3 b ( Ξ), the operator f → µ t , f is a linear operator with µ t (IR D ) = 1. Hence, from Riesz-Markov-Kakutani representation theorem [44,52] by assuming µ t ∈ M(IR D ) , existence of unique operator implies f → f, µ t implies the existence of the unique probability measure µ t . Now, we equip the measure space M(IR D ) with the following norm 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 (µ * ,1 t ) 0≤t≤T , (µ * ,2 t ) 0≤t≤T are two solutions satisfying the mean-field equations (116) with the initial distributions µ 1 0 , µ 2 0 ∈ M(IR D ), respectively. For any test function f ∈ C 3 b ( Ξ) we have that × ∇f (ξ)(∇ ξ (ϕ( x; ξ)ϕ(x; ξ))) T , µ * ,1 s − µ * ,2 s P ⊗2 x,y ((dz, d z) ds.
We bound the first term on the right side of Equation (153) as follows xlvi where used the definition of the norm · on the measure space M(IR D ) from Eq. (152).
A.3. Proof of Corollary 4.2.1. To establish the proof, we recall from Eq. (107) that for all f ∈ C b (IR p ).The total variation (TV) distance between two measures ν 1 , ν 2 ∈ M(IR D ) is defined Therefore, by taking the supremum of Eq. (165) over all the bounded functions f ∈ C b (IR p ), b = 1, we obtain that Based on the upper bound (204) on the remainder term, we have for some constant C 0 > 0. Moreover, from the concentration inequality (215), we also have that with the probability of at least 1 − δ, the following inequality holds Lastly, recall the definition of the drift term in Eq. (105a). By carrying out a similar bounding method leading to Eq. (135), it can be shown that By plugging Eq. (168), (169), and (170) into Eq. (167), we derive that with the probability of 1 − δ. We now leverage the following lemma: xlviii Lemma A.17. (Bounded Equivalence of the Wasserstein and Total Variation Distances, [47]) 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 where Diam(X ) def = sup x,y∈X d(x, y), and Sep(X ) def = inf x =y∈X d(x, y).
Consider the metric space (IR D , · 2 ). Note that the empirical measures µ N m , µ N 0 have a countable support Therefore, using the upper bounds in (172) of Lemma A.17 and 171, we conclude that when the step-size is of the order A.4. Proof of Theorem 4.6. The proof is built upon [53]. It suffices to show that for every integer ∈ IN, and for all the test functions Using the triangle inequality, we now have that For the first term on the right side of Eq. (175) we have where (a) is by Jensen's inequality, (b) is by Fatou's lemma, (c) follows from the basic N , as well as the fact that µ * t , f k ≤ b and µ N t , f k ≤ b for all k = 1, 2, · · · , N due to the boundedness of the test functions f 1 , · · · , f ∈ C 3 b (IR D ), and (d) follows from the weak convergence µ N t weakly → µ * t to the mean-field solution (116).
xlix Now, consider the second term on the right hand side of Eq. (175). Due to the exchangability of the initial states (ξ k 0 ) 1≤k≤N , the law of the random variables (ξ k 0 ) 1≤k≤N is also exchangable. Therefore, we obtain that where Π( , N ) is the set of all permutations of numbers selected from {1, 2, · · · , N }. Notice that the right hand side of Eq. (177) is the symmetrized version of the left hand side (178).
By Doob's Martingale inequality [12], the following inequality holds In particular, with the probability of at least 1 − ρ, we have