# A Continuous-time Stochastic Gradient Descent Method for Continuous Data

Kexin Jin<EMAIL_ADDRESS>
Department of Mathematics
Princeton University
Princeton, NJ 08544-1000, USA Jonas Latz<EMAIL_ADDRESS>
School of Mathematical and Computer Sciences
Heriot-Watt University
Edinburgh, EH14 4AS, United Kingdom Chenguang Liu<EMAIL_ADDRESS>
Delft Institute of Applied Mathematics
Technische Universiteit Delft
2628 Delft, The Netherlands Carola-Bibiane Schönlieb<EMAIL_ADDRESS>
Department of Applied Mathematics and Theoretical Physics
University of Cambridge
Cambridge, CB3 0WA, United Kingdom

###### Abstract

Optimization problems with continuous data appear in, e.g., robust machine
learning, functional data analysis, and variational inference. Here, the
target function is given as an integral over a family of (continuously)
indexed target functions – integrated with respect to a probability measure.
Such problems can often be solved by stochastic optimization methods:
performing optimization steps with respect to the indexed target function with
randomly switched indices.

In this work, we study a continuous-time variant of the stochastic gradient
descent algorithm for optimization problems with continuous data. This so-
called stochastic gradient process consists in a gradient flow minimizing an
indexed target function that is coupled with a continuous-time index process
determining the index. Index processes are, e.g., reflected diffusions, pure
jump processes, or other Lévy processes on compact spaces. Thus, we study
multiple sampling patterns for the continuous data space and allow for data
simulated or streamed at runtime of the algorithm. We analyze the
approximation properties of the stochastic gradient process and study its
longtime behavior and ergodicity under constant and decreasing learning rates.
We end with illustrating the applicability of the stochastic gradient process
in a polynomial regression problem with noisy functional data, as well as in a
physics-informed neural network.

Keywords: Stochastic optimization, functional data analysis, robust learning,
arbitrary data sources, Markov processes

##  1 Introduction

The training of a machine learning model is often represented through an
optimization problem. The goal is to calibrate the model’s parameters to
optimize its goodness-of-fit with respect to training data. The goodness-of-
fit is usually quantified through a loss function that sums up the losses from
the misrepresentation of every single training data set, see, e.g., Goodfellow
et al. (2016) or Hansen (2010); Cam (1990) for similar optimization problems
in imaging and statistics. In ‘big data’ settings, the sum of these loss
functions consists of thousands or millions of terms, making classical
optimization methods computationally infeasible. Thus, efficiently solving
optimization problems of this form has been a focus of machine learning and
optimization research in the past decades. Here, methods often build upon the
popular stochastic gradient descent method.

Originally, stochastic gradient descent was proposed by Robbins and Monro
(1951) to optimize not only sums of loss functions, but also expectations of
randomized functions.111Actually, the ‘stochastic approximation method’ of
Robbins and Monro (1951) aims at finding roots of functions that are given as
expectations of randomized functions. The method they construct resembles
stochastic gradient descent for a least squares loss function. Of course, a
normalized sum is just a special case of an expected value, making stochastic
gradient descent available for the kind of training problem described above.
Based on stochastic gradient descent ideas, improved algorithms have been
proposed for optimizing sums of loss functions, such as Chambolle et al.
(2018); Defazio et al. (2014); Duchi et al. (2011). Unfortunately, these
methods often specifically target sums of loss functions and are often
infeasible to optimize general expected values of loss functions.

The optimization of expected values of loss functions appears in the presence
of countably infinite and continuous data in functional data analysis and non-
parametric statistics (e.g., Sinova et al. (2018)), physics-informed deep
learning (e.g., Raissi et al. (2019)), inverse problems (e.g., Bredies and
Lorenz (2018)), and continuous data augmentation/adversarial robustness (e.g.,
Cohen et al. (2019); Shorten and Khoshgoftaar (2019); Pinto et al. (2017)).
Some of these problems are usually studied after discretising the data.
Algorithms for discrete data sometimes deteriorate at the continuum limit,
i.e. as the number of data sets goes to infinity. Thus, we prefer studying the
continuum case immediately. Finally, ‘continuous data’ can also refer to
general noise models. Here, expected values are minimized in robust
optimization (e.g., Nemirovski et al. (2009)), variational Bayesian inference
(e.g., Cherief-Abdellatif (2019)), and optimal control (e.g., May et al.
(2013)). Overall, the optimization of general expected values is a very
important task in modern data science, machine learning, and related fields.

In this work, we study stochastic gradient descent for general expected values
in a continuous-time framework. We now proceed with the formal introduction of
the optimization problem, the stochastic gradient descent algorithm, its
continuous-time limits, and current research in this area.

###  1.1 Problem setting and state of the art

We study optimization problems of the form

$\min_{\theta\in X}\Phi(\theta):=\int_{S}f(\theta,y)\pi(\mathrm{d}y),$ (1)

where $X:=\mathbb{R}^{K}$, $S$ is a Polish space, $f:S\times
X\rightarrow\mathbb{R}$ is a measurable function that is continuously
differentiable in the first variable, and $\pi$ is a probability measure on
$S$. Moreover, we assume that the integral above always exists. We refer to
$X$ as _parameter space_ , $S$ as _index set_ , $\Phi$ as _full target
function_ , and $f$ as _subsampled target function_. In these optimization
problems, it is usually impossible or intractable to evaluate the integral
$\Phi(\theta)$ or its gradient $\nabla\Phi(\theta)$ for $\theta\in X$. Hence,
traditional optimization algorithms, such as steepest gradient descent or
Newton methods are not applicable.

As mentioned above, it is possible to employ stochastic optimization methods,
such as the _stochastic gradient descent (SGD)_ method, see Kushner and Yin
(2003); Robbins and Monro (1951). The stochastic gradient descent method for
(1) proceeds through the following discrete-time dynamic that iterates over
$n\in\mathbb{N}:=\\{1,2,\ldots\\}$:

$\displaystyle\theta_{n}=\theta_{n-1}-\eta_{n}\nabla_{\theta}f(\theta_{n-1},y_{n}),$
(2)

where $y_{1},y_{2},\ldots\sim\pi$ independent and identically distributed
(i.i.d.), $(\eta_{n})_{n=1}^{\infty}\in(0,\infty)^{\mathbb{N}}$ is a non-
increasing sequence of _learning rates_ , and $\theta_{0}\in X$ is an
appropriate initial value. Hence, SGD is an iterative method that employs only
the gradient of the integrand $f$, but not $\Phi$. SGD converges to the
minimizer of $\Phi$, if $\eta_{n}\rightarrow 0$, as $n\rightarrow\infty$,
sufficiently slowly, and $f(\cdot,y)$ $(y\in S)$ is strongly convex; see,
e.g., Bubeck (2015). Additionally, SGD is used in practice also for non-convex
optimization problems and with constant learning rate. The constant learning
rate setting is popular especially due to its regularizing properties; see Ali
et al. (2020); Smith et al. (2021).

To understand, improve, and study discrete-time dynamical systems, it is
sometimes advantageous to represent them in continuous time, see, e.g. the
works by de Wiljes et al. (2018); Kovachki and Stuart (2021); Trillos and
Sanz-Alonso (2020). Continuous-time models allow us to concentrate on the
underlying dynamics and omit certain numerical considerations. Moreover, they
give us natural ways to construct new, efficient algorithms.

The discrete-time dynamic in (2) is sometimes represented through a
continuous-time diffusion process, see Ali et al. (2020); Li et al. (2019,
2017); Mandt et al. (2016, 2017); Wojtowytsch (2021):

$\mathrm{d}\theta_{t}=-\nabla\Phi(\theta_{t})\mathrm{d}t+\sqrt{\eta(t)}\Sigma(\theta_{t})^{1/2}\mathrm{d}W_{t},$

where
$\Sigma(\theta)=\int(\nabla_{\theta}f(\theta;y)-\nabla_{\theta}\Phi(\theta))\otimes(\nabla_{\theta}f(\theta;y)-\nabla_{\theta}\Phi(\theta))\pi(\mathrm{d}y)$,
$(W_{t})_{t\geq 0}$ is a $K$-dimensional Brownian motion, and
$(\eta(t))_{t\geq 0}$ is an interpolation of the learning rate sequence. While
this diffusion approach is suitable to describe the dynamic of the moments of
SGD, it does not immediately allow us to construct new stochastic optimization
algorithms, as the system depends on the inaccessible $\nabla\Phi$.

$-1$$1$$.6$$.2$$-.2$$-.6$$0$$10$$(\mathbf{i}(t))_{t\geq
0}$$(\theta^{\dagger}_{t})_{t\geq 0}$$0$$10$ Figure 1: Cartoon of the
stochastic gradient process $(\theta^{\dagger}_{t})_{t\geq 0}$ with index
process $(\mathbf{i}(t))_{t\geq 0}$ on the discrete index set
$S:=\\{-1,-0.6,\ldots,1\\}$. The index process is a Markov pure jump process
on $S$. The process $(\theta^{\dagger}_{t})_{t\geq 0}$ aims at optimizing the
$\mathrm{Unif}(S)$-integral of the subsampled target functional is
$f(\theta,y):=\frac{1}{2}(\theta-y^{2})^{2}$ $(\theta\in X:=\mathbb{R},y\in
S)$.

A continuous-time representation of stochastic gradient descent that does not
depend on $\Phi$ has recently been proposed by Latz (2021). This work only
considers the discrete data case, i.e., $S$ is finite and
$\pi:=\mathrm{Unif}(S)$. SGD is represented by the _stochastic gradient
process_ $(\theta^{\dagger}_{t})_{t\geq 0}$. It is defined through the coupled
dynamical system

$\mathrm{d}\theta_{t}^{\dagger}=-\nabla_{\theta}f(\theta_{t}^{\dagger};\mathbf{i}(t))\mathrm{d}t,$
(3)

where $(\mathbf{i}(t))_{t\geq 0}$ is a suitable continuous-time Markov process
on $S$, which we call _index process_. Hence, the process
$(\theta_{t}^{\dagger})_{t\geq 0}$ represents gradient flows with respect to
the subsampled target functions that are switched after random waiting times.
The random waiting times are controlled by the continuous-time Markov process
$(\mathbf{i}(t))_{t\geq 0}$. We show an example of the coupling of exemplary
process $(\mathbf{i}(t))_{t\geq 0}$ and $(\theta_{t}^{\dagger})_{t\geq 0}$ in
Figure 1. The setting is $S:=\\{-1,-0.6,\ldots,1\\}$, $\pi:=\mathrm{Unif}(S)$,
$X:=\mathbb{R}$, and $f(\theta,y):=\frac{1}{2}(\theta-y^{2})^{2}$ $(\theta\in
X,y\in S)$. There, we see that the sample path of
$(\theta_{t}^{\dagger})_{t\geq 0}$ is piecewise smooth, with non-smooth
behavior at the jump times of $(\mathbf{i}(t))_{t\geq 0}$.

If the process $(\mathbf{i}(t))_{t\geq 0}$ is homogeneous-in-time, the
dynamical system represents a constant learning rate. Inhomogeneous
$(\mathbf{i}(t))_{t\geq 0}$ with decreasing mean waiting times, on the other
hand, model a decreasing learning rate. Under certain assumptions, the process
$(\theta_{t}^{\dagger})_{t\geq 0}$ converges to a unique stationary measure
when the learning rate is constant or to the minimizer of $\Phi$ when the
learning rate decreases.

###  1.2 This work.

We now briefly introduce the continuous-time stochastic gradient descent
methods that we study throughout this work. Then, we summarize our main
contributions and give a short paper outline.

In the present work, we aim to generalize the dynamical system (3) to include
more general spaces $S$ and probability measures $\pi$ – studying the more
general optimization problems of type (1). We proceed as follows: We define a
stationary continuous-time Markov process $(V_{t})_{t\geq 0}$ on $S$ that is
geometrically ergodic and has $\pi$ as its stationary measure. This process
$(V_{t})_{t\geq 0}$ is now our _index process_. Similarly to (3), we then
couple $(V_{t})_{t\geq 0}$ with the following gradient flow:

$\mathrm{d}\theta_{t}=-\nabla_{\theta}f(\theta_{t},V_{t})\mathrm{d}t.$ (4)

Note that the index process $(V_{t})_{t\geq 0}$ can be considerably more
general than the Markov jump processes studied by Latz (2021); we discuss
examples below and in Section 2. As the dynamical system (4) contains the
discrete version (3) as a special case, we refer to $(\theta_{t})_{t\geq 0}$
also as _stochastic gradient process_.

We give an example for $(\theta_{t},V_{t})_{t\geq 0}$ in Figure 2. There, we
consider $S:=[-1,1]$, $\pi:=\mathrm{Unif}[-1,1]$, $X:=\mathbb{R}$, and
$f(\theta,y):=\frac{1}{2}(\theta-y^{2})^{2}$ $(\theta\in X,y\in S)$. A
suitable choice for $(V_{t})_{t\geq 0}$ is a reflected Brownian motion on
$[-1,1]$. Although it is coupled with $(V_{t})_{t\geq 0}$, the process
$(\theta_{t})_{t\geq 0}$ appears to be relatively smooth. This may be due to
the smoothness of the subsampled target function $f$. Moreover, we note that
the example in Figure 1 is a discretized data version of the example here in
Figure 2.

$-1$$1$$0$$0$$10$$(V_{t})_{t\geq 0}$$(\theta_{t})_{t\geq 0}$$0$$10$ Figure 2:
Cartoon of the stochastic gradient process $(\theta_{t})_{t\geq 0}$ with index
process $(V_{t})_{t\geq 0}$, with continuous index set $S:=[-1,1]$. The index
process is a reflected Brownian motion on $S$.

More similarly to the discrete data case (3), one could also choose
$(V_{t})_{t\geq 0}$ to be a Markov pure jump process on $S$ that has $\pi$ as
a stationary measure. Indeed, the reflected Brownian motion was constructed
rather artificially. Sampling from $\mathrm{Unif}[-1,1]$ is not actually
difficult in practice and we just needed a way to find a continuous-time
Markov process that is stationary with respect to $\mathrm{Unif}[-1,1]$.
However, there are cases, where one may not be able to sample independently
from $\pi$. For instance, $\pi$ could be the measure of interest in a
statistical physics simulation or Bayesian inference. In those cases, Markov
chain Monte Carlo methods are used to approximate $\pi$ through a Markov chain
stationary with respect to it, see, e.g., Robert and Casella (2004). In other
cases, the data might be time series data that is streamed at runtime of the
algorithm – a related problem has been studied by Sirignano and Spiliopoulos
(2017). Hence, in the work, we also discuss stochastic optimization in those
cases or – more generally – stochastic optimization with respect to data from
arbitrary sources.

As the index process $(V_{t})_{t\geq 0}$ is stationary, the stochastic
gradient process as defined above would, again, represent the situation of a
constant learning rate $(\eta_{n})_{n=1}^{\infty}$. However, as before we
usually cannot hope for convergence to a stationary point if there is not a
sense of a decreasing learning rate. Hence, we need to introduce an
inhomogeneous variant of $(V_{t})_{t\geq 0}$ that represents a decreasing
learning rate.

We start with the stochastic process $(V_{t}^{\rm dc})_{t\geq 0}$ that
represents the index process associated to the discrete-time stochastic
gradient descent dynamic (2) with constant learning rate parameter $\eta=1$.
This index process is given by

$V_{t}^{\rm
dc}=\sum_{n=1}^{\infty}y_{n}\mathbf{1}[t\in[n-1,n)]=\sum_{n=1}^{\infty}y_{n}\mathbf{1}[t-n+1\in[0,1)]\qquad(t\geq
0),$

where $y_{1},y_{2},\ldots\sim\pi$ i.i.d. and $\mathbf{1}[\cdot]$ is the
indicator function: $\mathbf{1}[{\rm true}]=1$ and $\mathbf{1}[{\rm
false}]=0$. We now want to turn the process $(V_{t}^{\rm dc})_{t\geq 0}$ into
the index process $(V_{t}^{\rm dd})_{t\geq 0}$ that represents a decreasing
learning rate $(\eta_{n})_{n=1}^{\infty}$. It is defined through:

$(V_{t}^{\rm dd})_{t\geq
0}=\sum_{n=1}^{\infty}y_{n}\mathbf{1}\left[t\in\left[H_{n-1},H_{n}\right)\right]=\sum_{n=1}^{\infty}y_{n}\mathbf{1}\left[\frac{t-H_{n-1}}{\eta_{n}}\in\left[0,1\right)\right],$

where we denote $H_{n}:=\sum_{m=1}^{n}\eta_{m}$. Hence, we can represent
$(V_{t}^{\rm dd})_{t\geq 0}:=(V_{\beta(t)}^{\rm dc})_{t\geq 0}$, where
$\beta:[0,\infty)\rightarrow[0,\infty)$ is given by

$\beta(t)=\sum_{n=1}^{\infty}\frac{t+n-1-H_{n-1}}{\eta_{n}}\mathbf{1}\left[t\in\left[H_{n-1},H_{n}\right)\right]\qquad(t\geq
0)$ (5)

is a piecewise linear, non-decreasing function with
$\beta(t)\rightarrow\infty,$ as $t\rightarrow\infty$.

Following this idea, we turn our homogeneous index process $(V_{t})_{t\geq 0}$
that represents a constant learning rate into an inhomogeneous process with
decreasing learning rate using a suitable rescaling function $\beta$. In that
case, we obtain a stochastic gradient process of type

$\mathrm{d}\xi_{t}=-\nabla_{\xi}f(\xi_{t},V_{\beta(t)})\mathrm{d}t,$

which we will use to represent the stochastic gradient descent algorithm with
decreasing learning rate. Note that while we require $\beta$ to satisfy
certain conditions that ensure the well-definedness of the dynamical system,
it is not strictly necessary for it to be of the form (5). Actually, we later
assume that $\beta$ is smooth.

The main contributions of this work are the following:

  * •

We study stochastic gradient processes for optimization problems of the form
(1) with finite, countably infinite, and continuous index sets $S$.

  * •

We give conditions under which the stochastic gradient process with constant
learning rate is well-defined and that it can approximate the full gradient
flow $\mathrm{d}\zeta_{t}=-\nabla\Phi(\zeta_{t})\mathrm{d}t$ at any accuracy.
In addition, we study the geometric ergodicity of the stochastic gradient
process and properties of its stationary measure.

  * •

We study the well-definedness of the stochastic gradient process with
decreasing learning rate and give conditions under which the process converges
to the minimizer of $\Phi$ in the optimization problem (1).

  * •

In numerical experiments, we show the suitability of our stochastic gradient
process for (convex) polynomial regression with continuous data and the (non-
convex) training of physics-informed neural networks with continuous sampling
of function-valued data.

This work is organized as follows. In Section 2, we study the index process
$(V_{t})_{t\geq 0}$ and give examples for various combinations of index spaces
$S$ and probability measures $\pi$. Then, in Sections 3 and 4, we analyze the
stochastic gradient process with constant and decreasing learning rate,
respectively. In Section 5, we review discretization techniques that allow us
to turn the continuous dynamical systems into practical optimization
algorithms. We employ these techniques in Section 6, where we present
numerical experiments regarding polynomial regression and the training of
physics-informed neural networks. We end with conclusions and outlook in
Section 7.

##  2 The index process: Feller processes and geometric ergodicity

Before we define the stochastic gradient flow, we introduce and study the
class of stochastic processes $(V_{t})_{t\geq 0}$ that can be used for the
data switching in (4). Moreover, we give an overview of appropriate processes
for various measures $\pi$. For more background material on (continuous-time)
stochastic processes, we refer the reader to the book by Revuz and Yor (2013),
the book by Liggett (2010), and other standard literature.

Let $\mathcal{S}=(S,m)$ be a compact Polish space and

$\displaystyle\Omega=\\{\omega:[0,\infty)\to S\ |\ \omega\ \text{is right
continuous with left limits.}\\}.$

We consider a filtered probability space
$(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},(\mathbb{P}_{x})_{x\in S})$,
where $\mathcal{F}$ is the smallest $\sigma$-algebra on $\Omega$ such that the
mapping $\omega\to\omega(t)$ is measurable for any $t\geq 0$ and the
filtration $\mathcal{F}_{t}$ is right continuous. Let $(V_{t})_{t\geq 0}$ be a
$(\mathcal{F}_{t})_{t\geq 0}$ adapted stochastic process from $\Omega$ to $S$.
We assume that $(V_{t})_{t\geq 0}$ is Feller with respect to
$(\mathcal{F}_{t})_{t\geq 0}.$ $(\mathbb{P}_{x})_{x\in S}$ is a collection of
probability measures on $\Omega$ such that $\mathbb{P}_{x}(V_{0}=x)=1.$ For
any probability measure $\mu$ on $S$, we define

$\mathbb{P}_{\mu}(\cdot):=\int_{S}\mathbb{P}_{y}(\cdot)\mu(\mathrm{d}y)$

and denote expectations with respect to $\mathbb{P}_{x}$ and
$\mathbb{P}_{\mu}$ by $\mathbb{E}_{x}$ and $\mathbb{E}_{\mu}$, respectively.

Below we give a set of assumptions on the process $(V_{t})_{t\geq 0}$. We need
those to ensure that a certain coupling property holds. We comment on these
assumptions after stating them.

###### Assumption 1

Let $(V_{t})_{t\geq 0}$ be a Feller process on
$(\Omega,\mathcal{F},(\mathcal{F}_{t}))_{t\geq 0},(\mathbb{P}_{x})_{x\in S})$.
We assume the following:

  * (i)

$(V_{t})_{t\geq 0}$ admits a unique invariant measure $\pi$.

  * (ii)

For any $x\in S$, there exist a family $(V^{x}_{t})_{t\geq 0}$ and a
stationary version $(V^{\pi}_{t})_{t\geq 0}$ defined on the same probability
space $(\tilde{\Omega},\mathcal{\tilde{F}},\tilde{\mathbb{P}})$ such that,
$(V^{x}_{t})_{t\geq 0}\stackrel{{\scriptstyle d}}{{=}}(V_{t})_{t\geq 0}$ in
$\mathbb{P}_{x}$ and $(V^{\pi}_{t})_{t\geq 0}\stackrel{{\scriptstyle
d}}{{=}}(V_{t})_{t\geq 0}$ in $\mathbb{P}_{\pi}$, i.e. for any $0\leq
t_{1}<\cdots<t_{n}$,

$\tilde{\mathbb{P}}(V_{t_{1}}^{x}\in A_{1},\cdots,V_{t_{n}}^{x}\in
A_{n})=\mathbb{P}_{x}(V_{t_{1}}\in A_{1},\cdots,V_{t_{n}}\in A_{n}),$
$\tilde{\mathbb{P}}(V_{t_{1}}^{\pi}\in A_{1},\cdots,V_{t_{n}}^{\pi}\in
A_{n})=\mathbb{P}_{\pi}(V_{t_{1}}\in A_{1},\cdots,V_{t_{n}}\in A_{n}),$

where $A_{1},\cdots,A_{n}\in\mathcal{B}(S)$.

  * (iii)

Let $T^{x}:=\inf{\\{t\geq 0\ |\ V^{x}_{t}=V^{\pi}_{t}\\}}$ be a stopping time.
There exist constants $C,\delta>0$ such that for any $t\geq 0$,

$\sup_{x\in S}\tilde{\mathbb{P}}(T^{x}\geq t)\leq C\exp({-\delta t}).$

First, we assume that $(V_{t})_{t\geq 0}$ has a stationary measure $\pi$.
Second, we assume that for the process $(V_{t})_{t\geq 0}$ that starts from
$x$ with probability $1$, we can find a coupled process $(V^{x}_{t})_{t\geq
0}$. Also, given that the process $(V_{t})_{t\geq 0}$ starts with its
invariant measure $\pi$, we can find a stationary version
$(V^{\pi}_{t})_{t\geq 0}$ of $(V_{t})_{t\geq 0}$. Here, the processes
$(V^{x}_{t})_{t\geq 0}$ and $(V^{\pi}_{t})_{t\geq 0}$ are defined on the same
probability space. Third, we assume that the processes $(V^{x}_{t})_{t\geq 0}$
and $(V^{\pi}_{t})_{t\geq 0}$ intersect exponentially fast. The exponential
rate can be chosen uniformly in $x$ since $S$ is compact. With Assumption 1,
we have the following lemma.

######  Lemma 1 (Geometric Ergodicity)

Under Assumption 1, there exist constants $C,\delta>0$ such that for any $x\in
S$ and $t\geq 0$,

$\displaystyle\sup_{A\in\mathcal{B}(S)}|\mathbb{P}_{x}(V_{t}\in A)-\pi(A)|\leq
C\exp({-\delta t}),$

where $\mathcal{B}(S)$ is the set of all Borel measurable sets of $S$.

Proof  For any given $x\in S$, we construct the following process by coupling
$(V^{x}_{t})_{t\geq 0}$ and $(V^{\pi}_{t})_{t\geq 0}$:

$\displaystyle\tilde{V}^{x}_{t}=\left\\{\begin{aligned} &V^{x}_{t},\ \ \ 0\leq
t\leq T^{x},\\\ &V^{\pi}_{t},\ \ \ t>T^{x}.\end{aligned}\right.$

By the strong Markov property, $(\tilde{V}^{x}_{t})_{t\geq
0}\stackrel{{\scriptstyle d}}{{=}}(V^{x}_{t})_{t\geq 0}$. For any
$A\in\mathcal{B}(S)$, notice that

$\displaystyle\left|\mathbb{P}_{x}(V_{t}\in A)-\pi(A)\right|$ $\displaystyle=$
$\displaystyle|\tilde{\mathbb{P}}(V^{x}_{t}\in
A)-\tilde{\mathbb{P}}(V^{\pi}_{t}\in A)|$ $\displaystyle=$
$\displaystyle|\tilde{\mathbb{P}}(\tilde{V}^{x}_{t}\in
A)-\tilde{\mathbb{P}}(V^{\pi}_{t}\in A)|$ $\displaystyle=$
$\displaystyle|\tilde{\mathbb{P}}(\tilde{V}^{x}_{t}\in A,\tilde{V}^{x}_{t}\neq
V^{\pi}_{t})+\tilde{\mathbb{P}}(\tilde{V}^{x}_{t}\in
A,\tilde{V}^{x}_{t}=V^{\pi}_{t})$ $\displaystyle\ \ \
-(\tilde{\mathbb{P}}(V^{\pi}_{t}\in A,\tilde{V}^{x}_{t}\neq
V^{\pi}_{t})+\tilde{\mathbb{P}}(V^{\pi}_{t}\in
A,\tilde{V}^{x}_{t}=V^{\pi}_{t}))|$ $\displaystyle=$
$\displaystyle|\tilde{\mathbb{P}}(\tilde{V}^{x}_{t}\in A,\tilde{V}^{x}_{t}\neq
V^{\pi}_{t})-\tilde{\mathbb{P}}(V^{\pi}_{t}\in A,\tilde{V}^{x}_{t}\neq
V^{\pi}_{t})|$ $\displaystyle\leq$ $\displaystyle
2\tilde{\mathbb{P}}(\tilde{V}^{x}_{t}\neq V^{\pi}_{t})$ $\displaystyle\leq$
$\displaystyle 2\tilde{\mathbb{P}}(T^{x}\geq t)\leq C\exp({-\delta t}).$

From the third assumption in Assumption 1, $C$ and $\delta$ are independent of
$x$ and this completes the proof.

In the lemma above, we have shown geometric ergodicity of $(V_{t})_{t\geq 0}$
in the total variation distance. Next, we show that the same rate of
convergence of $(V_{t})_{t\geq 0}$ holds in the weak topology.

###### Corollary 2

Under Assumption 1, there exist constants $C,\delta>0$ such that for any
$h\in\mathcal{C}(S),$ i.e. the set of all continuous function on $S$, we have

$\displaystyle\sup_{x\in
S}\left|\mathbb{E}_{x}[h(V_{t})]-\int_{S}h(y)\pi(\mathrm{d}y)\right|\leq
C\left\|h\right\|_{\infty}\exp({-\delta t})$

where $\left\|h\right\|_{\infty}:=\sup_{x\in S}|h(x)|$

Proof  Rewrite $\mathbb{E}_{x}[h(V_{t})]$ as

$\displaystyle\mathbb{E}_{x}[h(V_{t})]=\int_{S}h(y)\mathbb{P}_{x}(V_{t}\in\mathrm{d}y).$

Then we have

$\displaystyle\left|\mathbb{E}_{x}[h(V_{t})]-\int_{S}h(y)\pi(\mathrm{d}y)\right|$
$\displaystyle=\left|\int_{S}h(y)[\mathbb{P}_{x}(V_{t}\in\mathrm{d}y)-\pi(\mathrm{d}y)]\right|$
$\displaystyle\leq\left\|h\right\|_{\infty}\left|\mathbb{P}_{x}(V_{t})-\pi\right|(S),$

where $\frac{1}{2}\left|\mathbb{P}_{x}(V_{t})-\pi\right|(S)$ is the total
variation of the measure $\mathbb{P}_{x}(V_{t})-\pi.$ Notice that

$\left|\mathbb{P}_{x}(V_{t})-\pi\right|(S)=2\sup_{A\in\mathcal{B}(S)}|\mathbb{P}_{x}(V_{t}\in
A)-\pi(A)|.$

By Lemma 1, we have

$\displaystyle\sup_{x\in
S}\left|\mathbb{E}_{x}[h(V_{t})]-\int_{S}h(y)\pi(\mathrm{d}y)\right|\leq$
$\displaystyle\left\|h\right\|_{\infty}\sup_{x\in
S}\left|\mathbb{P}_{x}(V_{t})-\pi\right|(S)$ $\displaystyle\leq$
$\displaystyle 2\left\|h\right\|_{\infty}\sup_{x\in
S}\sup_{A\in\mathcal{B}(S)}|\mathbb{P}_{x}(V_{t}\in A)-\pi(A)|$
$\displaystyle\leq$ $\displaystyle C\left\|h\right\|_{\infty}\exp({-\delta
t}),$

which completes the proof.

We now study four examples for processes that satisfy our assumptions: Lévy
processes with two-sided reflections on a compact interval, continuous-time
Markov processes on finite and countably infinite spaces, and processes on
rectangular sets with independent coordinates.

###  2.1 Example 1: Lévy processes with two-sided reflection

For any $b>0,$ we say a triplet $((X_{t})_{t\geq 0},(L_{t})_{t\geq
0},(U_{t})_{t\geq 0})$ is a solution to the Skorokhod problem of the Lévy
process $(X_{t})_{t\geq 0}$ on the space $S:=[0,b]$ if for all $t\geq 0,$

$\displaystyle V_{t}=X_{t}+L_{t}-U_{t},$ (6)

where $(L_{t})_{t\geq 0},\ (U_{t})_{t\geq 0}$ are non-decreasing right
continuous processes such that

$\displaystyle\int_{0}^{\infty}V_{t}\mathrm{d}L_{t}=\int_{0}^{\infty}(b-V_{t})\mathrm{d}U_{t}=0.$

In other words, $(L_{t})_{t\geq 0}$ and $(U_{t})_{t\geq 0}$ can only increase
when $(V_{t})_{t\geq 0}$ is at the lower boundary $0$ or the upper boundary
$b$. From Andersen et al. (2015, Proposition 5.1), we immediately have that
the process $(V_{t})_{t\geq 0}$ in (6) satisfies Assumption 1. The geometric
ergodicity follows from Andersen et al. (2015, Remark 5.3). As an example, the
standard Brownian Motion (BM) reflected at 0 and 1 can be written as

$V_{t}=B_{t}+\tilde{L}_{t}^{0}-\tilde{L}_{t}^{1},$

where $(B_{t})_{t\geq 0}$ is a standard BM and $(\tilde{L}^{a}_{t})_{t\geq 0}$
is the symmetric local time of $(V_{t})_{t\geq 0}$ at $a$. Intuitively, a
local time describes the time spent at a given point of a continuous
stochastic process. The formal definition of symmetric local time of
continuous semimartingales can be found, for example, in Revuz and Yor (2013,
Chapter VI). For the optimization problem (1) with $S=[0,1]$ and $\pi$ being
the uniform measure on $S$, the corresponding stochastic process in (4) can be
chosen to be this Brownian Motion with two-sided refection since its invariant
measure is the uniform measure on $[0,1]$. To see this, for $x\in[0,1]$, from
Andersen et al. (2015, Theorem 5.4),

$\displaystyle\pi([x,1])=\mathbb{P}(B_{\tau_{x}\land\tau_{x-1}}=x)=\mathbb{P}(\tau_{x}<\tau_{x-1})=1-x,$

where $\tau_{a}=\inf\\{t\geq 0|B_{t}=a\\}$.

###  2.2 Example 2: Continuous-time Markov processes

We consider a continuous-time Markov process $V_{t}$ on state space
$I=\\{1,2,\ldots,N\\}$ with transition rate matrix

$\mathbf{A}_{N}=\mathbf{\Lambda}_{N}-N\lambda\mathbf{I}_{N},$

where $\lambda>0$, $\mathbf{\Lambda}_{N}$ is a $N\times N$ matrix whose entry
is always $\lambda$, and $\mathbf{I}_{N}$ is the identity matrix. From Latz
(2021), we know that the transition probability is given by

$\displaystyle\mathbb{P}(V_{t+s}=i|V_{s}=j)=\frac{1-\exp(-\lambda
Nt)}{N}+\exp(-\lambda Nt)\mathbf{1}[{i=j}].$

The invariant measure $\pi$ of $V_{t}$ is the uniform measure on $I$, i.e.
$\pi(i)=1/N$ for $i\in\\{1,...,N\\}$. To see that $(V_{t})_{t\geq 0}$
satisfies the rest of Assumption 1, consider a stationary version
$(\hat{V}_{t})_{t\geq 0}$ that is independent of $(V_{t})_{t\geq 0}$. Let
$V_{0}=1$. We define $T=\inf\\{t\geq 0,V_{t}=\hat{V}_{t}\\}$. Moreover, for
$i,\ j\in\mathbb{N}$, we denote the $i$-th and $j$-th jump time of
$(V_{t})_{t\geq 0}$ and $(\hat{V}_{t})_{t\geq 0}$ by $T_{i}$ and
$\hat{T}_{j}$, respectively. Then we have

$\mathbb{P}(T=0)=\mathbb{P}(V_{0}=1,\hat{V}_{0}=1)=\mathbb{P}(V_{0}=1)\mathbb{P}(\hat{V}_{0}=1)=\frac{1}{N}.$

For any $i,\ j\in\mathbb{N}$, since $T_{i}$ and $\hat{T}_{j}$ are independent,
$\mathbb{P}(T_{i}=\hat{T}_{j})=0$. Let $Y_{t}=(V_{t},\hat{V}_{t})$ be a Markov
process on $I\times I$ with transition probability:

$\displaystyle\mathbb{P}(Y_{t+s}=(i,j)|Y_{s}=(i_{0},j_{0}))$
$\displaystyle=\frac{1-\exp(-2\lambda
Nt)}{2N}\Big{(}\mathbf{1}[{i=i_{0}}]+\mathbf{1}[{j=j_{0}}]\Big{)}$
$\displaystyle\qquad+\exp(-2\lambda Nt)\mathbf{1}[{(i,j)=(i_{0},j_{0})}].$

Thus, $T$ is the first time when $Y_{t}$ hits $\\{(i,i)|i=1,...,N\\}.$ Let the
n-th jump time of $Y_{t}$ be $\tau_{n}$, for $t>0$, we have

$\displaystyle\mathbb{P}(T\geq t)=$ $\displaystyle\sum_{n\geq
1}\mathbb{P}(T=\tau_{n},\tau_{n}\geq t)$ $\displaystyle=$
$\displaystyle\sum_{n\geq 1}\exp({-2(N-1)n\lambda
t})\frac{N-1}{N}\frac{1}{N-1}\Big{(}\frac{N-2}{N-1}\Big{)}^{n-1}$
$\displaystyle\leq$ $\displaystyle C\exp({-2(N-1)\lambda t}),$

where the second equality follows from

$\mathbb{P}(T=\tau_{n})=\frac{N-1}{N}\frac{1}{N-1}\Big{(}\frac{N-2}{N-1}\Big{)}^{n-1}$

since there are $2N-4$ states available for the next jump. Thus,
$(V_{t})_{t\geq 0}$ satisfies Assumption 1.

###  2.3 Example 3: Continuous-time Markov processes with countable states

We consider a continuous-time Markov process $(V_{t})_{t\geq 0}$ on state
space $\mathbb{N}_{0}:=\mathbb{N}\cup\\{0\\}$ with exponential jump times. At
time $t$, if $V_{t}\in\mathbb{N},$ it jumps to $0$ with probability $1$ at the
next jump time. Otherwise, if $V_{t}=0$, it jumps to $i$ with probability
$1/2^{i}.$ It is easy to verify that the invariant measure $\pi$ of $V_{t}$ is
$\pi(\\{i\\})=1/2^{i+1}$. One may consider $\mathbb{N}$ as one state and view
$V_{t}$ as a Markov process with two states.

To verify that $(V_{t})_{t\geq 0}$ satisfies the rest of Assumption 1,
similarly to the previous example, we consider a stationary version
$(\hat{V}_{t})_{t\geq 0}$ that is independent of $(V_{t})_{t\geq 0}$. Let
$V_{0}=0$ and $T=\inf\\{t\geq 0,V_{t}=\hat{V}_{t}=0\\}$. For $i,\
j\in\mathbb{N}_{0}$, we denote the $i$-th and $j$-th jump time of
$(V_{t})_{t\geq 0}$ and $(\hat{V}_{t})_{t\geq 0}$ by $T_{i}$ and
$\hat{T}_{j}$, respectively. Then we have

$\mathbb{P}(T=0)=\mathbb{P}(V_{0}=0,\hat{V}_{0}=0)=\mathbb{P}(V_{0}=0)\mathbb{P}(\hat{V}_{0}=0)=\frac{1}{2}.$

For any $i,\ j\in\mathbb{N}_{0}$, since $T_{i}$ and $\hat{T}_{j}$ are
independent, $\mathbb{P}(T_{i}=\hat{T}_{j})=0$. Let
$Y_{t}=(V_{t},\hat{V}_{t})$ be a Markov process on
$\mathbb{N}_{0}\times\mathbb{N}_{0}$. Notice that $T$ is the first time when
$Y_{t}$ hits $(0,0).$ Let the n-th jump time of $Y_{t}$ be $\tau_{n}$, for
$t>0$, we have

$\displaystyle\mathbb{P}(T\geq t)=$ $\displaystyle\sum_{n\geq
1}\mathbb{P}(T=\tau_{n},\tau_{n}\geq t)$ $\displaystyle=$
$\displaystyle\sum_{k\geq 0}\mathbb{P}(T=\tau_{2k+1},\tau_{2k+1}\geq t)$
$\displaystyle=$ $\displaystyle\sum_{k\geq
0}\exp({-2(2k+1)t})\frac{1}{2^{k+2}}\leq\exp({-t}),$

where the second and the third equality follows from
$\mathbb{P}(T=\tau_{2k+1})=2^{-k-1}$ and $\mathbb{P}(T=\tau_{2k})=0$ for
$k\geq 1$. Since $\inf\\{t\geq 0,V_{t}=\hat{V}_{t}\\}$ is upper bounded by
$T$, $(V_{t})_{t\geq 0}$ satisfies Assumption 1.

###  2.4 Example 4: Multidimensional processes

For multidimensional processes, Assumption 1 is satisfied if each component
satisfies Assumption 1 and all components are mutually independent. We
illustrate this by discussing the 2-dimensional case – higher-dimensional
processes can be constructed inductively. Multidimensional processes arise,
e.g., when the underlying space $S$ is multidimensional. They also arise when
the $S$ is one-dimensional, but we run multiple processes in parallel to
obtain a mini-batch SGD instead of single-draw SGD.

Let $(S^{1},m^{1})$ and $(S^{2},m^{2})$ be two compact Polish spaces. We
consider the probability triples $(\Omega^{1},(\mathcal{F}^{1}_{t}))_{t\geq
0},(\mathbb{P}^{1}_{a})_{a\in S^{1}})$ and
$(\Omega^{2},(\mathcal{F}^{2}_{t}))_{t\geq 0},(\mathbb{P}^{2}_{b})_{b\in
S^{2}})$ with
$\mathbb{P}^{1}_{a}((V^{1}_{0}=a)=\mathbb{P}^{2}_{b}(V^{2}_{0}=b)=1$. Let
$(V^{1}_{t})_{t\geq 0}$ and $(V^{2}_{t})_{t\geq 0}$ be
$(\mathcal{F}^{1}_{t}))_{t\geq 0}$ and $(\mathcal{F}^{2}_{t}))_{t\geq 0}$
adapted and from $\Omega^{1}$ to $S^{1}$ and $\Omega^{2}$ to $S^{2}$
respectively. In the following proposition, we construct a 2-dimensional
process $(V^{1}_{t},V^{2}_{t})_{t\geq 0}$ from $\Omega^{1}\times\Omega^{2}$ to
$(S^{1}\times S^{2},m^{1}+m^{2})$ with a family of probability measures
$(\mathbb{P}_{(a,b)})_{(a,b)\in S^{1}\times S^{2}}$ such that
$\mathbb{P}_{(a,b)}(A\times B)=\mathbb{P}^{1}_{a}(A)\mathbb{P}^{2}_{b}(B)$ for
$A\in\mathcal{F}^{1}$ and $B\in\mathcal{F}^{2}$.

We now show that the joint process $(V^{1}_{t},V^{2}_{t})_{t\geq 0}$ is Feller
and satisfies Assumption 1, if the marginals do.

###### Proposition 1

Let $(V^{1}_{t})_{t\geq 0}$ and $(V^{2}_{t})_{t\geq 0}$ be càdlàg and Feller
with respect to $(\mathcal{F}^{1}_{t})_{t\geq 0}$ and
$(\mathcal{F}^{2}_{t})_{t\geq 0},$ respectively, and satisfy Assumption 1 with
probability $(\mathbb{P}^{1}_{a})_{a\in S^{2}}$ and
$(\mathbb{P}^{2}_{b})_{b\in S^{2}}$, respectively. Then
$(V^{1}_{t},V^{2}_{t})_{t\geq 0}$ is also càdlàg and Feller with respect to
$\sigma(\mathcal{F}^{1}_{t}\times\mathcal{F}^{2}_{t})_{t\geq 0}$ and satisfies
Assumption 1 with $(\mathbb{P}_{(a,b)})_{(a,b)\in S^{1}\times S^{2}}$.

Proof  It is obvious that the process $(V^{1}_{t},V^{2}_{t})_{t\geq 0}$ is
càdlàg and Markovian. To verify the Feller property, we show that for any
continuous function $F$ on $S^{1}\times S^{2}$,
$\mathbb{E}_{(x,y)}[F(V^{1}_{t},V^{2}_{t})]$ is continuous in $(x,y)$. We
shall prove this by showing this property for separable $F$ and approximate
general continuous functions using this special case. Let $f$ and $g$ be
continuous functions on $S^{1}$ and $S^{2}$ respectively, then we have

$\displaystyle\mathbb{E}_{(x,y)}[f(V^{1}_{t})g(V^{2}_{t})]=\mathbb{E}^{1}_{x}[f(V^{1}_{t})]\mathbb{E}^{2}_{y}[g(V^{2}_{t})],$
(7)

which implies $\mathbb{E}_{(x,y)}[f(V^{1}_{t})g(V^{2}_{t})]$ is continuous in
$(x,y)$ since $(V_{t}^{1})_{t\geq 0}$ and $(V_{t}^{2})_{t\geq 0}$ are Feller.
By the Stone–Weierstrass theorem, for any $k\geq 1$, any continuous function
$F$ on $S^{1}\times S^{2}$ can be approximated as the following,

$\displaystyle\sup_{(x,y)\in S^{1}\times
S^{2}}\left|F(x,y)-\sum_{i=1}^{n_{k}}f_{i}^{k}(x)g_{i}^{k}(y)\right|\leq\frac{1}{k}$

where $f_{i}^{k}$ and $g_{i}^{k}$ are continuous. From (7), this implies
$\mathbb{E}_{(x,y)}[F(V^{1}_{t},V^{2}_{t})]$ is continuous on $S^{1}\times
S^{2}$.

Next, we prove that $(V^{1}_{t},V^{2}_{t})_{t\geq 0}$ satisfies Assumption 1.
Let $\pi^{1}$ and $\pi^{2}$ be the invariant measures of $(V^{1}_{t})_{t\geq
0}$ and $(V_{t}^{2})_{t\geq 0}$, respectively. Then $\pi^{1}\times\pi^{2}$ is
the invariant measure of $(V^{1}_{t},V^{2}_{t})_{t\geq 0}$ since
$(V^{1}_{t})_{t\geq 0}$ and $(V^{2}_{t})_{t\geq 0}$ are independent. From
Assumption 1, we know there exist
$(\tilde{\Omega}^{1},\mathcal{\tilde{F}}^{1},\tilde{\mathbb{P}}^{1}),$
$(\tilde{\Omega}^{2},\mathcal{\tilde{F}}^{2},\tilde{\mathbb{P}}^{2}),$ such
that for any $a\in S^{1}$ and $b\in S^{2}$ , $(V^{1,a}_{t})_{t\geq
0}\stackrel{{\scriptstyle d}}{{=}}(V^{1}_{t})_{t\geq 0}$ in
$\mathbb{P}^{1}_{a}$ and $(V^{2,b}_{t})_{t\geq 0}\stackrel{{\scriptstyle
d}}{{=}}(V^{2}_{t})_{t\geq 0}$ in $\mathbb{P}^{2}_{b}$. We define
$\tilde{\mathbb{P}}$ on $\tilde{\Omega}^{1}\times\tilde{\Omega}^{2}$ such that

$\tilde{\mathbb{P}}(A\times
B)=\tilde{\mathbb{P}}^{1}(A)\tilde{\mathbb{P}}^{2}(B),\ \
(A\in\mathcal{\tilde{F}}^{1},\ B\in\mathcal{\tilde{F}}^{2}).$

Then we have that $((V^{1,a}_{t})_{t\geq 0})_{a\in S^{1}}$ and
$(V^{1,\pi^{1}}_{t})_{t\geq 0}$ are independent of $((V^{2,b}_{t})_{t\geq
0})_{b\in S^{2}}$ and $(V^{2,\pi^{2}}_{t})_{t\geq 0}$ under
$\tilde{\mathbb{P}}.$ Similar to the proof of Lemma 1, we construct the
following processes by the coupling method:

$\displaystyle\tilde{V}^{1,a}_{t}=\left\\{\begin{aligned} &V^{1,a}_{t},\ \ \
0\leq t\leq T^{1,a},\\\ &V^{1,\pi^{1}}_{t},\ \ \
t>T^{1,a},\end{aligned}\right.$

and

$\displaystyle\tilde{V}^{2,b}_{t}=\left\\{\begin{aligned} &V^{2,b}_{t},\ \ \
0\leq t\leq T^{2,b},\\\ &V^{2,\pi^{2}}_{t},\ \ \
t>T^{2,b}.\end{aligned}\right.$

Then the distribution of $(\tilde{V}^{1,a}_{t},\tilde{V}^{2,b}_{t})_{t\geq 0}$
under $\tilde{\mathbb{P}}$ is the same as the distribution of
$(V^{1}_{t},V^{2}_{t})_{t\geq 0}$ under $\mathbb{P}_{(a,b)};$ the distribution
of $(V^{1,\pi^{1}}_{t},V^{2,\pi^{2}}_{t})_{t\geq 0}$ under
$\tilde{\mathbb{P}}$ is the same as the distribution of
$(V^{1}_{t},V^{2}_{t})_{t\geq 0}$ under $\mathbb{P}_{\pi^{1}\times\pi^{2}}.$
Moreover, $(\tilde{V}^{1,a}_{t},\tilde{V}^{2,b}_{t})_{t\geq 0}$ intersects the
invariant state $(V^{1,\pi^{1}}_{t},V^{2,\pi^{2}}_{t})_{t\geq 0}$ at time
$T^{1,a}\lor T^{1,b}.$ For any $(a,b)\in S^{1}\times S^{2}$,

$\tilde{\mathbb{P}}(T^{1,a}\lor T^{1,b}\geq
t)\leq\tilde{\mathbb{P}}(T^{1,a}\geq t)+\tilde{\mathbb{P}}(T^{1,b}\geq t)\leq
C\exp({-\delta t}).$

We have now discussed various properties of potential index processes and move
on to study the stochastic gradient process.

##  3 Stochastic gradient processes with constant learning rate

We now define and study the stochastic gradient process with constant learning
rate. Here, the switching between data sets is performed in a homogeneous-in-
time way. Hence, it models the discrete-time stochastic gradient descent
algorithm when employed with a constant learning rate. Although, one can
usually not hope to converge to the minimizer of the target functional in this
case, this setting is popular in practice.

To obtain the stochastic gradient process with constant learning rate, we will
couple the gradient flow (4) with the an appropriate process
$(V_{t/\varepsilon})_{t\geq 0}$. Here, $(V_{t})_{t\geq 0}$ is a Feller process
introduced in Section 2 and $\varepsilon>0$ is a scaling parameter that allows
us to uniformly control a switching rate parameter. To define the stochastic
process associated with this stochastic gradient descent problem, we first
introduce the following assumptions that guarantee the existence and
uniqueness of the solution of the associated stochastic differential equation.
After its formal definition and the proof of well-definedness, we move on to
the analysis of the process. Indeed, we show that the process approximates the
full gradient flow (9), as $\varepsilon\downarrow 0$. Moreover, we show that
the process has a unique stationary measure to which it converges in the
longtime limit at geometric speed.

We commence with regularity properties of the subsampled target function $f$
that are necessary to show the well-definedness of the stochastic gradient
process.

###### Assumption 2

Let $f(\theta,y)\in\mathcal{C}^{2}(\mathbb{R}^{K}\times S,\mathbb{R})$.

1\. $\nabla_{\theta}f,$ $H_{\theta}f$ are continuous.

2\. $\nabla_{\theta}f(\theta,y)$ is Lipschitz in $x$ and the Lipschitz
constant is uniform for $y\in S.$

3\. For $\theta\in\mathbb{R}^{K}$, $f(\theta,\cdot)$ and $\nabla_{\theta}f$
are integrable w.r.t to the probability measure $\pi(\cdot)$.

Now, we move on to the formal definition of the stochastic gradient process.

###### Definition 1

For $\varepsilon>0$, the _stochastic gradient process with_ _constant learning
rate_ _(SGPC)_ is a solution of the following stochastic differential
equation,

$\left\\{\begin{array}[]{l}\mathrm{d}\theta^{\varepsilon}_{t}=-\nabla_{\theta}f(\theta^{\varepsilon}_{t},V_{t/\varepsilon})\mathrm{d}t,\\\
\theta_{0}^{\varepsilon}=\theta_{0},\end{array}\right.$ (8)

where f satisfies Assumption 2 and $(V_{t})_{t\geq 0}$ is a Feller process
that satisfies Assumption 1.

Given these two assumptions, we can indeed show that the SGPC is a well-
defined Markov process. Moreover, we show that the stochastic gradient process
is Markovian, a property it shares with the discrete-time stochastic gradient
descent method.

###### Proposition 2

Let Assumptions 1 and 2 hold. Then, equation (8) has a unique strong solution,
i.e. the solution $(\theta^{\varepsilon}_{t})_{t\geq 0}$ is measurable with
respect to $\mathcal{F}^{\varepsilon}_{t}:=\mathcal{F}_{t/\varepsilon}$ for
any $t\geq 0$. For $y\in S$,
$(\theta_{t}^{\varepsilon},V_{t/\varepsilon})_{t\geq 0}$ is a Markov process
under $\mathbb{P}_{y}$ with respect to $(\mathcal{F}^{\varepsilon}_{t})_{t\geq
0}$.

Proof  The existence and the uniqueness of the strong solution to the equation
(8) can be found in Kushner (1990, Chapter 2, Theorem 4.1). To prove the
Markov property, we define the operator $(Q^{\varepsilon}_{t})_{t\geq 0}$ such
that

$Q^{\varepsilon}_{t}h(x,y):=\mathbb{E}_{y}[h(\theta^{\varepsilon}_{t},V_{t/\varepsilon})|\theta^{\varepsilon}_{0}=x],$

for any function $h$ bounded and measurable on $\mathbb{R}^{K}\times S$. For
any $s,t\geq 0$, we want to show

$\displaystyle\mathbb{E}[h(\theta^{\varepsilon}_{t+s},V_{(t+s)/\varepsilon})|\mathcal{F}^{\varepsilon}_{s}]=Q^{\varepsilon}_{t}h(\theta^{\varepsilon}_{s},V_{s/\varepsilon}).$

We set $\hat{\theta}^{\varepsilon}_{t}:=\theta^{\varepsilon}_{t+s},\
\mathcal{\hat{F}}_{t}:=\mathcal{F}^{\varepsilon}_{t+s},\
\hat{V}_{t/\varepsilon}:=V_{(t+s)/\varepsilon}.$ Since

$\displaystyle\theta^{\varepsilon}_{t+s}=\theta^{\varepsilon}_{s}-\int_{s}^{t+s}\nabla_{\theta}f(\theta^{\varepsilon}_{m},V_{m/\varepsilon})\mathrm{d}m,$

we have

$\displaystyle\hat{\theta}^{\varepsilon}_{t}=\hat{\theta}^{\varepsilon}_{0}-\int_{0}^{t}\nabla_{\theta}f(\hat{\theta}^{\varepsilon}_{m},\hat{V}_{m/\varepsilon})\mathrm{d}m.$

Hence $\hat{\theta}^{\varepsilon}_{t}$ is the solution of equation (8) with
$\hat{\theta}^{\varepsilon}_{0}=\theta^{\varepsilon}_{s}$ and
$\hat{V}_{0}=V_{s/\varepsilon}.$ Moreover,

$\displaystyle\mathbb{E}[h(\theta^{\varepsilon}_{t+s},V_{(t+s)/\varepsilon})|\mathcal{F}^{\varepsilon}_{s}]$
$\displaystyle=\mathbb{E}[h(\hat{\theta}^{\varepsilon}_{t},\hat{V}_{t/\varepsilon})|\hat{\theta}^{\varepsilon}_{0}=\theta^{\varepsilon}_{s},\hat{V}_{0}=V_{s/\varepsilon}]$
$\displaystyle=\mathbb{E}_{\hat{V}_{0}}[h(\hat{\theta}^{\varepsilon}_{t},\hat{V}_{t/\varepsilon})|\hat{\theta}^{\varepsilon}_{0}=\theta^{\varepsilon}_{s}]$
$\displaystyle=Q^{\varepsilon}_{t}h(\theta^{\varepsilon}_{s},V_{s/\varepsilon}),$

where the second equality and third equality follow from the homogeneous
Markov property of $(V^{\varepsilon}_{t})_{t\geq 0}$.

###  3.1 Approximation of the full gradient flow

We now let $\varepsilon\to 0$ and study the limiting behavior of SGPC. Indeed,
we aim to show that here the SGPC converges to the _full gradient flow_

$\mathrm{d}\zeta_{t}=-\Big{[}\int_{S}\nabla_{\zeta}f(\zeta_{t},v)\pi(\mathrm{d}v)\Big{]}\mathrm{d}t.$
(9)

We study this topic for two reasons: First, we aim to understand the
interdependence of $(V_{t})_{t\geq 0}$ and $(\theta_{t}^{\varepsilon})_{t\geq
0}$. Second, we understand SGPC as an approximation to the full gradient flow
(9), as motivated in the introduction. Hence, we should show that SGPC can
approximate the full gradient flow at any accuracy.

We now denote
$g(\cdot):=\int_{S}\nabla_{\zeta}f(\cdot,v)\pi(\mathrm{d}v)\in\mathcal{C}^{1}(\mathbb{R}^{K},\mathbb{R}^{K}).$
Then, we can define $(\zeta_{t})_{t\geq 0}$ through the dynamical system
$\mathrm{d}\zeta_{t}=-g(\zeta_{t})\mathrm{d}t$. Moreover, let
$\mathcal{C}([0,\infty):\mathbb{R}^{K})$ be the space of continuous functions
from $[0,\infty)$ to $\mathbb{R}^{K}$ equipped with the distance

$\rho\Big{(}(\varphi_{t})_{t\geq 0},(\varphi^{\prime}_{t})_{t\geq
0}\Big{)}:=\int_{0}^{\infty}\exp({-t})(1\land\sup_{0\leq s\leq
t}\left\|\varphi_{s}-\varphi_{s}^{\prime}\right\|)\mathrm{d}t,$

where $(\varphi_{t})_{t\geq 0},(\varphi^{\prime}_{t})_{t\geq
0}\in\mathcal{C}([0,\infty):\mathbb{R}^{K})$. We study the weak limit of the
system (8) as $\varepsilon\to 0$. Similar problems have been discussed in, for
example, Kushner (1990) and Kushner (1984).

###### Theorem 3

Let $\theta^{\varepsilon}_{0}=\theta_{0}$ and $\zeta_{0}=\theta_{0}$.
Moreover, let $(\theta^{\varepsilon}_{t})_{t\geq 0}$ and $(\zeta_{t})_{t\geq
0}$ solve (8) and (9), respectively. Then $(\theta^{\varepsilon}_{t})_{t\geq
0}$ under $\mathbb{P}_{\pi}$ converges weakly to $(\zeta_{t})_{t\geq 0}$ in
$\mathcal{C}([0,\infty):\mathbb{R}^{K})$ as $\varepsilon\to 0$, i.e. for any
bounded continuous function $F$ on $\mathcal{C}([0,\infty):\mathbb{R}^{K})$,

$\mathbb{E}_{\pi}[F\big{(}(\theta^{\varepsilon}_{t})_{t\geq
0}\big{)}]\to\mathbb{E}_{\pi}[F\big{(}(\zeta_{t})_{t\geq
0}\big{)}]=F\big{(}(\zeta_{t})_{t\geq 0}\big{)}.$

Proof  We shall verify that $(\theta^{\varepsilon}_{t})_{t\geq 0}$ is tight by
checking:

$\displaystyle 1.$ $\displaystyle\
\sup_{0<\varepsilon<1}\left\|\theta^{\varepsilon}_{0}\right\|<+\infty;$
$\displaystyle 2.$ $\displaystyle\ \textit{For any fixed }T>0,\
\lim_{\delta\to
0}\sup_{0<\varepsilon<1}\sup_{s,t\in[0,T],|s-t|\leq\delta}\left\|\theta^{\varepsilon}_{t}-\theta^{\varepsilon}_{s}\right\|\to
0.$

The first condition follows from $\theta^{\varepsilon}_{0}=\theta_{0}$. For
the second condition, by Assumption 2, let $C_{0}=\sup_{y\in
S}\left\|\nabla_{\theta}f(0,y)\right\|$ and $L_{f}$ be the Lipschitz constant
of $\nabla_{\theta}f(\cdot,y)$, we have

$\displaystyle\frac{\mathrm{d}\left\|\theta^{\varepsilon}_{t}\right\|^{2}}{\mathrm{d}t}$
$\displaystyle=-2\left\langle\theta^{\varepsilon}_{t},\nabla_{\theta}f(\theta^{\varepsilon}_{t},V_{t/\varepsilon})\right\rangle$
$\displaystyle=-2\left\langle\theta^{\varepsilon}_{t},\nabla_{\theta}f(\theta^{\varepsilon}_{t},V_{t/\varepsilon})-\nabla_{\theta}f(0,V_{t/\varepsilon})\right\rangle-2\left\langle\theta^{\varepsilon}_{t},\nabla_{\theta}f(0,V_{t/\varepsilon})\right\rangle$
$\displaystyle\leq
2L_{f}\left\|\theta^{\varepsilon}_{t}\right\|^{2}+2C_{0}\left\|\theta^{\varepsilon}_{t}\right\|$
$\displaystyle\leq
2L_{f}\left\|\theta^{\varepsilon}_{t}\right\|^{2}+\left\|\theta^{\varepsilon}_{t}\right\|^{2}+C^{2}_{0}$
$\displaystyle=(2L_{f}+1)\left\|\theta^{\varepsilon}_{t}\right\|^{2}+C^{2}_{0}.$

By Grönwall’s inequality,

$\displaystyle\left\|\theta^{\varepsilon}_{t}\right\|^{2}\leq(\left\|\theta_{0}\right\|^{2}+C^{2}_{0})e^{(2L_{f}+1)t}.$
(10)

Therefore, $\theta^{\varepsilon}_{t}$ is bounded on any finite time interval.
For any fixed $T>0,$ let

$C_{T,f,\theta_{0}}=\sup_{\left\|x\right\|\leq(\left\|\theta_{0}\right\|^{2}+C^{2}_{0})e^{(2L_{f}+1)T},\
y\in S}\left\|\nabla_{\theta}f(x,y)\right\|.$

Then for any $s,t\in[0,T]$,

$\displaystyle\left\|\theta^{\varepsilon}_{t}-\theta^{\varepsilon}_{s}\right\|\leq\int_{s}^{t}\left\|\nabla_{\theta}f(\theta^{\varepsilon}_{m},V_{m/\varepsilon})\right\|\mathrm{d}m\leq
C_{T,f,\theta_{0}}|t-s|.$

Hence, $(\theta^{\varepsilon}_{t})_{t\geq 0}$ is tight in
$\mathcal{C}([0,\infty):\mathbb{R}^{K})$. By Prokhorov’s theorem, let
$(\theta_{t})_{t\geq 0}$ be a weak limit of $(\theta^{\varepsilon}_{t})_{t\geq
0}$. We shall verify that $(\theta_{t})_{t\geq 0}$ satisfies equation (9),
which is equivalent to show that for any bounded differentiable function
$\varphi$, $h$

$\displaystyle\mathbb{E}_{\pi}\Big{[}\Big{(}\varphi(\theta_{t})-\varphi(\theta_{s})+\int_{s}^{t}\left\langle\nabla_{\theta}\varphi(\theta_{m}),g(\theta_{m})\right\rangle\mathrm{d}m\Big{)}h\Big{(}(\theta_{t_{i}})_{i=1,...,n}\Big{)}\Big{]}=0,$

$\forall\ 0\leq t_{1}<\cdots<t_{n}\leq s.$ The case $t=0$ is obvious. Since
$(\theta^{\varepsilon}_{t})_{t\geq 0}$ is a strong solution to equation (8),
for any $0\leq s<t$,

$\displaystyle\varphi(\theta^{\varepsilon}_{t})=\varphi(\theta^{\varepsilon}_{s})-\int_{s}^{t}\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{m}),\nabla_{\theta}f(\theta^{\varepsilon}_{m},V_{m/\varepsilon})\right\rangle\mathrm{d}m.$
(11)

Hence, we have

$\displaystyle\mathbb{E}_{\pi}\Big{[}\Big{(}\varphi(\theta^{\varepsilon}_{t})-\varphi(\theta^{\varepsilon}_{s})+\int_{s}^{t}\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{s}),\nabla_{\theta}f(\theta^{\varepsilon}_{m},V_{m/\varepsilon})\right\rangle\mathrm{d}m\Big{)}h\Big{(}(\theta^{\varepsilon}_{t_{i}})_{i=1,...,n}\Big{)}\Big{]}=0,$

Moreover, when $\varepsilon\to 0,$

$\displaystyle\mathbb{E}_{\pi}\Big{[}\Big{(}\varphi(\theta^{\varepsilon}_{t})-\varphi(\theta^{\varepsilon}_{s})\Big{)}h\Big{(}(\theta^{\varepsilon}_{t_{i}})_{i=1,...,n}\Big{)}\Big{]}\to\mathbb{E}_{\pi}\Big{[}\Big{(}\varphi(\theta_{t})-\varphi(\theta_{s})\Big{)}h\Big{(}(\theta_{t_{i}})_{i=1,...,n}\Big{)}\Big{]}.$

Hence, all we need to show is the following

$\displaystyle\mathbb{E}_{\pi}\Big{[}\Big{(}\int_{s}^{t}\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{m}),\nabla_{\theta}f(\theta^{\varepsilon}_{m},V_{m/\varepsilon})\right\rangle\mathrm{d}m-\int_{s}^{t}\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{m}),g(\theta^{\varepsilon}_{m})\right\rangle\mathrm{d}m\Big{)}h\Big{(}(\theta^{\varepsilon}_{t_{i}})_{i=1,...,n}\Big{)}\Big{]}\to
0,$ (12)

which is equivalent to prove that

$\displaystyle\mathbb{E}_{\pi}\Big{[}\int_{s}^{t}\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{m}),\nabla_{\theta}f(\theta^{\varepsilon}_{m},V_{m/\varepsilon})\right\rangle\mathrm{d}m-\int_{s}^{t}\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{m}),g(\theta^{\varepsilon}_{m})\right\rangle\mathrm{d}m\Big{|}\mathcal{F}^{\varepsilon}_{s}\Big{]}\to
0.$ (13)

Let $\tilde{\varepsilon}:=1/[1/\sqrt{\varepsilon}]$, where $[x]$ is the
greatest integer less than or equal to $x$. Then we have the following
decomposition

$\displaystyle\mathbb{E}_{\pi}\Big{[}\int_{s}^{t}\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{m}),\nabla_{\theta}f(\theta^{\varepsilon}_{m},V_{m/\varepsilon})\right\rangle\mathrm{d}m-\int_{s}^{t}\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{m}),g(\theta^{\varepsilon}_{m})\right\rangle\mathrm{d}m\Big{|}\mathcal{F}^{\varepsilon}_{s}\Big{]}$
$\displaystyle=$
$\displaystyle\tilde{\varepsilon}\sum^{1/\tilde{\varepsilon}}_{i=0}\tilde{\varepsilon}^{-1}\mathbb{E}_{\pi}\Big{[}\int_{s+i(t-s)\tilde{\varepsilon}}^{s+(i+1)(t-s)\tilde{\varepsilon}}\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{m}),\nabla_{\theta}f(\theta^{\varepsilon}_{m},V_{m/\varepsilon})-g(\theta^{\varepsilon}_{m})\right\rangle\mathrm{d}m\Big{|}\mathcal{F}^{\varepsilon}_{s}\Big{]}$
$\displaystyle=$
$\displaystyle\tilde{\varepsilon}\sum^{1/\tilde{\varepsilon}}_{i=0}\mathbb{E}_{\pi}\Big{[}\tilde{\varepsilon}^{-1}\mathbb{E}_{\pi}\Big{[}\int_{s+i(t-s)\tilde{\varepsilon}}^{s+(i+1)(t-s)\tilde{\varepsilon}}\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{m}),\nabla_{\theta}f(\theta^{\varepsilon}_{m},V_{m/\varepsilon})-g(\theta^{\varepsilon}_{m})\right\rangle\mathrm{d}m\Big{|}\mathcal{F}^{\varepsilon}_{s+i(t-s)\tilde{\varepsilon}}\Big{]}\Big{|}\mathcal{F}^{\varepsilon}_{s}\Big{]}.$

We claim that as $\varepsilon\to 0$,

$\displaystyle\sup_{0\leq
r<t}\tilde{\varepsilon}^{-1}\mathbb{E}_{\pi}\Big{[}\int_{r}^{r+(t-s)\tilde{\varepsilon}}G(\theta^{\varepsilon}_{m},V_{m/\varepsilon})\mathrm{d}m\Big{|}\mathcal{F}^{\varepsilon}_{r}\Big{]}\to
0,$ (14)

where
$G(x,y):=\left\langle\nabla_{\theta}\varphi(x),\nabla_{\theta}f(x,y)-g(x)\right\rangle$.
Notice that for any fixed $t>0$, $(\theta^{\varepsilon}_{s})_{0\leq s\leq t}$
is uniformly equicontinuous. Hence, we have

$\displaystyle\sup_{0\leq r\leq t}\sup_{r\leq m\leq
r+(t-s)\tilde{\varepsilon}}\left\|\theta^{\varepsilon}_{m}-\theta^{\varepsilon}_{r}\right\|$
$\displaystyle=\sup_{0\leq r\leq t}\sup_{r\leq m\leq
r+(t-s)\tilde{\varepsilon}}\int_{r}^{r+(t-s)\tilde{\varepsilon}}\left\|\nabla_{\theta}f(\theta^{\varepsilon}_{m},V_{m/\varepsilon})\right\|\mathrm{d}m$
$\displaystyle\leq\tilde{\varepsilon}\sup_{0\leq m\leq
t}\left\|\nabla_{\theta}f(\theta^{\varepsilon}_{m},V_{m/\varepsilon})\right\|.$

Therefore, as $\varepsilon\to 0,$

$\displaystyle\sup_{0\leq
r<t}\left|\tilde{\varepsilon}^{-1}\mathbb{E}_{\pi}\Big{[}\int_{r}^{r+(t-s)\tilde{\varepsilon}}G(\theta^{\varepsilon}_{m},V_{m/\varepsilon})\mathrm{d}m\Big{|}\mathcal{F}^{\varepsilon}_{r}\Big{]}-\tilde{\varepsilon}^{-1}\mathbb{E}_{\pi}\Big{[}\int_{r}^{r+(t-s)\tilde{\varepsilon}}G(\theta^{\varepsilon}_{r},V_{m/\varepsilon})\mathrm{d}m\Big{|}\mathcal{F}^{\varepsilon}_{r}\Big{]}\right|\to
0.$

Hence, (14) is equivalent to

$\displaystyle\sup_{0\leq
r<t}\left|\tilde{\varepsilon}^{-1}\mathbb{E}_{\pi}\Big{[}\int_{r}^{r+(t-s)\tilde{\varepsilon}}G(\theta^{\varepsilon}_{r},V_{m/\varepsilon})\mathrm{d}m\Big{|}\mathcal{F}^{\varepsilon}_{r}\Big{]}\right|\to
0.$ (15)

By Corollary 2,

$\displaystyle\sup_{0\leq
r<t}\left|\tilde{\varepsilon}^{-1}\mathbb{E}_{\pi}\Big{[}\int_{r}^{r+(t-s)\tilde{\varepsilon}}G(\theta^{\varepsilon}_{r},V_{m/\varepsilon})\mathrm{d}m\Big{|}\mathcal{F}^{\varepsilon}_{r}\Big{]}\right|$
$\displaystyle=$ $\displaystyle\sup_{0\leq
r<t}\left|\tilde{\varepsilon}^{-1}\mathbb{E}_{V_{r/\varepsilon},x=\theta^{\varepsilon}_{r}}\Big{[}\int_{r}^{r+(t-s)\tilde{\varepsilon}}G(x,V_{(m-r)/\varepsilon})\mathrm{d}m\Big{]}\right|$
$\displaystyle=$ $\displaystyle\sup_{0\leq
r<t}\left|\tilde{\varepsilon}^{-1}\mathbb{E}_{V_{r/\varepsilon},x=\theta^{\varepsilon}_{r}}\Big{[}\int_{r}^{r+(t-s)\tilde{\varepsilon}}\left\langle\nabla_{\theta}\varphi(x),\nabla_{\theta}f(x,V_{(m-r)/\varepsilon})-g(x)\right\rangle\mathrm{d}m\Big{]}\right|$
$\displaystyle=$ $\displaystyle\sup_{0\leq
r<t}\left|\tilde{\varepsilon}^{-1}\int_{r}^{r+(t-s)\tilde{\varepsilon}}\mathbb{E}_{V_{r/\varepsilon},x=\theta^{\varepsilon}_{r}}\Big{[}\left\langle\nabla_{\theta}\varphi(x),\nabla_{\theta}f(x,V_{(m-r)/\varepsilon})-g(x)\right\rangle\Big{]}\mathrm{d}m\right|$
$\displaystyle\leq$ $\displaystyle\sup_{0\leq
r<t}\tilde{\varepsilon}^{-1}\int_{r}^{r+(t-s)\tilde{\varepsilon}}\left|\mathbb{E}_{V_{r/\varepsilon},x=\theta^{\varepsilon}_{r}}\Big{[}\left\langle\nabla_{\theta}\varphi(x),\nabla_{\theta}f(x,V_{(m-r)/\varepsilon})-g(x)\right\rangle\Big{]}\right|\mathrm{d}m$
$\displaystyle\leq$ $\displaystyle\sup_{0\leq
r<t}\tilde{\varepsilon}^{-1}\left\|\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{r}),\nabla_{\theta}f(\theta^{\varepsilon}_{r},\cdot)\right\rangle\right\|_{\infty}\int_{r}^{r+(t-s)\tilde{\varepsilon}}e^{-\delta(m-r)/\varepsilon}dm$
$\displaystyle=$ $\displaystyle\sup_{0\leq
r<t}\tilde{\varepsilon}^{-1}\left\|\left\langle\nabla_{\theta}\varphi(\theta^{\varepsilon}_{r}),\nabla_{\theta}f(\theta^{\varepsilon}_{r},\cdot)\right\rangle\right\|_{\infty}\int_{0}^{(t-s)\tilde{\varepsilon}}e^{-\delta
k/\varepsilon}dk$ $\displaystyle\leq$ $\displaystyle
C_{t,\varphi,f}\frac{\varepsilon}{\delta\tilde{\varepsilon}}\leq
C_{t,\varphi,f}\frac{\sqrt{\varepsilon}}{2\delta}\to 0.$

This completes the proof of (13). Hence, any weak limit of
$(\theta^{\varepsilon}_{t})_{t\geq 0}$ is a martingale solution to equation
(9). Since equation (9) is a deterministic ordinary differential equation and
$\theta^{\varepsilon}_{0}=\theta_{0}$ is independent of $\varepsilon$, we have
$(\theta^{\varepsilon}_{t})_{t\geq 0}$ converges weakly to $(\zeta_{t})_{t\geq
0}$ as $\varepsilon\to 0$.

Instead of looking at the full trajectories of the processes, we can also
study their distributions and show convergence in the Wasserstein distance. We
first need to introduce some notation.

Let $\nu$ and $\nu^{\prime}$ be two probability measures on
$(\mathbb{R}^{K},\mathcal{B}(\mathbb{R}^{K}))$. We define the Wasserstein
distance between those measures by

${{\mathcal{W}}}_{d}(\nu,\nu^{\prime})=\inf_{\Gamma\in\mathcal{H}(\nu,\nu^{\prime})}\int_{\mathbb{R}^{K}\times\mathbb{R}^{K}}d(y,y^{\prime})\Gamma(\mathrm{d}y,\mathrm{d}y^{\prime}),$

where $d(y,y^{\prime}):=1\land\left\|y-y^{\prime}\right\|$ and
$\mathcal{H}(\nu,\nu^{\prime})$ is the set of coupling between $\nu$ and
$\nu^{\prime}$, i.e.

$\mathcal{H}(\nu,\nu^{\prime})=\\{\Gamma\in\text{Pr}(\mathbb{R}^{K}\times\mathbb{R}^{K}):\Gamma(A\times\mathbb{R}^{K})=\nu(A),\Gamma(\mathbb{R}^{K}\times
B)=\nu^{\prime}(B),\forall A,B\in\mathcal{B}(\mathbb{R}^{K})\\}.$

To simplify the notation, for $B\in\mathcal{B}(\mathbb{R}^{K})$,
$\theta\in\mathbb{R}^{K}$, and $y\in S$, we denote

$C^{\varepsilon}_{t}(B|\theta,y):=\mathbb{P}_{y}(\theta^{\varepsilon}_{t}\in
B|\theta^{\varepsilon}_{0}=\theta),$
$C^{\varepsilon}_{t}(B|\theta,\pi):=\mathbb{P}_{\pi}(\theta^{\varepsilon}_{t}\in
B|\theta^{\varepsilon}_{0}=\theta),$

where $\pi$ is the invariant measure of $(V_{t})_{t\geq 0}$.

Now we study the approximation property of SGPC in the Wasserstein distance.
Indeed, the following corollary follows immediately from Theorem 3.

###### Corollary 4

There exists a function $\alpha:(0,1)\to[0,1],$ such that

${{\mathcal{W}}}_{d}(C^{\varepsilon}_{t}(\cdot|\theta_{0},\pi),\delta(\cdot-\zeta_{t}))\leq(\exp(t)\alpha(\varepsilon))\land
1.$

Moreover, $\lim_{\varepsilon\to 0}\alpha(\varepsilon)=0$.

Proof  By Theorem 3, we have $(\theta^{\varepsilon}_{t})_{t\geq
0}\Rightarrow(\zeta_{t})_{t\geq 0}$. By Skorokhod’s representation theorem,
there exists a sequence $(\tilde{\theta}^{\varepsilon}_{t})_{t\geq 0}$ such
that

$\displaystyle(\tilde{\theta}^{\varepsilon}_{t})_{t\geq
0}\overset{d}{=}(\theta^{\varepsilon}_{t})_{t\geq 0}\ \text{under}\
\mathbb{P}_{\pi},$
$\displaystyle\rho\Big{(}(\tilde{\theta}^{\varepsilon}_{t}-\zeta_{t})_{t\geq
0},0\Big{)}\to 0\ \text{almost\ surely\ in}\ \mathbb{P}_{\pi}.$

This implies

$\displaystyle\mathbb{E}_{\pi}[F((\tilde{\theta}^{\varepsilon}_{t}-\zeta_{t})_{t\geq
0})]\to F(0),$

for any bounded continuous function $F$ on
$\mathcal{C}([0,\infty):\mathbb{R}^{K})$. By taking

$F((\tilde{\theta}^{\varepsilon}_{t}-\zeta_{t})_{t\geq 0}))=\sup_{t\geq
0}\exp({-t})\left(1\land\sup_{0\leq s\leq
t}\left\|\tilde{\theta}^{\varepsilon}_{t}-\zeta_{t}\right\|\right)$

and

$\displaystyle\alpha(\varepsilon):=\mathbb{E}_{\pi}\Big{[}\sup_{t\geq
0}\exp({-t})\left(1\land\sup_{0\leq s\leq
t}\left\|\tilde{\theta}^{\varepsilon}_{t}-\zeta_{t}\right\|\right)\Big{]}\to
0,$

we have, for all $t\geq 0$,

$\displaystyle\mathbb{E}_{\pi}\Big{[}1\land\left\|\tilde{\theta}^{\varepsilon}_{t}-\zeta_{t}\right\|\Big{]}\leq\exp(t)\alpha(\varepsilon).$

Since $1\land\left\|\tilde{\theta}^{\varepsilon}_{t}-\zeta_{t}\right\|\leq 1$
and $\tilde{\theta}^{\varepsilon}_{t}\overset{d}{=}\theta^{\varepsilon}_{t}$,
denoting the distribution of $\tilde{\theta}^{\varepsilon}_{t}$ as
$F_{\tilde{\theta}^{\varepsilon}_{t}}$,

$\displaystyle{{\mathcal{W}}}_{d}(C^{\varepsilon}_{t}(\cdot|\theta_{0},\pi),\delta(\cdot-\zeta_{t}))\leq$
$\displaystyle{{\mathcal{W}}}_{d}(F_{\tilde{\theta}^{\varepsilon}_{t}}\
,C^{\varepsilon}_{t}(\cdot|\theta_{0},\pi))+\mathbb{E}_{\pi}\Big{[}1\land\left\|\tilde{\theta}^{\varepsilon}_{t}-\zeta_{t}\right\|\Big{]}$
$\displaystyle\leq$ $\displaystyle(\exp(t)\alpha(\varepsilon))\land 1.$

Finally in this section, we look at a technical result concerning the
asymptotic behavior of the full gradient flow $(\zeta_{t})_{t\geq 0}.$ First,
we will additionally assume that the subsampled target function $f(\cdot,y)$
in the optimization problem is strongly convex, with a convexity parameter
that does not depend on $y\in S$. We state this assumption below.

######  Assumption 3 (Strong Convexity)

For any $x_{1},x_{2}\in\mathbb{R}^{K},$

$\left\langle
x_{1}-x_{2},\nabla_{\theta}f(x_{1},y)-\nabla_{\theta}f(x_{2},y)\right\rangle\geq\kappa\left\|x_{1}-x_{2}\right\|^{2}$

where $\kappa>0$ and $\kappa$ is independent of $y\in S$.

Strong convexity implies, of course, that the full target function
$g:=\int_{S}\nabla_{\theta}f(\cdot,y)\pi(\mathrm{d}y)$ has a unique minimizer
$\theta^{*}$. It also implies that the associated full gradient flow
$(\zeta_{t})_{t\geq 0}$ converges at exponential speed to this unique
minimizer. We give a short proof of this statement below.

###### Lemma 5

Let $(\zeta_{t})_{t\geq 0}$ be the process that solves (9) with initial data
$\theta_{0}$. Under Assumption 3, we have

$\displaystyle\left\|\zeta_{t}-\theta_{*}\right\|^{2}\leq\left\|\theta_{0}-\theta_{*}\right\|^{2}\exp({-\kappa
t}),$

where $\theta_{*}$ is a stationary solution of (9).

Proof  Since $\theta_{*}$ is a stationary solution,

$g(\theta_{*})=0\ \ \text{and}\ \
\mathrm{d}(\zeta_{t}-\theta_{*})=-(g(\zeta_{t})-g(\theta_{*}))\mathrm{d}t.$

Therefore,

$\displaystyle\frac{\mathrm{d}\left\|\zeta_{t}-\theta_{*}\right\|^{2}}{\mathrm{d}t}=2\left\langle\zeta_{t}-\theta_{*},\frac{\mathrm{d}(\zeta_{t}-\theta_{*})}{\mathrm{d}t}\right\rangle=-2\left\langle\zeta_{t}-\theta_{*},g(\zeta_{t})-g(\theta_{*})\right\rangle\leq-2\kappa\left\|\zeta_{t}-\theta_{*}\right\|^{2}.$

By Grönwall’s inequality,

$\displaystyle\left\|\zeta_{t}-\theta_{*}\right\|^{2}\leq\left\|\theta_{0}-\theta_{*}\right\|^{2}\exp({-\kappa
t}).$

###  3.2 Longtime behavior and ergodicity

We now study the longtime behavior of SGPC, i.e. the behavior and distribution
of $(\theta_{t}^{\varepsilon},V_{t/\varepsilon})$ for $t\gg 0$ large. Indeed,
the main result of this section will be the geometric ergodicity of this
coupled process and a study of its stationary measure. Initially, we study
stability of the stochastic gradient process
$(\theta^{\varepsilon}_{t})_{t\geq 0}.$

###### Lemma 6

Under Assumption 3, we have

$\left\|\theta^{\varepsilon}_{t}\right\|^{2}\leq\left\|\theta^{\varepsilon}_{0}\right\|^{2}\exp({-\kappa
t})+\frac{8K_{f}^{2}}{\kappa^{2}},$

where $K_{f}:=\sup_{y\in S}\left\|\nabla_{\theta}f(0,y)\right\|$.

Proof  By Itô’s formula, we have

$\displaystyle\frac{\mathrm{d}\left\|\theta^{\varepsilon}_{t}\right\|^{2}}{\mathrm{d}t}=2\left\langle\theta^{\varepsilon}_{t},d\theta^{\varepsilon}_{t}/dt\right\rangle=-2\left\langle\theta^{\varepsilon}_{t},\nabla_{\theta}f(\theta^{\varepsilon}_{t},V_{t/\varepsilon})\right\rangle.$
(16)

By Assumption 3,

$\displaystyle\left\langle\theta^{\varepsilon}_{t},\nabla_{\theta}f(\theta^{\varepsilon}_{t},V_{t/\varepsilon})\right\rangle=$
$\displaystyle\left\langle\theta^{\varepsilon}_{t}-0,\nabla_{\theta}f(\theta^{\varepsilon}_{t},V_{t/\varepsilon})-\nabla_{\theta}f(0,V_{t/\varepsilon})\right\rangle+\left\langle\theta^{\varepsilon}_{t},\nabla_{\theta}f(0,V_{t/\varepsilon})\right\rangle$
$\displaystyle\geq$
$\displaystyle\kappa\left\|\theta^{\varepsilon}_{t}\right\|^{2}-\left\|\theta^{\varepsilon}_{t}\right\|\left\|\nabla_{\theta}f(0,V_{t/\varepsilon})\right\|$
$\displaystyle\geq$
$\displaystyle\frac{\kappa}{2}\left\|\theta^{\varepsilon}_{t}\right\|^{2}-\frac{4}{\kappa}\left\|\nabla_{\theta}f(0,V_{t/\varepsilon})\right\|^{2}$
$\displaystyle\geq$
$\displaystyle\frac{\kappa}{2}\left\|\theta^{\varepsilon}_{t}\right\|^{2}-\frac{4K_{f}^{2}}{\kappa}.$

Hence (16) implies

$\displaystyle\frac{\mathrm{d}\left\|\theta^{\varepsilon}_{t}\right\|^{2}}{\mathrm{d}t}\leq-\kappa\left\|\theta^{\varepsilon}_{t}\right\|^{2}+\frac{8K_{f}^{2}}{\kappa}.$
(17)

Multiplying $\exp({\kappa t})$ on both sides of (17), we get

$\displaystyle\frac{\mathrm{d}(\left\|\theta^{\varepsilon}_{t}\right\|^{2}\exp({\kappa
t}))}{\mathrm{d}t}\leq\frac{8K_{f}^{2}\exp({\kappa t})}{\kappa},$

that is

$\displaystyle\left\|\theta^{\varepsilon}_{t}\right\|^{2}\exp({\kappa
t})-\left\|\theta^{\varepsilon}_{0}\right\|^{2}\leq\frac{8K_{f}^{2}(\exp({\kappa
t})-1)}{\kappa^{2}}\leq\frac{8K_{f}^{2}\exp({\kappa t})}{\kappa^{2}}.$

Therefore,

$\left\|\theta^{\varepsilon}_{t}\right\|^{2}\leq\left\|\theta^{\varepsilon}_{0}\right\|^{2}\exp({-\kappa
t})+\frac{8K_{f}^{2}}{\kappa^{2}}.$

Using this lemma, we are now able to prove the first main result of this
section, showing geometric ergodicity of
$(\theta^{\varepsilon}_{t},V_{t/\varepsilon})_{t\geq 0}.$ First, we introduce
a Wasserstein distance, on the space on which
$(\theta^{\varepsilon}_{t},V_{t/\varepsilon})_{t\geq 0}$ lives. Let $\Pi$ and
$\Pi^{\prime}$ be two probability measures on $(\mathbb{R}^{K}\times
S,\mathcal{B}(\mathbb{R}^{K}\times S))$. We define the Wasserstein distance
between those measures by

$\widetilde{{\mathcal{W}}}_{\tilde{d}}(\Pi,\Pi^{\prime})=\inf_{\widetilde{\Gamma}\in\mathcal{H}(\Pi,\Pi^{\prime})}\int_{(\mathbb{R}^{K}\times
S)\times(\mathbb{R}^{K}\times
S)}\tilde{d}((u,v),(u^{\prime},v^{\prime}))\widetilde{\Gamma}(\mathrm{d}u\mathrm{d}v,\mathrm{d}u^{\prime}\mathrm{d}v^{\prime}),$

where $\tilde{d}((u,v),(u^{\prime},v^{\prime})):=\boldsymbol{1}_{v\neq
v^{\prime}}+(1\land\left\|u-u^{\prime}\right\|)\boldsymbol{1}_{v=v^{\prime}}$.
For $a\in S$ and $m\in\mathbb{R}^{K}$, let $H^{\varepsilon}_{t}(\cdot|m,a)$ be
the distribution of $(\theta^{\varepsilon}_{t},V_{t/\varepsilon})$ under
$\mathbb{P}_{a}$ with $\theta^{\varepsilon}_{0}=m$. Moreover, recall that
$(V_{t})_{t\geq 0}$ is a Feller process that satisfies Assumption 1. More
specifically, it satisfies Assumption 1 (iii) with a constant $\delta$:

$\sup_{x\in S}\tilde{\mathbb{P}}(T^{x}\geq t)\leq C\exp({-\delta t}),$

where $T^{x}:=\inf{\\{t\geq 0\ |\ V^{x}_{t}=V^{\pi}_{t}\\}}$. With the
constant $\delta$ defined this way, we have the following theorem.

###### Theorem 7

Under Assumption 3, for any $0<\varepsilon\leq 1\land({\delta}/{2\kappa})$,
the (coupled) process $(\theta_{t}^{\varepsilon},V_{t/\varepsilon})_{t\geq 0}$
admits an unique stationary measure $\Pi^{\varepsilon}$ on
$(\mathbb{R}^{K}\times S,\mathcal{B}(\mathbb{R}^{K}\times S)).$ Moreover,

$\displaystyle\widetilde{{\mathcal{W}}}^{2}_{\tilde{d}}(H^{\varepsilon}_{t}(\cdot|m,a),\Pi^{\varepsilon})\leq
C_{f}\exp({-\kappa
t})\int_{\mathbb{R}^{K}}(1+\left\|x-m\right\|^{2})\Pi^{\varepsilon}(\mathrm{d}x,S),$
(18)
$\displaystyle\widetilde{{\mathcal{W}}}^{2}_{\tilde{d}}(H^{\varepsilon}_{t}(\cdot|m,\pi),\Pi^{\varepsilon})\leq
C_{f}\exp({-\kappa
t})\int_{\mathbb{R}^{K}}\left\|x-m\right\|^{2}\Pi^{\varepsilon}(\mathrm{d}x,S),$
(19)

where the constant $C_{f}$ only depends on $f$.

Proof  To obtain the existence of the invariant measure, we apply the weak
form of Harris’ Theorem in Cloez and Hairer (2015, Theorem 3.7) by verifying
the Lyapunov condition, the $\tilde{d}$-contracting condition, and the
$\tilde{d}$-small condition.

(i) Lyapunov condition: Let $V(x,y)=\left\|x\right\|^{2}$. To verify that it
satisfies (2.1) in Cloez and Hairer (2015, Definition 2.1), we take
$x=\theta^{\varepsilon}_{0}$ and $(P_{t})_{t\geq 0}$ to be the semi-group
associated with $(\theta^{\varepsilon}_{t})_{t\geq 0}$. Then Lemma 6 yields
that $V(x,y)$ is a Lyapunov function. The existence of the Lyapunov function
can be understood as that the coupled process does not go to infinity.

(ii) $\tilde{d}$-contracting condition: $\tilde{d}$-contracting states that
there exists $t^{*}>0$ such that for any $t>t^{*}$, there exists some
$\alpha<1$ such that

$\widetilde{{\mathcal{W}}}_{\tilde{d}}(\delta_{(m,a)}P^{\varepsilon}_{t},\delta_{(n,b)}P^{\varepsilon}_{t})\leq\alpha\tilde{d}((m,a),(n,b))$

for any $(m,a),(n,b)\in\mathbb{R}^{K}\times S$ such that
$\tilde{d}((m,a),(n,b))<1$. Here, $(P^{\varepsilon}_{t})_{t\geq 0}$ is the
semi-group operator associated with (8). Notice that
$\tilde{d}((m,a),(n,b))<1$ implies $a=b$. Let $(\theta^{(m,a)}_{t})_{t\geq
0},\ (\theta^{(n,a)}_{t})_{t\geq 0}$ solve the following equations:

$\displaystyle\theta^{(m,a)}_{t}=m-\int_{0}^{t}\nabla_{\theta}f(\theta^{(m,a)}_{t},\tilde{V}^{a}_{s/\varepsilon})\mathrm{d}s,$
$\displaystyle\theta^{(n,a)}_{t}=n-\int_{0}^{t}\nabla_{\theta}f(\theta^{(n,a)}_{t},\tilde{V}^{a}_{s/\varepsilon})\mathrm{d}s,$

where $\tilde{V}^{a}_{t}=V^{a}_{t}$ for $t\leq T^{a}$ and
$\tilde{V}^{a}_{t}=V^{\pi}_{t}$ for $t>T^{a}.$ Then by Itô’s formula and
Assumption 3,

$\displaystyle\mathrm{d}\left\|\theta^{(m,a)}_{t}-\theta^{(n,a)}_{t}\right\|^{2}/\mathrm{d}t=$
$\displaystyle-2\left\langle\theta^{(m,a)}_{t}-\theta^{(n,a)}_{t},\nabla_{\theta}f(\theta^{(m,a)}_{t},\tilde{V}^{a}_{s/\varepsilon})-\nabla_{\theta}f(\theta^{(n,a)}_{t},\tilde{V}^{a}_{s/\varepsilon})\right\rangle$
$\displaystyle\leq$
$\displaystyle-\kappa\left\|\theta^{(m,a)}_{t}-\theta^{(n,a)}_{t}\right\|^{2}.$

By Grönwall’s inequality,

$\left\|\theta^{(m,a)}_{t}-\theta^{(n,a)}_{t}\right\|^{2}\leq\exp({-\kappa
t})\left\|m-n\right\|^{2}.$ (20)

Noticing that $\tilde{d}((m,a),(n,b))<1$ implies $\left\|m-n\right\|^{2}<1$,
by choosing $t\geq\frac{1}{\kappa},$ we obtain

$\displaystyle\left\|\theta^{(m,a)}_{t}-\theta^{(n,a)}_{t}\right\|^{2}$
$\displaystyle\leq\exp({-1})\left\|m-n\right\|^{2}$
$\displaystyle=\exp({-1})(\left\|m-n\right\|^{2}\land
1)=\exp({-1})\tilde{d}^{2}((m,a),(n,a)).$

Therefore, with $t^{*}=1/\kappa$,

$\displaystyle\widetilde{{\mathcal{W}}}_{\tilde{d}}(H^{\varepsilon}_{t}(\cdot|m,a),H^{\varepsilon}_{t}(\cdot|n,b))\leq\tilde{\mathbb{E}}\Big{[}\left\|\theta^{(m,a)}_{t}-\theta^{(n,a)}_{t}\right\|\Big{]}\leq\exp({-1})\tilde{d}((m,a),(n,b)).$

(iii) $\tilde{d}$-small condition: We shall verify that there exists $t_{*}>0$
such that for any $t>t_{*}$, the sublevel set
$\mathscr{V}:=\\{(x,y)\in\mathbb{R}^{K}\times S\ |\
V(x,y)\leq{32K_{f}^{2}}/{\kappa^{2}}\\}$ is $\tilde{d}$-small for
$(P^{\varepsilon}_{t})_{t\geq 0}$, meaning that there exists a constant
$\zeta$ such that

$\widetilde{{\mathcal{W}}}_{\tilde{d}}(\delta_{(m,a)}P^{\varepsilon}_{t},\delta_{(n,b)}P^{\varepsilon}_{t})\leq
1-\zeta,$

for all $(m,a),(n,b)\in\mathscr{V}$. Let $(\theta^{(m,a)}_{t})_{t\geq 0},\
(\theta^{(n,b)}_{t})_{t\geq 0}$ solve the following equations:

$\displaystyle\theta^{(m,a)}_{t}=m-\int_{0}^{t}\nabla_{\theta}f(\theta^{(m,a)}_{t},\tilde{V}^{a}_{s/\varepsilon})\mathrm{d}s,$
$\displaystyle\theta^{(n,b)}_{t}=n-\int_{0}^{t}\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{b}_{s/\varepsilon})\mathrm{d}s,$

where $\tilde{V}^{a}_{t}=V^{a}_{t}$ for $t\leq T^{a}$ and
$\tilde{V}^{a}_{t}=V^{\pi}_{t}$ for $t>T^{a}$; $\tilde{V}^{b}_{t}=V^{b}_{t}$
for $t\leq T^{b}$ and $\tilde{V}^{b}_{t}=V^{\pi}_{t}$ for $t>T^{b}$. By Itô’s
formula, Assumption 3, and the $\varepsilon$-Young inequality,

$\displaystyle\mathrm{d}\left\|\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t}\right\|^{2}/\mathrm{d}t$
$\displaystyle=$
$\displaystyle-2\left\langle\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t},\nabla_{\theta}f(\theta^{(m,a)}_{t},\tilde{V}^{a}_{t/\varepsilon})-\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{b}_{t/\varepsilon})\right\rangle$
$\displaystyle=$
$\displaystyle-2\left\langle\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t},\nabla_{\theta}f(\theta^{(m,a)}_{t},\tilde{V}^{a}_{t/\varepsilon})-\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{a}_{t/\varepsilon})\right\rangle$
$\displaystyle\ \
-2\left\langle\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t},\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{a}_{t/\varepsilon})-\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{b}_{t/\varepsilon})\right\rangle$
$\displaystyle\leq$
$\displaystyle-2\kappa\left\|\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t}\right\|^{2}+2\left|\left\langle\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t},\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{a}_{t/\varepsilon})-\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{b}_{t/\varepsilon})\right\rangle\right|$
$\displaystyle\leq$
$\displaystyle-\kappa\left\|\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t}\right\|^{2}+\frac{4}{\kappa}\left\|\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{a}_{t/\varepsilon})-\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{b}_{t/\varepsilon})\right\|^{2}.$

Multiplying $\exp({\kappa t})$ on both sides, we obtain

$\displaystyle\mathrm{d}\Big{(}\exp({\kappa
t})\left\|\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t}\right\|^{2}\Big{)}/\mathrm{d}t\leq\frac{4\exp({\kappa
t})}{\kappa}\left\|\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{a}_{t/\varepsilon})-\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{b}_{t/\varepsilon})\right\|^{2},$

that is

$\displaystyle\exp({\kappa
t})\left\|\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t}\right\|^{2}$
$\displaystyle\leq\left\|m-n\right\|^{2}$
$\displaystyle\qquad+\frac{4}{\kappa}\int_{0}^{t}\exp({\kappa
s})\left\|\nabla_{\theta}f(\theta^{(n,b)}_{s},\tilde{V}^{a}_{t/\varepsilon})-\nabla_{\theta}f(\theta^{(n,b)}_{s},\tilde{V}^{b}_{t/\varepsilon})\right\|^{2}\mathrm{d}s.$

Notice that $\tilde{V}^{a}_{t}=\tilde{V}^{b}_{t}$ if $t>T^{a}\lor T^{b}$.
Hence, we have

$\displaystyle\exp({\kappa
t})\left\|\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t}\right\|^{2}$
$\displaystyle\qquad\leq\left\|m-n\right\|^{2}+\frac{4}{\kappa}\int_{0}^{t\land\varepsilon(T^{a}\lor
T^{b})}\exp({\kappa
s})\left\|\nabla_{\theta}f(\theta^{(n,b)}_{s},\tilde{V}^{a}_{t/\varepsilon})-\nabla_{\theta}f(\theta^{(n,b)}_{s},\tilde{V}^{b}_{t/\varepsilon})\right\|^{2}\mathrm{d}s.$

For $(m,a),(n,b)\in\mathscr{V}\times S$, by Lemma 6, we have that

$\left\|\nabla_{\theta}f(\theta^{(m,a)}_{t},\tilde{V}^{a}_{t/\varepsilon})\right\|\text{
and
}\left\|\nabla_{\theta}f(\theta^{(n,b)}_{t},\tilde{V}^{b}_{t/\varepsilon})\right\|$

are bounded by some constant $C_{f}.$ This implies

$\displaystyle\left\|\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t}\right\|^{2}\leq\exp({-\kappa
t})\left\|m-n\right\|^{2}+4C_{f}\exp({-\kappa
t})\exp({\kappa\varepsilon(T^{a}\lor T^{b})}).$ (21)

Recall $0<\varepsilon\leq 1\land\frac{\delta}{2\kappa}$. By (iii), Assumption
1,

$\displaystyle\tilde{\mathbb{E}}[\exp({\kappa\varepsilon(T^{a}\lor T^{b})})]=$
$\displaystyle\kappa\varepsilon\int_{0}^{\infty}\exp({\kappa\varepsilon
x})\tilde{\mathbb{P}}(T^{a}\lor T^{b}\geq x)\mathrm{d}x$ $\displaystyle\leq$
$\displaystyle C\kappa\varepsilon\int_{0}^{\infty}\exp({\kappa\varepsilon
x})\exp({-\delta x})\mathrm{d}x$ $\displaystyle\leq$
$\displaystyle\frac{C\delta}{2}\int_{0}^{\infty}\exp\left({-\frac{\delta
x}{2}}\right)\mathrm{d}x=C.$

Therefore,

$\displaystyle\tilde{\mathbb{E}}\Big{[}\left\|\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t}\right\|^{2}\Big{]}\leq\exp({-\kappa
t})\left\|m-n\right\|^{2}+C_{f}\exp({-\kappa t}).$

Moreover,

$\tilde{\mathbb{P}}(\tilde{V}^{a}_{t/\varepsilon}\neq\tilde{V}^{b}_{t/\varepsilon})=\tilde{\mathbb{P}}(T^{a}\lor
T^{b}\geq{t}/{\varepsilon})\leq C\exp({-\delta{t}/{\varepsilon}})\leq
C\exp({-2\kappa t}).$

Hence for any $(m,a),(n,b)\in\mathscr{V}\times S$, by taking
$t\geq\frac{1}{\kappa}[\log(8C_{f})+\log({512K_{f}^{2}}/{\kappa^{2}})+\frac{1}{2}\log(8C)],$
we have

$\widetilde{{\mathcal{W}}}^{2}_{\tilde{d}}(H^{\varepsilon}_{t}(\cdot|m,a),H^{\varepsilon}_{t}(\cdot|n,b))\leq\tilde{\mathbb{P}}(\tilde{V}^{a}_{t/\varepsilon}\neq\tilde{V}^{b}_{t/\varepsilon})+\tilde{\mathbb{E}}\Big{[}\left\|\theta^{(m,a)}_{t}-\theta^{(n,b)}_{t}\right\|^{2}\Big{]}\leq\frac{1}{2}.$

We have verified all three conditions from Cloez and Hairer (2015, Theorem
3.7) and hence conclude the existence and uniqueness of the invariant measure
of $(\theta_{t}^{\varepsilon},V_{t/\varepsilon})_{t\geq 0}$ denoted as
$\Pi^{\varepsilon}$.

Next, we are going to prove (18) and (19). Define $\Theta^{\varepsilon}$ such
that $(\Theta^{\varepsilon},\tilde{V}^{\pi}_{0})\sim\Pi^{\varepsilon}$ in
$\tilde{\mathbb{P}}.$ Let $\theta^{\Pi^{\varepsilon}}_{t}$ and
$\theta^{(m,\pi)}_{t}$ solve following equations

$\displaystyle\theta^{\Pi^{\varepsilon}}_{t}=\Theta^{\varepsilon}-\int_{0}^{t}\nabla_{\theta}f(\theta^{\Pi^{\varepsilon}}_{t},\tilde{V}^{\pi}_{t/\varepsilon})\mathrm{d}s,$
$\displaystyle\theta^{(m,\pi)}_{t}=m-\int_{0}^{t}\nabla_{\theta}f(\theta^{(m,\pi)}_{t},\tilde{V}^{\pi}_{t/\varepsilon})\mathrm{d}s.$

Recall that from (21) and (20), we have

$\displaystyle\left\|\theta^{(m,a)}_{t}-\theta^{\Pi^{\varepsilon}}_{t}\right\|^{2}\leq$
$\displaystyle\exp({-\kappa
t})\left\|m-\Theta^{\varepsilon}\right\|^{2}+4C_{f}\exp({-\kappa
t})\exp({\kappa\varepsilon T^{a}}),$
$\displaystyle\left\|\theta^{(m,\pi)}_{t}-\theta^{\Pi^{\varepsilon}}_{t}\right\|^{2}\leq$
$\displaystyle\exp({-\kappa t})\left\|m-\Theta^{\varepsilon}\right\|^{2}.$

Therefore, by (iii), Assumption 1 and $0<\varepsilon\leq
1\land({\delta}/{2\kappa})$,

$\displaystyle\widetilde{{\mathcal{W}}}^{2}_{\tilde{d}}(H^{\varepsilon}_{t}(\cdot|m,a),\Pi^{\varepsilon})\leq$
$\displaystyle\tilde{\mathbb{P}}(\tilde{V}^{a}_{t/\varepsilon}\neq\tilde{V}^{\pi}_{t/\varepsilon})+\tilde{\mathbb{E}}\Big{[}\left\|\theta^{(m,a)}_{t}-\theta^{\Pi^{\varepsilon}}_{t}\right\|^{2}\Big{]}$
$\displaystyle\leq$
$\displaystyle\tilde{\mathbb{P}}(T^{a}\geq\frac{t}{\varepsilon})+\exp({-\kappa
t})\tilde{\mathbb{E}}[\left\|m-\Theta^{\varepsilon}\right\|^{2}]+4C_{f}\exp({-\kappa
t})\tilde{\mathbb{E}}[\exp({\kappa\varepsilon T^{a}})]$ $\displaystyle\leq$
$\displaystyle C_{f}\exp({-\kappa
t})\int_{\mathbb{R}^{K}}(1+\left\|x-m\right\|^{2})\Pi^{\varepsilon}(\mathrm{d}x,S),$
$\displaystyle\widetilde{{\mathcal{W}}}^{2}_{\tilde{d}}(H^{\varepsilon}_{t}(\cdot|m,\pi),\Pi^{\varepsilon})\leq$
$\displaystyle\tilde{\mathbb{E}}\Big{[}\left\|\theta^{(m,\pi)}_{t}-\theta^{\Pi^{\varepsilon}}_{t}\right\|^{2}\Big{]}$
$\displaystyle\leq$ $\displaystyle\exp({-\kappa
t})\tilde{\mathbb{E}}[\left\|m-\Theta^{\varepsilon}\right\|^{2}]$
$\displaystyle\leq$ $\displaystyle C_{f}\exp({-\kappa
t})\int_{\mathbb{R}^{K}}\left\|x-m\right\|^{2}\Pi^{\varepsilon}(\mathrm{d}x,S).$

In the following corollaries, we study the integrals on the right-hand side of
the inequalities in Theorem 7. Moreover, we show that the result above
immediately implies not only geometric ergodicity of the coupled process
$(\theta^{\varepsilon}_{t},V_{t/\varepsilon})_{t\geq 0}$, but also of its
marginal, the stochastic gradient process $(\theta^{\varepsilon}_{t})_{t\geq
0}$.

###### Corollary 8

Under the same assumptions as Theorem 7, there exists a constant $C_{f,m}$
that depends only on $f$ and the initial value $m=\theta_{0}^{\varepsilon}$,
such that

$\displaystyle\widetilde{{\mathcal{W}}}_{\tilde{d}}(H^{\varepsilon}_{t}(\cdot|m,a),\Pi^{\varepsilon})\leq
C_{f,m}\exp\left({-\frac{\kappa t}{2}}\right),$ (22)
$\displaystyle\widetilde{{\mathcal{W}}}_{\tilde{d}}(H^{\varepsilon}_{t}(\cdot|m,\pi),\Pi^{\varepsilon})\leq
C_{f,m}\exp\left({-\frac{\kappa t}{2}}\right),$ (23)
$\displaystyle{{\mathcal{W}}}_{d}(C^{\varepsilon}_{t}(\cdot|m,a),\Pi^{\varepsilon}(\cdot,S))\leq
C_{f,m}\exp\left({-\frac{\kappa t}{2}}\right),$ (24)
$\displaystyle{{\mathcal{W}}}_{d}(C^{\varepsilon}_{t}(\cdot|m,\pi),\Pi^{\varepsilon}(\cdot,S))\leq
C_{f,m}\exp\left({-\frac{\kappa t}{2}}\right).$ (25)

Proof  By Lemma 6,

$\int_{\left\|x\right\|^{2}\geq\frac{8K^{2}}{\kappa^{2}}+\left\|m\right\|^{2}+1}\Pi^{\varepsilon}(\mathrm{d}x,S)=\lim_{t\to\infty}\tilde{\mathbb{P}}(\left\|\theta^{(m,a)}_{t}\right\|\geq\frac{8K_{f}^{2}}{\kappa^{2}}+\left\|m\right\|^{2}+1)=0.$

Let $C_{f,m}=C_{f}(2\left\|m\right\|+\frac{4K_{f}}{\kappa}+2)$. From (18), we
have

$\displaystyle\tilde{{\mathcal{W}}}_{\tilde{d}}(H^{\varepsilon}_{t}(\cdot|m,a),\Pi^{\varepsilon})\leq$
$\displaystyle C_{f}\exp\left({\frac{-\kappa
t}{2}}\right)\Big{(}\int_{\mathbb{R}^{K}}(1+\left\|x-m\right\|^{2})\Pi^{\varepsilon}(\mathrm{d}x,S)\Big{)}^{1/2}$
$\displaystyle=$ $\displaystyle C_{f}\exp\left({\frac{-\kappa
t}{2}}\right)\Big{(}\int_{\left\|x\right\|^{2}\leq\frac{8K_{f}^{2}}{\kappa^{2}}+\left\|m\right\|^{2}+1}(1+\left\|x-m\right\|^{2})\Pi^{\varepsilon}(\mathrm{d}x,S)\Big{)}^{1/2}$
$\displaystyle\leq$ $\displaystyle C_{f,m}\exp\left({\frac{-\kappa
t}{2}}\right).$

(23) can be derived similarly from (19). Moreover, notice that

$\displaystyle{{\mathcal{W}}}_{d}(C^{\varepsilon}_{t}(\cdot|m,a),\Pi^{\varepsilon}(\cdot,S))\leq$
$\displaystyle\tilde{\mathbb{E}}\Big{[}\left\|\theta^{(m,a)}_{t}-\theta^{\Pi^{\varepsilon}}_{t}\right\|\Big{]},$
(26)
$\displaystyle{{\mathcal{W}}}_{d}(C^{\varepsilon}_{t}(\cdot|m,\pi),\Pi^{\varepsilon}(\cdot,S))\leq$
$\displaystyle\tilde{\mathbb{E}}\Big{[}\left\|\theta^{(m,\pi)}_{t}-\theta^{\Pi^{\varepsilon}}_{t}\right\|\Big{]}.$
(27)

(24) and (25) can be derived similarly from (26) and (27).
Combining (24) and (25), we immediately have:

###### Corollary 9

Under the same assumptions as Theorem 7,

$\displaystyle{{\mathcal{W}}}_{d}(C^{\varepsilon}_{t}(\cdot|m,a),C^{\varepsilon}_{t}(\cdot|m,\pi))\leq
C_{f,m}\exp\left({-\frac{\kappa t}{2}}\right).$

So far, we have shown that the stochastic gradient process
$(\theta^{\varepsilon}_{t})_{t\geq 0}$ converges to a unique stationary
measure $\Pi^{\varepsilon}(\cdot,S)$. It is often not possible to determine
this stationary measure. However, we can comment on its asymptotic behavior as
$\varepsilon\rightarrow 0$. Indeed, we will show that
$\Pi^{\varepsilon}(\cdot,S)$ concentrates around the minimizer $\theta_{*}$ of
the full target function.

###### Proposition 3

Under Assumption 3, the measure $\Pi^{\varepsilon}(\cdot,S)$ on
$(\mathbb{R}^{K},\mathcal{B}(\mathbb{R}^{K}))$ approximates
$\delta(\cdot-\theta_{*}).$ In other words, we have

$\displaystyle{{\mathcal{W}}}_{d}(\Pi^{\varepsilon}(\cdot,S),\delta(\cdot-\theta_{*}))\leq\rho(\varepsilon)$

where $\rho:(0,1)\rightarrow[0,1]$ and $\lim_{\varepsilon\to
0}\rho(\varepsilon)=0.$

Proof  By the triangle inequality,

$\displaystyle{{\mathcal{W}}}_{d}(\Pi^{\varepsilon}(\cdot,S),\delta(\cdot-\theta_{*}))$
$\displaystyle\leq$
$\displaystyle{{\mathcal{W}}}_{d}(\Pi^{\varepsilon}(\cdot,S),C^{\varepsilon}_{t}(\cdot|m,\pi))+{{\mathcal{W}}}_{d}(C^{\varepsilon}_{t}(\cdot|m,\pi),\delta(\cdot-\zeta_{t}))+{{\mathcal{W}}}_{d}(\delta(\cdot-\zeta_{t}),\delta(\cdot-\theta_{*})).$

Let $\theta_{0}=\theta^{\varepsilon}_{0}=\theta_{*}.$ Then by Lemma 5, we have
the last term

${{\mathcal{W}}}_{d}(\delta(\cdot-\zeta_{t}),\delta(\cdot-\theta_{*}))\leq\left\|\theta_{*}-\theta_{*}\right\|\exp({-\kappa
t})=0.$

By (22) and Corollary 4, for any $t\geq 0,$

$\displaystyle{{\mathcal{W}}}_{d}(\Pi^{\varepsilon}(\cdot,S),C^{\varepsilon}_{t}(\cdot|m,\pi))+{{\mathcal{W}}}_{d}(C^{\varepsilon}_{t}(\cdot|m,\pi),\delta(\cdot-\zeta_{t}))\leq
C_{f,\theta_{*}}\exp\left({-\frac{\kappa
t}{2}}\right)+(\exp(t)\alpha(\varepsilon))\land 1.$

By choosing $t=-\log({{1\land\alpha(\varepsilon)}})/2,$ we get

$\displaystyle{{\mathcal{W}}}_{d}(\Pi^{\varepsilon}(\cdot,S),\delta(\cdot-\theta_{*}))\leq$
$\displaystyle C_{f,\theta_{*}}\exp\left({-\frac{\kappa
t}{2}}\right)+(\exp(t)\alpha(\varepsilon))\land 1$ $\displaystyle\leq$
$\displaystyle
C_{f,\theta_{*}}(\alpha(\varepsilon))^{\frac{\kappa}{4}}+(\alpha(\varepsilon))^{1/2}.$

Taking
$\rho(\varepsilon):=\big{(}C_{f,\theta_{*}}(\alpha(\varepsilon))^{\frac{\kappa}{4}}+(\alpha(\varepsilon))^{1/2}\big{)}\land
1$ completes the proof.

##  4 Stochastic gradient processes with decreasing learning rate

Constant learning rates are popular in some practical situations, but the
associated stochastic gradient process usually does not converge to the
minimizer of $\Phi$. This is also true for the discrete-time stochastic
gradient descent algorithm. However, SGD can converge to the minimizer if the
learning rate is decreased over time. In the following, we discuss a
decreasing learning version of the stochastic gradient process and show that
this dynamical system indeed converges to the minimizer of $\Phi$.

As discussed in Section 1.2, we obtain the stochastic gradient process with
decreasing learning rate by non-linearly rescaling the time in the constant-
learning-rate index process $(V_{t})_{t\geq 0}$. Indeed, we choose a function
$\beta:[0,\infty)\rightarrow[0,\infty)$ and then define the decreasing
learning rate index process by $(V_{\beta(t)})_{t\geq 0}$. We have discussed
an intuitive way to construct a rescaling function $\beta$ also in Section
1.2.

In the following, we define $\beta$ through an integral
$\beta(t)=\int_{0}^{t}\mu(s)\mathrm{d}s,$ $t\geq 0$. We commence this section
with necessary growth conditions on $\mu$ which allow us to then give the
formal definition of the stochastic gradient process with decreasing learning
rate. Then, we study the longtime behavior of this process.

###### Assumption 4

Let $\mu:[0,\infty)\to(0,\infty)$ be a non-decreasing continuously
differentiable function with $\lim_{t\to\infty}\mu(t)=\infty$ and

$\displaystyle\lim_{t\to\infty}\frac{\mu^{\prime}(t)t}{\mu(t)}=0.$

Assumption 4 implies that $\mu$ goes to infinity, but at a very slow pace.
Indeed it says that $\lim_{t\to\infty}\mu(t)/t^{\gamma}=0$, $\gamma>0$, that
is $\mu(t)$ grows slower than any polynomial.

###### Definition 2

The _stochastic gradient process with decreasing learning rate (SGPD)_ is a
solution of the following stochastic differential equation,

$\left\\{\begin{array}[]{l}\mathrm{d}\xi_{t}=-\nabla_{\xi}f(\xi_{t},V_{\beta{(t)}})\mathrm{d}t,\\\
\xi_{0}=\theta_{0},\end{array}\right.$ (28)

where f satisfies Assumption 2, $(V_{t})_{t\geq 0}$ is a Feller process that
satisfies Assumption 1, and $\beta(t)=\int_{0}^{t}\mu(s)\mathrm{d}s$ with
$\mu$ satisfying Assumption 4.

To see that $(\xi_{t})_{t\geq 0}$ is well-defined, consider the following:
$t\mapsto\beta(t)$ is an increasing continuous function. Thus,
$(V_{\beta(t)})_{t\geq 0}$ is càdlàg and Feller with respect to
$(\mathcal{F}_{\beta(t)})_{t\geq 0}$. We then obtain well-definedness of
$(\xi_{t})_{t\geq 0}$ by replacing $(V_{t/\varepsilon})_{t\geq 0}$ by
$(V_{\beta(t)})_{t\geq 0}$ in the proof of Proposition 2.

We now move on to studying the longtime behavior of the SGPD $(\xi_{t})_{t\geq
0}$. In a first technical result, we establish a connection between SGPD
$(\xi_{t})_{t\geq 0}$ and a time-rescaled version of SGPC
$(\theta_{t}^{\varepsilon})_{t\geq 0}$. To this end, note that
$\dot{\beta}(t)=\mu(t)>0$, $\beta(t)$ is strictly increasing. Hence, the
inverse function of $\beta(t)$ exists and

$(\beta)^{-1}(t)=\int_{0}^{t}\frac{1}{\mu((\beta)^{-1}(s)))}\mathrm{d}s.$

This gives us the following inequality.

###### Proposition 4

For any $0<\varepsilon<1$,

$\displaystyle\left\|\xi_{t}-\theta^{\varepsilon}_{\varepsilon\beta(t)}\right\|^{2}\leq
C_{f,\theta_{0},\mu}\left[\frac{\exp({-2\varepsilon\kappa(\beta(t)-\beta(\frac{t}{2}))})}{\varepsilon}+\frac{1}{\varepsilon}\Big{(}\left|\frac{1}{\mu(t)}-\varepsilon\right|+\left|\frac{1}{\mu(\frac{t}{2})}-\varepsilon\right|\Big{)}\right]$

almost surely, where the constant $C_{f,\theta_{0},\mu}$ depending only on
$f$, the initial data $\theta_{0}$, and $\mu$.

Proof  From (28) and (8),

$\displaystyle\xi_{t}=$
$\displaystyle\theta_{0}-\int_{0}^{\beta(t)}\nabla_{\xi}f(\xi_{(\beta)^{-1}(s)},V_{s})\mathrm{d}(\beta)^{-1}(s),$
$\displaystyle\theta^{\varepsilon}_{t}=$
$\displaystyle\theta_{0}-\varepsilon\int_{0}^{\frac{t}{\varepsilon}}\nabla_{\theta}f(\theta^{\varepsilon}_{\varepsilon
s},V_{s})\mathrm{d}s.$

Let $b_{t}:=\mathrm{d}(\beta)^{-1}(t)/\mathrm{d}t=1/\mu((\beta)^{-1}(t)))>0,$
we have

$\displaystyle\xi_{(\beta)^{-1}(t)}=$
$\displaystyle\theta_{0}-\int_{0}^{t}\nabla_{\xi}f(\xi_{(\beta)^{-1}(s)},V_{s})b_{s}\mathrm{d}s$
$\displaystyle\theta^{\varepsilon}_{\varepsilon t}=$
$\displaystyle\theta_{0}-\varepsilon\int_{0}^{t}\nabla_{\theta}f(\theta^{\varepsilon}_{\varepsilon
s},V_{s})\mathrm{d}s$

Therefore, by Itô’s formula and Assumption 3,

$\displaystyle\mathrm{d}\left\|\theta^{\varepsilon}_{\varepsilon
t}-\xi_{(\beta)^{-1}(t)}\right\|^{2}/\mathrm{d}t=$
$\displaystyle-2\left\langle\theta^{\varepsilon}_{\varepsilon
t}-\xi_{(\beta)^{-1}(t)},\varepsilon\nabla_{\theta}f(\theta^{\varepsilon}_{\varepsilon
t},V_{t})-\varepsilon\nabla_{\xi}f(\xi_{(\beta)^{-1}(t)},V_{t})\right\rangle$
$\displaystyle-2(\varepsilon-
b_{t})\left\langle\theta^{\varepsilon}_{\varepsilon
t}-\xi_{(\beta)^{-1}(t)},\nabla_{\xi}f(\xi_{(\beta)^{-1}(t)},V_{t})\right\rangle$
$\displaystyle\leq$
$\displaystyle-2\varepsilon\kappa\left\|\theta^{\varepsilon}_{\varepsilon
t}-\xi_{(\beta)^{-1}(t)}\right\|^{2}+C_{f,\theta_{0}}\left|b_{t}-\varepsilon\right|,$

where the last step follows from the boundedness of
$\theta^{\varepsilon}_{\varepsilon t}$, $\xi_{(\beta)^{-1}(t)}$, and
$\nabla_{\xi}f(\xi_{(\beta)^{-1}(t)},V_{t})$. $\xi_{(\beta)^{-1}(t)}$ is
bounded can be showed similarly to Lemma 6. Multiplying
$\exp({2\varepsilon\kappa t})$ on both sides, we obtain

$\displaystyle\mathrm{d}\Big{(}\exp({2\varepsilon\kappa
t})\left\|\theta^{\varepsilon}_{\varepsilon
t}-\xi_{(\beta)^{-1}(t)}\right\|^{2}\Big{)}/\mathrm{d}t\leq$ $\displaystyle
C_{f,\theta_{0}}\left|b_{t}-\varepsilon\right|\exp({2\varepsilon\kappa t}),$

which implies

$\displaystyle\left\|\theta^{\varepsilon}_{\varepsilon
t}-\xi_{(\beta)^{-1}(t)}\right\|^{2}\leq$ $\displaystyle
C_{f,\theta_{0}}\exp({-2\varepsilon\kappa
t})\int_{0}^{t}\left|b_{s}-\varepsilon\right|\exp({2\varepsilon\kappa
s})\mathrm{d}s.$

Notice that $b_{s}$ is bounded and non-increasing, hence we have

$\displaystyle\left\|\xi_{t}-\theta^{\varepsilon}_{\varepsilon\beta(t)}\right\|^{2}\leq$
$\displaystyle
C_{f,\theta_{0}}\exp({-2\varepsilon\kappa\beta(t)})\int_{0}^{\beta(t)}\left|b_{s}-\varepsilon\right|\exp({2\varepsilon\kappa
s})\mathrm{d}s$ $\displaystyle=$ $\displaystyle
C_{f,\theta_{0}}\exp({-2\varepsilon\kappa\beta(t)})\Big{(}\int_{0}^{\beta(\frac{t}{2})}+\int^{\beta(t)}_{\beta(\frac{t}{2})}\Big{)}\left|b_{s}-\varepsilon\right|\exp({2\varepsilon\kappa
s})\mathrm{d}s$ $\displaystyle\leq$ $\displaystyle
C_{f,\theta_{0},\mu}\frac{\exp({-2\varepsilon\kappa(\beta(t)-\beta(\frac{t}{2}))})}{\varepsilon}$
$\displaystyle\qquad+C_{f,\theta_{0}}\exp({-2\varepsilon\kappa\beta(t)})\int^{\beta(t)}_{\beta(\frac{t}{2})}\left|b_{s}-\varepsilon\right|\exp({2\varepsilon\kappa
s})\mathrm{d}s$ $\displaystyle\leq$ $\displaystyle
C_{f,\theta_{0},\mu}\left[\frac{\exp({-2\varepsilon\kappa(\beta(t)-\beta(\frac{t}{2}))})}{\varepsilon}+\frac{1}{\varepsilon}\Big{(}\left|\frac{1}{\mu(t)}-\varepsilon\right|+\left|\frac{1}{\mu(\frac{t}{2})}-\varepsilon\right|\Big{)}\right].$

Now, we get to the main result of this section, where we show the convergence
of $(\xi_{t})_{t\geq 0}$ to the minimizer $\theta_{*}$ of $\Phi$. In the
following, we denote

$\displaystyle D_{t}(B|\theta_{0},a)$
$\displaystyle:=\mathbb{P}_{a}(\xi_{t}\in B|\xi_{0}=\theta_{0}),$
$\displaystyle D_{t}(B|\theta_{0},\pi)$
$\displaystyle:=\mathbb{P}_{\pi}(\xi_{t}\in
B|\xi_{0}=\theta_{0})\qquad\qquad(B\in\mathcal{B}(\mathbb{R}^{K}),\theta_{0}\in\mathbb{R}^{K}),$

where $a\in S$ and $\pi$ is the invariant measure of $(V_{t})_{t\geq 0}$,
respectively.

###### Theorem 10

Under Assumption 3, given $\theta_{0}\in\mathbb{R}^{K}$ and $a\in S$, there
exists $T>0$ such that for any $t>T$,

$\displaystyle{{\mathcal{W}}}_{d}(D_{t}(\cdot|\theta_{0},\pi),\delta(\cdot-\theta_{*}))\leq
C_{f,\theta_{0},\mu}A(t),$ (29)
$\displaystyle{{\mathcal{W}}}_{d}(D_{t}(\cdot|\theta_{0},a),\delta(\cdot-\theta_{*}))\leq
C_{f,\theta_{0},\mu}A(t),$ (30)

where

$A(t):=\exp\left({\frac{-\kappa
t}{8}}\right)+\left[\frac{\mu(t)-\mu(\frac{t}{2})}{\mu(t)}\right]^{1/2}+\rho\left(\frac{1}{\mu(\frac{t}{2})}\right)$

and $\lim_{t\to\infty}A(t)=0.$

Proof  To prove (29), by the triangle inequality,

$\displaystyle{{\mathcal{W}}}_{d}(D_{t}(\cdot|\theta_{0},\pi),\delta(\cdot-\theta_{*}))$
$\displaystyle\leq$
$\displaystyle{{\mathcal{W}}}_{d}(D_{t}(\cdot|\theta_{0},\pi),C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},\pi))+{{\mathcal{W}}}_{d}(C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},\pi),\Pi^{\varepsilon}(\cdot,S))+{{\mathcal{W}}}_{d}(\Pi^{\varepsilon}(\cdot,S),\delta(\cdot-\theta_{*})).$

For the last two terms, by (25) and Proposition 3,

$\displaystyle{{\mathcal{W}}}_{d}(C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},\pi),\Pi^{\varepsilon}(\cdot,S))+{{\mathcal{W}}}_{d}(\Pi^{\varepsilon}(\cdot,S),\delta(\cdot-\theta_{*}))\leq
C_{f,m}\exp({{-\kappa\varepsilon\beta(t)}/{2}})+\rho(\varepsilon).$

For the first term, by Proposition 4,

$\displaystyle{{\mathcal{W}}}_{d}(D_{t}(\cdot|\theta_{0},\pi),C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},\pi))$
$\displaystyle\qquad\leq
C_{f,\theta_{0},\mu}\left[\frac{\exp({-2\varepsilon\kappa(\beta(t)-\beta(\frac{t}{2}))})}{\varepsilon}+\frac{1}{\varepsilon}\Big{(}\left|\frac{1}{\mu(t)}-\varepsilon\right|+\left|\frac{1}{\mu(\frac{t}{2})}-\varepsilon\right|\Big{)}\right]^{1/2}.$

Since $\lim_{t\to\infty}\mu(t)=\infty$, there exists $T>0$ such that
$1/\mu(\frac{T}{2})<\frac{\delta}{2\kappa}$. Let
$\varepsilon=1/\mu(\frac{t}{2})$, $t>T$, we have

$\displaystyle\exp\left({\frac{-\kappa\varepsilon\beta(t)}{2}}\right)=\exp\left({\frac{-\kappa\int_{0}^{t}\mu(s)\mathrm{d}s}{2\mu(\frac{t}{2})}}\right)\leq\exp\left({\frac{-\kappa
t}{8}}\right)$

and

$\displaystyle\frac{\exp\left({-2\varepsilon\kappa(\beta(t)-\beta(\frac{t}{2}))}\right)}{\varepsilon}$
$\displaystyle=\mu\left(\frac{t}{2}\right)\exp\left({\frac{-\kappa\int_{\frac{t}{2}}^{t}\mu(s)\mathrm{d}s}{\mu(\frac{t}{2})}}\right)$
$\displaystyle\leq\mu\left(\frac{t}{2}\right)\exp\left({\frac{-\kappa
t}{2}}\right)\leq C\exp\left({\frac{-\kappa t}{8}}\right).$

Therefore,

$\displaystyle{{\mathcal{W}}}_{d}(D_{t}(\cdot|\theta_{0},\pi),\delta(\cdot-\theta_{*}))\leq
C_{f,\theta_{0},\mu}\Big{[}\exp\left({\frac{-\kappa
t}{8}}\right)+\left(\frac{\mu(t)-\mu(\frac{t}{2})}{\mu(t)}\right)^{1/2}+\rho\left(\frac{1}{\mu\left(\frac{t}{2}\right)}\right)\Big{]}$

From Assumption 4, by the mean value theorem,

$\displaystyle\frac{\mu(t)-\mu(\frac{t}{2})}{\mu(t)}=\frac{t\mu^{\prime}(\tau_{t})}{2\mu(t)}=\frac{\tau_{t}\mu^{\prime}(\tau_{t})}{\mu(\tau_{t})}\frac{t}{2\tau_{t}}\frac{\mu(\tau_{t})}{\mu(t)}\leq\frac{\tau_{t}\mu^{\prime}(\tau_{t})}{\mu(\tau_{t})}\to
0$

where $\tau_{t}\in[\frac{t}{2},t].$ Thus (29) is obtained by taking

$A(t):=\exp\left({\frac{-\kappa
t}{8}}\right)+\left[\frac{\mu(t)-\mu(\frac{t}{2})}{\mu(t)}\right]^{1/2}+\rho\left(\frac{1}{\mu(\frac{t}{2})}\right).$

To prove (30), by the triangle inequality,

$\displaystyle{{\mathcal{W}}}_{d}(D_{t}(\cdot|\theta_{0},a),\delta(\cdot-\theta_{*}))$
$\displaystyle\leq{{\mathcal{W}}}_{d}(D_{t}(\cdot|\theta_{0},a),C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},a))$
$\displaystyle\ \ \ \ \
+{{\mathcal{W}}}_{d}(C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},a),C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},\pi))$
$\displaystyle\ \ \ \ \
+{{\mathcal{W}}}_{d}(C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},\pi),\Pi^{\varepsilon}(\cdot,S))$
$\displaystyle\ \ \ \ \
+{{\mathcal{W}}}_{d}(\Pi^{\varepsilon}(\cdot,S),\delta(\cdot-\theta_{*})).$

By Corollary 9, we have

${{\mathcal{W}}}_{d}(C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},a),C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},\pi))\leq
C_{f,m}\exp\left({-\frac{\kappa\varepsilon\beta(t)}{2}}\right).$

Notice that $C_{f,m}\exp({-\frac{\kappa t}{4}})\leq C_{f,m}A(t)$ when
$\varepsilon=1/\mu(\frac{t}{2}).$ Similar to the proof of (29), we have

$\displaystyle{{\mathcal{W}}}_{d}(D_{t}(\cdot|\theta_{0},a),C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},a))+{{\mathcal{W}}}_{d}(C^{\varepsilon}_{\varepsilon\beta(t)}(\cdot|\theta_{0},\pi),\Pi^{\varepsilon}(\cdot,S))+$
$\displaystyle{{\mathcal{W}}}_{d}(\Pi^{\varepsilon}(\cdot,S),\delta(\cdot-\theta_{*}))$
$\displaystyle\qquad\qquad\qquad\leq C_{f,\theta_{0},\mu}A(t),$

which completes the proof.

Thus, we have shown that the distribution of $(\xi_{t})_{t\geq 0}$ converges
in Wasserstein distance to the Dirac measure concentrated in the minimizer
$\theta_{*}$ of $\Phi$. This result is independent of whether we initialize
the index process $(V_{\beta(t)})_{t\geq 0}$ with its stationary measure or
with any deterministic value.

##  5 From continuous dynamics to practical optimization.

So far, we have discussed the stochastic gradient process as a continuous-time
coupling of an ODE and a stochastic process. In order to apply the stochastic
gradient process in practice, we need to discretize ODE and stochastic process
with appropriate time-stepping schemes. That means, for a given increasing
sequence $(t(k))_{k=0}^{\infty}$, with $t(0):=0$ and
$\lim_{k\rightarrow\infty}t(k)=\infty,$ we seek discrete-time stochastic
processes $(\widehat{V}_{k},\widehat{\theta}_{k})_{k=0}^{\infty}$, such that
$(\widehat{V}_{k},\widehat{\theta}_{k})_{k=0}^{\infty}\approx(V_{t(k)},\theta_{t(k)}^{\varepsilon})_{k=0}^{\infty}$
and analogous discretizations for $(V_{\beta(t)},\xi_{t})_{t\geq 0}$.

In the following, we propose and discuss time stepping strategies and the
algorithms arising from them. We discuss the index process and gradient flow
separately, which we consider sufficient as the coupling is only one-sided.

###  5.1 Discretization of the index process

We have defined the stochastic gradient process for a huge range of potential
index processes $(V_{t})_{t\geq 0}$. The discretization of such processes has
been the topic of several works, see, e.g., Gillespie (1977); Lord et al.
(2014). In the following, we focus on one case and refer to those previous
works for other settings and details.

Indeed, we study the setting $S:=[-1,1]$ and $\pi:=\mathrm{Unif}[-1,1]$ and
discuss the discretization of $(V_{t})_{t\geq 0}$ as a Markov pure jump
process and as a reflected Brownian motion.

#### Markov pure jump process.

A suitable Markov pure jump process is a piecewise constant càdlàg process
$(V_{t})_{t\geq 0}$ with Markov transition kernel

$\displaystyle\mathbb{P}_{x}(V_{t}\in\cdot)=\exp(-\lambda
t)\delta(\cdot-x)+(1-\exp(-\lambda t))\mathrm{Unif}[-1,1]\qquad\qquad(t\geq
0),$

where $\lambda>0$ is a rate parameter. We can now discretize the process
$(V_{t})_{t\geq 0}$ just through sampling from this Markov kernel for our
discrete time points. We describe this in Algorithm 1.

Algorithm 1 Discretized Markov pure jump process

1: initialize $\widehat{V}_{0}$, $\lambda>0$, and a sequence of points
$(t(k))_{k=0}^{\infty}$

2: for $k=1,2,\ldots$ do

3: sample $U\sim\mathrm{Unif}[0,1]$

4: if $U\leq\exp(-\lambda(t(k)-t(k-1))$ then

5: $\widehat{V}_{k}\leftarrow\widehat{V}_{k-1}$ {process stays at its current
position}

6: else

7: sample $\widehat{V}_{k}\sim\mathrm{Unif}[-1,1]$ {process jumps to a new
position}

8: end if

9: end for

10: return $(\widehat{V}_{k})_{k=0}^{\infty}$

#### Reflected Brownian motion.

We have defined the reflected Brownian Motion on a non-empty compact interval
through the Skorohod problem in Subsection 2.1.

Let $\sigma>0$ and $(W_{t})_{t\geq 0}$ be a standard Brownian motion. Probably
the easiest way to sample a reflected Brownian motion is by discretizing the
rescaled Brownian motion $(\sigma\cdot W_{t})_{t\geq 0}$ using the
Euler–Maruyama scheme and projecting back to $S$, whenever the sequence leaves
$S$. This scheme has been studied by Pettersson (1995). We describe the full
scheme in Algorithm 2.

Pettersson (1995) shows that this scheme converges at a rather slow rate. As
we usually assume that the domain on which we move is rather low-dimensional
and the sampling is rather cheap, we can afford small discretization stepsizes
$t(k)-t(k-1)$, for $k\in\mathbb{N}$. Thus, the slow rate of convergence is
manageable. Other schemes for the discretization of reflected Brownian motions
have been discussed by, e.g., Blanchet and Murthy (2018); Liu (1995).

Algorithm 2 Discretized Reflected Brownian motion on $S$

1: initialize $\widehat{V}_{0}$, $\sigma>0$, a sequence of points
$(t(k))_{k=0}^{\infty}$, and the projection operator $\mathrm{proj}_{S}$
mapping onto $S$

2: for $k=1,2,\ldots$ do

3: $V^{\prime}\leftarrow V_{k-1}+\sigma\sqrt{t(k)-t(k-1)}\psi$,
$\psi\sim\mathrm{N}(0,1^{2})$ {Euler-Maruyama update}

4: if $V^{\prime}\not\in S$ then

5: $\widehat{V}_{k}\leftarrow\mathrm{proj}_{S}V^{\prime}$ {project back}

6: else

7: $\widehat{V}_{k}\leftarrow V^{\prime}$ {accept Euler–Maruyama update}

8: end if

9: end for

10: return $(\widehat{V}_{k})_{k=0}^{\infty}$

###  5.2 Discretization of the gradient flow

We now briefly discuss the discretization of the gradient flow in the
stochastic gradient process. Based on these ideas, we will conduct numerical
experiments in Section 6.

#### Stochastic gradient descent

In stochastic gradient descent, the gradient flow is discretized with a
forward Euler method. This method leads to an accurate discretization of the
respective gradient flow if the stepsize/learning rates are sufficiently
small. In the presence of rather large stepsizes and stiff vector fields,
however, the forward Euler method may be inaccurate and unstable, see, e.g.,
Quarteroni et al. (2007).

#### Stability

Several ideas have been proposed to mitigate this problem. The stochastic
proximal point method, for instance, uses the backward Euler method to
discretize the gradient flow; see Bianchi (2015). Unfortunately, such implicit
ODE integrators require us to invert a possibly highly non-linear and complex
vector field. In convex stochastic optimization this inversion can be replaced
by evaluating a proximal operator. For strongly convex optimization, on the
other hand, Eftekhari et al. (2021) proposes stable explicit methods.

#### Efficient optimizers

Plenty of highly efficient methods for stochastic optimization methods are
nowadays available, especially in machine learning. Those have often been
proposed without necessarily thinking of the stable and accurate
discretization of a gradient flow: such are adaptive methods Kingma and Ba
(2015), variance reduced methods (e.g., Defazio et al. (2014)), or momentum
methods (e.g., Kovachki and Stuart (2021) for an overview), which have been
shown in multiple works to be highly efficient; partially also in non-convex
optimization. We could understand those methods also as certain
discretizations of the gradient flow. Thus, we may also consider the
combination of a feasible index process $(V_{t})_{t\geq 0}$ with the discrete
dynamical system in, e.g., the Adam method (Kingma and Ba (2015)).

##  6 Applications

We now study two fields of application of the stochastic gradient process for
continous data. In the first example, we consider regularized polynomial
regression with noisy functional data. In this case, we can easily show that
the necessary assumptions for our analysis hold. Thus, we use it to illustrate
our analytical results and especially to learn about the implicit
regularization that is put in place due to different index proccesses.

In the second example, we study so-called physics-informed neural networks. In
these continuous-data machine learning problems, a deep neural network is used
to approximate the solution of a partial differential equation. The associated
optimization problem is usually non-convex. Our analysis does not hold in this
case: We study it to get more insights in the behavior of the stochastic
gradient process in state-of-the-art deep learning problems.

Figure 3: True function $\Theta$ (red) and noisy observation $g$ (grey) in the
polynomial regression example.

###  6.1 Polynomial regression with functional data

We begin with a simple polynomial regression problem with noisy functional
data. We observe the function $g:[-1,1]\rightarrow\mathbb{R}$ which is given
through

$g(y):=\Theta(y)+\Xi(y)\qquad(y\in[-1,1]),$

where $\Theta:[-1,1]\rightarrow\mathbb{R}$ is a smooth function and $\Xi$ is a
Gaussian process with highly oscillating, continuous realizations. We aim at
identifying the unknown function $\Theta$ subject to the observational noise
$\Xi$. Here, we represent the function $\Theta$ on a basis consisting of a
finite number of Legendre polynomials on $[-1,1]$. We denote this basis of
Legendre polynomials by $(\ell_{k})_{k=1}^{K}$. To estimate the prefactors of
the polynomials, we minimize the potential

$\Phi(\theta):=\frac{1}{2}\int_{[-1,1]}\left(g(y)-\sum_{k=1}^{K}\theta_{k}\ell_{k}(y)\right)^{2}\mathrm{d}y+\frac{\alpha}{2}\|\theta\|^{2}_{2}\qquad\qquad(\theta\in
X),$ (31)

where $\alpha>0$ is a regularization parameter. This can be understood as a
maximum-a-posteriori estimation of the unknown $\theta$ with Gaussian prior
under the (misspecified) assumption that the data is perturbed with Gaussian
white noise. We employ the following associated subsampled potentials:

$f(\theta,y):=\frac{1}{2}\left(g(y)-\sum_{k=1}^{K}\theta_{k}\ell_{k}(y)\right)^{2}+\frac{\alpha}{2}\|\theta\|^{2}_{2}\qquad\qquad(\theta\in
X,y\in[-1,1]).$ (32)

Those subsampled potentials satisfy the strong convexity assumption, i.e.,
Assumption 3.

Figure 4: Estimation results of the polynomial regression problem using
stochastic gradient descent with constant learning rate $\eta=0.1$ (top row)
and a version of stochastic gradient descent that uses the implicit midpoint
rule (bottom row). The figures depict the mean over 100 runs (black solid
line), mean $\pm$ standard deviation (black dotted line). Left column:
trajectory of the rel_err over time; centre column: comparison of $\Theta$
(solid red line) and estimated polynomial; right column: estimation error in
terms of abs_err.

#### Setup

In particular, we have produced artificial data $g$, by setting
$\Theta:=\sin(\pi\cdot)$ and choosing

$\Xi(x)=\sum_{j=1}^{200}\frac{10}{1000+(\pi j)^{3/2}}\sin(2\pi
j(x-0.5))\Xi_{j}\qquad\qquad(x\in[-1,1])$

and i.i.d. random variables $\Xi_{1},\ldots,\Xi_{200}\sim\mathrm{N}(0,1^{2})$.
Note that $\Xi$ is a Gaussian random field given through the truncated
Karhunen-Loève expansion of a covariance operator that is related to the
Matérn family, see, e.g., Lindgren et al. (2011).

We show $\Theta$ and $g$ in Figure 3. For our estimation, we set
$\alpha:=10^{-4}$ and use the $K=9$ Legendre polynomials with degrees
$0,\ldots,8$. We employ the stochastic gradient process with constant learning
rate, using either a reflected diffusion process or a pure Markov jump process
for the index process $(V_{t})_{t\geq 0}$. We discretize the gradient flow
using the implicit midpoint rule: an ODE $z^{\prime}=q(z),z(0)=z_{0}$ is then
discretized with stepsize $h>0$ by successively solving the implicit formula

$z_{k}=z_{k-1}+\frac{h}{2}q(z_{k})+\frac{h}{2}q(z_{k-1})\qquad(k\in\mathbb{N}).$

In our experiments, we choose $h=0.1$. We use Algorithms 1 and 2 to discretize
the index processes with constant stepsize $t(\cdot)-t(\cdot-1)=10^{-2}.$ We
perform $J:=100$ repeated runs for each of the considered settings for
$N:=5\cdot 10^{4}$ time steps and thus, obtain a family of trajectories
$(\theta^{(j,n)})_{n=1,\ldots,N,j=1,\ldots,J}$. In each case, we choose the
initial values $V(0):=0$ and the $\theta^{(j,0)}:=(0.5,\ldots,0.5).$

We study the distance of the estimated polynomial to the true function
$\Theta$ by the relative error:

$\mathrm{rel\\_err}_{n,j}:=\frac{\sum_{l=1}^{L}\left(\Theta(x_{l})-\sum_{k=1}^{K}\theta_{k}^{(j,n)}\ell_{k}(x_{l})\right)^{2}}{\sum_{l^{\prime}=1}^{L}\Theta(x_{l^{\prime}})^{2}},$

for trajectory $j\in\\{1,\ldots,J\\}$ and time step $n\in\\{1,\ldots,N\\}$.
Here $(x_{l})_{l=1}^{L}$ are $L:=10^{3}$ equispaced points in $[-1,1]$.
Moreover, we compare the estimated polynomial to the true function $\Theta$ by

$\mathrm{abs\\_err}_{j,x}:=\left\lvert\Theta(x)-\sum_{k=1}^{K}\theta_{k}^{(j,N)}\ell_{k}(x)\right\rvert$

for trajectory $j\in\\{1,\ldots,J\\}$ at position $x\in[-1,1].$ In each case,
we study mean and standard deviation (StD) computed over the $100$ runs.

Figure 5: Estimation results of the polynomial regression problem using the
stochastic gradient process with reflected Brownian motion process with
$\sigma=5$ (top row), $\sigma=0.5$ (centre row), and $\sigma=0.05$ (bottom
row). The figures depict the mean over 100 runs (black solid line), mean $\pm$
standard deviation (black dotted line). Left column: trajectory of the rel_err
over time; centre column: comparison of $\Theta$ (solid red line) and
estimated polynomial; right column: estimation error in terms of abs_err.

#### Results and discussion

For the polynomial regression problem we now study:

  * •

stochastic gradient descent, as given in (2), with constant learning rate
$\eta_{(\cdot)}=h=0.1$ (Figure 4 top row),

  * •

stochastic gradient descent algorithm, for which the forward Euler update is
replaced by an implicit midpoint rule update, with constant learning rate
$\eta_{(\cdot)}=h=0.1$ (Figure 4 bottom row),

  * •

the stochastic gradient process with reflected Brownian motion as an index
process with standard deviation $\sigma\in\\{5,0.5,0.05\\}$ (Figure 5), and

  * •

the stochastic gradient process with Markov pure jump process as an index
process with rate parameter $\lambda\in\\{10,1,0.1,0.01\\}$ (Figure 6).

In addition to those plots, we give means and standard deviations of the
relative errors at the terminal state of the iterations in Table 1. To compare
the convergence behavior of the different methods, we plot the rel_err within
the first 2000 discrete time steps in Figure 7.

Figure 6: Estimation results of the polynomial regression problem using the stochastic gradient process with pure jump index process with $\lambda=10$ (first row), $\lambda=1$ (second row), $\lambda=0.1$ (third row), and $\lambda=0.01$ (fourth row). The figures depict the mean over 100 runs (black solid line), mean $\pm$ standard deviation (black dotted line). Left column: trajectory of the rel_err over time; centre column: comparison of $\Theta$ (solid red line) and estimated polynomial; right column: estimation error in terms of abs_err. Method | Parameters | Mean of $\mathrm{rel\\_err}_{N,(\cdot)}$ |  $\pm$ StD
---|---|---|---
SGD | $\eta_{(\cdot)}=0.1$ | $1.844\cdot 10^{-2}$ | $\pm 4.012\cdot 10^{-3}$
SGD implicit | $\eta_{(\cdot)}=0.1$ | $1.719\cdot 10^{-2}$ | $\pm 3.939\cdot 10^{-3}$
SGPC with reflected diffusion index process | $\sigma=5$ | $1.586\cdot 10^{-2}$ | $\pm 4.038\cdot 10^{-3}$
$\sigma=0.5$ | $1.587\cdot 10^{-2}$ | $\pm 2.979\cdot 10^{-3}$
$\sigma=0.05$ | $4.637\cdot 10^{-2}$ | $\pm 8.776\cdot 10^{-2}$
SGPC with Markov pure jump index process | $\lambda=10$ | $2.100\cdot 10^{-2}$ | $\pm 6.049\cdot 10^{-3}$
$\lambda=1$ | $3.427\cdot 10^{-2}$ | $\pm 1.105\cdot 10^{-2}$
$\lambda=0.1$ | $3.866\cdot 10^{-2}$ | $\pm 1.142\cdot 10^{-2}$
$\lambda=0.01$ | $3.178\cdot 10^{-1}$ | $\pm 2.124\cdot 10^{-1}$
Table 1: Accuracy of the estimation in the polynomial regression model. Mean
and standard deviation of the relative error of the methods at the final point
of their trajectory. In particular, sample mean and sample standard deviation
of $j\mapsto\mathrm{rel\\_err}_{N,j}$, with $N=5\cdot 10^{4}$, computed over
$100$ independent runs.

We learn several things from these results. Unsurprisingly, the index
processes with a strong autocorrelation $(\lambda=0.01,\sigma=0.05)$ lead to
larger errors in the reconstruction: the processes move too slowly to capture
the index spaces appropriately. In the other cases, we can assume that the
processes have reached their stationary regime. Thus, in the figures and
table, we should learn about the implicit regularization that is implicated by
the different subsampling schemes, see Ali et al. (2020); Smith et al. (2021).
We especially see that the mean errors are reduced as $\sigma$ respectively
$\lambda$ increases, which illustrates the approximation of the full gradient
flow as shown in Theorem 3. Although, we should note that we compute the error
to the truth $\Theta$, which is likely not the true minimizer of the full
optimization problem 31.

Figure 7: Comparison of the mean rel_err of the stochastic methods for time
$t\leq 200$.

It appears that the stochastic gradient processes with reflected diffusion
index process $\sigma\in\\{0.5,5\\}$ returns the best results. Looking at the
error plots in the right column of Figure 5, we see that SGPC especially
outperforms the other algorithms close to the boundary. For $\sigma=5$ this
could be seen as a numerical artefact due to the time step $t-t(\cdot-1)$
being too large. This though is likely not the case for $\sigma=0.5$, where we
see a similar effect, albeit a bit weaker.

In the convergence plot, Figure 7, we see for different methods different
speeds of convergence to their respective stationary regime. Those speeds
again depend on the autocorrelation of the processes. Interestingly, the SGPC
with reflected diffusion index process and $\sigma=5$ appears to be the best
of the algorithms.

###  6.2 Solving partial differential equations using neural networks (NN)

Partial differential equations (PDEs) are used in science and engineering to
model systems and processes, such as: turbulent flow, biological growth, or
elasticity. Due to the implicit nature of a PDE and its complexity, the model
they represent usually needs to be approximated (‘solved’) numerically. Finite
differences, elements, and volumes have been the state of the art for solving
PDEs for the last decades. Recently, deep learning approaches have gained
popularity for the approximation of PDE solutions. Here, deep learning is
particularly successful in high-dimensional settings, where classical methods
suffer from the curse of dimensionality. See for example Raissi et al. (2019);
Lu et al. (2021) for physics-informed neural networks (PINN). Integrated
PyTorch-based packages are available for example see Chen et al. (2020); Pedro
et al. (2019). More recently, see Li et al. (2021) for a state-of-the-art
performance based on the Fourier neural operator.

Physics-informed neural networks are a very natural field of application of
deep learning with continuous data. Below we introduce PINNs, the associated
continuous-data optimization problem, and the state-of-the-art in the training
of PINNs. Then we consider a particular PDE, showcase the applicability of
SGP, and compare its performance with the standard SGD-type algorithm.

The basic idea of PINNs consists in representing the PDE solution by a deep
neural network where the parameters of the network are chosen such that the
PDE is optimally satisfied. Thus, the problem is reduced to an optimization
problem with the loss function formulated from differential equations,
boundary conditions, and initial conditions. More precisely, for PDE problems
of Dirichlet type, we aim to solve a system of equations of type

$\left\\{\begin{array}[]{rlll}\mathcal{L}(u(t,x))&=&s(t,x)&(t\in[0,\infty),\
x\in D)\\\ u(0,x)&=&u_{0}(x)&(x\in D)\\\ u(t,x)&=&b(t,x)&(t\in[0,\infty),\
x\in D)\end{array}\right.$ (33)

where $D\subset\mathbb{R}^{d}$ is an open, connected, and bounded set and
$\mathcal{L}$ is a differential operator defined on a function space $V$ (e.g.
$H^{1}(D)$). The unknown is $u:\bar{D}\to\mathbb{R}^{n}$. Functions $s(t,x)$,
$b(t,x)$, and $u_{0}(x)$ are given. In numerical practice, we need to replace
the infinite-dimensional space $V$ by a – in some sense – discrete
representation. Traditionally, one employs a finite-dimensional subspace of
$V$, say $\mathrm{span}\\{\psi_{1},\ldots,\psi_{K}\\}$, where
$\psi_{1},\ldots,\psi_{K}$ are basis functions in a finite element method. To
take advantage of the recent development of machine learning, one could solve
the problem on a set of deep neural networks contained in $V$, say

$\displaystyle{\Big{\\{}}\psi(\cdot;\theta):\ $
$\displaystyle\psi(x;\theta)=(W^{(K)}\sigma(\cdot)+b^{(K)})\circ\cdots\circ(W^{(1)}\sigma(x)+b^{(1)}),x\in[0,\infty)\times
D,$
$\displaystyle\theta=\left((W^{(K)},b^{(K)}),\ldots,(W^{(1)},b^{(1)})\right)\in\prod_{k=1}^{K}\left(\mathbb{R}^{n_{k}\times
n_{k-1}}\times\mathbb{R}^{n_{k}}\right)=:X{\Big{\\}}},$

where $\sigma:\mathbb{R}\rightarrow\mathbb{R}$ is an activation function,
applied component-wise, $n_{0}=d+1$ and $n_{K}=1$ to match input and output of
the PDE’s solution space, and $n_{1},\ldots,n_{K-1}$ determine the network’s
architecture.

In simpler terms, let $u(\cdot;\theta)\in V$ be the output of a feedforward
neural network (FNN) with parameters (biases/weights) denoted by $\theta\in
X$. The parameters can be learned by minimizing the mean squared error (MSE)
loss

$\displaystyle\Phi(\theta;\mathcal{L},s,u_{0},b):=$
$\displaystyle\int_{0}^{\infty}w(t)\int_{D}\left(\mathcal{L}(u(t,x;\theta))-s(x)\right)^{2}\mathrm{d}x\mathrm{d}t+\int_{\partial
D}\left(u(0,x;\theta)-u_{0}(x)\right)^{2}\mathrm{d}x$ $\displaystyle\
+\int_{0}^{\infty}w(t)\int_{\partial
D}\left(u(t,x;\theta)-b(t,x)\right)^{2}\mathrm{d}x\mathrm{d}t,$

where the first term is the $L^{2}$ norm of the PDE residual, the second term
is the $L^{2}$ norm of the residual for the initial condition, the third term
is the $L^{2}$ norm of the residual for the boundary conditions, and
$w:[0,\infty)\rightarrow[0,\infty)$ is an appropriate weight function. The FNN
then represents the solution via solving the following minimization problem

$\min_{\theta\in X}\Phi(\theta;\mathcal{L},s,u_{0},b).$ (34)

Note that in physics-informed neural networks, differential operators w.r.t.
the input $x$ and the gradient w.r.t the parameter $\theta$ are both obtained
using automatic differentiation.

#### Training of physics-informed neural networks

In practice, the optimization problem (34) is often replaced by an
optimization problem with discrete potential

$\displaystyle\widehat{\Phi}(\theta;\mathcal{L},s,u_{0},b):=$
$\displaystyle\sum_{k=1}^{K}\left(\mathcal{L}(u(t_{k},x_{k};\theta))-s(x_{k})\right)^{2}+\sum_{k^{\prime}=1}^{K^{\prime}}\left(u(0,x^{\prime}_{k^{\prime}};\theta)-u_{0}(x^{\prime}_{k^{\prime}})\right)^{2}$
$\displaystyle\
+\sum_{k^{\prime\prime}=1}^{K^{\prime\prime}}\left(u(t^{\prime\prime}_{k^{\prime\prime}},x^{\prime\prime}_{k^{\prime\prime}};\theta)-b(t_{k^{\prime\prime}}^{\prime\prime},x_{k^{\prime\prime}}^{\prime\prime})\right)^{2},$

for appropriate continuous indices

$(x_{k},t_{k})_{k=1}^{K}\in[0,\infty)^{K}\times
D^{K},(x_{k^{\prime}}^{\prime})_{k^{\prime}=1}^{K^{\prime}}\in\partial
D^{K},(x_{k^{\prime\prime}}^{\prime\prime},t_{k^{\prime\prime}}^{\prime\prime})_{k^{\prime\prime}=1}^{K^{\prime\prime}}\in[0,\infty)^{K}\times\partial
D^{K}$

that may be chosen deterministically or randomly, see for example Pedro et al.
(2019); Lu et al. (2021).

Focusing the training on a fixed set of samples can be problematic: fixing a
set of random samples might be unreliable; a reliable cover of the domain will
likely only be reached through tight meshing, which scales badly. Sirignano
and Spiliopoulos (2018) propose to use SGD on the continuous data space. They
employ the discrete dynamic in (2). Naturally, we would like to follow
Sirignano and Spiliopoulos (2018) and employ the SGP dynamic on the continuous
index set.

To train the PINNs with SGP, we again choose the reflected Brownian motion as
an index process, which we discretize with the Euler–Maruyama scheme in
Algorithm 2. In addition, we employ mini-batching to reduce the variance in
the estimator: We sample $M\in\mathbb{N}$ independent index processes
$(V_{t}^{(1)})_{t\geq 0},\ldots,(V_{t}^{(M)})_{t\geq 0}$ and then employ the
dynamical system

$\mathrm{d}\theta_{t}=-\frac{1}{M}\sum_{m=1}^{M}\nabla_{\theta}f(\theta_{t},V_{t}^{(m)})\mathrm{d}t.$

Hence, rather than optimizing with respect to a single data set, we optimize
with respect to $M$ different data sets in each iteration. While we only
briefly mention the mini-batching throughout our analysis, one can easily see
that it is fully contained in our framework.

In preliminary experiments, we noticed that the Brownian motion for the
sampling on the boundary is not very effective: possibly due to its localizing
effect. Hence, we obtain training data on the boundary by sampling uniformly,
which we consider justified as a mesh on the boundary scales more slowly as a
mesh in the interior and as the boundary behavior of the considered PDE is
rather predictable.

#### PDE and results

We now describe the partial differential equation that we aim to solve with
our PINN model. After introducing the PDE we immediately outline the PINN’s
architectures and show our estimation results. We the train networks on Google
Colab Pro using GPUs (often T4 and P100, sometimes K80). We are certain that a
more efficient PDE solution could be obtained by classical methods, e.g., the
finite element method. We do not compare the deep learning methods with
classical methods, as we are mainly interested in SGP and SGD in non-convex
continuous-data settings. Other methods that could approximate the PDE
solution are not our focus.

The PDE we study is a transport equation; which is a linear first order, time-
dependent model. One of the main advantages of studying this particular model
is that we know an analytical solution that allows us to compute a precise
test error.

######  Example 1 (1D Transport equation)

We solve the one-dimensional transport equation on the space $[0,1]$ with
periodic boundary condition:

$\left\\{\begin{array}[]{l}u_{t}+u_{x}=0,\ \ t\in[0,\infty),\ \ x\in[0,1]\\\
u(t=0)=\sin(2\pi x),\\\ u(t,0)=u(t,1).\end{array}\right.$ (35)

The neural network approximation of this PDE has already been studied by Pedro
et al. (2019), our experiments partially use the code associated to this work.
The network architecture is defined by a three-layer deep neural network with
128 neurons per layer and a Rectified Linear Unit (ReLu) activation function.
While theoretically the solution exists globally in time, we restrict $t$ to a
compact domain and w.l.o.g, we assume $t\in[0,1]$. From the interior of the
domain of time and space variables, i.e. $(0,1)\times(0,1)$, we use Algorithm
2 with $\sigma=0.5$ to sample the train set of size $3\cdot 10^{4}$ for SGPC
and SGPD and we uniformly sample $600$ points for the train set of SGD. In
addition, as a part of the train set for all three methods, we sample
uniformly $20$ and $60$ points for the initial condition and periodic boundary
condition, respectively.

The learning rate for SGD and SGPC is $0.01$. The learning rate for SGPD is
defined as

$\eta(t)=\frac{0.01}{\log(t+2)^{0.3}},$

which is chosen such that the associated $\mu:=1/\eta$ satisfies Assumption 4.
For all three methods, we use Adam (see Kingma and Ba, 2015) as the optimizer
to speed up the convergence; we use an $L^{2}$ regularizer with weight $0.1$
to avoid overfitting. Each model is trained over $600$ iterations with batch
size $50$. The training process for SGPC and SGPD contains only one epoch,
while we train $50$ epochs in the SGD case. We evaluate the models by testing
on a uniformly sampled test set of size $2\cdot 10^{3}$ and compare the
predicted values with the theoretical solution

$u(t,x)=\sin(2\pi(x-t)).$

We obtain the losses, the predicted solutions, and the test errors by
averaging over $30$ random experiments, i.e. $30$ independent runs of SGD,
SGPC, and SGPD, respectively. We give the results in Figures 8, 9, and 10.
Note that the timings are very similar for each of the algorithms, the fact
that SGPC and SGPD require us to first sample reflected Brownian motions is
negligible.

Figure 8: The plots of the loss vs iteration and its log scale for SGD, SGPC,
and SGPD. The losses are obtained by averaging over $30$ random experiments.

Figure 9: The plots of the solutions at time $t=0.1,0.5,0.9$. We evaluate the
models at $30$ uniformly sampled points. For each method, the predicted values
are taken by averaging over the predicted values from the best models (the
model that achieves the lowest training loss within the 600 iteration steps)
in $30$ random experiments. The black curve is the theoretical solution.

Figure 10: The plots of the test error at time $t=0.1,0.5,0.9$. We evaluate
the models at $2000$ uniformly sampled points. For each method, the predicted
values are taken by averaging over the predicted values from the best models
(the model that achieves the lowest training loss within the 600 iteration
steps) in $30$ random experiments. At each point $x$, the error calculated by
taking the absolute value of the difference between the predicted value and
the true solution.

From Figure 8, we notice that while SGD and SGPC behave similarly, SGPD does
converge faster. Here, Assumption 4 provides a way of designing a non-constant
learning rate in practice. On the test set, the mean squared errors for SGD,
SGPC, and SGPD are $4.5\cdot 10^{-4}$, $3.5\cdot 10^{-4}$, and $2.8\cdot
10^{-4}$. These test errors refer to the averaged model output of the 30
models from independent experiments. Combined with Figure 9 and Figure 10, we
observe that SGPC and SGPD generalize at least slightly better on the test
set. This improved generalization error might be due to the additional test
data generated by the Brownian motion, as compared to the fixed training set
used in PINNs. The combination with the reduction of the learning rate in
SGPD, appears to be especially effective.

##  7 Conclusions and outlook

In this work we have proposed and analyzed a continuous-time stochastic
gradient descent method for optimization with respect to continuous data. Our
framework is very flexible: it allows for a whole range of random sampling
patterns on the continuous data space, which is particularly useful when the
data is streamed or simulated. Our analysis shows ergodicity of the dynamical
system under convexity assumptions – converging to a stationary measure when
the learning rate is constant and to the minimizer when the learning rate
decreases. In experiments we see the suitability of the method and the effect
of different sampling patterns on its implicit regularization.

We end this work by now briefly listing some interesting problems for future
research in this area. First, we would like to learn how the SGP sampling
patterns perform in large-scale (adversarially-)robust machine learning and in
other applications we have mentioned but not studied here. Moreover, from a
both practical and analytical perspective, it would be interesting to also
consider non-compact index spaces $S$. Those appear especially in robust
optimal control and variational Bayes. Finally, we consider the following
generalization of the optimization problem (1) to be of high interest:

$\min_{\theta\in X}\int_{S}f(\theta,y)\Pi(\mathrm{d}y|\theta),$

where $\Pi$ is now a Markov kernel from $X$ to $S$. Hence, in this case the
probability distribution and the sampling pattern itself depend on the
parameter $\theta$. Optimization problems of this form appear in the optimal
control of random systems (e.g., Deqing and Sergei (2015)) and empirical Bayes
(e.g., Casella (2001)) but also in reinforcement learning (e.g., Sutton and
Barto (2018)).

Acknowledgments

JL and CBS acknowledge support from the EPSRC through grant EP/S026045/1
“PET++: Improving Localisation, Diagnosis and Quantification in Clinical and
Medical PET Imaging with Randomised Optimisation”.

## References

  * Ali et al. (2020) Alnur Ali, Edgar Dobriban, and Ryan Tibshirani.  The implicit regularization of stochastic gradient flow for least squares.  In _Proceedings of the 37th International Conference on Machine Learning (ICML 2020)_ , pages 233–244, 2020.
  * Andersen et al. (2015) Lars Nørvang Andersen, Søren Asmussen, Peter W. Glynn, and Mats Pihlsgård.  Lévy processes with two-sided reflection.  In _Lévy Matters V: Functionals of Lévy Processes_ , pages 67–182. Springer International Publishing, 2015.
  * Bianchi (2015) Pascal Bianchi.  A stochastic proximal point algorithm: convergence and application to convex optimization.  In _2015 IEEE 6th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)_ , pages 1–4, 2015.
  * Blanchet and Murthy (2018) Jose Blanchet and Karthyek Murthy.  Exact simulation of multidimensional reflected brownian motion.  _Journal of Applied Probability_ , 55(1):137–156, 2018.