# Model reduction of Brownian oscillators: quantification of errors and long-
time behaviour

Matteo Colangeli<EMAIL_ADDRESS>Department of Information
Engineering, Computer Science and Mathematics, University of L’Aquila, Italy.
Manh Hong Duong<EMAIL_ADDRESS>School of Mathematics, University of
Birmingham, UK.  Adrian Muntean<EMAIL_ADDRESS>Department of
Mathematics and Computer Science & Centre for Societal Risk Research (CSR),
Karlstad University, Sweden.

###### Abstract

A procedure for model reduction of stochastic ordinary differential equations
with additive noise was recently introduced in [CDM22], based on the Invariant
Manifold method and on the Fluctuation-Dissipation relation. A general
question thus arises as to whether one can rigorously quantify the error
entailed by the use of the reduced dynamics in place of the original one. In
this work we provide explicit formulae and estimates of the error in terms of
the Wasserstein distance, both in the presence or in the absence of a sharp
time-scale separation between the variables to be retained or eliminated from
the description, as well as in the long-time behaviour.

Keywords: Model reduction, Wasserstein distance, error estimates, coupled
Brownian oscillators, invariant manifold, Fluctuation-Dissipation relation.

##  1 Introduction

The notion of scale separation is largely invoked in multiscale modelling and
homogeneization methods (including model reduction and operator splitting
techniques) [GKS04, PS08], and has also found far-reaching applications in
different areas of science and engineering, e.g. in climate dynamics [GL20],
biochemical systems [SS17], chemical reaction networks [KK13], smoldering
combustion [IOMF14], and so on. A neat illustration of this notion can be
traced in the preface of Haken’s seminal book on Synergetics [Hak04], where
the author writes: “In large classes of systems that are originally described
by many variables, the behavior of a system is described and determined by
only few variables, the order parameters. They fix the behavior of the
individual parts via the slaving principle”. A physical rationale behind the
slaving principle amounts to the assumption of decomposition of motions: there
exists a short time-scale during which the slow variable does not change
significantly, while the fast variable rapidly settles on a value determined
by the slow one. The evolution of the latter, in turn, takes place on a much
longer scale. A specific form of such principle is realized through the method
of adiabatic elimination of fast variables, which underlies the derivation of
the Smoluchowski equation from the underdamped Langevin equation. A sharp
distinction between slow and fast variables is also a prerequisite for
application of the Mori-Zwanzing method [Zwa01] in the derivation of reduced
equations from higher dimensional stochastic dynamics, where the Markovian
structure of the original process is preserved in the reduced description by
stipulating a perfect time-scale separation. The same guiding principle
underpins, in kinetic theory, the Grad moment method [Gra49, CKK07], and has
also been exploited in the derivation of linear hydrodynamics from the
Boltzmann equation using the framework of the Invariant Manifold [GK05,
CKO09]. The latter method has also been exploited in [CDM22] to characterize
the deterministic component of the contracted description in a system of two
coupled (underdamped) Brownian harmonic oscillators. The structure of the
noise term of the Markovian reduced dynamics, in turn, was determined via the
Fluctuation-Dissipation relation. A general question, then, concerns the
derivation of a quantitative estimate of the error stemming from the use of
the reduced dynamics in place of the original one. A first attempt, in this
direction, was proposed in [CM22], and it was based on the study of the
equilibrium correlation functions in the reduced and in the original
processes. A uniform-in-time type of convergence of the correlations evaluated
in the two processes was proven to hold in the so-called overdamped limit,
where the friction parameter diverges.

In this work we take a step further, and compute explicitly the Wasserstein
distance between the laws of the original and reduced processes. This paves
the way to explicitly quantify the error inherent to the contracted
description. We focus on two classical models thoroughly studied in
statistical physics and molecular dynamics, namely the underdamped Brownian
harmonic oscillator and a system of two coupled overdamped Brownian harmonic
oscillators. In the more traditional approach based on the slow-fast
decomposition of motions, a reduced description can be achieved by passing the
parameter to a certain limit, thus establishing a perfect time-scale
separation, see e.g. [Zwa01, GLCG21]. In the present work, instead, we derive
the reduced dynamics in a regime characterized by a finite time-scale
separation, which is controlled, in the two considered models, by either the
friction parameter or the coupling parameter. We show that the reduced and
original dynamics are exponentially close at any time, and they coincide if we
pass the parameter to the corresponding limit. We also prove that the two
dynamics have the same equilibrium measure and, furthermore, they
exponentially converge to the equilibrium measure with the same rate. This
notable property is a direct consequence of the proposed reduction scheme, in
particular of the selection of solutions to the invariance equation obtained
from the Invariant Manifold method. As a consequence of this, the spectrum of
the reduced drift matrix is a subset of the spectrum of the original drift
matrix. The models and precise statements of the results are presented in
Section 3 and Section 4.

The work is structured as follows. In Sec. 2 we review the definition of the
Wasserstein distance between two probability measures and introduce the basic
notation used throughout the manuscript. In Sec. 3 we compute our error
estimate based on the Wasserstein distance for a Brownian harmonic oscillator,
for which the laws of the original and the contracted descriptions are
analytically known. In Sec. 4 we apply our method to a slightly more involved
model, constituted by a pair of coupled overdamped Brownian harmonic
oscillators. Conclusions and a final outlook are finally drawn in Sec. 5.

##  2 Preliminaries

In this Section we introduce the Wasserstein distance between two probability
measures and also fix the notation used throughout the manuscript.

###  2.1 Wasserstein distance

In this section we recall the definition of the Wasserstein distance between
two probability measures and its explicit formula when the two probability
measures are Gaussian distributions. The Wasserstein metric plays an central
role in many research fields such as optimal transport, partial differential
equations and data science. For a detailed account of the topics, we refer the
reader to Villani’s monograph [Vil03].

Let $P_{2}(\mathbb{R}^{d})$ be the space of probability measures $\mu$ on
$\mathbb{R}^{d}$ with finite second moment, namely

$\int_{\mathbb{R}^{d}}|x|^{2}\mu(dx)<\infty.$

Let $\mu$ and $\nu$ be two probability measures belonging to
$P_{2}(\mathbb{R}^{d})$. The $L^{2}$-Wasserstein distance, $W_{2}(\mu,\nu)$,
between $\mu$ and $\nu$ is defined via

$W^{2}_{2}(\mu,\nu):=\inf_{\gamma\in\Gamma(\mu,\nu)}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}|x-y|^{2}\,\gamma(dx,dy),$
(1)

where $\Gamma(\mu,\nu)$ denotes the set of all couplings between $\mu$ and
$\nu$, i.e., the set of all probability measures on
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ having $\mu$ and $\nu$ as the first and
the second marginals respectively. More precisely,

$\Gamma(\mu,\nu):=\\{\gamma\in
P(\mathbb{R}^{d}\times\mathbb{R}^{d}):\gamma(A\times\mathbb{R}^{d})=\mu(A)~{}\text{and}~{}\gamma(\mathbb{R}^{d}\times
A)=\nu(A)\\},$

for all Borel measurable sets $A\subset\mathbb{R}^{d}$.

In particular, the Wasserstein distance between two Gaussian measures can be
computed explicitly in terms of the means and covariance matrices [GS84], see
also e.g., [Tak12]

$W_{2}(\mathcal{N}(u,U),\mathcal{N}(v,V))^{2}=|u-v|^{2}+\mathrm{tr}U+\mathrm{tr}V-2\mathrm{tr}\sqrt{V^{\frac{1}{2}}UV^{\frac{1}{2}}},$
(2)

where $u,v$ are the means and $U,V$ are the covariance matrices. In a one
dimensional space, the above formula reduces to

$W_{2}(\mathcal{N}(u_{1},\sigma_{1}^{2}),\mathcal{N}(u_{2},\sigma_{2}^{2})^{2}=(u_{1}-u_{2})^{2}+(\sigma_{1}-\sigma_{2})^{2}.$
(3)

###  2.2 Linear drift-diffusion equations

We recall here a well-known result concerning the explicit solution of a
general linear drift-diffusion where the initial data is a Gaussian
distribution. In the subsequent sections, we will apply this result to our
models of (coupled) Brownian oscillators.

To set the stage, we consider the following general linear drift-diffusion
equation

$\partial_{t}\rho=-\mathop{\mathrm{div}}\nolimits(Cx\rho)+\mathop{\mathrm{div}}\nolimits(D\nabla\rho),\quad\rho(0)=\rho_{0}.$
(4)

In the above equation, the unknown is a probability measure $\rho=\rho(t,x)$
with $(t,x)\in(0,\infty)\times\mathbb{R}^{d}$; $C$ and $D$ are two constant
matrices of order $d$ representing the drift and diffusion matrices; the
initial data $\rho_{0}$ is a probability measure on $\mathbb{R}^{d}$.

The following lemma provides the explicit formula for the solution of (4) when
the initial data is a Gaussian distribution, see for instance [GP18].

######  Lemma 2.1.

Suppose the initial data is a Gaussian,
$\rho_{0}\sim\mathcal{N}(\mu(0),\Sigma(0))$, then the solution to (4) is given
by

$\rho(t,x)=\frac{1}{\sqrt{(2\pi)^{d}\det\Sigma(t)}}\exp\Big{[}-\frac{1}{2}(x-\mu(t))^{T}\Sigma^{-1}(t)(x-\mu(t))\Big{]}$
(5)

where $\mu(t)$ and $\Sigma(t)$ are given by

$\mu(t):=e^{tC}\mu(0),\quad\Sigma(t):=e^{tC}\Sigma(0)e^{tC^{T}}+2\int_{0}^{t}e^{sC}De^{sC^{T}}\,ds.$
(6)

Under suitable conditions on $C$ and $K$, we have $\mu(t)\rightarrow 0$ and
$\Sigma(t)\rightarrow\Sigma_{\infty}$ where

$\Sigma_{\infty}:=2\int_{0}^{\infty}e^{sC}De^{sC^{T}}\,ds.$

Note that $\Sigma_{\infty}$ satisfies the so-called Lyapunov equation

$2D=C\Sigma_{\infty}+\Sigma_{\infty}C^{T}.$

###  2.3 Exponential of a $2\times 2$ matrix

Lemma 2.1 provides the explicit form of the unique solution to the linear
drift-diffusion equation (4) when the initial data is a Gaussian. However, in
general the formula (6) is analytically hard to compute since it involves
exponential of matrices. The following lemma provides an explicit formula for
the exponential of a $2\times 2$ matrix, which will be used in the subsequent
analysis.

######  Lemma 2.2.

Let $a,b,c,d\in\mathbb{R}$ be taken arbitrarily with
$a^{2}+b^{2}+c^{2}+d^{2}>0$. The following identity holds

$\exp\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\frac{1}{\Delta}\begin{pmatrix}m_{11}&m_{12}\\\
m_{21}&m_{22}\end{pmatrix},$ (7)

where $\Delta:=\sqrt{(a-d)^{2}+4bc}$ and

$\displaystyle m_{11}$
$\displaystyle:=e^{(a+d)/2}\Big{[}\Delta\cosh{\frac{1}{2}\Delta}+(a-d)\sinh{\frac{1}{2}\Delta}\Big{]},$
$\displaystyle m_{12}$ $\displaystyle:=2be^{(a+d)/2}\sinh{\frac{1}{2}\Delta},$
$\displaystyle m_{21}$ $\displaystyle:=2ce^{(a+d)/2}\sinh{\frac{1}{2}\Delta},$
$\displaystyle m_{22}$
$\displaystyle:=e^{(a+d)/2}\Big{[}\Delta\cosh{\frac{1}{2}\Delta}+(d-a)\sinh{\frac{1}{2}\Delta}\Big{]}.$

###### Proof.

We refer the reader to [BS93] for a justification of the formula (7). ∎

##  3 Model reduction of a Brownian oscillator

To start off the discussion, we begin with the investigation of a  simple
model of an underdamped Brownian oscillator considered in [CM22], which is
amenable to an explicit analytical solution. The original dynamics reads as
follows:

$\displaystyle dx(t)$ $\displaystyle=v(t)\,dt$ $\displaystyle dv(t)$
$\displaystyle=-\omega^{2}x(t)\,dt-\gamma
v(t)\,dt+\sqrt{2\gamma\beta^{-1}}\,dW(t),$ $\displaystyle(x(0),v(0))$
$\displaystyle=(x_{0},v_{0})$

Exploiting the Invariant Manifold method and the Fluctuation-Dissipation
relation (for a short summary of the method, see Section 4 below, where the
same reduction procedure is applied to a system of coupled overdamped Brownian
harmonic oscillators), the reduced dynamics attains the form:

$d\bar{x}(t)=-\alpha\bar{x}(t)\,dt+\sqrt{2D_{r}}\,dW(t),\quad\bar{x}(0)=x_{0},$

where

$\alpha=\frac{\gamma-\sqrt{\gamma^{2}-4\omega^{2}}}{2},\quad
D_{r}=\frac{\alpha}{\omega^{2}\beta}.$

The reader is referred to [CM22] to see the details of the calculations. The
main result of this section is the following theorem.

######  Theorem 3.1.

  1. (i)

(exact solutions of the original and the reduced dynamics) $\mu_{t}$ and
$\bar{\mu}_{t}$ are Gaussian measures

$\mu_{t}=\mathcal{N}(m(t),\sigma(t)),\quad\bar{\mu}_{t}=\mathcal{N}(\bar{m}_{t},\bar{\sigma}(t)),$
(8)

where

$\displaystyle m(t)$
$\displaystyle=\frac{\lambda_{1}e^{-\lambda_{2}t}-\lambda_{2}e^{-\lambda_{1}t}}{\lambda_{1}-\lambda_{2}}x_{0}+\frac{e^{-\lambda_{2}t}-e^{-\lambda_{1}t}}{\lambda_{1}-\lambda_{2}}v_{0},$
$\displaystyle\sigma(t)$
$\displaystyle=\frac{\gamma\beta^{-1}}{(\lambda_{1}-\lambda_{2})^{2}}\Big{[}\frac{\lambda_{1}+\lambda_{2}}{\lambda_{1}\lambda_{2}}+\frac{4}{\lambda_{1}+\lambda_{2}}(e^{-(\lambda_{1}+\lambda_{2})t}-1)-\frac{1}{\lambda_{1}}e^{-2\lambda_{1}t}-\frac{1}{\lambda_{2}}e^{-2\lambda_{2}t}\Big{]},$
$\displaystyle\bar{m}(t)$ $\displaystyle=e^{-\lambda_{2}t}\bar{x}_{0},$
$\displaystyle\bar{\sigma}(t)$
$\displaystyle=\frac{1}{\omega^{2}\beta}(1-e^{-2\lambda_{2}t})$

where

$\lambda_{1}=\frac{\gamma+\sqrt{\gamma^{2}-4\omega^{2}}}{2},\quad\lambda_{2}=\frac{\gamma-\sqrt{\gamma^{2}-4\omega^{2}}}{2}=\frac{2\omega^{2}}{\gamma+\sqrt{\gamma^{2}-4\omega^{2}}}.$
(9)

  2. (ii)

(Exact Wasserstein distance between the laws of the original and reduced
dynamics) The Wasserstein distance between $\mu_{t}$ and $\bar{\mu}_{t}$ can
be computed explicitly via

$W_{2}^{2}(\mu_{t},\bar{\mu}_{t})=(m(t)-\bar{m}(t))^{2}+\Big{(}\sqrt{\sigma_{xx}(t)}-\sqrt{\bar{\sigma}(t)}\Big{)}^{2}.$
(10)

  3. (iii)

(explicit rate of convergence in the high-friction limit) It holds that

$W_{2}^{2}(\mu_{t},\bar{\mu}_{t})\leq\frac{4}{\gamma^{2}-4\omega^{2}}\Big{[}(\omega|x_{0}|+|v_{0}|)^{2}+\frac{4}{\beta}\Big{]}\quad\forall
t>0.$ (11)

As a consequence,

$\lim_{\gamma\rightarrow+\infty}W_{2}^{2}(\mu_{t},\bar{\mu}_{t})=0.$

Note that (11) is a much stronger statement providing an explicit rate of
convergence.

  4. (iv)

(Common rates of convergence to equilibrium) There exists a constant $C>0$,
which can be found explicitly, such that

$W_{2}(\mu_{t},\mu_{\infty}),~{}W_{2}(\bar{\mu}_{t},\bar{\mu}_{\infty})\leq
Ce^{-\lambda_{2}t},$

where

$\mu_{\infty}=\bar{\mu}_{\infty}=\mathcal{N}\Big{(}0,\frac{1}{\beta\omega^{2}}\Big{)}.$

This result shows that the original dynamics and the reduced one not only
share the same equilibrium, they have the same rates of convergence to
equilibrium in the Wasserstein distance.

  5. (v)

(long-time behaviour) It holds that

$W_{2}^{2}(\mu_{t},\bar{\mu}_{t})\leq\Big{[}\frac{\omega|x_{0}|+|v_{0}|}{\sqrt{\gamma^{2}-4\omega^{2}}}+\frac{10}{\beta(\gamma^{2}-4\omega^{2})}\Big{]}e^{-\lambda_{2}t}.$
(12)

As a consequence of this, we also have

$\lim_{t\rightarrow+\infty}W_{2}(\mu_{t},\bar{\mu}_{t})=0,$

which is already obtained in the previous part. Estimate (12) is a stronger
statement, showing that the two dynamics are exponentially close at any time
$t>0$.

  6. (vi)

Suppose that the initial data $x_{0}$ is randomly distributed according to an
even probability measure $\rho_{0}\in L^{1}(\mathbb{R})$ then the estimates in
parts $(iii)$ and $(iv)$ still hold true.

###### Proof.

$(i)$. The law $\rho_{t}$ of $z(t)=\begin{pmatrix}x(t)\\\ v(t)\end{pmatrix}$
satisfies the kinetic Fokker Planck equation

$\partial_{t}\rho_{t}=\mathscr{L}^{*}\rho_{t},\quad\rho|_{t=0}=\delta_{(x_{0},v_{0})},$

where
$\mathscr{L}^{*}\rho:=-v\partial_{x}\rho+\omega^{2}x\partial_{v}\rho+\gamma\big{[}\partial_{v}(v\rho)+\beta^{-1}\partial^{2}_{vv}\rho\big{]}$.

According to [Risken, Section 10.2] $\mu_{t}$ is a bivariate Gaussian measure
with mean $M(t)\in\mathbb{R}^{2}$ and covariane matrix
$\Sigma(t)\in\mathbb{R}^{2\times 2}$. They are $t$ dependent objects given by

$M(t)=\begin{pmatrix}m_{x}(x)\\\
m_{v}(t)\end{pmatrix},\quad\Sigma^{-1}(t)=\begin{pmatrix}[\sigma_{xx}(t)]^{-1}&[\sigma_{xv}(t)]^{-1}\\\
[\sigma_{vx}(t)]^{-1}&[\sigma_{vv}(t)]^{-1}\end{pmatrix},$

where

$\displaystyle m_{x}(t)$
$\displaystyle=\frac{\lambda_{1}e^{-\lambda_{2}t}-\lambda_{2}e^{-\lambda_{1}t}}{\lambda_{1}-\lambda_{2}}x_{0}+\frac{e^{-\lambda_{2}t}-e^{-\lambda_{1}t}}{\lambda_{1}-\lambda_{2}}v_{0},$
$\displaystyle m_{v}(t)$
$\displaystyle=\omega^{2}\frac{e^{-\lambda_{1}t}-e^{-\lambda_{2}t}}{\lambda_{1}-\lambda_{2}}x_{0}+\frac{\lambda_{1}e^{-\lambda_{1}t}-\lambda_{2}e^{-\lambda_{2}t}}{\lambda_{1}-\lambda_{2}}v_{0},$
$\displaystyle\sigma_{xx}(t)$
$\displaystyle=\frac{\gamma\beta^{-1}}{(\lambda_{1}-\lambda_{2})^{2}}\Big{[}\frac{\lambda_{1}+\lambda_{2}}{\lambda_{1}\lambda_{2}}+\frac{4}{\lambda_{1}+\lambda_{2}}(e^{-(\lambda_{1}+\lambda_{2})t}-1)-\frac{1}{\lambda_{1}}e^{-2\lambda_{1}t}-\frac{1}{\lambda_{2}}e^{-2\lambda_{2}t}\Big{]},$
$\displaystyle\sigma_{xv}(t)$
$\displaystyle=\frac{\gamma\beta^{-1}}{(\lambda_{1}-\lambda_{2})^{2}}(e^{-\lambda_{1}t}-e^{-\lambda_{2}t})^{2},$
$\displaystyle\sigma_{vv}(t)$
$\displaystyle=\frac{\gamma\beta^{-1}}{(\lambda_{1}-\lambda_{2})^{2}}\Big{[}\lambda_{1}+\lambda_{2}+\frac{4\lambda_{1}\lambda_{2}}{\lambda_{1}+\lambda_{2}}(e^{-(\lambda_{1}+\lambda_{2})t}-1)-\lambda_{1}e^{-2\lambda_{1}t}-\lambda_{2}e^{-2\lambda_{2}t}\Big{]},$

where

$\lambda_{1}=\frac{\gamma+\sqrt{\gamma^{2}-4\omega^{2}}}{2},\quad\lambda_{2}=\frac{\gamma-\sqrt{\gamma^{2}-4\omega^{2}}}{2},\quad\text{thus}\quad\lambda_{1}+\lambda_{2}=\gamma,\quad\lambda_{1}\lambda_{2}=\omega^{2},\quad\lambda_{1}-\lambda_{2}=\sqrt{\gamma^{2}-4\omega^{2}}.$
(13)

Note that, since in the overdamped regime $\gamma\geq 2\omega$, we have

$\lambda_{2}=\frac{\gamma-\sqrt{\gamma^{2}-4\omega^{2}}}{2}=\frac{4\omega^{2}}{2(\gamma+\sqrt{\gamma^{2}-4\omega^{2}})}\leq\frac{4\omega^{2}}{4\omega}=\omega.$

Since $z(t)$ is a bivariate Gaussian, it follows that the law of $x(t)$, which
is the first marginal of $z(t)$, is a univariate Gaussian measure,
$\mu_{t}=\mathcal{N}(m(t),\sigma(t))$, with mean $m(t)=m_{x}(t)$ and variance
$\sigma(t)=\sigma_{xx}(t)$, where $m_{x}(t)$ and $\sigma_{xx}(t)$ are defined
above. Using (13) we can re-write $m(t)$ and $\sigma(t)$ as follows

$\displaystyle m(t)$
$\displaystyle=e^{-\lambda_{2}t}x_{0}+\frac{e^{-\lambda_{2}t}-e^{-\lambda_{1}t}}{\lambda_{1}-\lambda_{2}}(\lambda_{2}x_{0}+v_{0}),$
(14) $\displaystyle\sigma(t)$
$\displaystyle=\frac{\gamma\beta^{-1}}{(\gamma^{2}-4\omega^{2})}\Big{[}\frac{\gamma}{\omega^{2}}+\frac{4}{\gamma}(e^{-\gamma
t}-1)-\frac{1}{\lambda_{1}}e^{-2\lambda_{1}t}-\frac{1}{\lambda_{2}}e^{-2\lambda_{2}t}\Big{]}$
(15)
$\displaystyle=\frac{1}{\beta\omega^{2}\,[1-4(\omega/\gamma)^{2}]}\big{(}1-e^{-2\lambda_{2}t}\big{)}+\frac{\gamma\beta^{-1}}{(\gamma^{2}-4\omega^{2})}\Big{[}\frac{4}{\gamma}(e^{-\gamma
t}-1)-\frac{e^{-2\lambda_{1}t}-e^{-2\lambda_{2}t}}{\lambda_{1}}\Big{]},$ (16)

where in the last equality we have used the following equality

$\displaystyle\frac{1}{\lambda_{1}}e^{-2\lambda_{1}t}+\frac{1}{\lambda_{2}}e^{-2\lambda_{2}t}$
$\displaystyle=\frac{(\lambda_{1}+\lambda_{2})e^{-2\lambda_{2}t}}{\lambda_{1}\lambda_{2}}+\frac{(e^{-2\lambda_{1}t}-e^{-2\lambda_{2}t})}{\lambda_{1}}$
$\displaystyle=\frac{\gamma
e^{-2\lambda_{2}t}}{\omega^{2}}+\frac{(e^{-2\lambda_{1}t}-e^{-2\lambda_{2}t})}{\lambda_{1}}$

The reduced dynamics is an Ornstein-Uhlenbeck process, therefore its law is a
Gaussian measure, $\bar{\mu}_{t}=\mathcal{N}(\bar{m}(t),\bar{\sigma}^{2}(t))$,
with mean

$\bar{m}(t)=e^{-\alpha t}x_{0}=e^{-\lambda_{2}t}x_{0},$ (17)

and variance

$\bar{\sigma}(t)=\frac{D_{r}}{\alpha}(1-e^{-2\alpha
t})=\frac{1}{\omega^{2}\beta}(1-e^{-2\lambda_{2}t}).$ (18)

$(ii)$ Using the general explicit formula for the Wasserstein distance between
two univariate Gaussian measures, we obtain the Wasserstein distance between
the original dynamics and the reduced dynamics,
$W^{2}_{2}(\mu_{t},\bar{\mu}_{t})$, as follows

$W_{2}^{2}(\mu_{t},\bar{\mu}_{t})^{2}=\Big{(}m_{x}(t)-\bar{m}(t)\Big{)}^{2}+\Big{(}\sqrt{\sigma_{xx}(t)}-\sqrt{\bar{\sigma}(t)}\Big{)}^{2},$
(19)

$(iii)$ We now provide estimate for $W_{2}^{2}(\mu_{t},\bar{\mu}_{t})$ in the
high-friction regime, which corresponds to a large time-scale separation,
since the difference $\lambda_{1}-\lambda_{2}=\sqrt{\gamma^{2}-4\omega^{2}}$
grows with $\gamma$ for fixed $\omega$. We have

$m(t)-\bar{m}(t)=\frac{e^{-\lambda_{2}t}-e^{-\lambda_{1}t}}{\lambda_{1}-\lambda_{2}}(\lambda_{2}x_{0}+v_{0}).$
(20)

Therefore, since
$|e^{-\lambda_{2}t}-e^{-\lambda_{1}t}\leq|e^{-\lambda_{2}t}|+|e^{-\lambda_{1}t}|\leq
2$,

$|m(t)-\bar{m}(t)|\leq\frac{2}{\sqrt{\gamma^{2}-4\omega^{2}}}\big{(}\lambda_{2}|x_{0}|+|v_{0}|\big{)}\leq\frac{2}{\sqrt{\gamma^{2}-4\omega^{2}}}\big{(}\omega|x_{0}|+|v_{0}|\big{)}$

Next we estimate $|\sigma(t)-\bar{\sigma}_{t}|$. Since

$|e^{-\gamma t}-1|\leq e^{-\gamma t}+1\leq
2,\quad|e^{-2\lambda_{1}t}-e^{-2\lambda_{2}t}|\leq
e^{-2\lambda_{1}t}+e^{-2\lambda_{2}t}\leq
2,\quad\gamma_{1}\geq\frac{\gamma}{2}$

we have

$\Big{|}\frac{4}{\gamma}(e^{-\gamma
t}-1)-\frac{e^{-2\lambda_{1}t}-e^{-2\lambda_{2}t}}{\lambda_{1}}\Big{|}\leq\frac{4}{\gamma}|e^{-\gamma
t}-1|+\frac{|e^{-2\lambda_{1}t}-e^{-2\lambda_{2}t}|}{\lambda_{1}}\leq\frac{12}{\gamma}.$

Therefore,

$\displaystyle|\sigma(t)-\bar{\sigma}(t)|$
$\displaystyle=\left|\frac{1-e^{-2\lambda_{2}t}}{\beta\omega^{2}}\Big{[}\frac{1}{1-4(\omega/\gamma)^{2}}-1\Big{]}+\frac{\gamma\beta^{-1}}{(\gamma^{2}-4\omega^{2})}\Big{[}\frac{4}{\gamma}(e^{-\gamma
t}-1)-\frac{e^{-2\lambda_{1}t}-e^{-2\lambda_{2}t}}{\lambda_{1}}\Big{]}\right|$
(21)
$\displaystyle=\left|\frac{1-e^{-2\lambda_{2}t}}{\beta}\frac{4(1/\gamma)^{2}}{1-4(\omega/\gamma)^{2}}+\frac{\gamma\beta^{-1}}{(\gamma^{2}-4\omega^{2})}\Big{[}\frac{4}{\gamma}(e^{-\gamma
t}-1)-\frac{e^{-2\lambda_{1}t}-e^{-2\lambda_{2}t}}{\lambda_{1}}\Big{]}\right|$
$\displaystyle\leq\frac{4}{\beta(\gamma^{2}-4\omega^{2})}+\frac{12}{\beta(\gamma^{2}-4\omega^{2})}=\frac{16}{\beta(\gamma^{2}-4\omega^{2})}.$

It follows that

$\displaystyle W_{2}^{2}(\mu,\bar{\mu})$
$\displaystyle=\big{(}m(t)-\bar{m}(t)\big{)}^{2}+\Big{(}\sqrt{\sigma(t)}-\sqrt{\bar{\sigma}(t)}\Big{)}^{2}$
$\displaystyle\leq\big{(}m_{x}(t)-\bar{m}(t)\big{)}^{2}+|\sigma_{xx}(t)-\bar{\sigma}(t)|$
$\displaystyle\leq\frac{4}{\gamma^{2}-4\omega^{2}}(\omega|x_{0}|+|v_{0}|)^{2}+\frac{16}{\beta(\gamma^{2}-4\omega^{2})}=\frac{4}{\gamma^{2}-4\omega^{2}}\Big{[}(\omega|x_{0}|+|v_{0}|)^{2}+\frac{4}{\beta}\Big{]},$

where to obtain the second line from the first line, we have used the
inequality $(a-b)^{2}\leq|a^{2}-b^{2}|$ for $a,b\geq 0$.

$(iv)$ We have

$\lim_{t\rightarrow\infty}m(t)=\lim_{t\rightarrow\infty}\bar{m}(t)=0\quad\forall
x_{0},v_{0};\quad\lim_{t\rightarrow\infty}\sigma(t)=\frac{1}{\beta\omega^{2}}=:\sigma_{\infty};\quad\lim_{t\rightarrow\infty}\bar{\sigma}(t)=\frac{1}{\beta\omega^{2}}=:\bar{\sigma}_{\infty}=\sigma_{\infty}.$

Thus the original dynamics and the reduced one share the same equilibrium
measure

$\mu_{\infty}=\bar{\mu}_{\infty}=\mathcal{N}(0,\sigma_{\infty}).$

Furthermore, we compute the rates of convergence explicitly

$\displaystyle W_{2}(\mu_{t},\mu_{\infty})^{2}$
$\displaystyle=(m(t)-m_{\infty})^{2}+(\sqrt{\sigma(t)}-\sqrt{\sigma_{\infty}})^{2}$
$\displaystyle\leq m(t)^{2}+|\sigma(t)-\sigma_{\infty}|$
$\displaystyle=\Big{(}e^{-\lambda_{2}t}x_{0}+\frac{e^{-\lambda_{2}t}-e^{-\lambda_{1}t}}{\lambda_{1}-\lambda_{2}}(\lambda_{2}x_{0}+v_{0})\Big{)}^{2}+\frac{\gamma\beta^{-1}}{(\gamma^{2}-4\omega^{2})}\Big{|}\frac{4}{\gamma}e^{-\gamma
t}-\frac{1}{\lambda_{1}}e^{-2\lambda_{1}t}-\frac{1}{\lambda_{2}}e^{-2\lambda_{2}}\Big{|}$
$\displaystyle=e^{-2\lambda_{2}t}\Big{(}x_{0}+\frac{1-e^{-(\lambda_{1}-\lambda_{2})t}}{\lambda_{1}-\lambda_{2}}(\lambda_{2}x_{0}+v_{0})\Big{)}^{2}+\frac{\gamma\beta^{-1}}{(\gamma^{2}-4\omega^{2})}e^{-2\lambda_{2}t}\Big{|}\frac{4}{\gamma}e^{-2\lambda_{1}t}-\frac{1}{\lambda_{1}}e^{-2(\lambda_{1}-\lambda_{2})t}-\frac{1}{\lambda_{2}}\Big{|}$
$\displaystyle\leq Ce^{-2\lambda_{2}t},$

for some constant $C$, which can be computed explicitly (but it is not the
focus of this part), where we have used the fact that
$\lambda_{1}>\lambda_{2}>0$. Thus

$W_{2}(\mu_{t},\mu_{\infty})\leq Ce^{-\lambda_{2}t}.$

Similarly

$\displaystyle W_{2}(\bar{\mu}_{t},\bar{\mu}_{\infty})^{2}$
$\displaystyle=(\bar{m}(t)-\bar{m}_{\infty})^{2}+(\sqrt{\bar{\sigma}(t)}-\sqrt{\bar{\sigma}_{\infty}})^{2}$
$\displaystyle\leq\bar{m}(t)^{2}+|\overline{\sigma}(t)-\overline{\sigma}_{\infty}|$
$\displaystyle=e^{-2\lambda_{2}t}\Big{[}x_{0}^{2}+\frac{1}{\beta\omega^{2}}\Big{]}.$

Thus we also obtain

$W_{2}(\bar{\mu}_{t},\bar{\mu}_{\infty})\leq Ce^{-\lambda_{2}t}.$

$(v)$ Now we estimate $W_{2}^{2}(\mu_{t},\bar{\mu}_{t})$ in the large time
regime. We only need to estimate the difference between the variances
$|\sigma(t)-\bar{\sigma}(t)|$. According to (21), we have

$\displaystyle\sigma(t)-\bar{\sigma}(t)$
$\displaystyle=\frac{1-e^{-2\lambda_{2}t}}{\beta\omega^{2}}\Big{[}\frac{1}{1-4(\omega/\gamma)^{2}}-1\Big{]}+\frac{\gamma\beta^{-1}}{(\gamma^{2}-4\omega^{2})}\Big{[}\frac{4}{\gamma}(e^{-\gamma
t}-1)-\frac{e^{-2\lambda_{1}t}-e^{-2\lambda_{2}t}}{\lambda_{1}}\Big{]}$
$\displaystyle=-\frac{e^{-2\lambda_{2}t}}{\beta\omega^{2}}\Big{[}\frac{1}{1-4(\omega/\gamma)^{2}}-1\Big{]}+\frac{\gamma\beta^{-1}}{(\gamma^{2}-4\omega^{2})}\Big{[}\frac{4}{\gamma}e^{-\gamma
t}-\frac{e^{-2\lambda_{1}t}-e^{-2\lambda_{2}t}}{\lambda_{1}}\Big{]}$

where, to obtain the second line, we have used the following cancellation

$\frac{1}{\beta\omega^{2}}\Big{[}\frac{1}{1-4(\omega/\gamma)^{2}}-1\Big{]}-\frac{\gamma\beta^{-1}}{(\gamma^{2}-4\omega^{2})}\frac{4}{\gamma}=0.$

Therefore,  it holds

$|\sigma(t)-\bar{\sigma}(t)|\leq\frac{e^{-2\lambda_{2}t}}{\beta\omega^{2}}\Big{[}\frac{1}{1-4(\omega/\gamma)^{2}}-1\Big{]}+\frac{\gamma\beta^{-1}}{(\gamma^{2}-4\omega^{2})}\Big{[}\frac{4}{\gamma}e^{-\gamma
t}+\frac{e^{-2\lambda_{2}t}-e^{-2\lambda_{1}t}}{\lambda_{1}}\Big{]}.$

$(vi)$ Now, we can estimate the Wasserstein distance
$W_{2}^{2}(\mu_{t},\bar{\mu}_{t})$ to explore  the long time behaviour, viz.

$\displaystyle W_{2}^{2}(\mu_{t},\bar{\mu}_{t})$
$\displaystyle=\big{(}m(t)-\bar{m}(t)\big{)}^{2}+\Big{(}\sqrt{\sigma(t)}-\sqrt{\bar{\sigma}(t)}\Big{)}^{2}$
$\displaystyle\leq\big{(}m_{x}(t)-\bar{m}(t)\big{)}^{2}+|\sigma_{xx}(t)-\bar{\sigma}(t)|$
$\displaystyle\leq\frac{e^{-\lambda_{2}t}-e^{-\lambda_{1}t}}{\lambda_{1}-\lambda_{2}}(\omega|x_{0}|+|v_{0}|)+\frac{e^{-2\lambda_{2}t}}{\beta\omega^{2}}\Big{[}\frac{1}{1-4(\omega/\gamma)^{2}}-1\Big{]}+\frac{\gamma\beta^{-1}}{(\gamma^{2}-4\omega^{2})}\Big{[}\frac{4}{\gamma}e^{-\gamma
t}+\frac{e^{-2\lambda_{2}t}-e^{-2\lambda_{1}t}}{\lambda_{1}}\Big{]}$
$\displaystyle=\frac{4}{\beta(\gamma^{2}-4\omega^{2})}e^{-\gamma
t}+\Big{[}\frac{\omega|x_{0}|+|v_{0}|}{\sqrt{\gamma^{2}-4\omega^{2}}}+\frac{4}{\beta(\gamma^{2}-4\omega^{2})}+\frac{\gamma}{\beta\lambda_{1}(\gamma^{2}-4\omega^{2})}\Big{]}e^{-\lambda_{2}t}$
$\displaystyle\qquad-\Big{[}\frac{\omega|x_{0}|+|v_{0}|}{\sqrt{\gamma^{2}-4\omega^{2}}}+\frac{\gamma}{\beta\lambda_{1}(\gamma^{2}-4\omega^{2})}\Big{]}e^{-\lambda_{1}t}$
$\displaystyle\leq\Big{[}\frac{\omega|x_{0}|+|v_{0}|}{\sqrt{\gamma^{2}-4\omega^{2}}}+\frac{8}{\beta(\gamma^{2}-4\omega^{2})}+\frac{\gamma}{\beta\lambda_{1}(\gamma^{2}-4\omega^{2})}\Big{]}e^{-\lambda_{2}t}$
$\displaystyle\leq\Big{[}\frac{\omega|x_{0}|+|v_{0}|}{\sqrt{\gamma^{2}-4\omega^{2}}}+\frac{10}{\beta(\gamma^{2}-4\omega^{2})}\Big{]}e^{-\lambda_{2}t}.$

Here we have used the fact that $\gamma\geq\lambda_{2}$ and
$\frac{\gamma}{\lambda_{1}}=\frac{2\gamma}{\gamma+\sqrt{\gamma^{2}-4\omega^{2}}}\leq
2$.

$(vi)$. Suppose that $x_{0}$ is randomly distributed following an even
distribution $\rho_{0}$. Then the laws of $x(t)$ and $\bar{x}(t)$ are given by

$\mu_{t}=\mathcal{N}(m(t),\sigma(t))\ast\rho_{0},\quad\mu_{t}=\mathcal{N}(\bar{m}(t),\bar{\sigma}(t))\ast\rho_{0}.$

Since
$\mathcal{N}(m(t),\sigma(t)),\mathcal{N}(\bar{m}(t),\bar{\sigma}(t))\in\mathcal{P}_{2}(\mathbb{R})$,
according to [San15, Lemma 5.2] we have

$W_{2}^{2}(\mu_{t},\bar{\mu}_{t})=W_{2}^{2}(\mathcal{N}(m(t),\sigma(t))\ast\rho_{0},\mathcal{N}(\bar{m}(t),\bar{\sigma}(t))\ast\rho_{0})\leq
W^{2}(\mathcal{N}(m(t),\sigma(t)),\mathcal{N}(\bar{m}(t),\bar{\sigma}(t)),$

thus the upper bound estimates in the two previous parts are still true. ∎

##  4 Model reduction of two coupled underdamped Brownian oscillators

We now proceed with the computation of the Wasserstein distance for a slightly
more elaborate model, corresponding to a system of two coupled overdamped
Brownian harmonic oscillators. The dynamics of the model can conveniently be
written as follows:

$\displaystyle\dot{x}_{1}$
$\displaystyle=ax_{1}+k(x_{2}-x_{1})+\sigma_{1}\dot{W}_{1}$ (22a)
$\displaystyle\dot{x}_{2}$
$\displaystyle=-k(x_{2}-x_{1})+dx_{2}+\sigma_{2}\dot{W}_{2},$ (22b)

where $\dot{W}$ denotes the formal derivative of a Wiener process,
corresponding to a white noise, $a,d<0$ are parameters characteristic of the
individual oscillator (without loss of generality we also assume $a\geq d$),
$\sigma_{1},\sigma_{2}>0$ denote the noise strenghts, and finally, $k>0$ is
the coupling parameter.

The system (22) represents the overdamped version of the coupled underdamped
Langevin dynamics of the two oscillators. A contracted description for the
deterministic case (i.e., with $\sigma_{1}=\sigma_{2}=0$) under a suitable
assumption of scale separation is studied, with applications to relaxation
dynamics in proteins, in [SMR11]. We can derive a reduced system by
eliminating the variable $x_{2}$, in (22), using the procedure introduced in
[CM22, CDM22]. This consists of two distinct steps: (i) the deterministic
component of the dynamics is obtained using the Invariant Manifold method,
then (ii) the diffusion terms are determined via fulfilling the Fluctuation-
Dissipation relation.

###  4.1 Deterministic evolution

Let $\langle\mathcal{O}\rangle$ denote the average over noise of the variable
$\mathcal{O}$. The original dynamics can be written as

$\dot{\mathbf{z}}=\mathbf{Q}~{}\mathbf{z}\;,$ (23)

where $\mathbf{z}=(\langle x_{1}\rangle,\langle x_{2}\rangle)$ and

$\mathbf{Q}=\mathbf{Q}(k)=\begin{pmatrix}a-k&k\\\ k&-k+d\end{pmatrix}$ (24)

The characteristic polynomial of $\mathbf{Q}$ is

$\lambda^{2}-(a+d-2k)\lambda+(ad-ak-dk)=0.$

Thus $\mathbf{Q}$ has two real negative eigenvalues:

$\lambda_{\pm}=\lambda_{\pm}(k):=\frac{(a+d-2k)\pm\sqrt{(a-d)^{2}+4k^{2}}}{2}\;,$
(25)

In this model, the time-scale separation is encoded in the difference
$\lambda_{+}-\lambda_{-}=\sqrt{(a-d)^{2}+4k^{2}}$, which grows with increasing
$k$, for fixed parameters $a,d$. We seek a closure of the form $\langle
x_{2}\rangle=\alpha\langle x_{1}\rangle$, hence, following [CDM22], we define
a macroscopic time derivative of $\langle x_{2}\rangle$ via the chain rule:

$\displaystyle\partial_{t}^{macro}\langle x_{2}\rangle$
$\displaystyle:=\frac{\partial\langle x_{2}\rangle}{\partial\langle
x_{1}\rangle}\langle\dot{x}_{1}\rangle$
$\displaystyle=(\alpha(a-k)+\alpha^{2}k)\langle x_{1}\rangle\;,$

which expresses the slaving principle mentioned in Sec. 1. Furthermore, we
also define the microscopic time derivative of $\langle x_{2}\rangle$ in terms
of the vector field given in Eq. (23), where $\langle x_{2}\rangle$ is
expressed through the aforementioned closure. We thus set:

$\displaystyle\partial_{t}^{micro}\langle x_{2}\rangle$
$\displaystyle:=k\langle x_{1}\rangle+(d-k)\langle x_{2}\rangle$
$\displaystyle=(k+\alpha(d-k))\langle x_{1}\rangle\;.$

The Invariant Manifold method requires that microscopic and macroscopic time
derivatives of $\langle x_{2}\rangle$ coincide, independently of the values of
the observable $\langle x_{1}\rangle$. Thus, we obtain the following
invariance equation

$\alpha(a-k)+\alpha^{2}k=k+\alpha(d-k)\quad\Longleftrightarrow\quad
k\alpha^{2}+(a-d)\alpha+k=0\;,$ (26)

which has two solutions

$\alpha_{\pm}=\alpha_{\pm}(k):=\frac{-(a-d)\pm\sqrt{(a-d)^{2}+4k^{2}}}{2k}\;.$

The reduced dynamics for the deterministic part is

$\langle\dot{x}_{1}\rangle=(a-k+k\hat{\alpha})\langle x_{1}\rangle\;,$ (27)

where $\hat{\alpha}\in\\{\alpha_{+},\alpha_{-}\\}$ which will be specified
later. It is noticeable that

$a-k+k\alpha_{\pm}=\frac{(a+d-2k)\pm\sqrt{(a-d)^{2}+4k^{2}}}{2}\equiv\lambda_{\pm}\;.$

Looking at (27),  we notice that  the coefficient multiplying $\langle
x_{1}\rangle$ coincides with one of the eigenvalues of the matrix
$\mathbf{Q}$. To pick up the right eigenvalue, we use the following criterion.
We select $\hat{\alpha}$ from solutions $\alpha_{\pm}$ to the invariance
equation (26) that satisfies $a-k+k\hat{\alpha}\rightarrow a<0$ as
$k\rightarrow 0$, that is $k\hat{\alpha}\rightarrow 0$ as $k\rightarrow 0$.
Since we assume that $a\geq d$, we take

$\hat{\alpha}=\alpha_{+}=\frac{-(a-d)+\sqrt{(a-d)^{2}-4k^{2}}}{2k}\;.$

###  4.2 Incorporating the noise

To characterize the noise term, we employ the methodology proposed in
[CDM22]. Therefore, we first define the diffusion matrix $\boldsymbol{D}$ as

$\boldsymbol{D}=\begin{pmatrix}\sigma_{1}&0\\\ 0&\sigma_{2}\end{pmatrix}\;,$
(28)

and we also denote

$\dot{\mathbf{W}}=(\dot{W}_{1},\dot{W}_{2})\;.$

The solution of Eqs. (22) reads:

$\mathbf{z}(t)=e^{\mathbf{Q}t}\mathbf{z}_{0}+\int_{0}^{t}e^{\mathbf{Q}(t-s)}\mathbf{D}~{}\dot{\mathbf{W}}ds\;.$
(29)

We thus find

$\displaystyle\lim_{t\rightarrow\infty}\mathbb{E}[z_{1}(t)^{2}]\equiv\overline{\Sigma}_{11}=-\frac{1}{2}\Big{(}\frac{1}{a-2k}+\frac{1}{a}\Big{)}.$

The full reduced system takes hence the form

$d\hat{x}(t)=\lambda_{+}\hat{x}(t)\,dt+\hat{D}dW_{t}\;,$ (30)

where the drift coefficient $\lambda_{+}$ is defined in (25) and the diffusion
coefficient $\hat{D}$ is given by

$\hat{D}=-\lambda_{+}\overline{\Sigma}_{11}\;.$ (31)

###  4.3 Quantification of errors and the long-time behaviour

In this section we will compute explicitly the Wasserstein distance between
the laws of the original dynamics of $x_{1}$ and of the reduced dynamics (30)
and study their long-time behaviour. The Fokker Planck equation associated to
the full original dynamics (22) is given by the following linear-drift
diffusion equation

$\partial_{t}\rho=-\mathop{\mathrm{div}}\nolimits(\mathbf{Q}\rho)+\mathop{\mathrm{div}}\nolimits(\boldsymbol{D}\nabla\rho),$
(32)

where $\rho=\rho(t,x_{1},x_{2})$ is the joint probability density of
$(x_{1},x_{2})$, the drift matrix $\boldsymbol{Q}$ and the diffusion matrix
$\boldsymbol{D}$ are given in (24) and (28) respectively. Note that the above
system is a special case of the general drift-diffusion equation introduced in
Section 2.2.

Since we are focusing on the role of the coupling parameter, for simplicity of
presentation, we consider identical oscillator, that is $a=d<0$ and
normalising $\sigma_{1}=\sigma_{2}=1$, so that

$\mathbf{Q}=\begin{pmatrix}a-k&k\\\
k&a-k\end{pmatrix},\quad\text{and}\quad\boldsymbol{D}=I.$

The main result of this section is the following theorem.

######  Theorem 4.1.

Let $\rho_{1}(t)$ be the distribution of $x_{1}(t)$ of the original coupled
dynamics (30) starting at a deterministic initial data
$(x_{1},x_{2})(0)=(x_{1},x_{2})$, and $\hat{\rho}_{1}(t)$ be the distribution
of the reduced dynamics (30) starting from $x_{1}$. Then there exists a
constant $C>0$ such that the following statements hold

  1. (i)

$W_{2}(\rho_{1}(t),\hat{\rho}_{1}(t))^{2}\leq Ck$.

  2. (ii)

$\max\\{W_{2}(\rho_{1}(t),\rho_{\infty}),W_{2}(\hat{\rho}_{1}(t),\rho_{\infty})\\}\leq
Ce^{at}$, where $\rho_{\infty}=\mathcal{N}(0,\overline{\Sigma}_{11})$

  3. (iii)

$W_{2}(\rho_{1}(t),\hat{\rho}_{1}(t)\leq Ce^{at}$.

###### Proof.

According to Lemma 2.1, the solution to (32) is given by
$\rho(t,x_{1},x_{2})=\mathcal{N}(\mu(t),\Sigma(t)),$ where

$\mu(t)=e^{t\mathbf{Q}}\begin{pmatrix}x_{1}\\\
x_{2}\end{pmatrix},\quad\Sigma(t)=2\int_{0}^{t}e^{s\mathbf{Q}}e^{s\mathbf{Q}^{T}}\,ds.$

Since $\mathbf{Q}\mathbf{Q}^{T}=\mathbf{Q}^{T}\mathbf{Q}$ and
$\mathbf{Q}=\mathbf{Q}^{T}$, we have

$e^{s\mathbf{Q}}e^{s\mathbf{Q}^{T}}=e^{s(\mathbf{Q}+\mathbf{Q}^{T})}=e^{2s\mathbf{Q}}.$

Thus, we can simplify $\Sigma(t)$ as

$\Sigma(t)=2\int_{0}^{t}e^{2s\mathbf{Q}}\,ds.$

Applying lemma 2.2, we compute

$\displaystyle e^{t\mathbf{Q}}$
$\displaystyle=\frac{1}{\Delta}\begin{pmatrix}m_{11}&m_{12}\\\
m_{21}&m_{22}\end{pmatrix},\quad\Delta=2kt,$ $\displaystyle m_{11}$
$\displaystyle=m_{22}=e^{(a-k)t}\Delta\cosh{\frac{1}{2}\Delta}=\frac{1}{2}e^{(a-k)t}\Delta(e^{kt}+e^{-kt})=\frac{1}{2}\Delta(e^{(a-2k)t}+e^{at})$
$\displaystyle m_{12}$
$\displaystyle=m_{21}=2kte^{(a-k)t}\sinh{\frac{1}{2}\Delta}=\frac{1}{2}\Delta
e^{(a-k)t}(e^{kt}-e^{-kt}).=\frac{1}{2}\Delta(e^{at}-e^{(a-2k)t}).$

Thus

$e^{t\mathbf{Q}}=\frac{1}{2}\begin{pmatrix}e^{(a-2k)t}+e^{at}&e^{at}-e^{(a-2k)t}\\\
e^{at}-e^{(a-2k)t}&e^{(a-2k)t}+e^{at}\end{pmatrix}.$

Similarly

$e^{2t\mathbf{Q}}=\frac{1}{2}\begin{pmatrix}e^{2(a-2k)t}+e^{2at}&e^{2at}-e^{2(a-2k)t}\\\
e^{2at}-e^{2(a-2k)t}&e^{2(a-2k)t}+e^{2at}\end{pmatrix}.$

Therefore,

$\Sigma(t)=2\int_{0}^{t}e^{2sQ}\,ds=\frac{1}{2}\begin{pmatrix}\frac{e^{2(a-2k)t}-1}{a-2k}+\frac{e^{2at}-1}{a}&\frac{e^{2at}-1}{a}-\frac{e^{2(a-2k)t}-1}{a-2k}\\\
\frac{e^{2at}-1}{a}-\frac{e^{2(a-2k)t}-1}{a-2k}&\frac{e^{2(a-2k)t}-1}{a-2k}+\frac{e^{2at}-1}{a}\end{pmatrix}.$

It follows that

$\rho_{1}(t)=\mathcal{N}(\mu_{1}(t),\Sigma_{11}(t))=\mathcal{N}\Bigg{(}\frac{1}{2}\Big{(}(e^{(a-2k)t}+e^{at})x_{1}+(e^{at}-e^{(a-2k)t})x_{2}\Big{)},\frac{1}{2}\Big{(}\frac{e^{2(a-2k)t}-1}{a-2k}+\frac{e^{2at}-1}{a}\Big{)}\Bigg{)},$

Since $\hat{x}$ is an OU process, we obtain

$\hat{\rho}_{1}(t)=\mathcal{N}(\hat{\mu}_{1}(t),\hat{\Sigma}_{1}(t))=\mathcal{N}\Big{(}e^{\lambda_{+}t}x_{1},-\frac{\hat{D}}{\lambda_{+}}(1-e^{2\lambda_{+}t})\Big{)}=\mathcal{N}\Big{(}e^{\lambda_{+}t}x_{1},\overline{\Sigma}_{11}(1-e^{2\lambda_{+}t})\Big{)},$

recalling that, with $a=d$

$\lambda_{+}=\frac{(a+d-2k)+\sqrt{(a-d)^{2}+4k^{2}}}{2}=a,\quad\overline{\Sigma}_{11}=-\frac{1}{2}\Big{(}\frac{1}{a-2k}+\frac{1}{a}\Big{)}.$

The Wasserstein distance between $\rho_{1}$ and $\hat{\rho}_{1}$ is given by

$W_{2}(\rho_{1}(t),\hat{\rho}_{1}(t))^{2}=(\mu_{1}(t)-\hat{\mu}_{1}(t))^{2}+\Big{(}\sqrt{\Sigma_{11}}(t)-\sqrt{\hat{\Sigma}_{1}}(t)\Big{)}^{2}$
(33)

$(i)$ We compute

$\displaystyle|\mu_{1}(t)-\hat{\mu}_{1}(t)|$
$\displaystyle=\frac{1}{2}\Big{|}\Big{(}e^{(a-2k)t}+e^{at})x_{1}+(e^{at}-e^{(a-2k)t})x_{2}\Big{)}-e^{at}x_{1}\Big{|}$
$\displaystyle=\frac{1}{2}|(e^{at}-e^{(a-2k)t})(x_{2}-x_{1})|$
$\displaystyle=\frac{1}{2}e^{at}|x_{2}-x_{1}|(1-e^{-2kt})$ (34)
$\displaystyle\leq|x_{2}-x_{1}|ke^{at}t$ $\displaystyle\leq
k\,|x_{2}-x_{1}|\frac{1}{|a|e},$

where in the first inequality we have used the elementary inequality
$1-e^{-x}\leq x$ for all $x>0$, and in the last inequality we have used
(noting that $a<0$)

$\max_{t>0}te^{at}=\frac{1}{|a|e}.$ (35)

We also estimate

$\displaystyle\Sigma_{11}(t)-\hat{\Sigma}_{1}(t)$
$\displaystyle=\frac{1}{2}\Bigg{[}\frac{1}{a-2k}\Big{(}e^{2(a-2k)t}-e^{2\lambda_{+}t}\Big{)}+\frac{1}{a}\Big{(}e^{2at}-e^{2\lambda_{+}t}\Big{)}\Bigg{]}$
$\displaystyle=\frac{1}{2}\frac{1}{a-2k}\Big{(}e^{2(a-2k)t}-e^{2at}\Big{)}$
$\displaystyle=\frac{1}{2}\frac{1}{2k-a}e^{2at}\Big{(}1-e^{-4kt}\Big{)}$ (36)
$\displaystyle\leq\frac{2k}{2k-a}e^{2at}t$
$\displaystyle\leq\frac{k}{a^{2}e},$

where to go from (36) to the next line, we have used $1-e^{-4kt}\leq 4kt$ and
(35) again (with $a$ replaced by $2a$). Therefore, we have

$\displaystyle W_{2}(\rho_{1}(t),\hat{\rho}_{1}(t))^{2}$
$\displaystyle\leq(\mu_{1}(t)-\hat{\mu}_{1}(t))^{2}+\Big{|}\Sigma_{11}(t)-\overline{\Sigma}_{1}(t)\Big{|}$
$\displaystyle\leq
k^{2}\,|x_{2}-x_{1}|^{2}\frac{1}{|a|^{2}e^{2}}+\frac{k}{a^{2}e}\leq Ck,$

for any bounded $k$.

$(iii)$ Since $a<0$,

$\lim_{t\rightarrow\infty}\mu_{1}(t)=\lim_{t\rightarrow\infty}\hat{\mu}_{1}(t)=0,\quad\lim_{t\rightarrow\infty}\Sigma_{11}(t)=-\frac{1}{2}\Big{(}\frac{1}{a-2k}+\frac{1}{a}\Big{)}=\overline{\Sigma}_{11}.$

it implies that

$\lim_{t\rightarrow\infty}\rho_{1}(t)=\lim_{t\rightarrow\infty}\hat{\rho}_{1}(t)=\rho_{\infty}=\mathcal{N}(0,\overline{\Sigma}_{11}).$

We can also compute explicitly the rates of convergence of these limits in the
Wasserstein distance. We have

$\displaystyle
W_{2}(\rho_{1}(t),\rho_{\infty})^{2}=\mu_{1}(t)^{2}+\Big{(}\sqrt{\Sigma_{11}}(t)-\sqrt{\overline{\Sigma}}\Big{)}^{2}\leq\mu_{1}(t)^{2}+\Big{|}\Sigma_{11}(t)-\overline{\Sigma}_{11}\Big{|}.$
(37)

We estimate each term on the right hand side of (37). For the first term, we
get

$\displaystyle\mu_{1}(t)$
$\displaystyle=\frac{1}{2}\Big{(}(e^{(a-2k)t}+e^{at})x_{1}+(e^{at}-e^{(a-2k)t})x_{2}\Big{)}$
$\displaystyle=\frac{1}{2}e^{at}\Big{(}(1+e^{-2kt})x_{1}+(1-e^{-2kt})x_{2}\Big{)}$
$\displaystyle\leq Ce^{at}.$ (38)

For the second term, we have

$\displaystyle|\Sigma_{11}(t)-\overline{\Sigma}_{11}|=\frac{1}{2}\Big{|}\frac{e^{2(a-2k)t}}{a-2k}+\frac{e^{2at}}{a}\Big{|}=\frac{1}{2}e^{2at}\Big{|}\frac{1}{a}+\frac{e^{-4kt}}{a-2k}\Big{|}\leq
Ce^{2at}.$ (39)

Substituting (4.3) and (39) to (37), we obtain

$W_{2}(\rho_{1}(t),\rho_{\infty})\leq Ce^{at},$

thus $\rho_{1}$ exponentially converges, with a rate $a$, to $\rho_{\infty}$.
Similarly,

$\displaystyle W_{2}(\hat{\rho}_{1}(t),\rho_{\infty})^{2}$
$\displaystyle=\hat{\rho}_{1}(t)^{2}+\Big{(}\sqrt{\hat{\Sigma}_{1}}(t)-\sqrt{\overline{\Sigma}}_{11}\Big{)}^{2}$
$\displaystyle\leq\hat{\rho}_{1}(t)^{2}+\Big{|}\hat{\Sigma}_{1}(t)-\overline{\Sigma}_{11}\Big{|}$
$\displaystyle=(x_{1}^{2}+\overline{\Sigma}_{11})e^{2at}.$

Hence $\hat{\rho}_{1}$ exponentially converges with the same rate $a$ to
$\rho_{\infty}$.

$(iv)$ According to (34) and (36) we have

$\displaystyle|\mu_{1}(t)-\hat{\mu}_{1}(t)|=\frac{1}{2}e^{at}|x_{2}-x_{1}|(1-e^{-2kt})\leq
Ce^{at},$
$\displaystyle|\Sigma_{11}(t)-\hat{\Sigma}_{1}(t)|=\frac{1}{2}\frac{1}{|2k-a|}e^{2at}\Big{(}1-e^{-4kt}\Big{)}\leq
Ce^{2at}.$

Thus

$W_{2}(\rho_{1}(t),\hat{\rho}_{1}(t))^{2}\leq(\mu_{1}(t)-\hat{\mu}_{1}(t))^{2}+|\Sigma_{11}(t)-\hat{\Sigma}_{1}(t)|\leq
Ce^{2at},$

that is

$W_{2}(\rho_{1}(t),\hat{\rho}_{1}(t))\leq Ce^{at}.$

This completes the proof of this theorem. We remark that we have assumed
deterministic initial data, but the theorem can also be extended to the case
where the initial data follow symmetric distributions as in Section 3. ∎

##  5 Summary and outlook

In this work we have employed the reduction scheme recently introduced in
[CM22, CDM22], which suitably combines the Invariant Manifold method  with the
Fluctuation-Dissipation relation, to derive a contracted description for two
classical models of statistical physics, namely the underdamped Brownian
harmonic oscillator and a system of two coupled overdamped Brownian harmonic
oscillators. The present work significantly extends the previous results: we
succeeded here to quantify explicitly the error between the original and the
reduced dynamics, as well as their rates of convergence to equilibrium. The
technical tool we used is the Wasserstein distance, which is widely employed
in the theory of optimal transport. We have thus shown that the two dynamics
are exponentially close at any time, share the same equilibrium measure, and
exponentially converge to the same equilibrium measure with the same rate.
Furthermore, the two dynamics are also found to coincide if the relevant
parameter controlling the time-scale separation of the original model is sent
to infinity. The linearity of the considered models has clearly played an
important role in the analysis of this work, enabling the explicit
computations of their solutions and of the involved Wasserstein distances. A
key challenge for future developments is to generalize our analysis in order
to deal with non-linear models, where explicit solutions and computations are
not accessible. Another direction of research points toward the investigations
of systems with a large numbers of degrees of freedom, e.g. models relevant to
climate dynamics [HAK23], or small systems of interest in modern
nanotechnologies, such as biomolecular motors [WKST16].

## Acknowledgements

MC’s research was performed under the auspices of Italian National Group of
Mathematical Physics (GNFM) of INdAM. MHD research was supported by EPSRC
grants EP/W008041/1 and EP/V038516/1.

## References

  * [BS93] D.S. Bernstein and W. So.  Some explicit formulas for the matrix exponential.  IEEE Transactions on Automatic Control, 38(8):1228–1232, 1993.
  * [CDM22] M. Colangeli, M. H. Duong, and A. Muntean.  A reduction scheme for coupled Brownian harmonic oscillators.  Journal of Physics A: Mathematical and Theoretical, 55:505002, 2022\.
  * [CKK07] M. Colangeli, I. V. Karlin, and M. Kröger.  Hyperbolicity of exact hydrodynamics for three-dimensional linearized Grad’s equations.  Phys. Rev. E, 76:022201, 2007.
  * [CKO09] M. Colangeli, M. Kröger, and H. C. Öttinger.  Boltzmann equation and hydrodynamic fluctuations.  Phys. Rev. E, 80:051202, 2009.
  * [CM22] M. Colangeli and A. Muntean.  Reduced markovian descriptions of brownian dynamics: Toward an exact theory.  Front. Phys. 10:903030. doi: 10.3389/fphy.2022.903030, 10:903030, 2022.
  * [GK05] A. N. Gorban and I. V. Karlin.  Invariant Manifolds for Physical and Chemical Kinetics, volume 660 of Lect. Notes Phys. Springer-Verlag, Berlin, 2005.
  * [GKS04] D. Givon, R. Kupferman, and A. Stuart.  Extracting macroscopic dynamics: model problems and algorithms.  Nonlinearity, 17(6):R55, 2004.
  * [GL20] M. Ghil and V. Lucarini.  The physics of climate variability and climate change.  Reviews of Modern Physics, 92:035002, 2020.
  * [GLCG21] M. S. Gutiérrez, V. Lucarini, M. D. Chekroun, and M. Ghil.  Reduced-order models for coupled dynamical systems: Data-driven methods and the Koopman operator.  Chaos, 31:053116, 2021.
  * [GP18] S. N. Gomes and G. A. Pavliotis.  Mean field limits for interacting diffusions in a two-scale potential.  Journal of nonlinear science, 28(3):905–941, 2018.
  * [Gra49] H. Grad.  On the kinetic theory of rarefied gases.  Comm. Pure App. Math., 2:331–407, 1949.
  * [GS84] C. R. Givens and R. M. Shortt.  A class of Wasserstein metrics for probability distributions.  Michigan Mathematical Journal, 31(2):231–240, 1984.
  * [Hak04] H. Haken.  Synergetics. Introduction and Advanced Topics.  Graduate studies in mathematics. Springer-Verlag Berlin Heidelberg, 2004\.
  * [HAK23] F. Hummel, P. Ashwin, and C. Kuehn.  Reduction methods in climate dynamics—a brief review.  Physica D: Nonlinear Phenomena, page 133678, 2023.
  * [IOMF14] E.R. Ijioma, T. Ogawa, A. Muntean, and T. Fatima.  Homogenization and dimension reduction of filtration combustion in heterogeneous thin layers.  Networks & Heterogeneous Media, 9(4):709–737, 2014.
  * [KK13] H.-W. Kang and T.G̃. Kurtz.  Separation of Time-Scales and Model Reduction for Stochastic Reaction Networks.  The Annals of Applied Probability, 23(2):529–583, 2013.
  * [PS08] G. Pavliotis and A. Stuart.  Multiscale methods: averaging and homogenization.  Springer Science & Business Media, 2008.
  * [San15] F. Santambrogio.  Optimal Transport for Applied Mathematicians.  Birkäuser, NY, 55(58-63):94, 2015.
  * [SMR11] R. Soheilifard, D. E. Makarov, and G. J. Rodin.  Rigorous coarse-graining for the dynamics of linear systems with applications to relaxation dynamics in proteins.  The Journal of Chemical Physics, 135(5):054107, 2011.
  * [SS17] S. K. Schoffner and S. Schnell.  The physics of climate variability and climate change.  Mathematical Biosciences, 287:122–129, 2017.
  * [Tak12] A. Takatsu.  Wasserstein geometry of porous medium equation.  In Annales de l’Institut Henri Poincaré C, Analyse non linéaire, volume 29, pages 217–232. Elsevier, 2012.
  * [Vil03] C. Villani.  Topics in Optimal Transportation.  Graduate Studies in Mathematics. American Mathematical Society, 2003\.
  * [WKST16] S.-W. Wang, K. Kawaguchi, S. Sasa, and K.-H. Tang.  Entropy Production of Nanosystems with Time Scale Separation.  Phys. Rev. Lett., 117:070601, 2016.
  * [Zwa01] R. Zwanzig.  Nonequilibrium Statistical Mechanics.  Oxford University Press, 2001.


# Deep Learning for Moving Blockage Prediction using Real Millimeter Wave
Measurements

Shunyao Wu, Muhammad Alrabeiah, Andrew Hredzak, Chaitali Chakrabarti, and
Ahmed Alkhateeb
School of Electrical, Computer and Energy Engineering
Arizona State University Tempe, AZ 85287
Email: {vincentw, malrabei, ahredzak, chaitali<EMAIL_ADDRESS>
###### Abstract

Millimeter wave (mmWave) communication is a key component of 5G and beyond.
Harvesting the gains of the large bandwidth and low latency at mmWave systems,
however, is challenged by the sensitivity of mmWave signals to blockages; a
sudden blockage in the line of sight (LOS) link leads to abrupt disconnection,
which affects the reliability of the network. In addition, searching for an
alternative base station to re-establish the link could result in needless
latency overhead. In this paper, we address these challenges collectively by
utilizing machine learning to anticipate dynamic blockages proactively. The
proposed approach sees a machine learning algorithm learning to predict future
blockages by observing what we refer to as the pre-blockage signature. To
evaluate our proposed approach, we build a mmWave communication setup with a
moving blockage and collect a dataset of received power sequences. Simulation
results on a real dataset show that blockage occurrence could be predicted
with more than 85% accuracy and the exact time instance of blockage occurrence
can be obtained with low error. This highlights the potential of the proposed
solution for dynamic blockage prediction and proactive hand-off, which
enhances the reliability and latency of future wireless networks.

###### Index Terms:

Millimeter wave, machine learning, blockage prediction, handover

##  I Introduction

Communication at the mmWave frequency range offers high bandwidth and higher
data rate demands required by 5G and beyond cellular systems [1]. However,
mmWave systems struggle in the presence of objects that block the LOS
connection between a base station and its users [2]. This is majorly rooted
into the poor penetration and reflection properties of mmWave signals [1, 3].
This struggle with blockages or Non-LOS (NLOS) connections ultimately affects
the reliability of the system. In addition to that, it introduces a latency
burden resulting from the need for user hand-off [2]. A key approach to
address that struggle lies in equipping the mmWave system with the capability
to predict possible blockage proactively. Such successful prediction allows a
base station to take mitigation measures, e.g., user hand-off, before the link
is blocked, thereby resolving the reliability and latency problems.

Many recently published studies have used machine learning to address problems
arising from link blockages in MIMO and mmWave communications [4, 5, 2, 6, 7].
The work in [4, 5, 6] collectively demonstrates the ability of a machine
learning algorithm (whether deep [5] or shallow [4]) to identify or
differentiate LOS and NLOS links. This identification task could be an
interesting ability to a system operating in the sub-6 GHz range. However, the
requirements are more stringent for a mmWave system and demand a proactive
approach. A step towards doing so is presented in [2], where a proactive
solution is proposed to predict stationary blockages. This solution depends on
beam sequences alone, and as such, it cannot handle dynamic blockages. Another
direction addressing blockage prediction relies on Vision-Aided Wireless
Communications [8, 7]. In [8], a camera feed along with sub-6 GHz channels are
used to identify currently blocked mmWave links. This work is the seed to that
presented in [7], where proactive blockage prediction using images and mmWave
beams is attempted using bimodal deep learning algorithm. Despite the appeal,
that work require extra sensory data (images), and it does not take full
advantage of the wireless data.

In this paper, we propose a machine learning algorithm that addresses
proactive dynamic-blockage prediction. The algorithm uses sequences of
received powers to predict whether a blockage is incoming or not. The basic
idea behind our algorithm is the ability to recognize what we have dubbed pre-
blockage signature. We argue that signature could serve as a important clue
for incoming blockages. Our contribution is summarized as follows:

  1. 1.

We propose a recurrent neural network (RNN) architecture based on Gated-
recurrent units (GRUs) to predict incoming blockages. The architecture is
designed to learn one of two tasks. The first is to predict whether a blockage
is incoming or not, and the other is to pinpoint the time instance at which
the blockage occurs.

  2. 2.

We develop a mmWave communication setup with a moving blockage. We use that
setup to build a dataset of received power sequences and their corresponding
future link statuses. The dataset is utilized to the proposed algorithm and
its ability to predict incoming blockages.

The rest of this paper is organized as follows. Section II presents the system
and channel models adopted to study the dynamic-blockage prediction. Section
III presents a formulation of the prediction problem. The proposed machine
learning model is presented in Section IV. The data collection scenario and
setup is introduced in Section V. Proposed algorithm evaluation and simulation
results are presented in Section VI, and, finally, the paper is concluded in
Section VII.

##  II System and Channel Models

The following two subsections will introduce the system and channel models
adopted in this paper.

Figure 1: An illustration of the considered system model.

System model: The communication system considered in this work is described in
Fig.1. It assumes a base station serving a static user who is in the vicinity
of a possible moving blockage. The base station employs an $M$-element Uniform
Linear Array (ULA) antenna operating at 60GHz carrier frequency with
Orthogonal Frequency Division Multiplexing (OFDM). It also assumes a fully
analog architecture with a predefined beam-steering codebook
$\mathcal{F}=\\{\mathbf{f}_{w}\\}_{w=1}^{W}$, where
$\mathbf{f}_{w}\in\mathbb{C}^{M\times 1}$ is given by:

$\mathbf{f}_{w}=\frac{1}{\sqrt{M}}\left[1,e^{j\frac{2\pi}{\lambda}d\sin(\phi_{w})},\dots,e^{j(M-1)\frac{2\pi}{\lambda}d\cos(\phi_{w})}\right]^{T},$
(1)

where $d$ is the spacing between the ULA elements, $\lambda$ is the
wavelength, and $\phi_{w}\in\\{\frac{2\pi w}{W}\\}_{w=0}^{W-1}$ is a uniformly
quantized azimuth angle with a step of $1/W$. At any time instance $t$, the
downlink received signal is expressed as follows:

$r_{k}[t]=\mathbf{h}_{k}[t]^{T}\mathbf{f}_{w}s_{k}[t]+n_{k}$ (2)

where $\mathbf{h}_{k}[t]\in\mathbb{C}^{M\times 1}$ is the downlink channel at
the $k$th subcarrier, $s_{k}[t]$ is the symbol transmitted on the $k$th
subcarrier, and, finally, $n_{k}$ is a complex Gaussian noise sample,
$\sim\mathcal{CN}(0,\sigma^{2})$ at the $k$th subcarrier.

Channel model: This work adopts the geometric (physical) channel model, which
captures the physical characteristics of signal propagation including the
dependence on the environment geometry, materials, frequency band, etc., [1].
With this model, the channel can be expressed as:

$\mathbf{h}_{k}=\sum_{d=0}^{D-1}\sum_{\ell=1}^{L}\alpha_{\ell}e^{-\j\frac{2\pi
k}{K}d}p\left(dT_{\mathrm{S}}-\tau_{\ell}\right){\mathbf{a}}\left(\theta_{\ell},\phi_{\ell}\right),$
(3)

where $L$ is number of channel paths, $\alpha_{\ell}$ is the path gain
(including the path-loss), $\tau_{\ell}$ is the delay, $\theta_{\ell}$ is the
azimuth angle of arrival, and $\phi_{\ell}$ is the elevation angle, of the
$\ell$th channel path. $T_{\mathrm{S}}$ represents the sampling time while $D$
denotes the cyclic prefix length (assuming that the maximum delay is less than
$DT_{\mathrm{S}}$).

##  III Problem Formulation

Proactively identifying Line of Sight (LOS) link status has significant
advantages both at the physical and network layer levels. In this paper, we
focus on two specific problems: (i) How to use the received mmWave signal
power information to predict whether there is a blockage in the near future or
not, and (ii) in case there is a blockage, how to use the received mmWave
signal power information to predict when that blockage will occur.

Problem 1: To formulate the presence of a blockage in the near future, let
$t\in\mathbb{Z}$ be the index of the discrete time instance, $x[t]$ be the
link status at the $t$th time instance, and let $S_{ob}=\\{\lvert
r[t+n]\rvert^{2}\\}_{n=-T_{o}+1}^{0}$ be the sequence of received signal power
for the observation interval of $T_{o}$ instances. Note that for simplicity of
expression, the subcarrier index $k$ is dropped from $r[t+n]$. Given a signal
power based observation sequence, we want to predict the occurrence of
blockage within a future time interval extending over $T_{P}$ instances. We
use $b_{T_{P}}$ to indicate whether there is a blockage occurrence within that
interval or not. More formally, $b_{T_{p}}$ is defined as follows:

$b_{T_{P}}=\begin{cases}0,&x[t+n^{\prime}]=0\quad\forall
n^{\prime}\in\\{1,\dots,T_{P}\\}\\\ 1,&\text{otherwise}\end{cases}$ (4)

where $1$ indicates the occurrence of a blockage and $0$ is the absence of
blockage. The goal of this problem is to predict $b_{T_{p}}$ with high
accuracy, i.e., high success probability
$\mathbb{P}(\hat{b}_{T_{p}}=b_{T_{p}}|S_{ob})$ where $\hat{b}_{T_{p}}$ is the
predicted link status. To that end, a prediction function $f_{\Theta}(S_{ob})$
parameterized by a set of parameters $\Theta$ could to be learned using a
machine learning algorithm such that it maximizes
$\mathbb{P}(\hat{b}_{T_{p}}=b_{T_{p}}|S_{ob})$.

Figure 2: Example of received signal power vs time. It also shows photos of
the location of the transmitter and blockage (from the receiver perspective).

Problem 2: Given signal power based observation sequence $S_{ob}$ and the
knowledge that there is a blockage in the future $T_{p}$ instances, the goal
is to predict $n^{\prime}$ at which $x[t+n^{\prime}]=1$. This represents the
exact instance at which the blockage occurs within a window of $T_{p}$
instances. Similar to Problem 1, the future instance is predicted by a
parameterized function $g_{\Theta}(\mathcal{S}_{ob})$ that could be learned
using a ML algorithm. The aim of the ML algorithm is to maximize the
prediction accuracy
$\mathbb{P}(\hat{n}^{\prime}=n^{\prime}|\mathcal{S}_{ob},b_{T_{p}}=1)$.

##  IV Moving Blockage Prediction using Recurrent Neural Networks

Deep neural networks have recently established themselves as a powerful
learning algorithm for many applications [9][10][11], and as such, they will
be the center of the proposed solution for the future-blockage prediction
problem.

###  IV-A Key Idea

Moving objects in a wireless communication environment contribute to changes
in the wireless channels, resulting in obvious fluctuations in received signal
power. This fluctuation pattern in received signal power is referred to as the
Pre-Blockage Signature [12]. Such pattern is illustrated in Fig. 2. It shows a
captured sequence of received power versus time instances, and the
corresponding photos show how far the blockages (the metal object in the
photos) is from the transmitter (the object circled with red in the photos).
The received power starts with smooth fluctuations (between the 1-st and 30-th
instances in Fig. 2), for the blockage is far from both the receiver and the
transmitter. However, as that blockage advances, its effect becomes clearer in
the received sequence, which could be seen in the red-shaded region of Fig. 2.
The sharper fluctuations is what we refer to as the pre-blockage signature,
and it could be utilized to identify incoming blockages.

Figure 3: The overall RNN architecture composed of a recurrent and a
prediction component

###  IV-B Deep Learning Model

Neural Network Architecture: Learning the pre-blockage signature from a
sequence of observed received signal power requires a neural network that
could process input data samples over time such as recurrent neural networks.
We design a Gated Recurrent Unit (GRU) network of $Q$-layers that takes in a
sequence of observed received signal powers (i.e., $S_{ob}$) and learns to
predict the link status $b_{T_{p}}$. Fig.3 depicts the schematic of such a
network. Each layer in the network comprises a sequence of GRUs equal to
$T_{o}$. The output of the last GRU of the last layer is fed to a Fully
Connected (FC) layer followed by either a classifier for problem 1 or a
regressor for problem 2. The classifier outputs a probability vector
($\hat{\mathbf{p}}$) of whether the link status is blocked or not in $T_{P}$
future time instances. For problem 2, the regressor outputs the predicted time
instance $\hat{n}^{\prime}$ indicating the time instance at when the blockage
will occur.

Pre-Processing: Before we input the data into our network, we need to pre-
process it to make it suitable for our model to learn, see [13] for more
information. We choose to standardize the inputs by subtracting the mean $\mu$
of the dataset and dividing by its standard deviation $\sigma$. Let
$\mathbf{A}\in\mathbb{R}^{U\times N}$ be the dataset matrix, where each row
represents a data point with a total number of data points of $U$. Data
standardization is done by computing:

$\mathbf{\hat{A}}_{u,n}=\frac{\mathbf{A}_{u,n}-\mu}{\sigma},$ (5)

where:

$\mu=\frac{1}{N\times U}\sum_{u=1}^{U}\sum_{n=1}^{N}\mathbf{A}_{u,n}$ (6)

$\sigma=\sqrt{\frac{1}{N\times
U}\sum_{u=1}^{U}\sum_{n=1}^{N}(\mathbf{A}_{u,n}-\mu})^{2}$ (7) $\forall
u\in\\{1,\dots,U\\}\text{and }n\in\\{1,\dots,N\\},$ (8)

Training loss: For problem 1 the future link-status prediction is posed as a
binary classification problem, in which the classifier attempts to determine
whether the link is blocked or not within the future time interval. As such,
the network training is performed with a cross entropy loss function [11]
computed over the outputs of the network:

$l_{\text{CH}}=\sum_{c=1}^{2}p_{c}\log{\hat{p}_{c}},$ (9)

where $\mathbf{p}=[p_{1},p_{2}]^{T}$ is the one-hot vector where a one hot
vector is a representation of categorical variables as binary vector, the
category with highest probability is encoded as 1 others are encoded as 0’s.
It takes one of two values: $[1,0]^{T}$ for the case when $b_{T_{p}}=0$ and
$[0,1]^{T}$ for the case when $b_{T_{p}}=1$, and $l_{\text{CH}}$ is the
training loss computed for one data point.

For problem 2, we pose the problem of predicting the blockage instance as a
regression problem. Our model tries to determine the exact time instance at
which the blockage occurs. We use Mean Square Error (MSE) loss as training
function. In formal terms, we aim to minimize the difference between the
predicted instance and groundtruth instance:

$l_{\text{MSE}}=(n^{\prime(u)}-\hat{n}^{\prime(u)})^{2}$ (10)

where $n^{\prime(u)}$ and $\hat{n}^{\prime(u)}$ are ground truth time instance
and predicted time instance, respectively.

##  V Experimental Setup and Evaluation Dataset

###  V-A Communication Scenario & Testbed

Figure 4: Experimental Scenario. The blockage moves along a trajectory between the transmitter (TX) and receiver (RX). TABLE I: Parameters for mmWave Communication Systems Name | Value
---|---
Carrier Frequency | 60GHz
Signal Bandwidth | 20MHz
number of subcarriers | 64
Horn Antenna Gain | 20dBi
Transmit Power | 30dBm

We build a mmWave communication system comprising of a transmitter with an
omnidirectional antenna communicating with a receiver with a 10-degree
beamwidth horn antenna. The operating parameters for our mmWave communication
system are shown in Table I. To simulate a moving blockage, we use a metal
cylinder with a height of 1 m, which can completely block the LOS link between
the transmitter and receiver. Then, we mount that cylinder onto a programmed
robot that moves along a pre-defined trajectory to simulate the moving
blockage. Fig. 4 depicts our system.

(a)

(b)

Figure 5: Top-view of the experimental scenario. (a) shows the trajectory of
the moving blockage and the relative positions of transmitter (TX), receiver
(RX) and moving blockage. (b) shows the rotated TX-RX setup.

We consider an indoor scenario where the transmitter and receiver are placed 8
m apart from each other and the robot moves between them in different
trajectories. To illustrate this further, let’s consider the coordinate system
in Fig. 4, where the x-axis extends along the LOS between the transmitter and
receiver and y-axis extends perpendicular to the LOS between the transmitter
and receiver. z-axis is perpendicular to the ground. The robot moves back and
forth on the y-axis to create multiple back and forward trajectories with a
spacing of 1 m. Then, the whole testbed is rotated with a small angle on the
x-y plane such that the background is slightly changed, and the robot is
programmed to do another back and forth cycle; see Fig. 5(a) and Fig. 5(b) for
an example of the robot motion in the rotated testbed. By continuously moving
and rotating, we collect a dataset of power sequences that encodes different
motion and propagation patterns.

###  V-B Dataset Generation

Figure 6: Sequence generation using sliding window. The top image shows the
received power for a raw sequence. The bottom image shows the correspond link
status.

Every robot trajectory provides us with a single received-power sequence,
which is manually annotated to create the link status labels. We use 1 to
indicate the time instances at which the LOS link is blocked and 0 otherwise.
We call the pair of received power and link status sequence a raw sequence
pair. Furthermore, since the number of raw sequence pairs we collected from
the experiment is limited, data augmentation is used to increase the dataset
size. Originally, we conducted 158 experiments and generated 158 sequence
pairs, i.e.
$\mathcal{S}_{1}=\\{(S_{d1},x_{d1})^{(u)}\\}_{u=1}^{U_{d}},U_{d}=158$, $d_{1}$
means the original input. We generate additional pairs by dropping samples at
rates 2, 3, and 4, which results in
$\mathcal{S}_{2}=\\{(S_{d2},x_{d2})^{(u)}\\}_{u=1}^{U_{d}}$,
$\mathcal{S}_{3}=\\{(S_{d3},x_{d3})^{(u)}\\}_{u=1}^{U_{d}}$,
$\mathcal{S}_{4}=\\{(S_{d4},x_{d4})^{(u)}\\}_{u=1}^{U_{d}}$, respectively. In
reality, this procedure means the blockage moves along the same trajectory at
2,3 or 4 times the original speed. Next we concatenate these sequences
together as
$\mathcal{S}=\mathcal{S}_{1}\bigcup\mathcal{S}_{2}\bigcup\mathcal{S}_{3}\bigcup\mathcal{S}_{4}=\\{(S,x)^{(u)}\\}_{u=1}^{U}$
where $U=4U_{d}$.

The method to generate the data points for Problem 1 and Problem 2 are as
follows.

Problem 1: A data point consists of an input received signal power sequence
$S_{ob}$ and the input label $b_{T_{p}}$. To generate the received signal
power sequence, we use a sliding window method, shown in Fig. 6. For example,
assume that current time is $t$, we generate $S_{ob}$ by extracting the
received power sequence from time instance $t-T_{o}+1$ to $t$ shown as red box
in the top subplot in Fig. 6. For input label generation, we first extract the
label sequence from time instance $t+1$ to $t+T_{P}$ shown as green box in
bottom subplot in Fig. 6, this gives us sequence
$\\{x[t+n]\\}_{n=t}^{t+T_{P}}$. Then we generate input labels $b_{t_{P}}$
using eqn. 4. Finally, we pair the received signal power sequences with input
labels, expressed as $S_{P1}=\\{(S_{ob},b_{T_{P}})^{(u)}\\}_{u=1}^{U_{P1}}$,
where $U_{P1}$ is the total samples we input to our model. In this problem,
$S_{P1}$ are mixed with two categories, non-transition sequence pairs where
$b_{T_{P}}=0$ and transition sequence pairs where $b_{T_{P}}=1$. To eliminate
the dataset bias, we keep the ratio of transition and non-transition sequence
pairs to 1:1.

Problem 2: A data point is represented by input received signal power sequence
$S_{ob}$ with the ground-truth time instance $n^{\prime}$ instead of the
label. So $S_{P2}=\\{(S_{ob},n^{\prime})^{(u)}\\}_{u=1}^{U_{P2}}$, where
$U_{P2}$ is the total number of samples that are input to our model. In this
problem, we only select transition-sequence pairs to be our input dataset.

In our experiments, we pick input sequence length ($T_{o}$) as 10 and
prediction interval ($T_{p}$) from 1 to 40. For each $T_{p}$, each raw
sequence pair generates one transition-sequence pair and multiple non-
transition sequence pairs. Totally, we get 632 transition-sequence pairs for
each $T_{p}$.

##  VI Experimental Results

###  VI-A Evaluation Metrics

Since problem 1 is a classification problem, we use Top-1 accuracy as our
evaluation metric. It is defined as:

$\text{Acc}_{\text{top-1}}=\frac{1}{U_{v1}}\sum_{u=1}^{U_{v1}}\mathbbm{1}(b_{T_{p}}^{(u)}=\hat{b}_{T_{p}}^{(u)}),$
(11)

where $\mathbbm{1}$ is the indicator function, $U_{v1}$ is total samples of
the validation set in problem 1, $b_{T_{p}}^{(u)}$ and $\hat{b}_{T_{p}}^{(u)}$
are, respectively, the target and predicted link status for a future interval
of $T_{P}$ instances.

Problem 2 is posed as a regression problem, and so we use Mean Absolute Error
(MAE) and its standard deviation to evaluate the quality of our model
predictions. The MAE is the mean absolute error between ground-truth value and
predicted value. For each prediction interval $T_{P}$, we calculate MAE across
the samples and standard deviation of these absolute errors.

$e^{(u)}_{T_{P}}=\lvert n^{\prime(u)}-\hat{n}^{\prime(u)}\rvert,\quad\forall
u\in\\{1,\dots,U_{v2}\\},$ (12)
$\bar{e}_{T_{P}}=\frac{1}{U_{v2}}\sum_{u=1}^{U_{v2}}\lvert
n^{\prime(u)}-\hat{n}^{\prime(u)}\rvert,$ (13)

$\text{std}_{T_{P}}={\sqrt{{\frac{1}{U_{v2}}}\sum_{i=1}^{U_{v2}}\left(e^{(u)}_{T_{P}}-{\bar{e}_{T_{P}}}\right)^{2}}},$
(14)

where, $e^{(u)}_{T_{P}}$ is the absolute error for $u$th sample,
$\bar{e}_{T_{P}}$ is the MAE, $\text{std}_{T_{P}}$ is the standard deviation
of absolute error, $U_{v2}$ is the total numbers of samples in validation set,
$n^{\prime(u)}$ and $\hat{n}^{\prime(u)}$ are target and predicted time
instances between current time and the time of blockage occurrence, assuming
prediction interval is $T_{P}$.

###  VI-B Network Training

We build the deep learning model described in Section IV using Pytorch. We
input 10 successive time-instances of received signal power and the
corresponding label for training. Our model consists of 1 RNN layer with 20
hidden states with a dropout rate = 0.2. The number of epochs is 1000. These
parameters are chosen based on empirical experiments. The detailed parameters
of our RNN for problem 1 and problem 2 are shown in Table II.

TABLE II: Parameters for Deep Learning Model | Value
---|---
Name | Problem 1 | Problem 2
Input sequence length | 10 | 10
Predicted future time steps | 1-40 | 1-40
Hidden state of RNN | 20 | 20
Output dimension | 2 | 1
Number of RNN layer | 1 | 1
Dropout rate | 0.2 | 0.2
Epoch | 1000 | 1000

###  VI-C Performance Evaluation

Problem 1: Fig. 7 plots the Top-1 accuracy as a function of the prediction
interval. We see that as the prediction interval increases the accuracy
decreases; the Top-1 accuracy decreases sharply at first, and then flattens
out. Our model achieves high accuracy when predicting the occurrence of
blockage in near future. i.e. we can achieve above 80% accuracy when
predicting the occurrence of blockage in the future 6 time instances. This is
because when the prediction interval is small, our training dataset contains a
large number of sequences with a pre-blockage signature, resulting in a high
ratio of these sequences in the training dataset. The prediction accuracy is
good until the prediction interval goes beyond 15 time instances. It finally
converges to the “random guess” performance as the prediction interval
approaches 40 time instances, i.e., approaching an accuracy of $50\%$. Such
general trend emphasizes value of pre-blockage signatures, which is commonly
clear when the blockage is very close to the transmitter and receiver.

Figure 7: Top-1 blockage prediction accuracy for different future prediction
intervals. Figure 8: Mean absolute error between the target (the exact time
instance where the blockage happens) and the prediction of this blockage time
instance for different future prediction intervals.

Problem 2: Fig. 8 plots the mean absolute error between prediction and ground-
truth as a function of prediction interval. For each prediction instance, we
show the standard deviation as an error bar. Given prediction interval below
10 time instance, our model can predict all the blockage transitions within 10
time instances under 1.9 with low volatility ($\pm$1.5). For a shorter
prediction intervals, our observation window more likely falls into a pre-
blockage signature resulting in lower prediction errors. However, as the
prediction interval increases, our observation window captures more sequences
without pre-blockage signature resulting in weaker prediction. However, when
the prediction interval is 40 time instance, we can still predict the exact
time of the blockage occurrence with the mean absolute error of under 8 time
instances.

##  VII Conclusion

In this paper, we proposed a deep-learning-based solution for the moving
blockage prediction in mmWave communication systems. Specifically, we
developed an RNN model to predict both the occurrence of the moving blockage
and the exact time when the LOS link is blocked. Simulation results on real
data showed that our model can achieve good performance for the moving
blockage prediction, which allows the user to be proactively handed over to
another base station without disconnecting the session with a high success
probability.

## References

  * [1] R. W. Heath, N. Gonzalez-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, “An overview of signal processing techniques for millimeter wave mimo systems,” _IEEE journal of selected topics in signal processing_ , vol. 10, no. 3, pp. 436–453, 2016.
  * [2] A. Alkhateeb, I. Beltagy, and S. Alex, “Machine learning for reliable mmwave systems: Blockage prediction and proactive handoff,” in _2018 IEEE Global Conference on Signal and Information Processing (GlobalSIP)_. IEEE, 2018, pp. 1055–1059.
  * [3] T. S. Rappaport, Y. Xing _et al._ , “Overview of millimeter wave communications for fifth-generation (5g) wireless networks—with a focus on propagation models,” _IEEE Transactions on Antennas and Propagation_ , vol. 65, no. 12, pp. 6213–6230, 2017.
  * [4] J.-S. Choi, W.-H. Lee, J.-H. Lee, J.-H. Lee, and S.-C. Kim, “Deep learning based nlos identification with commodity wlan devices,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 4, pp. 3295–3303, 2017.
  * [5] C. Huang, A. F. Molisch _et al._ , “Machine learning-enabled los/nlos identification for mimo system in dynamic environment,” _IEEE Transactions on Wireless Communications_ , 2020.
  * [6] M. Alrabeiah and A. Alkhateeb, “Deep learning for mmwave beam and blockage prediction using sub-6 ghz channels,” _IEEE Transactions on Communications_ , vol. 68, no. 9, pp. 5504–5518, 2020.
  * [7] G. Charan, M. Alrabeiah, and A. Alkhateeb, “Vision-aided dynamic blockage prediction for 6g wireless communication networks,” _arXiv preprint arXiv:2006.09902_ , 2020.
  * [8] M. Alrabeiah, A. Hredzak, and A. Alkhateeb, “Millimeter wave base stations with cameras: Vision-aided beam and blockage prediction,” in _2020 IEEE 91st Vehicular Technology Conference (VTC2020-Spring)_ , 2020, pp. 1–5.
  * [9] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [10] G. Huang, Z. Liu, and K. Q. Weinberger, “Densely connected convolutional networks,” _CoRR_ , vol. abs/1608.06993, 2016. [Online]. Available: http://arxiv.org/abs/1608.06993
  * [11] I. Goodfellow, Y. Bengio, and A. Courville, “Deep learning,” 2016, book in preparation for MIT Press. [Online]. Available: http://www.deeplearningbook.org
  * [12] A. Alkhateeb, S. Alex _et al._ , “Deep learning coordinated beamforming for highly-mobile millimeter wave systems,” _IEEE Access_ , vol. 6, pp. 37 328–37 348, 2018.
  * [13] Y. A. LeCun, L. Bottou, G. B. Orr, and K.-R. Müller, “Efficient backprop,” in _Neural networks: Tricks of the trade_. Springer, 2012, pp. 9–48.