A Deep Neural Network Optimization Method Via A Traffic Flow Model

A Deep Neural Network Optimization Method Via A Traffic Flow Model Adeyemi Damilare Adeoye African Master’s in Machine Intelligence African Institute for Mathematical Sciences (AIMS) Kigali, Rwanda [email protected] Philipp Petersen Faculty of Mathematics University of Vienna Vienna, Austria [email protected] Abstract We present, via the solution of nonlinear parabolic partial differential equations (PDEs), a continuous-time formulation for stochastic optimization algorithms used for training deep neural networks. Using continuous-time formulation of stochastic differential equations (SDEs), relaxation approaches like the stochastic gradient descent (SGD) method are interpreted as the solution of nonlinear PDEs that arise from modeling physical problems. We reinterpret, through homogenization of SDEs, the modified SGD algorithm as the solution of the viscous Burgers’ equation that models a highway traffic flow. 1 Introduction Deep neural networks (DNNs) have achieved massive success in several areas including image classification, speech recognition, and natural language processing, particularly when there is a nonlinear relationship between the given data and the labels [1–4]. The depths of these networks allow complex data-label relationships to be expressed since each layer nonlinearly transforms the features and therefore effectively filters the information content. This makes them superior to the traditional techniques. Neural networks (NNs) are considered a class of parametric functions each of which explains some data-label relation. In order to train a given network architecture, we need to find the parameters of these functions that generalize well to new unlabeled data. This leads to solving an inverse problem, also called the learning problem, which involves minimizing a non-convex function. With the presence of several hidden layers, the use of deeper network architectures introduces increased network capacity, and thus a high computational complexity, of the parameter estimation problem. Therefore, despite being used for decades, deep learning has only recently revolutionized many applications incited by advances in computational hardware and availability of large datasets. When considering very large datasets, notable sources of limitations encountered in deep learning are the dimensionality and non-convexity of the relative optimization problem. The stochastic gradient descent (SGD) method [5, 6] and its variants [7–10] have been extensively used as a traditional learning algorithm for finding optimal solutions. In these SGD-based algorithms, iterative procedures for parameter estimation are performed using gradient information computed via backpropagation [11–13]. While they are efficient to implement and allow the estimation to scale to massive datasets, the incremental updates to the parameter tend to be slow, especially in the initial stages of the training. Additionally, difficulties stemming from numerical instabilities of the neural network model, and statistical inefficiency (as the algorithms do not use all the information in the dataset) have been reported about these algorithms [14, 15]. A related problem is the observation of vanishing or exploding gradients [16]. Since the gradient represents the sensitivity of the output with respect to a perturbation in the input, a vanishing gradient implies that the output is insensitive with respect to the input, while an exploding gradient implies that the output is unstable with respect to the input. These observations lead to the question of whether there are other methods for training deep neural networks. Attempts to answer this question result to the development of alternative training methods that outperform and eliminate some of the problems associated with the parameter estimation methods based on SGD. In recent times, researchers have developed frameworks inspired by differential equations for training algorithms. The continuous dynamical system approach to deep learning has been explored, where DNNs are being idealized as a discretization of an ordinary differential equation (ODE) [17–23], or as a discretization of a partial differential equation (PDE) [[24, 25]]. Motivated by previous work on local entropy that was introduced to replace the modified loss function of a non-convex optimization problem, [26] reinterpretes relaxation techniques arising in statistical physics as solutions of a visous Hamilton-Jacobi PDE through a stochastic differential equation (SDE) homogenization problem [27]. The authors proved, through a stochastic control interpretation, that the modified algorithm performs better than the SGD. In this work, we reinterprete the modified SGD as the solution of the viscous Burgers’ equation that models a highway traffic flow. The rest of this work is organized as follows: In section 2, we introduce the concept of DNNs, and study how optimization problems arise in machine learning. In particular, we present the use of DNNs in perceptual tasks such as speech or image recognition categorized as supervised learning tasks that involve large-scale, highly nonlinear, and non-convex optimization problems. We also provide a review of the SGD algorithm for optimizing DNNs, and a motivation for seeking modified SGD algorithms such as the one presented in this work. In section 3, we derive the viscous Burgers’ equation as a model for a highway traffic flow, from which we obtain the stochastic differential equation (SDE) satisfied by the parameter updates of the neural network. In section 4, we solve the viscous Burgers’ equation and show that the solution exactly gives the gradient of the local entropy function interpreted as the objective function of the optimization problem that results from training a DNN. In section 5, we introduce the concept of homogenization for SDEs [28], and obtain the continuous-time dynamics, which is then sampled numerically to obtain the discrete-time solution. 2 Continuous-time stochastic optimization in deep neural networks (DNNs) Suppose we are given a set of data points {(θ1, y1), . . . , (θn, yn)}, where θi∈Rdθ, yi∈Rdy, i= 1, . . . , n represent input-output pairs, and nis the number of samples. Here, the input instance θimay represent a word vector, or the feature vector of an object or image, etc, while the output instance yimay represent a real-valued vector for a regression problem, or an integer-valued vector for a classification problem. A deep neural network (DNN) is a family H(called the hypothesis space) of prediction functions defined by H:={y(·;x):x∈W},(2.1) where Wis a set of parameters specified below, for some given y(·;·):Rdθ×W→Rdywhose value, the predicted output vector ˆy∈Rdy, is computed by applying iterative transformations to a given input vector θ∈Rdθthrough the recursion formula z0=θ;zl=σ(Wlzl−1+bl), l = 1, . . . , L, (2.2) where Lis the number of hidden layers (or depth) of the neural network, zlis the hidden state at layer l, the weight matrix Wl∈Rdl×dl−1together with the bias vector bl∈Rdlconstitute the trainable parameters in the l-th layer, and x= (Wl, bl)∈Wrepresents the collection of all the trainable parameters with W=Rdl×dl−1×Rdl. The function σ:R→Rknown as the neuron activation function introduces nonlinearity to the network, and is applied element-wise. Hence, the output zLof the entire network is given by the function y:Rdθ×W→Rdydefined by y:Rdθ×W→Rdy(2.3) (θ, x)7→ y(θ;x):=TLσ(TL···σ(T2σ(T1(θ)))),(2.4) where Tl(x) = Wlθ+bl, and l= 1, . . . , L. 2.1 Non-convex optimization problems in DNNs Let ℓ:Rdy×Rdy→Rbe a given distance metric called the loss function. The optimization problem we want to solve involves choosing the parameter of the neural network that minimizes 2 the distance ℓ(yi,ˆyi)between the predicted output ˆyi=yx(θi)and the true output yi. Assume that the input-output space Rdθ×Rdyis endowed with the unknown joint probability distribution P:Rdθ×Rdy→[0,1], and consider the following optimization problem min x E[ℓ(y, yx(θ))]; E[ℓ(y, yx(θ))] :=∫Rdθ×Rdy ℓ(y, yθ(x))dP(θ, y).(2.5) Since Pis unknown, and that in supervised learning we have access to a set of N∈Nindependently drawn input-output samples {(θi, yi)}N i=1 ⊆Rdθ×Rdyfrom P, we rather attempt to minimize the empirical risk F:W→Rgiven as F(x):=1 N N ∑ i=1 ℓ(yi, yx(θi)).(2.6) Therefore, the optimization problem we wish to solve is the empirical risk minimization problem min xF(x),(2.7) and we write x∗= argmin x F(x),(2.8) as the value of xthat minimizes F(x). The choice of the loss function ℓ(y, ˆy)mostly depends on the type of problem we want to solve. For example, a popular choice of ℓ(y, ˆy)for a regression problem is the quadratic loss function ℓ(y, ˆy) = ky−ˆyk2. For binary classification problems, ℓ(y, ˆy)is usually chosen to be ℓ(y, ˆy) = log(1 + exp(−yˆy)). The optimization problem (2.7) is highly nonlinear and nonconvex in general, making it difficult to solve to global optimality [29, 30]. Gradient-based methods have however been successfully used to compute approximate solutions, where the gradient of the objective in (2.7) is computed by the chain rule using back-propagation. A typical DNN contains very large numbers of parameters in their layers, each of which has to be optimized. This can take several hundreds of hours of training, and weeks of computation on a high-end computer. Popular gradient-based methods that have been used extensively to optimize DNNs are the first-order algorithms such as the gradient descent (GD), and the stochastic gradient descent (SGD) and its variants. Developments in this work are based on continuous-time SDE interpretation of the SGD. 2.2 Gradient-based optimization The GD [31, 32] used for solving optimization problems of the form (2.8) performs iteration of the form xk+1 =xk−δt∇F(xk),(2.9) where δt>0is the step-size (or learning rate) and ∇F(xk)is the gradient of the loss function at the k-th iteration (k∈N). In practice, it is not possible to load all samples into the memory of a single GPU or CPU for computing the entire gradient 1 N∑N i=1 ∇ℓ(yi, Fx(θi)). A more practical choice is SGD or any of its variants. Theoretically, SGD works as follows: at the k-th iteration, randomly select iand update the parameter by xk+1 =xk−δt∇Fi(xk),(2.10) where Fi(x):=ℓ(yi, Fx(θi)). The stochastic nature of SGD arises from the approximation of the gradient over a subset of the data points. Also common is “mini-batch" SGD, where we choose a random subset Ik⊆1, . . . , N, |Ik|=BNto write (2.7) as a finite-sum optimization problem: min xFmb(x); Fmb(x):=1 B B ∑ i=1 Fi(x),(2.11) 3 where each Fi(x)represents the sum of training loss for a “mini-batch" of training samples, and B is the total number of mini-batches (smaller than the total number of training samples N) [5]. We then perform the iteration xk+1 =xk−δt∇Fmb(xk).(2.12) One reason for the use of SGD instead of GD is the faster convergence in a “memory-constraint” system where it is difficult to find a parallel way to load and process all samples in a single machine [33]. In such a system, batch GD has to scan through the entire training set during one iteration before taking the next iteration, whereas the SGD starts to converge right away from the first sample, hence a noticeable faster convergence speed in the early iterations. 2.3 Continuous-time stochastic optimization While it converges very fast in the early iterations, the SGD algorithm slows down in later iterations and thus struggles to reach an accurate solution. This is because although the stochastic gradient is an unbiased estimator of the gradient: E[∇Fmb(x)] = ∇F(x),(2.13) it may end up with a high variance. In order to get the best of both the fast initial convergence of SGD and the steady linear convergence of GD to a great extent, we further assume a bounded variance for the stochastic gradient, that is, for all iterates xkof the SGD, there exists νmb >0such that [34] E[k∇Fmb(x)−∇F(x)k2]≤2νmb.(2.14) Here, the "mb" in the notation νmb only indicates the noise is coming from the mini-batch gradients. Henceforth, we will simply write ν. We start with the standard observation that the discrete-time dynamics in (2.10) can be viewed as a discretization of the continuous-time Langevin dynamics described by the Itˆ o stochastic differential equation (SDE) [35, 36] dx(t) = −∇F(x(t)) dt+√2νdW(t), t ≥0, x(0) = x0,(2.15) where W(t)∈Rmis the standard m-dimensional Wiener process. Under suitable assumptions on F, it can be shown that the unique invariant distribution of (2.15) is the Gibbs distribution ρ∞(ν;x)∝exp (−1 2νF(x)),(2.16) and that using the idea of simulated annealing to evade metastability1, the distributions of x(t) are made to converge rapidly to ρ∞as t→ ∞ [37–39]. Moreover from (2.16), ρ∞concentrates around the minimizers of F(x)for small enough values of ν[40, 41]. Consequently, as the Gibbs distribution turns out to be an approximate minimizer for the empirical risk (2.6), we will observe that the resulting modified SGD algorithm tracks the Langevin dynamics in a suitable sense. 2.4 Regularization of the loss function As discussed, training DNNs involves modeling a high-dimensional nonconvex loss function. In order to improve the computational efficiency with a large number of training data, the stochastic estimation of the gradients of these loss functions are carried out using the SGD algorithm or any of its variants. However, to improve performance of the optimization algorithm and the generalization of the model, it is generally required to introduce a regularization term in the loss function. Several regularization techniques that show outstanding results have been developed for optimizing nonconvex empirical risks of deep networks. These regularization techniques are categorized as either explicit, e.g. weight decay, sparsity constraint, and entropy minimization or implicit, e.g. noise injection, dropout, learning rate decay, model ensemble, and batch normalization [42]. 1Metastability is a phenomenon that describes a set of states of a dynamic system that can persist for a long time. In other words, a dynamic system is said to exhibit metastability if the convergence of the probability density of its states takes an exponentially long time. 4 Our focus in this work is to derive, through homogenization of SDEs, a form of the regularized SGD presented by [26], from the solution of the viscous Burgers’ equation that models a highway traffic flow. [26] models the loss function of a class of deep networks as a spin glass Hamiltonian system, and with results that show good empirical performance, they introduce a smoothing technique that stems from imposing PDE regularity which contributes to the analysis of the geometric property of the energy landscape of the networks. By combining tools from the theory of PDEs and tools from modeling stochastic first-order algorithms as continuous-time stochastic processes, we are able to provide insights on how new approaches for solving optimization problems that arise in DNNs could be developed starting with mathematical models that describe real life situations. As shown, this process can often lead to improvements in state-of-the-art DNN algorithms. 3 The viscous Burgers’ equation as the model for a highway traffic flow Consider the viscous Burgers’ equation ut+u·∇u=ν∆ufor (x, t)∈Rd×[0,∞),(3.1) where ν∈R+is a positive coefficient, and ∇= n ∑ i=1 ∂ ∂xiand ∆ = n ∑ i=1 ∂2 ∂x2 i are respectively the gradient and Laplacian in the space variables. We will express the model for a highway traffic flow as the Burgers’ equation (3.1). The traffic flow model will be derived using a continuous, deterministic approach. By continuous, we mean that the traffic will be observed from a distance far enough that the size and behaviour of a car becomes negligible, and will rather deal with average concepts, such as the number of cars per mile and the average speed of the car. By deterministic, we mean to consider the situation where the change of the traffic with time is exclusively specified by the model. Hence, in this setting, each car entering and exiting a highway will be replaced by a continuous density function ρ:Rd×(0,∞)→R (x, t)7→ ρ(x, t).(3.2) With the additional constraint that follows from the realization that the typical car speed depends on the level of traffic congestion (the car density), that is q=Q(ρ), the differential form of the equation for car (nonlinear) conservation is given by the closed system ρt+∇q= 0 (3.3) q=Q(ρ)(3.4) where the traffic density ρ(x, t)∈Rdmeasures the number of cars per unit length of the road at position xand time t, and the traffic flow Q(x, t)∈Rd, also known as the flux, represents the number of cars crossing position xper unit time at time t[43]. A third fundamental (dependent) quantity that helps to describe the traffic flow is the velocity u∈Rdof a car at any point on the highway, and is assumed to depend only on the traffic density, that is, u=U(ρ). These three fundamental variables are related by the identity Q=ρU. (3.5) The prescribed function Q(ρ)depends on the characteristics of the road, cars and drivers, such as the car speed, speed limits, weather conditions, and the road conditions. Equations (3.3) and (3.4) together gives the first-order, non-linear, partial differential equation for traffic flow in the form of a closed equation for the single variable ρ: ρt+Q′(ρ)∇ρ= 0 for (x, t)∈Rd×(0,∞).(3.6) As the traffic flow model is based on the first-order approximation, we could define a better approximation for qby assuming it is a function of the density ρas well as its gradient ρx, and define q=Q(ρ)−ν∇ρ. (3.7) Consider a linear relationship between velocity and traffic density [43]: u(ρ) = um(1−ρ ρm),(3.8) 5 where umis the maximum velocity, and ρmis the maximum density. We have (from (3.5)) Q′(ρ) = d dρ(ρu) =d dρ [ρum(1−ρ ρm)] =um(1−2ρ ρm).(3.9) Also, from (3.7) ∇q=Q′(ρ)∇ρ−ν∆ρ, or (using (3.3)) ρt+Q′(ρ)∇ρ−ν∆ρ= 0.(3.10) Assume w, ρ ∈C∞(Rd×(0,∞)), and define w(ρ) = Q′(ρ) = um(1−2ρ ρm),(3.11) then ρ=ρm 2um (1 −w), ρt=−ρmwt 2um ,∇ρ=−ρm∇w 2um ,and ∆ρ=ρm∆w 2um .(3.12) Substituting (3.12) into (3.10) and multiplying the resulting equation by w′(ρ) = −2um ρm, we get the viscous Burgers’ equation for the highway traffic flow: wt+w∇w=ν∆wfor (x, t)∈Rd×(0,∞).(3.13) 4 Nonlinear parabolic PDEs and the local entropy Introduced for studying the energy landscape of neural network learning algorithms, the local entropy is a modified form of the objective function F(x)of the neural network parameters [see for example, 44]. While extending the notion of local entropy to continuous variables, [27] replaced the objective function F(x), defined for x∈Rd, with the function Fγ(x)given by Fγ(x) = 2νlog (G2νγ ∗exp (−1 2vF(x)))(x∈Rd),(4.1) where Gγ(x) = (2πγ)d/2exp (−|x|2 2γ)is the Green’s function (or heat kernel), and f∗gdefines the convolution of fand g: (f∗g)(x):=∫∞ −∞ f(t)g(x−t)dt. (4.2) Assuming F∈C∞(Rd), we will derive a closed form solution w(x, t)of the initial-value problem {wt+w∇w=ν∆win (x, t)∈Rd×(0,∞), w(x, 0) = F(x)on (x, t)∈Rd×{t= 0},(4.3) for the viscous Burgers’ equation (3.13), by a non-linear change of variables, due to the transformation of Cole [45] and Hopf [46], which turns it into an initial-value problem for a linear diffusion equation (the heat equation). We will then show that this solution w(x, t)recovers the gradient of Fγ(x). 4.1 Proposition. The local entropy function Fγ(x)defined by (4.1) is the solution of the initialvalue problem for the viscous Hamilton-Jacobi equation ϕt+(∇ϕ)2 2=ν∆ϕfor (x, t)∈Rd×(0,∞),(4.4) with initial values ∇ϕ(x, 0) = F(x)[26]. Moreover, the gradient ∇Fγ(x)solves the viscous Burgers’ equation (4.3). 6 Proof. We start by writing the viscous Burgers’ equation (3.13) as wt+∇(w2 2−ν∇w)= 0.(4.5) As it is not yet specified, let us assume for the moment that wis a smooth solution of (4.5), we set ψ:=ϕ(w), where ϕ:Rd×(0,∞)→Ris a smooth function. We want to choose ϕso that ψsolves a linear equation. Let wsatisfy the irrotational condition ∇×W= d ∑ i,j=1 (∂wj ∂xi−∂wi ∂xj)(ei∧ej)= 0,(4.6) where ei, i ≤1≤dare basis of Rd, and for vectors u, v ∈Rd,u∧vdefines the wedge/exterior product given by u∧v=u⊗v−v⊗u. (4.7) We define ϕ(x, t):=∫x −∞ w(z, t)dz, =⇒w(x, t) = ∇ϕ(x, t),(4.8) and write the Hopf-Cole transformation [46, 45] ϕ=−2νlog(ψ),(4.9) where log(v) = (log vi)1≤i≤dfor a vector v∈Rd. We have                    ϕt=−2ν ψψt, ∇ϕ=−2ν ψ∇ψ, ∆ϕ=−2ν ψ∆ψ+ 2ν(∇ψ ψ)2 . (4.10) Putting (4.10) in (4.8) and (4.5), we get the linear d-dimensional heat equation ψt=ν∆ψ(x∈Rd, t > 0),(4.11) with the initial condition ψ(x, 0) = g1(x). In terms of the initial value w(x, 0) = F(x), and from (4.10), the initial data for ψbecomes w(x, 0) = F(x) = −2ν∇ψ(x, 0) ψ(x, 0) .(4.12) Integrating, we get ψ(x, 0) = g1(x) = exp (−1 2ν∫x 0 F(z)dz).(4.13) The general solution to the resulting initial-value problem for ψis given as convolution of the initialvalue with the heat kernel H(x, t):            ψ(x, t) = ∫∞ −∞ H(x−y, t)ψ(y, 0)dy for (x, t)∈Rd×(0,∞), H(x, t) = 1 √4πνt exp (−x2 4νt). (4.14) By normalizing ϕ(x, t)by ϕ(0, t) = 0, we obtain from ψ(x, t), the solution for the initial-value problem for the viscous Hamilton-Jacobi equation (4.4): ϕ(x, t) = −2νlog   ∫∞ −∞ 1 √4πνt exp [−(x−y)2 4νt −1 2ν∫y 0 F(z)dz]dy   =:Ft(x).(4.15) 7 Finally, the solution of the viscous Burgers’ equation (3.13) is w(x, t) = ∇ϕ(x, t) = ∫∞ −∞ x−y texp (−(x−y)2 4νt −1 2ν∫y 0F(z)dz)dy ∫∞ −∞ exp (−(x−y)2 4νt −1 2ν∫y 0F(z)dz)dy ,(4.16) where x∈Rd, t > 0.■ 4.2 Remark. The solution ψ(x, t)of the heat equation (4.11), given by (4.14), can be written in a more convenient form as ψ(x, t) = 1 2√πνt ∫∞ −∞ exp (−G 2ν)dy, (4.17) where G(x, y;t) = ∫y 0 F(z)dz +(x−y)2 2t(x, y ∈Rd, t > 0),(4.18) and the gradient is given by ∇ψ(x, t) = −1 4ν√πνt ∫∞ −∞ x−y texp (−G 2ν)dy. (4.19) 4.3 Remark. We can equivalently write (4.16) as w(x, t) = ∇ϕ(x, t) = ∫∞ −∞ y−x tρ∞ 1(dy;x) = ∫∞ −∞ ∇F(x−y)ρ∞ 2(dy;x),(4.20) where x∈Rd, t > 0. The functions ρ∞ 1(y;x), ρ∞ 2(y;x)are invariant distributions given by              ρ∞ 1(y;x) = 1 Z1 exp (−1 2νF(y)−|x−y|2 4νt ), ρ∞ 2(y;x) = 1 Z2 exp (−1 2νF(x−y)−|y|2 4νt ), (4.21) and              Z1(x) = −∫∞ −∞ exp (−1 2νF(y)−|x−y|2 4νt )dy, Z2(x) = −∫∞ −∞ exp (−1 2νF(x−y)−|y|2 4νt )dy, (4.22) are normalizing constants. 5 The local entropy via homogenization for SDEs We assume a smooth, steady, periodic, and incompressible (or divergence-free) flow velocity, that is, w(x, t) = −h(x), where his a smooth periodic vector field with ∇ · h(x)=0. Consider F(x) = gϵ:=g(ϵx), with ϵ > 0, so that initially the velocity is slowly varying in space. We expect that the velocity will not vary significantly on small length and time scales, hence to see the effective behaviour of (3.13), we need to look at large length and time scales. Furthermore, if haverages over the unit cube (not necessarily to zero), and if we consider the rescaling x→x/ϵ and t→t/ϵp(p > 0), so as to bring out order-one variations in velocity on large length and time scales, homogenization techniques enable us to show that the rescaled velocity field wϵ(x, t):= w(x/ϵ, t/ϵp)converges, for ϵ1, to the solution w(x, t)of the viscous Burgers’ equation (4.3). In this setting, equation (4.3) becomes      1 ϵ2−pwϵ t−1 ϵhϵwϵ x=νwϵ xx in (x, t)∈Rd×(0,∞), wϵ(x, 0) = g(ϵx)on (x, t)∈Rd×{t= 0}, (5.1) 8 where hϵ=h(x/ϵ). As to what time scale to consider (choice of p) in order to observe interesting dynamical behaviour, it is intuitive to study (5.1) for different values of pin cases when haverages to zero and when it does not average to zero. For the purpose of this work, with v > 0, we are considering the case when haverages to zero, hence the choice p= 2, and equation (5.1) becomes      wϵ t−1 ϵhϵwϵ x=νwϵ xx in (x, t)∈Rd×(0,∞), wϵ(x, 0) = g(ϵx)on (x, t)∈Rd×{t= 0}, (5.2) Let W(t)∈Rmdenote an independent standard m-dimensional Wiener process, h:Z=Rn→ Rna smooth periodic vector-valued function, and β:Z=Rn→Rn×ma smooth matrix-valued or scalar-valued function. Consider the Itˆ o SDE dz dt =h(z) + β(z)dW dt , z(0) = z0,(5.3) which is interpreted as an integral equation for z(t)∈C(R+,Z): z(t) = z0+∫t 0 h(z(s))ds +∫t 0 β(z(s))dW(s).(5.4) We introduce the auxiliary dynamics y=x/ϵ and consider the case where z= (xT, yT)T, with x∈ X =Rd,y∈ Y =Rn−d. We then write the system of SDEs      dx(s) = h(x, y)ds, x(0) = x0, dy(s) = 1 ϵg(x, y)ds +1 √ϵβ(x, y)dW(s), y(0) = y0,(5.5) where we have scaled time to s=ϵt, and gis a sufficiently smooth function. Both x(s)and y(s)contain fast dynamics, but with the homogenization parameter ϵ > 0, the dynamics in y(s)is provided with a faster time-scale than that in x(s). The system (5.5) corresponds to the backward Kolmogorov equation    ∂w ∂t =1 ϵL0win (x, y, t)∈ X ×Y ×(0,∞), w(x, 0) = g(ϵx)on (x, y, t)∈ X ×Y ×{t= 0}, (5.6) where the generator L0is defined as L0w=g(x, y)·∇yw+1 2B∆yw, (5.7) with B(x, y) = β(x, y)β(x, y)T. The corresponding Fokker-Planck (forward Kolmogorov) equation is    ∂w ∂t =1 ϵL∗ 0win (x, y, t)∈ X ×Y ×(0,∞), w(x, 0) = g(ϵx)on (x, y, t)∈ X ×Y ×{t= 0}, (5.8) where the L2(Rd)-adjoint operator L∗ 0is defined as L∗ 0w=−∇y(g(x, y)w) + 1 2∆y(B(x, y)w).(5.9) Next, we eliminate the y(s)dependence in the Kolmogorov equation (5.6) to identify only a simplified equation for the dynamics of x(s). To do this, we make an ergodicity assumption: we asume that for each fixed x,L0has a one-dimensional null space, and L∗ 0has a null space spanned by the invariant distribution ρ∞ 1(y;x). That is, L0satisfies {L01(y) = 0, L∗ 0ρ∞ 1(y;x) = 0,(5.10) where 1(y)denotes constants (equal to 1a.e.) in y(s)and ρ∞ 1(y;x)is the density of an ergodic measure µx(dy) = ρ∞ 1(y;x)dy. 9