Understanding Mass Transfer Directions via Data-Driven Models with Application to Mobile Phone Data

TL;DR

This paper presents a data-driven method combining Dynamic Mode Decomposition and Wasserstein distance to infer unknown mass transfer velocity fields from snapshots of mass distribution, applied to mobile phone data to analyze human travel flows.

Contribution

It introduces a novel algorithm that integrates DMD and Wasserstein distance to solve inverse mass transfer problems using real-world mobile phone data.

Findings

Successfully reconstructed human travel flows from mobile phone density data.

Demonstrated the effectiveness of the method in large populated areas.

Provided insights into mass transfer dynamics using real telecommunication data.

Abstract

The aim of this paper is to solve an inverse problem which regards a mass moving in a bounded domain. We assume that the mass moves following an unknown velocity field and that the evolution of the mass density can be described by partial differential equations (PDEs), which is also unknown. The input data of the problems are given by some snapshots of the mass distribution at certain times, while the sought output is the velocity field that drives the mass along its displacement. To this aim, we put in place an algorithm based on the combination of two methods: first, we use the Dynamic Mode Decomposition to create a mathematical model describing the mass transfer; second, we use the notion of Wasserstein distance (also known as earth mover's distance) to reconstruct the underlying velocity field that is responsible for the displacement. Finally, we consider a real-life application: the algorithm is employed to study the travel flows of people in large populated areas using, as input data, density profiles (i.e. the spatial distribution) of people in given areas at different time instances. This kind of data are provided by the Italian telecommunication company TIM and are derived by mobile phone usage.

Figures

Tables

Table 1. Table 1: Comparison of the errors defined in ( 3.5 ).

$N$	$\Delta x$	$\Delta t$	$E^G$	$E^L$
20	0.20	0.100	0.185	0.030
30	0.13	0.067	0.159	0.028
40	0.10	0.050	0.122	0.020

Table 2. Table 2: Computational time.

$N$	Global Exact	Global DMD	Local Exact	Local DMD
20	18.10 s	18.89 s	0.28 s	0.36 s
30	11 min	11 min	0.93 s	1.41 s
40	3 h 6 min	3 h 7 min	2.40 s	4.85 s

Table 3. Table 3: Comparison between dimensions of matrices and vectors for the two methods.

Algoritm	Vectors Dimension	Matrix Dimension
Global	$N^2$	$2N^3$
Local

Equations

$$\frac{d\hat{\mathbf{y}}}{dt}=\hat{\mathbf{A}}\hat{\mathbf{y}}$$

$$\hat{\mathbf{y}}(t)=\sum _{i=1}^r{\beta }_i{\mathbf{\psi }}_i\exp ({\omega }_{i}t),$$

$$\mathbf{Y}=\left[\begin{array}{cccc}\text{vline}& \text{vline}& & \text{vline}\\ \mathbf{y}(t_0)& \mathbf{y}(t_1)& \cdots & \mathbf{y}(t_{n_{\text{f}}-1})\\ \text{vline}& \text{vline}& & \text{vline}\end{array}\right],{\mathbf{Y}}^{\prime }=\left[\begin{array}{cccc}\text{vline}& \text{vline}& & \text{vline}\\ \mathbf{y}(t_1)& \mathbf{y}(t_2)& \cdots & \mathbf{y}(t_{n_{\text{f}}})\\ \text{vline}& \text{vline}& & \text{vline}\end{array}\right]$$

$${\mathbf{Y}}^{\prime }=\mathbf{AY},$$

$$\mathbf{A}={\mathbf{Y}}^{\prime }{\mathbf{Y}}^{†}.$$

$$\underset{{\mathbf{\phi }}_1,\dots ,{\mathbf{\phi }}_r\in {\mathbb{R}}^n}{\min }\sum _{j=0}^{n_{\text{f}}-1}{\Vert \mathbf{y}(t_j)-\sum _{i=1}^r⟨\mathbf{y}(t_j),{\mathbf{\phi }}_i⟩{\mathbf{\phi }}_i\Vert }^2\text{mbox}suchthat⟨{\mathbf{\phi }}_i,{\mathbf{\phi }}_j⟩={\delta }_{ij},$$

$$\mathbf{Y}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T.$$

$$\hat{\mathbf{A}}={\mathbf{U}}^T\mathbf{AU}={\mathbf{U}}^T{\mathbf{Y}}^{\prime }\mathbf{V}{\mathbf{\Sigma }}^{-1},$$

$$\hat{\mathbf{A}}\mathbf{W}=\mathbf{W\Lambda },$$

$${\mathbf{\Psi }}^{\text{mbox}DMD}={\mathbf{Y}}^{\prime }\mathbf{V}{\mathbf{\Sigma }}^{-1}\mathbf{W}.$$

$$W_{p}(\rho^{\textup{i}},\rho^{\textup{f}})=\bigg{(}\min_{T\in\mathcal{T}}\int_{\Omega}c(\bm{\xi},T(\bm{\xi}))^{p}\,\rho^{\textup{i}}(\bm{\xi})d\bm{\xi}\bigg{)}^{\frac{1}{p}}$$

$$\mathcal{T}:=\Biggr{\{}T\colon\Omega\to\Omega\,:\,\int\displaylimits_{B}\rho^{\textup{f}}(\bm{\xi})d\bm{\xi}=\int\displaylimits_{\{\bm{\xi}:T(\bm{\xi})\in B\}}\rho^{\textup{i}}(\bm{\xi})d\bm{\xi},\quad\forall\,B\subset\Omega\Biggr{\}}$$

$$\text{minimize }\mathcal{H}:=\sum _{j,k=1}^N{c}_{jk}{x}_{jk}$$

$$\sum _k{x}_{jk}=m_j^{\text{i}}\forall j,\sum _j{x}_{jk}=m_k^{\text{f}}\forall k\text{and}{x}_{jk}\ge 0.$$

$$\mathbf{x}=(x_{11},x_{12},\dots ,x_{1N},x_{21},\dots ,x_{2N},\dots ,x_{N1},\dots ,x_{NN}{)}^{\text{T}},$$

$$\mathbf{c}=(c_{11},c_{12},\dots ,c_{1N},c_{21},\dots ,c_{2N},\dots ,c_{N1},\dots ,c_{NN}{)}^{\text{T}},$$

$$\mathbf{b}=(m_1^{\text{i}},\dots ,m_N^{\text{i}},m_1^{\text{f}},\dots ,m_N^{\text{f}}{)}^{\text{T}},$$

$$\mathbf{M}=\left[\begin{array}{ccccc}{\text{1}}_N& 0& 0& \dots & 0\\ 0& {\text{1}}_N& 0& \dots & 0\\ 0& 0& {\text{1}}_N& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & {\text{1}}_N\\ {\mathbf{I}}_N& {\mathbf{I}}_N& {\mathbf{I}}_N& {\mathbf{I}}_N& {\mathbf{I}}_N\end{array}\right],$$

$$\min$$

$$\dim \mathbf{x}=\dim \mathbf{c}=N^2,\dim \mathbf{b}=2N,\dim \mathbf{M}=2N^3.$$

$$V_{\text{max}}\delta t\le \Delta x$$

$$\tilde{\mathbf{M}}=\left[\begin{array}{ccccc}{\text{1}}_{r_1}& 0& 0& \dots & 0\\ 0& {\text{1}}_{r_2}& 0& \dots & 0\\ 0& 0& {\text{1}}_{r_3}& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & {\text{1}}_{r_N}\\ & & \tilde{\mathbf{D}}& & \end{array}\right],$$

$$\min$$

$$\dim \tilde{\mathbf{x}}=\dim \tilde{\mathbf{c}}\le dN,\dim \tilde{\mathbf{b}}=2N,\dim \tilde{\mathbf{M}}\le 2d{N}^2.$$

$$\{\begin{array}{ll}{\partial }_{t}u(\mathbf{x},t)+\mathbf{v}\cdot \nabla u(\mathbf{x},t)=0& \mathbf{x}\in \Omega ,t\in [0,T]\\ u(\mathbf{x},t)=0& \mathbf{x}\in \partial \Omega \\ u(\mathbf{x},0)=u_0(\mathbf{x})& \mathbf{x}\in \Omega ,\end{array}$$

$$c({\mathbf{\xi }}_1,{\mathbf{\xi }}_2)=\Vert {\mathbf{\xi }}_1-{\mathbf{\xi }}_2{\Vert }_{{\mathbb{R}}^2}^{1+\epsilon },$$

$$E^G:=\frac{{\Vert {\mathbf{X}}_E^{\ast }-{\mathbf{X}}_D^{\ast }\Vert }_F}{{\Vert {\mathbf{X}}_E^{\ast }\Vert }_F},E^L:=\frac{{\Vert {\tilde{\mathbf{X}}}_E^{\ast }-{\tilde{\mathbf{X}}}_D^{\ast }\Vert }_F}{{\Vert {\tilde{\mathbf{X}}}_E^{\ast }\Vert }_F},$$

$$E=\frac{{\Vert {\mathbf{P}}^{\mathrm{data}}-{\tilde{\mathbf{P}}}^{\mathrm{DMD}}\Vert }_F}{{\Vert {\mathbf{P}}^{\mathrm{data}}\Vert }_F}.$$

$$\tilde{N}=\underset{\text{corners}}{\underbrace{4\cdot 4}}+\underset{\text{boundaries}}{\underbrace{6\cdot 2(N_R+N_C-4)}}+\underset{\text{internal nodes}}{\underbrace{9(N-4-2(N_R+N_C-4))}}<9N.$$

$$c_{jk}=\{\begin{array}{ll}0& \text{if j=k or j and k are not adjacent}\\ {\ell }_x& \text{if j and k are horizontally adjacent}\\ {\ell }_y& \text{if j and k are vertically adjacent}\\ \sqrt{{\ell }_x^2+{\ell }_y^2}& \text{if j and k are diagonally adjacent.}\end{array}$$

$$\{\begin{array}{ll}{\partial }_{t}y-\epsilon \Delta y+y({\partial }_{x_1}y+{\partial }_{x_2}y)=0& \mathbf{x}\in \Omega ,t\in [0,T]\\ y(\mathbf{x},t)=0& \mathbf{x}\in \partial \Omega \\ y(\mathbf{x},0)=\sin (\pi x_1)\sin (\pi x_2)& \mathbf{x}\in \Omega ,\end{array}$$

Full text

Keywords.

Data-driven methods, dynamic mode decomposition, Wasserstein distance, earth mover’s distance, cellular data, presence data.

Mathematics Subject Classification.

37C10, 35R30, 76D55.

Introduction

In this paper we aim at solving an inverse problem which regards a mass moving in a bounded domain with finite velocity. We assume that the mass moves following an unknown velocity field and that the evolution of the mass density can be described by an unknown PDE. The input data of the problem are given by some snapshots of the mass distribution at certain times, while the sought information is the velocity field that drives the mass along its displacement.

The basic idea can be summarized as follows: given two snapshots of the mass distribution at two instants of time, we want to understand where each portion of the mass (which is assumed to be conserved from one instant of time to the other) is transported from/to, i.e. how the first spatial distribution is rearranged in the second one. The goal is pursued by computing a numerical approximation of the Wasserstein distance (also known as earth mover’s distance or Mallows distance) between the two consecutive density profiles, specifying a suitable cost function which measures the “energy” consumed by the system for moving forward. The computation of the Wasserstein distance gives, as a by-product, the minimum-cost mass flow from the first to the second configuration, i.e. how the mass distributes in space and time.

Despite the good preliminary results obtained by the above described Wasserstein-based approach [2, 19], the algorithm is found to be excessively expensive both in terms of CPU time and memory requirements. This fact strongly restricts the applicability of the method. To fix this, in this paper we propose to couple the method with the Dynamic Mode Decomposition (DMD): a data-driven technique that takes in input the snapshots of the mass distribution and returns an analytical approximation of the dynamics underlying the mass transfer. More precisely, it provides a system of ODEs which describes the evolution of the mass in any point of the domain. Solving the ODEs, we are able to recover the mass distribution at any time, thus increasing at will the number of available snapshots or, analogously, decreasing at will the time frame between them. Controlling the time frame between two consecutive snapshots is the key to simplify the computation of the Wasserstein distance and makes the computation of the flows feasible even on large domains.

Finally, a real-world application of the proposed methodology is illustrated. We are interested in inferring activity-based human mobility flows from mobile phone data. We assume that mobile devices are not singularly tracked, but their logs are aggregated in order to obtain the total number of users in a given area. In this way we get the density profiles (i.e. the spatial distribution) of people in a given area at various instants of time. The dataset we have at our disposal is provided by the Italian telecommunication company TIM. The time frame between two consecutive snapshots is 15 minutes.

As before, the goal is to “assign a direction” to the presence data. In fact, the mere representation of time varying density of people clearly differentiate attractive from repulsive or neutral areas but does not provide any information about the directions of flows of people. In other words, we are interested in a “where-from-where-to” type of information, which reveals travel flows and patterns of people providing a sort of origin-destination matrix.

Paper organization

The paper is organized as follows. In Section 2 we present the DMD method and the Wasserstein distance. In Section 3 we show how to couple those methods to obtain an efficient algorithm. In Section 4 we apply the proposed approach to real-life data. In Section 5 we draw our conclusions. Finally, in Appendix A we apply the DMD method to the viscous Burgers’ equation.

Mathematical background

In this section we recall the building blocks of the methodology proposed in the paper, namely the DMD method used to build a data-driven model, and the Wasserstein distance to determine the transport map driving the moving mass.

DMD method

DMD is a data-driven method capable of providing accurate assessments of the spatio-temporal coherent structures in a given complex system, or short-time future estimates of such a system. Although a complete list of references for DMD goes beyond the scopes of this work, we would like to mention [3, 5, 6, 13, 15]. In the current work we use the DMD algorithm as in [16].

To begin, we suppose to have a set of data $\mathcal{X}=\{{\bf y}(t_{0}),\ldots,{\bf y}(t_{n_{\textup{f}}})\}$ for some time instances $\{t_{j}\}_{j=0}^{n_{\textup{f}}}$ with ${\bf y}(t_{j})\in\mathbb{R}^{N},j=0,\ldots,n_{\textup{f}}$ and $\Delta t=t_{j+1}-t_{j}$ for $j=0,\dots,n_{\textup{f}}-1$. The goal of the method is to build a mathematical model upon the dataset $\mathcal{X}\in\mathbb{R}^{N\times(n_{\textup{f}}+1)}$. The DMD procedure thus constructs the approximate linear evolution $\widehat{\bf y}(t)$ for the dataset $\mathcal{X}$ exploiting its low-rank structure:

[TABLE]

where $\widehat{\bf A}\in\mathbb{R}^{N\times N}$ is unknown, $\widehat{\bf y}(0)=\widehat{\bf y}_{0}$, and the solution has the form

[TABLE]

where $r<N,{\boldsymbol{\psi}}_{i}$ and $\omega_{i}$ are the eigenvectors and eigenvalues of the unknown matrix $\widehat{\bf A}$. The coefficients $\beta_{i}$ of the vector ${\boldsymbol{\beta}}$ can be determined from the initial data. For example, at $t=t_{0}$ we have ${\bf y}(t_{0})={\bf y}_{0}$ so that (2.2) gives ${\boldsymbol{\beta}}={\bf\Psi}^{\dagger}{\bf y}_{0}$, where ${\bf\Psi}$ is a matrix comprised of the DMD modes $\boldsymbol{\psi}_{i}$ and ${\dagger}$ denotes the Moore-Penrose pseudo-inverse. To compute the matrix $\widehat{\bf A}$, we first split the dataset into two snapshot matrices

[TABLE]

and suppose the following linear relation hold true:

[TABLE]

where ${\bf A}:=\exp{(\widehat{\bf A}\Delta t}$).

Specifically, we assume that ${\bf y}(t_{j})$ is an initial condition to obtain ${\bf y}(t_{j+1})$, i.e. its corresponding output after some prescribed evolution time $\Delta t>0$. Thus, the DMD method computes the best linear operator ${\bf A}$ relating to the matrices above:

[TABLE]

We refer to ${\bf Y}$ and ${\bf Y}^{\prime}$ as input and output snapshot matrices respectively.

The DMD algorithm aims at optimally constructing the matrix ${\bf A}$ so that the error between the true and approximate solution is small in a least-square sense, i.e. $\|{\bf y}(t)-\widehat{\bf y}(t)\|\ll 1$. Of course, the optimality of the approximation holds only over the sampling window where ${\bf A}$ is constructed, but the approximate solution can be used to make future state predictions, and to decompose the dynamics into various time-scales.

The matrix ${\bf A}$ is often highly ill-conditioned and when the state dimension $n$ is large, the aforementioned matrix may be even intractable to analyze directly. Instead, DMD circumvents the eigen-decomposition of ${\bf A}$ by considering a low rank representation in terms of a matrix $\widehat{\bf A}$ projected with the Proper Orthogonal Decomposition (POD). Although the description of the POD method goes beyond the scopes of this paper, we recall that the POD projection solves the following optimization problem

[TABLE]

where $\{{\boldsymbol{\varphi}}_{i}\}_{i=1}^{r}$ are the POD projectors. The solution of the optimization problem (2.6) is obtained by means of a Singular Value Decomposition (SVD) of the dataset ${\bf Y}$, where the first $r\ll N$ columns of the left singular eigenvectors are the required POD basis. We refer to [18] for a complete description of the POD method. We also mention that POD method is equivalent to Principal Component Analysis (PCA) or Karhunen-Loéve expansion in other contexts (see e.g. [8, 9]).

The exact DMD algorithm proceeds as follows [16]: first, we collect data ${\bf Y,Y^{\prime}}$ as in (2.3) and compute the reduced, or economy, singular value decomposition of Y

[TABLE]

We note that the use of the economy SVD is suggested since the matrices ${\bf Y},{\bf Y}^{\prime}\in\mathbb{R}^{n\times n_{\textup{f}}}$ with $N\gg n_{\textup{f}}$ and that the economy ${\bf U}\in\mathbb{R}^{N\times r}$ is sufficient to provide the same approximation of the regular SVD given the limited amount of snapshots.

Furthermore, the DMD typically makes use of low-rank structure so that the total number of modes, $r\ll N$, allows for dimensionality reduction of the dynamical system. Then, we compute the least-squares fit ${\bf A}$ that satisfies ${\bf Y}^{\prime}={\bf AY}$ and project onto the POD modes ${\bf U}$. Specifically, the Moore-Penrose pseudo-inverse of ${\bf Y}$ allows us to compute ${\bf A}={\bf Y}^{\prime}{\bf Y}^{\dagger}$, where the Moore-Penrose algorithm provides the least-square fitting procedure. In terms of its low-rank projection, this yields

[TABLE]

and then, we compute the eigen-decomposition of ${\bf\widehat{A}}\in\mathbb{R}^{r\times r}$

[TABLE]

where ${\bf\Lambda}$ are the DMD eigenvalues. Finally, the DMD modes ${\bf\Psi}^{\mbox{\tiny DMD}}$ are given by

[TABLE]

The algorithm is summarized in Algorithm 1.

In Appendix A we provide a numerical experiment to show the effectiveness of the DMD method. We compute the data from a nonlinear PDE, i.e. the viscous Burgers’ equation.

Wasserstein distance and optimal mass transfer problem

The notion of Wasserstein distance is strictly related to the Monge–Kantorovich optimal mass transfer problem [17], which can be easily explained as follows: given a sandpile with mass distribution $\rho^{\textup{i}}$ and a pit with equal volume and mass distribution $\rho^{\textup{f}}$, find a way to minimize the cost of transporting sand into the pit. The cost for moving mass depends on both the distance from the point of origin to the point of arrival and the amount of mass is moved along that path. We are interested in minimizing this cost by finding the optimal path to transport the mass from the initial to the final configuration.

Given two density functions $\rho^{\textup{i}},\rho^{\textup{f}}:\Omega\to\mathbb{R}$ for some bounded $\Omega\subset\mathbb{R}^{n}$ such that $\int_{\mathbb{R}^{n}}\rho^{\textup{i}}=\int_{\mathbb{R}^{n}}\rho^{\textup{f}}$, we define the $L^{p}$-Wasserstein distance between $\rho^{\textup{i}}$ and $\rho^{\textup{f}}$ as

[TABLE]

where

[TABLE]

and $c:\Omega\times\Omega\to\mathbb{R}$ is a given cost function, which defines the cost of transferring a unit mass between any two points in $\Omega$. Note that $\mathcal{T}$ is the set of all possible maps which transfer the mass from one configuration to the other.

It is important to note here that we are not really interested in the actual value of the Wasserstein distance $W_{p}$, instead we look for the optimal map $T^{*}$ which realizes the $\arg\min$ in (2.8), and represents the paths along which the mass is transferred.

Numerical approximation of the Wasserstein distance

A direct numerical approximation of definition (2.8) is unfeasible, but a discrete approach is still possible. Indeed, we can resort to classical problems (see Hitchcock’s paper [7]) and methods (see e.g., [12, Sec. 6.4.1] and [14, Chap. 19]) to recast the original mass transfer problem in the framework of linear programming (LP). We also refer to [4] for a recent application of this methodology to traffic flow problems.

The idea consists in approximating the set $\Omega$ with a structured grid with $N$ square cells $C_{1},\ldots,C_{N}$, as it is commonly done for the numerical approximation of PDEs. We denote by $\Delta x$ the length of each side of the cells. Then, we define a graph $\mathcal{G}$ whose nodes coincide with the centers of the $N$ cells. Graph’s edges are defined in such a way that each node is directly connected with each other, including itself.

Introducing a numerical error (controlled by the parameter $\Delta x$), we are allowed to assume that $\forall j=1,\ldots,N$ all the mass distributed in the cell $C_{j}$ is concentrated in its center, i.e. in a node of the graph. We come up with an initial mass $m^{\textup{i}}_{j}:=\int_{C_{j}}\rho^{\textup{i}}dx$ and a final mass $m^{\textup{f}}_{j}:=\int_{C_{j}}\rho^{\textup{f}}dx$, for $j=1,\dots,N$, distributed on the graph nodes. Now, we simply aim at optimally rearranging the first mass into the second one moving it among the graph’s nodes.

We denote by $c_{jk}$ the cost to transfer a unit mass from node $j$ to node $k$, and by $x_{jk}$ the (unknown) mass moving from node $j$ to node $k$. The problem is then formulated as

[TABLE]

subject to

[TABLE]

Defining

[TABLE]

and the $2N\times N^{2}$ matrix

[TABLE]

where ${\bf I}_{N}$ is the $N\times N$ identity matrix and ${\boldsymbol{\text{1}}_{N}}=\displaystyle\underbrace{(1\ 1\dots\ 1)}_{N\text{ times}}$, our problem is written as a standard LP problem:

[TABLE]

The result of the algorithm is a vector ${\bf x}^{*}:=\arg\min{\bf c}^{\textup{T}}{\bf x}$ whose elements $x^{*}_{jk}$ represent how much mass moves from node $j$ to node $k$ employing the minimum-cost mass rearrangement.

Remark 2.1**.**

The dimension of the LP problem (2.9) is given by the dimensions of the matrix/vectors involved:

[TABLE]

Remark 2.2**.**

Hereafter we will refer to problem (2.9) as global, in order to stress the fact that it is possible to move the mass from and to any node of the graph.

Coupling DMD and optimal mass transfer problem

In this section we describe how we can drastically reduce the size of the LP problem (2.9) by using the DMD method. The resulting algorithm will allow us to study the mass transfer problem on large domains.

Main idea

The large size of the LP problem (2.9), stressed in Remark 2.1, mainly comes from the fact that the mass is allowed to be transferred from any graph node to any other graph node. However physical constraints prevent this to happening: assuming that the maximal velocity of the mass is $V_{\textup{max}}$ and denoting by $\Delta t$ the time frame from one snapshot to the following one, the maximal distance travelled is $V_{\textup{max}}\Delta t$.

In order to reduce the size of the LP problem (2.9), we restrict the set of all possible movements trajectories. The ideal time frame $\delta t$ of the snapshots would be the one which guarantees that the CFL-like condition (see, e.g., [11, Chapter 14])

[TABLE]

holds true. Indeed, under condition (3.1) the mass is allowed to move only towards the adjacent cells or not move at all.

Consider now a generic set of mass distributions $\mathcal{X}=\{{\bf m}(t_{0}),\ldots,{\bf m}(t_{n_{\textup{f}}})\}$, where $t_{j}=j\Delta t$, $j=0,\dots,n_{\textup{f}}$, and $\Delta t$ is the time frame of the snapshots of the set $\mathcal{X}$. A priori, the time step $\Delta t$ does not necessarily satisfy condition (3.1). In particular, a large distance in time between the snapshots means that there are not enough information for reducing the set of possible movements. However, with the DMD we are able to reconstruct the state of the system at any time instance, even if it is not provided in the original dataset. The coupling of DMD and of optimal mass transfer problem is done in this order:

1. i)

We first fix $\delta t<\Delta t$ such that it satisfies condition (3.1) and then we reconstruct the solution for each $\widetilde{t}_{j}=j\delta t$ via DMD, using Algorithm 1. The new set of snapshots we work with is $\widetilde{\mathcal{X}}=\{{\bf m}(\widetilde{t}_{0}),\ldots,{\bf m}(\widetilde{t}_{\widetilde{n_{\textup{f}}}})\}$, where ${\bf m}(\widetilde{t}_{j})$ is computed from (2.2), $j=0,\dots,\widetilde{n_{\textup{f}}}$. We observe that $\widetilde{n_{\textup{f}}}>n_{\textup{f}}$. 2. ii)

We recover the flows from the new set $\widetilde{\mathcal{X}}$ by means of an approximation of Wasserstein distance similar to the one done in Section 2.3, but with reduced size, as described in detail below. Note that, since $\widetilde{n_{\textup{f}}}>n_{\textup{f}}$, we have to solve more LP problems with respect to the global approach, but, despite this, we will get advantages by using this new approach.

Let us denote by $d$ the maximum number of neighbors per cell. The new unknown $\widetilde{x}_{jk}$, corresponding to the mass to be moved from node $j$ to node $k$, is defined only if $j$ and $k$ are adjacent or $j$ is equal to $k$. Analogously we define the cost function $\widetilde{c}_{jk}$. We denote by $\bf\widetilde{x}$ the vector of the unknowns and by $\bf\widetilde{c}$ the vector associated to the cost function. We introduce the vector ${\bf s}$ whose components are indexes of nodes. The first components are the indexes of the nodes adjacent to the node 1, then those adjacent to the node 2 and so until the node $N$. Vectors ${\bf\widetilde{x}}$ and ${\bf\widetilde{c}}$ are ordered similarly to ${\bf s}$. Since the mass can move only towards a maximum of $d$ nodes, the dimensions of $\bf\widetilde{x}$ and $\bf\widetilde{c}$ are lower or equal than $dN$. We introduce the matrix

[TABLE]

where ${\bf\text{1}}_{r_{i}}=\displaystyle\underbrace{(1\ 1\dots\ 1)}_{r_{i}\text{ times}}$, with $r_{i}$ the number of adjacent nodes for each node $i$, $r_{i}\leq d$, and ${\bf\widetilde{D}}_{ij}=1$ if the $j$-th element of ${\bf s}$ is equal to $i$, otherwise ${\bf\widetilde{D}}_{ij}=0$. Defining the vector ${\bf\widetilde{b}}=(m^{\textup{i}}_{1},\dots,m^{\textup{i}}_{N},m^{\textup{f}}_{1},\dots,m^{\textup{f}}_{N})^{\textup{T}}$ the LP problem becomes

[TABLE]

Remark 3.1**.**

The dimension of the LP problem (3.2) is given by the dimensions of the matrix/vectors involved:

[TABLE]

Remark 3.2**.**

Hereafter we will refer to problem (3.2) as local, in order to stress the fact that it is possible to move the mass from and to adjacent nodes of the graph only.

A toy example for the complete algorithm: the advection equation

In this test we propose an example for the complete algorithm described in Section 3.1. Let us consider the advection equation:

[TABLE]

with ${\bf x}=(x_{1},x_{2})$, $\Omega=[-2,2]\times[-2,2]$, $u_{0}({\bf x})=\max(0.5-x_{1}^{2}-x_{2}^{2},0)$ and constant velocity ${\bf v}=(0.5,0.5)$. It is well known that the analytical solution of (3.3) is $u({\bf x},t)=u_{0}({\bf x}-{\bf v}t)$ provided that we set $T$ small enough to have inactive boundary conditions, as it is the case here. Hereafter we denote by reference solution the analytical solution $u$ of (3.3).

Starting from some snapshots of the solution, we aim at reconstructing the velocity field driving the dynamics, i.e. the vector ${\bf v}=(0.5,0.5)$. In doing this, we will also compare the global and the local problem in terms of CPU time, see (2.9) and (3.2) respectively.

We choose a time step $\delta t$, and we collect snapshots with a larger temporal step size $\Delta t=\kappa\delta t$ with $\kappa>1$ and we reconstruct the solution via the DMD method. We note that here the snapshots are computed from the analytical solution of (3.3). In particular, we work on a grid $40\times 40$ and we choose $T=2$, $\Delta x=0.1$, $\delta t=0.05$ and $\Delta t=2\,\delta t$. We observe that, since $V_{\textup{max}}=\left\lVert\bf v\right\rVert=\sqrt{2}/2$, this choice of $\delta t$ fullfills the condition (3.1). The number of snapshots is 40 and the rank $r$ in Algorithm 1 for the DMD reconstruction is 20. Moreover we identify the nodes of the graph $\mathcal{G}$ of the numerical approximation of Wasserstein distance with the cells of the grid.

Choice of the cost function

Since the cost function is used to figure out the “price to pay” for moving the mass from a node of the graph to another one, the most intuitive definition of $c$ is the Euclidean distance in $\mathbb{R}^{2}$. This choice was proved to be unsuitable for the global problem, see [2]. Indeed, since the global algorithm allows any movements between the nodes of the graph, using the Euclidean distance for the function $c$ we loose the uniqueness of the optimal transfer map ${T}^{*}$. To see this let us assume that we have to move three unit masses one to the right. In the picture on the left of Figure 2 we move the three unit masses of one to the right while in the picture on the right we move only the first mass of three to the right. The Wasserstein distance between the two configuration is clearly the same and equal to three.

The solution proposed in [2] to fix this issue was to force the minimization algorithm to select primarly the small movements penalizing the large ones. To get this, the cost function was defined as

[TABLE]

where $\boldsymbol{\xi}_{1}$ and $\boldsymbol{\xi}_{2}$ are the coordinates of the nodes of the graph and $\varepsilon>0$ is a small parameter which accounts for the penalization.

In Figure 3 we show the level sets of the solution to (3.3) together with some arrows indicating the reconstructed velocity field $\bf v$. More precisely, on the left column of figure we show the results obtained with the local algorithm (3.2), on the center with the global algorithm with Euclidean distance and on the right with the global algorithm with the cost function defined in (3.4) (with $\varepsilon=0.1$). Similarly, the panels on the top show the results obtained with the reference solution whereas the panels on the bottom those obtained with the DMD solution.

The local algorithm and the global one with the correction in (3.4) are able to reconstruct the velocity field ${\bf v}=(0.5,0.5)$ with accurate approximation. The global algorithm with Euclidean distance, instead, is less precise, since the optimal flow does not correspond to the velocity field. Moreover, the algorithm based on the DMD reconstruction of the solution introduces small oscillations in the solution to (3.3). This is expected in such hyperbolic problems since the decay of the singular values of the dataset is very slow and the initial condition is non-smooth. However, such oscillations do not have an effect in the reconstruction of the flow.

To further validate our approach we compute the numerical error of the proposed algorithm. Since the results obtained with the global algorithm (2.9) with the Euclidean distance as cost function are the least accurate, we focus only on the other two approaches. For each time step $t_{n}=n\delta t$ we define ${\bf x}^{*}_{E}(t_{n})$ and ${\bf\widetilde{x}}^{*}_{E}(t_{n})$ as the solutions to the LP problems (2.9) and (3.2) respectively at time $t_{n}$, when the known terms ${\bf b}(t_{n})$ and ${\bf\widetilde{b}}(t_{n})$ are chosen as the reference solution to (3.3) at time $t_{n}$. Analogously, ${\bf x}^{*}_{D}$ and ${\bf\widetilde{x}}^{*}_{D}$ are the solution when the known terms are obtained with the DMD method. The vectors ${\bf x}^{*}_{E}(t_{n})$, ${\bf\widetilde{x}}^{*}_{E}(t^{n})$, ${\bf x}^{*}_{D}(t^{n})$ and ${\bf\widetilde{x}}^{*}_{D}(t_{n})$, for $n=1,\dots,\lfloor\frac{T}{\delta t}\rfloor+1$, are finally collected as the columns of the matrices ${\bf X}^{*}_{E}$, ${\bf\widetilde{X}}^{*}_{E}$, ${\bf X}^{*}_{D}$ and ${\bf\widetilde{X}}^{*}_{D}$ respectively. We define the errors

[TABLE]

where $\left\lVert\cdot\right\rVert_{F}$ is the Frobenius norm. In Table 1 we compare the errors defined in (3.5) obtained from the simulations on a grid with $N\times N$ nodes, for $N=20,\,30$ and $40$. As we can see from the table, increasing the number of nodes we reduce the space step $\Delta x=4/N$ and the time step $\Delta t=\Delta x/2$ and we obtain the decrease of the error between the reference solution and the DMD reconstruction. Moreover, the error obtained with the local algorithm is significantly smaller than the one obtained with the global approach.

Finally, in Table 2 we compare the computational time between the global approach (2.9), with the cost function as in (3.4), and the local approach (3.2) with respect to the nodes of the graph. We observe that the local algorithm is always faster than the global one. The difference between the two approaches becomes more relevant when we refine the grid by increasing the number of nodes, and thus the number of time steps. Specifically, for a grid $40\times 40$, the local algorithm required a few seconds whereas the global one more than three hours.

Application to real mobile phone data

In this section we focus on a specific application of the proposed approach. The real dataset gives information about the spatial distribution of people in a large populated area. The goal is to understand the travel flows of people, focusing in particular on recurring patterns and daily flows of commuters.

Dataset

The Italian telecommunication company TIM provides estimates of mobile phones presence in a given area in raster form: the area under analysis is split into a number of elementary territory units (ETUs) of the same size (about 130 $\times$ 140 m2 in urban areas). The estimation algorithm does not singularly recognize users and does not track them using GPS. It simply counts the number of phone attached to network nodes and, knowing the location and radio coverage of the nodes, estimates the number of TIM users within each ETU at any time. TIM has now a market share of 30% with about 29.4 million mobile lines in Italy (AGCOM, Osservatorio sulle comunicazioni 2/2017).

The data we considered refer to the area of the province of Milan (Italy), which is divided in 198,779 ETUs, distributed in a rectangular grid 389 $\times$ 511. Data span six months (February, March and April 2016 and 2017). The entire period is divided into time intervals of 15 minutes, therefore we have 96 data per day per ETU in total. In Figure 4 we graphically represent presence data at a fixed time. We observe that the peak of presence is located in correspondence of Milan city area.

Figure 5 shows the presences in the province of Milan in a typical working day in the left panel. The curve in the image decreases during the night, it increases in the day-time and decreases again in the evening. These variations are due to two main reasons: first, the arrival to and departure from Milan’s province of visitors and commuters. Second, the fact that when a mobile phone is switched off or is not communicating for more than six hours, its localization is lost. The presence value that most represents the population of the province is observed around 9 pm., when an equilibrium between traveling and phone usage is reached. This value changes between working days and weekends, but it is always in the order of $1.3\times 10^{6}$.

On the right panel of Figure 5 we show the trend of presence data during April 2017. We can observe a cyclical behavior: in the working days the number of presences in the province is significantly higher than during the weekends. It is interesting to note the presence of two low-density periods on April 15-18 and on April 22-26, 2017, determined respectively by the Easter and the long weekend for the Italy’s Liberation Day holiday.

DMD approach on TIM data

As explained in Section 2.1, we can reconstruct the presence data in each cell at any time. We denote by ${\bf m}(t^{\textup{i}})\in\mathbb{R}^{N}$ the vector containing the number of people present in the $N$ cells at a certain quarter of an hour and by ${\bf m}(t^{\textup{i}}+15\min)\in\mathbb{R}^{N}$ the same quantities at the consecutive time step. Applying the DMD method we are able to calculate ${\bf m}(t)\in\mathbb{R}^{N}$ for any $t\in[t^{\textup{i}},t^{\textup{i}}+15\min]$, see (2.2). We also note that DMD can be applied to higher dimensional dataset through randomized methods as explained in [1].

To validate the DMD approach on the TIM dataset we study the daily error using only half of the data at our disposal in the DMD reconstruction. More precisely, we denote by ${\bf P^{\mathrm{data}}}$ the matrix whose columns correspond to the real data of presences stored every $15$ minutes. Then, we reconstruct the data of presences every minute with the DMD algorithm, using snapshots every 30 minutes. In other words, we exploit only one column out of two of ${\bf P^{\mathrm{data}}}$ to build the matrix ${\bf{\widetilde{P}}}^{DMD}$, which collects the reconstructed data every minute. Since a day contains 96 quarters of an hour and 1440 minutes, ${\bf P^{\mathrm{data}}}$ is a matrix $96\times N$, while ${\bf{\widetilde{P}}}^{DMD}$ is a matrix $1440\times N$. Finally, to compare the real data with the DMD reconstruction, we extract from ${\bf{\widetilde{P}}}^{DMD}$ the rows corresponding to the original interval of $15$ minutes into the matrix ${\bf\widetilde{P}}^{\mathrm{DMD}}$, of dimensions $96\times N$, and then we define the error as:

[TABLE]

In Figure 6 we show the daily error (4.1) for an entire month of data in the whole area of the province of Milan. The daily error is of order $10^{-2}$, which certifies the accuracy of the DMD method.

Understanding human mobility flows

Following the approach described in Section 3, we define the graph $\mathcal{G}$ by exploiting the subdivision of the area of the province of Milan (Italy) into ETUs; we identify the $N$ nodes of the graph with the corresponding center of the ETUs, ordered from the left to the right and from the top to the bottom. The result is a rectangular graph $\mathcal{G}$ divided into $N_{R}$ rows and $N_{C}$ columns ($N=N_{R}\times N_{C}$). The mass $m_{j}(t_{n})$ is defined as the average number of presences in the node $j$ at time $t_{n}=n\,15\min$.

Let us assume that $V_{\textup{max}}$ in (3.1) is equal to 50 km/h. Since the dimensions of the ETUs is around 150 m, to apply the DMD we fix the new time step $\delta t$ equal to 10 seconds. With this choice we assume that people can move only towards the eight adjacent nodes of the rectangular graph, or not move at all. We observe that the mass in the nodes on the four corners of the graph can move only towards four directions (adjacent nodes or no movement), while the mass along the boundaries can move only towards six directions. In this way, the total number of possible movements between the cells is given by

[TABLE]

The vectors $\bf\widetilde{x}$ and $\bf\widetilde{c}$ associated to the unknown moving mass and to the cost function have length $\widetilde{N}$, while the matrix $\bf\widetilde{M}$ of the LP problem (3.2) has dimension $2N\times\widetilde{N}$. In Table 3, we compare the dimensions of the vectors and the matrices of the two different LP problems: it is clear that the computational time to solve problem (3.2) is significantly reduced respect to problem (2.9).

Choice of the cost function $c$

Since the ETUs are rectangles of length $\ell_{x}=130$ m and $\ell_{y}=140$ m, we define the cost function $c_{jk}$ for the local algorithm as follows:

[TABLE]

For the global approach, we use the cost function defined in (3.4) with $\varepsilon=0.1$.

To sum up, in order to solve the mass transfer problem for a whole day using snapshots taken every 15 minutes (real data) we have to solve 96 global LP-problems (2.9), whereas with the DMD algorithm we have to solve 8640 local LP-problems (3.2). As we will see in the following section, despite the larger number of LP problems, the local approach is more convenient than the global one.

Remark 4.1**.**

The Wasserstein distance is defined between two distributions of equal mass (see (2.8)). In our case the conservation of mass between two consecutive snapshots is not guaranteed. Let us consider a couple of snapshots with different total mass $\sum_{j}m_{j}^{0}\neq\sum_{j}m_{j}^{1}$. To correctly apply the algorithms for the identification of flows, we compute the mass in excess between the two snapshots and then we uniformly distribute it in all the nodes of the graph with lower mass. A more sophisticated approach could be the one suggested in [10], where a definition of Wasserstein distance between two distributions with different mass is proposed.

Numerical results

In this section we show the results obtained with the local algorithm (3.2) to study the flows of commuters and the influence of great events on human mobility. In both cases the rank $r$ in step 2 of Algorithm 1, used for the construction of the DMD solution, is 95. The flows are represented by arrows; we draw only those which correspond to the most significant movements, and we associate a darker color to the arrows corresponding to a larger movement of people. For graphical purposes, in the following plots we aggregate 6 time steps $\delta t$ to show 1-minute mass transport.

4.4.1 Flows of commuters

In the following simulations we consider the area of the Province of Milan during a generic working day. Milan is one of the biggest Italian city and it attracts many workers from outside. The city of Milan is located in the right part of the Province, therefore we mainly see movements from/to the left part of the analyzed area. In the top panel of Figure 7 we show the morning flows of a working day: we clearly see that the arrows are directed towards the city of Milan. In particular, we zoom on the arrows which overlap the roads heading to Milan. In the bottom panel of Figure 7 we show the opposite phenomenon: in the evening people go away from work to come back home inside the Province of Milan.

CPU times

Considering data for a 6 hours frame on an area of $144\times 240$ ETUs the local algorithm requires 144 hours of CPU time and works with 360 snapshots. The global approach (2.9) is not able to analyze such an area, since the matrix $\bf M$ in (2.9) has a computationally unmanageable dimension.

4.4.2 Flows influenced by a great event

In this test we show how the algorithm is able to capture the way a big event influences human flows. The event we have considered is the exhibition of the Salone del Mobile, held every April at Fiera Milano exhibition center in Rho, near Milan. We analyze a square area of $31\times 31$ ETUs centered around Fiera Milano. In the left panel of Figure 8 we show the morning fluxes directed to the exhibition area whereas in the right panel we show the evening flows directed from the exhibition area to the outside.

CPU times

For a simulation of 18 hours, from 06:00 to 23:45, on an area of $31\times 31$ ETUs, the local algorithm requires 6 minutes of CPU time while the global approach requires 30 minutes. Furthermore, we observe that the local algorithm works with $1065$ snapshots of data, whereas the global one with $266$ snapshots.

Conclusions

In this paper we have proposed an efficient method to solve an inverse mass transfer problem, consisting in recovering the dynamics underlying the mass displacement. The proposed algorithm prevents to handle large displacements of the mass, thus saving CPU time and memory allocation with respect to other recently proposed solutions.

The application of the methodology to a real dataset describing the movement of people was also investigated. It is useful to note here that the applicability of the Wasserstein-based methodology was not at all obvious, since it is based on assumptions which are not totally fullfilled. Indeed, the choice of the cost function $c$ does not take into account the fact that people move mainly along roads, are stopped by obstacles, buildings, etc., and in general are not free to move in all directions. Moreover, and most important, the computation of the Wasserstein distance stems from a global optimization problem in which the mass is considered as a whole. In other words, the optimal transport map $T^{*}$ is found by minimizing the cost of the displacement of the whole crowd, without any distinction among people, and with no regards about individual targets. Despite this, the results we have obtained (see especially those in Figure 7 and 8) are exactly as one can expect, meaning that the method is, overall, robust enough to work well even if the constitutive assumptions are not totally fullfilled.

Appendix A The DMD for the viscous Burgers’ equation

Here we propose a numerical experiment using data generated by the 2D viscous Burgers’ equation:

[TABLE]

with ${\bf x}=(x_{1},x_{2})$, $\Omega=[0,1]\times[0,1]$, $T=1$ and $\varepsilon=0.01$. The numerical approximation of equation (A.1) has been carried out by a finite difference method on a square grid with $N=1600$ nodes, using an adaptive Runge-Kutta method for the time integration with temporal step size $\Delta t>0$. We refer to [11] for further details on the numerical method. In what follows, since we do not know the analytical solution of equation, we consider as reference solution of (A.1) its numerical approximation with $\Delta t=0.0125$.

To show how DMD reconstructs the model we consider as dataset $\mathcal{X}$ the numerical approximation of (A.1) for the following temporal step sizes $\Delta t=\{0.1,0.05,0.025,0.0125\}$ with snapshots matrices of dimension $N\times 11,N\times 21,N\times 41,N\times 81,$ respectively. Once the snapshots have been computed we apply Algorithm 1 to reconstruct the solution corresponding to a fixed $\delta t=0.0125.$ Here, we aim at testing the capability of the DMD method to approximate the dataset for the missing information.

In Figure 9, we show the reference solution of (A.1) in the top row and the absolute difference between the reference solution and the reconstructed one with DMD in the bottom row for two time instances. We note that the dataset for the DMD reconstruction corresponds to the solution of (A.1) with $\Delta t=0.1$.

Then, in Figure 10 we show the relative error between the DMD approximation and the reference solution of (A.1) as follows:

[TABLE]

where ${\bf y}^{R}(t)$ is the reference solution and ${\bf y}^{DMD}(t)$ is its DMD approximation. The notation ${\bf Y}^{R},{\bf Y}^{DMD}$ refers to the matrices where each column contains the reference solution and its DMD approximation, respectively, for different time instances.

It is clear that the more information on the system, e.g. amount of snapshots, the better error behavior for the DMD solution (see left panel in Figure 10). However, we can already obtain accurate solutions with a poor dataset. In the middle panel of Figure 10 we show the error of the model for each time instance for the dataset corresponding to $\Delta t=\{0.1,0.05,0.025\}$ using $\widehat{E}^{DMD}(t)$ defined in (A.2). We note that DMD behaves very well for time instances provided in the snapshot set. This is to show that DMD is able to reconstruct accurately the solution for missing information. In this example the rank $r$ in the DMD method is 9. Finally, we would like to emphasize that, although the DMD model is linear, it works even for nonlinear problem as, e.g., (A.1). For the sake of completeness we also show in the right panel of Figure 10 the decay of error $E^{DMD}$ in (A.2) modifying the rank $r$ in the DMD approximation. As expected the error decays as the rank increases. The singular values of snapshots set are also in the picture.