Estimating ensemble flows on a hidden Markov chain

TL;DR

This paper introduces a new method to estimate the evolution of an ensemble of indistinguishable agents on a hidden Markov chain using aggregate data, extending optimal transport concepts to finite state spaces.

Contribution

It develops a convex maximum likelihood framework and a fast algorithm for estimating ensemble flows on hidden Markov chains from aggregate data.

Findings

Effective tracking of ensemble flows demonstrated in numerical examples.

Framework extends optimal transport to finite state hidden Markov models.

Convex formulation enables efficient computation.

Abstract

We propose a new framework to estimate the evolution of an ensemble of indistinguishable agents on a hidden Markov chain using only aggregate output data. This work can be viewed as an extension of the recent developments in optimal mass transport and Schr\"odinger bridges to the finite state space hidden Markov chain setting. The flow of the ensemble is estimated by solving a maximum likelihood problem, which has a convex formulation at the infinite-particle limit, and we develop a fast numerical algorithm for it. We illustrate in two numerical examples how this framework can be used to track the flow of identical and indistinguishable dynamical systems.

Equations

$${\rho }_1(x_1)={\int }_{{\mathbb{R}}^n}{q}_{ϵ}(0,x_0,1,x_1){\rho }_0(x_0)d{x}_0,$$

$$q_{ϵ}(s,x,t,y)=\frac{1}{{\left(2\pi (t-s)ϵ\right)}^{n/2}}\exp \left(-\frac{\Vert x-y{\Vert }^2}{2(t-s)ϵ}\right),$$

$$H(\mathcal{P}\mid \mathcal{W})=\int \log \left(\frac{d\mathcal{P}}{d\mathcal{W}}\right)d\mathcal{P}$$

$$H(p|q):=\sum _i{p}_i\log \left(\frac{p_i}{q_i}\right)$$

$$\underset{M\in {\mathbb{R}}^{n\times n}}{\mathrm{minimize}}$$

$$M1={\mu }_0,M^{T}1={\mu }_1,$$

$$P_{{\mu }_0,A}(M)=\prod _{i=1}^n\left(\left(\genfrac{}{}{0ex}{}{({\mu }_0{)}_i}{m_{i1},m_{i2},\dots ,m_{in}}\right)\prod _{j=1}^n{a}_{ij}^{m_{ij}}\right),$$

$$\begin{split}&\left|\log\left(P_{\mu_{0}^{(N)}\!,\,A}(M^{(N)})\right)+H\left(M^{(N)}\,\big{|}\,{\rm diag}(\mu_{0}^{(N)})A\right)\right|\\ &\phantom{\leq}\leq C\log(N).\end{split}$$

$$\frac{1}{N}{\mu }_0^{(N)}\to {\bar{\mu }}_0\text{mbox}and\frac{1}{N}{M}^{(N)}\to \bar{M}$$

$$\frac{1}{N}\log\left(P_{\mu_{0}^{(N)}\!,\,A}(M^{(N)})\right)\to-H\left(\bar{M}\,\big{|}\,{\rm diag}(\bar{\mu}_{0})A\right)$$

$$P_{{\mu }_0,A}(M)\sim e^{-H\left(M|\mathrm{diag}({\mu }_0)A\right)}.$$

$$\underset{M\in {\mathbb{R}}^{n\times n}}{\mathrm{minimize}}$$

$$M1={\mu }_0,M^{T}1={\mu }_1.$$

$$C=-ϵ\log (\mathrm{diag}({\mu }_0)A),$$

$$H(M\mid A)-H({\mu }_0\mid 1)$$

$$\underset{M_{[1:T]},{\mu }_{[1:T-1]}}{\mathrm{minimize}}$$

$$M_{t}1={\mu }_{t-1},M_t^{T}1={\mu }_t,$$

$$\text{for }t=1,\dots ,T.$$

$$\underset{{\bar{M}}_{[1:T]},{\mu }_{[1:T-1]}}{\mathrm{minimize}}$$

$${\bar{M}}_{t}1=1,{\mu }_t={\bar{M}}_t^T{\mu }_{t-1},$$

$$\text{for }t=1,\dots ,T.$$

$$M_{t}1={\mu }_{t-1},M_t^{T}1={\mu }_t,$$

$$D_{t}1={\mu }_t,D_t^{T}1={\Phi }_t,\text{for }t=1,\dots ,T.$$

$$\underset{M_{[1:T]},D_{[1:T]},{\mu }_{[1:T]}}{\mathrm{maximize}}\prod _{t=1}^T{P}_{{\mu }_{t-1},A}(M_t)P_{{\mu }_t,B}(D_t)$$

$$\log \left(\prod _{t=1}^T{P}_{{\mu }_{t-1},A}(M_t)P_{{\mu }_t,B}(D_t)\right)$$

$$-\sum_{t=1}^{T}\Big{(}H(M_{t}\mid{\rm diag}(\mu_{t-1})A)+H(D_{t}\mid{\rm diag}(\mu_{t})B)\Big{)}$$

$$\underset{M_{[1:T]},D_{[1:T]},{\mu }_{[1:T]}}{\mathrm{minimize}}$$

$$+\sum _{t=1}^{T}H(D_t|\mathrm{diag}({\mu }_t)B)$$

$$M_{t}1={\mu }_{t-1},M_t^{T}1={\mu }_t$$

$$D_{t}1={\mu }_t,D_t^{T}1={\Phi }_t$$

$$\text{for }t=1,\dots ,T.$$

$$y_1=A^T({\mu }_0\odot u_1),$$

$$y_t=A^T\left(y_{t-1}\odot (B{v}_{t-1})\right),t=2,\dots ,T$$

$$w_T=B{v}_T$$

$$w_t=\left(B{v}_t\right)\odot (A{w}_{t+1}),t=1,\dots ,T-1.$$

$${\mu }_t=\mathrm{diag}(w_t)A^T\left({\mu }_{t-1}./(A{w}_t)\right),t=1,\dots ,T.$$

$$M_t$$

Full text

I Introduction

State tracking of a set of agents is an important issue in many areas, e.g., target tracking, see [4] and references therein. In this case, one is often interested in tracking one single or a set of multiple distinct targets. However, in many applications information for each agent may not be available, e.g., if the population is too large to track every single agent, as in many biological systems, or due to data privacy [24]. In this work, we thus consider tracking the evolution of a finite ensemble of indistinguishable agents. Based on reduced and incomplete measurements of the whole population at different time points, we aim to recover an estimate of the discrete-time flow of the ensemble. Related state estimation problems for a continuum of agents and in continuous time have been considered in [10, 35] (see also [6]).

In this work, we use a hidden Markov model (HMM) to describe the particle flows and aggregate observations, similar to [2], and seek the most likely paths that the agents have taken. These paths are found by maximizing the log-likelihood function of the flow, subject to the constraint that the flow matches the given measurements. This gives rise to a convex maximum entropy type optimization problem, and we derive an efficient algorithm for solving it.

The problem of finding the most likely path for the evolution of a distribution is related to a discrete Schrödinger bridge problem [28].

Schrödinger’s thought experiment [31] has indirectly given rise to the concept of reciprocal processes [22, 25, 5], which connects this work to tracking of moving objects using reciprocal processes [18, 34, 32]. However, as mentioned before, we consider estimating the flow of an ensemble rather than single target tracking.

The outline of the paper is as follows: Section II presents background material, in particular on HMMs, Schrödinger bridges, and optimal mass transport. In section III we derive the maximum likelihood problem for a Markov chain with a known initial and final distribution, and relate it to prior work on the Schrödinger bridge problem [28]. In Section IV, which contains the main contribution, we extend this maximum entropy framework to HMMs with indirect and noisy observations. Moreover, we derive the corresponding maximum likelihood problem, and develop a fast iterative algorithm for solving it. The method is demonstrated on two examples in section V, and section VI contains conclusions and future directions. Some proofs are deferred to the appendix for improved readability.

II Background

II-A Notation

By $./$, $\odot$, $\log(\cdot)$, and $\exp(\cdot)$ we denote elementwise division, multiplication, logarithm, and exponential of matrices and vectors. Moreover, by $\rm supp(\cdot)$ we denote the support of a matrix, i.e., the non-zero elements.

II-B Hidden Markov chains

In this work, we consider hidden Markov models for stochastic modeling of a group of indistinguishable agents/particles. For an introduction to HMMs, see, e.g., [29, 21]. An HMM is a structure that consists of two stochastic processes. The first part is a Markov chain that evolves over a hidden set of states $X=\left\{X_{1},X_{2},\dots,X_{n}\right\}$ and is used to model the unobserved, underlying state of the system. We denote the state at time $t$ by $q_{t}$. The stochastic state transitions are encoded in the state transition matrix $A=\left[a_{ij}\right]_{i,j=1}^{n}$, where $a_{ij}=P(q_{t+1}=X_{j}|q_{t}=X_{i})$. The second part is an observation process providing partial and noisy information of the underlying process; here we use the observation symbols $Y=\left\{Y_{1},Y_{2},\dots,Y_{m}\right\}$. Moreover, the observation process is also Markovian with respect to the underlying state in the hidden Markov chain, i.e., the observation probabilities can be summarized in a matrix $B\in{\mathbb{R}}^{n\times m}$ with elements $b_{jk}=P(Y_{k}\text{ at }t|q_{t}=X_{j})$.

II-C Schrödinger bridges and large deviations

In the early 1930s, Schrödinger discussed the problem of determining the evolution of particles between two observed distributions [31]. Assuming a cloud of independent Brownian particles is observed at time instance $t=0$, the expected distribution at $t=1$ would be described by

[TABLE]

where $q_{\varepsilon}$ is the Brownian transition probability kernel

[TABLE]

and where the parameter $\epsilon$ denotes a diffusion coefficient. Schrödinger studied the problem where the observed particle distribution differs from the expected distribution (1). The most likely particle evolution connecting, hence bridging, the two marginals is called the Schrödinger bridge.

The Schrödinger bridge problem was later formulated in the context of large deviation theory [19, Sec. II.1.3], the study of rare events in the sense of deviations from the law of large numbers [16, 15]. As the number of trials (or particles) goes to infinity, the probability of such rare events approaches zero. Large deviation theory studies the rate of this decay, which can often be characterized by the exponential of a so called rate function.

Modeling the particle evolutions as independent identically distributed random variables on path space, a Schrödinger bridge is a probability measure ${\mathcal{P}}$ on path space that is most likely to describe the rare event of observing the two particle distributions. Such a measure is obtained by minimizing the corresponding rate function, which turns out to be the relative entropy with respect to the underlying probability law of the Brownian motion. In other words, ${\mathcal{P}}$ is the measure that is “most similar” to the Wiener measure ${\mathcal{W}}$ in the sense that it minimizes the relative entropy [26]

[TABLE]

over all probability measures that are absolutely continuous with respect to ${\mathcal{W}}$ and have the given particle distributions as marginals. The Schrödinger bridge can be constructed from the solution to a certain system of equations, called the Schrödinger system. A space and time discrete Schrödinger bridge problem for Markov chains is analysed in [28, 20, 8].

II-D Optimal mass transport

Another recently established connection of Schrödinger bridges is to the problem of optimal mass transport (OMT) [9, 7, 26, 27]. As the diffusion coefficient $\epsilon$ in (1) approaches [math], the solution to the Schrödinger bridge problem tends to the solution to a corresponding optimal mass transport problem[26]. Moreover, the Schrödinger bridge formulation is a regularization of OMT, as it is strictly convex and therefore guarantees a unique solution.

We introduce a discretized formulation of the OMT problem. For an extensive discussion of OMT see, e.g., [33]. Consider a discretization $\left\{x_{1},\dots,x_{n}\right\}$ of a compact space $X$ and two distributions $\mu_{0},\mu_{1}\in{\mathbb{R}}^{n}$ defined on this discretization. Given a cost matrix $C=[c_{ij}]_{i,j=1}^{n}$, where $c_{ij}$ denotes the cost of transporting a unit mass from point $x_{i}$ to $x_{j}$, we seek a transport plan $M=[m_{ij}]_{i,j=1}^{n}$, where $m_{ij}$ denotes the amount of mass being transported from $x_{i}$ to $x_{j}$, that minimizes the total transportation cost $\text{tr}\left(C^{T}M\right)$ between the two distributions, i.e., the transport plan is required to satisfy $M{\bf 1}=\mu_{0}$ and $M^{T}{\bf 1}=\mu_{1}$, where ${\bf 1}$ denotes an $n\times 1$-vector of ones.

Solving this linear program is computationally expensive for large $n$. It was therefore proposed to regularize the problem by introducing a Kullback-Leibler divergence term (sometimes called entropy term) to the objective [12].

Definition 1

Let $p$ and $q$ be two nonnegative vectors or matrices of the same dimension. The Kullback-Leibler (KL) divergence between $p$ from $q$ is defined as

[TABLE]

where $0\log 0$ is defined to be [math]. Note that $H(p|q)$ is jointly convex over $p,q$. See, e.g., [11] for more properties and interpretation of the KL divergence.

The discretized and regularized OMT problem then reads

[TABLE]

where $\epsilon>0$ is a regularization parameter and ${\bf 1}_{n\times n}$ denotes an $n\times n$-matrix of ones. The solution to this problem may be found by Sinkhorn iterations, which correspond to the fixed point iteration for the Schrödinger system in [7].

III Particle dynamics over a Markov chain

Consider a cloud of $N$ particles, where each particle is evolving according to a Markov chain as described in Section II-B. Let the vectors $\mu_{t}\in{\mathbb{N}}^{n}$ describe the particle distributions at time $t\in\{0,1\}$, where the $i$-th element $(\mu_{t})_{i}$ denotes the number of particles in state $X_{i}$ at time $t$. In analogy to the OMT framework, we define the mass transfer matrix $M=[m_{ij}]_{i,j=1}^{n}$, where $m_{ij}$ denotes the number of particles that transit from state $X_{i}$ to state $X_{j}$. Note that the mass transport matrix satisfies $M{\bf 1}=\mu_{0}$ and $M^{T}{\bf 1}=\mu_{1}$.

The state transition matrix $A=[a_{ij}]_{i,j=1}^{n}$ contains the particle transition probabilities. Thus, given the initial state $\mu_{0}$, the probability of a mass transfer matrix $M$ is

[TABLE]

where $\binom{\cdot}{\cdot,\ldots,\cdot}$ denotes a multinomial coefficient. The expected distribution at time $t=1$ is then given by $E(\mu_{1}|\mu_{0})=A^{T}\mu_{0}$. If $\mu_{1}$ is observed to be different from $A^{T}\mu_{0}$ a discrete version of the Schrödinger bridge problem can be solved (see Section II-C). That is, to find the matrix $M$ that maximizes $P_{\mu_{0},A}(M)$ subject to that the constraint $M^{T}{\bf 1}=\mu_{1}$ on the final marginal is satisfied.

If the number of particles is large, then the log-likelihood of (4) can be approximated in terms of a KL divergence.

Proposition 1

Given $A$, let $\mu^{(N)}_{0}\in{\mathbb{N}}^{n}$ be a sequence of distributions with $N$ particles, and $M^{(N)}\in{\mathbb{N}}^{n\times n}$ be a sequence of mass transfer matrices such that $M^{(N)}\mathbf{1}=\mu^{(N)}_{0}$ and $\rm supp(M^{(N)})\subseteq\rm supp({\rm diag}(\mu_{0}^{(N)})A)$. Then there exists a constant $C>0$ such that for all $N$ it holds that

[TABLE]

Proof:

See appendix. ∎

Proposition 1 implies that for sequences $\mu^{(N)}_{0}$ and $M^{(N)}$ satisfying the assumptions, if

[TABLE]

as $N\to\infty$, then

[TABLE]

as $N\to\infty$. This means that the KL divergence approximates the log-likelihood of $P_{\mu_{0},A}(M)$ with increasing accuracy as the number of particles increases. We write this informally as

[TABLE]

In terms of large deviation theory, we thus interpret $H(\cdot|{\rm diag}(\mu_{0})A)$ as the rate function for $P_{\mu_{0},A}(\cdot)$. In fact, Proposition 1 can also be derived from a large deviation principle (see, e.g., [15, Ch. 2.1.1]).

For systems with many particles, we may therefore formulate the problem of finding the most likely mass transfer matrix $M$ between distributions $\mu_{0}$ and $\mu_{1}$ with underlying state transition matrix $A$ as the convex optimization problem

[TABLE]

Remark 1

Let $A$ and $\mu_{0}$ be strictly positive. With the cost matrix

[TABLE]

the entropy regularized OMT problem (3) is equivalent to problem (5). Note that entropy regularized OMT problems have previously been solved by formulating them in terms of KL-projection problems [1].

III-A Connection to Schrödinger bridges

We note that given the prior distribution $\mu_{0}$, the objective in (5) may be written as

[TABLE]

where the second term is constant. Hence, if we associate $A$ and $M$ with the measures $d{\mathcal{W}}$ and $d{\mathcal{P}}$ in (2), the problem in Proposition 1 relates to a discretized Schrödinger bridge. Our problem formulation indeed corresponds to the time and space discrete Schrödinger bridge from [28]. To see this, consider a Markov chain of length $T$. Using Proposition 1, knowing the marginals $\mu_{0}$ and $\mu_{T}$, we can find the most likely evolution of the particles between them as the solution to

[TABLE]

Note that for a nonnegative matrix $M_{t}$ and strictly positive marginal $\mu_{t-1}$, the first constraint asserts that there is a row-stochastic matrix $\bar{M}_{t}$ such that $M_{t}={\rm diag}(\mu_{t-1})\bar{M}_{t}$. Plugging this expression for the matrices $M_{t}$ into (6) gives

[TABLE]

Here $A_{i\cdot}$ denotes the $i$-th row of $A$. This is precisely the formulation of a Schrödinger bridge over a Markov chain from [28, eq. (24)] with time invariant transition probabilities. In [28] it is shown that a unique solution to a corresponding Schrödinger system exists if $\mu_{T}$ is a strictly positive distribution and all elements are strictly positive in the matrix $A$ raised to the power $T$. The solution to the Schrödinger system may be obtained by a fixed point iteration [20], which is linked to the Sinkhorn iterations for entropy regularized OMT problems, cf. Section II-D.

We note that the optimization problem (7) is non-convex and will thus work with the formulation (6) in the remaining part of this article.

IV Particle dynamics over hidden Markov chain

In this section, we extend our framework to the setting of a hidden Markov chain. The initial marginal distribution is assumed to be known. In case the hidden states are linked to the observations by a deterministic linear mapping, they may be estimated in a similar fashion as in [17]. Here instead, we consider the non-deterministic case where the available observations emerge from the hidden distributions through an observation probability matrix $B\in{\mathbb{R}}^{n\times m}$.

Equivalently to the mass transfer plans $M$, define the observation matrix $D\in{\mathbb{N}}^{n\times m}$, with entries $d_{jk}$ denoting the number of particles that are in hidden state $X_{j}$ and observed in state $Y_{k}$. Given a hidden state $\mu$, the probability for any observation matrix $D$ is given by $P_{\mu,B}(D)$ as defined in (4). Hence, the large deviation result in Proposition 1 holds for $D$ with rate function $H(\cdot|{\rm diag}(\mu)B)$.

Given an initial distribution $\mu_{0}\in{\mathbb{N}}^{n}$ and a set of measurements $\Phi_{1},\dots,\Phi_{T}\in{\mathbb{N}}^{m}$, we seek the most likely set of matrices $M_{1},\dots,M_{T}$ and $D_{1},\dots,D_{T}$ such that for some set of hidden distributions $\mu_{1},\dots,\mu_{T}$ it holds that

[TABLE]

This model is illustrated in Figure 1. The maximum likelihood solution is obtained by solving the optimization problem

[TABLE]

subject to (8). From Proposition 1 it follows that

[TABLE]

can be approximated by

[TABLE]

when the number of particles is large. We thus estimate the matrices $M_{1},\dots,M_{T}$, $D_{1},\dots,D_{T}$ and the hidden states $\mu_{1},\dots,\mu_{T}$ by maximizing (9) subject to the constraints (8), i.e., by solving

[TABLE]

Remark 2

The modeling assumptions leading to this optimization problem require knowledge of the initial distribution as well as the transition and observation probabilities for the hidden Markov model. However, in practice these are typically not known exactly. In the examples in Section V, we illustrate that the estimation is accurate even when there are significant model errors. The generalization of the proposed method to the case where the initial distribution is not known will be discussed in a forthcoming paper.

IV-A Computational method

In this section we develop a numerical method to solve problem (10). To this end, recall from Remark 1 the connection between KL-minimization problems and entropy regularized OMT problems. The latter can be efficiently solved by Sinkhorn iterations[12], which in turn are equivalent to a block coordinate ascent in a dual problem[23]. Motivated by this, we choose to follow a similar approach.

Proposition 2

Let $u_{1}\in{\mathbb{R}}^{n}$ and $v_{t}\in{\mathbb{R}}^{m}$, for $t=1,\ldots,T$, be positive initial values, and iterate the following steps:

1. (1)

$u_{1}=e{\bf 1}./(Aw_{1})$** 2. (2)

$v_{t}=e\Phi_{t}./\left(B^{T}(y_{t}\odot(Aw_{t+1}))\right)$* for $t=1,\dots,T$,*

where in each step and for each $t$ in step (2), the vectors $y_{t}$ and $w_{t}$ are recursively defined as

[TABLE]

and

[TABLE]

In the limit point of the iteration, the estimates for the hidden marginals are then recursively constructed, starting from the known $\mu_{0}$, as

[TABLE]

Furthermore, the corresponding mass transfer matrices are given by

[TABLE]

where

[TABLE]

for $t=1,\dots,T.$

Proof:

See appendix. ∎

It is worth noting that intermediate results of the vectors $y_{t}$ and $w_{t}$ may be stored, such that the update of $u_{1}$ requires only one matrix-vector multiplications with $A$, and the update of $v_{t}$, for any $t=1,\dots,T$, involves two multiplications with $B$ and one with $A$. One iteration sweep, i.e. one update of $u_{1}$ and the set $v_{t}$, for $t=1,\dots,T$, thus requires $\mathcal{O}(Tn\max(n,m))$ operations.

V Simulations

V-A Particle dynamics

Consider a cloud of 1000 particles evolving from an initial distribution $\mu_{0}\in{\mathbb{R}}^{n}$ with $n=100$ states. The particles transition matrix is given by $\tilde{A}\in{\mathbb{R}}^{n\times n}$ with elements

[TABLE]

which corresponds to a discretization of a normal distribution ${\mathcal{N}}(1,0.5)$, and thus induces a drift on the dynamics of the cloud. The true dynamics of the particles are assumed to be unknown and instead modeled by a transition matrix $A\in{\mathbb{R}}^{n\times n}$ with elements

[TABLE]

At each time instance, the particles are observed in $m=5$ bins, where the observation probability matrix $B\in{\mathbb{R}}^{n\times m}$ has elements

[TABLE]

We estimate the flow of the particles and hidden particle distributions solving problem (10) for $T=50$ time instances with the method proposed in Proposition 2. One estimate is formed using the true initial distribution $\mu_{0}$ as a prior distribution, and for a second estimate we use a uniform prior.

Figure 2 shows the true hidden particle cloud, the corresponding observations, and the two estimates. With full information of the initial states available, the proposed method provides a good estimate of the hidden states despite discrepancies between the true and assumed transition matrices $\tilde{A}$ and $A$. In the case of no prior information, i.e., the prior distribution is set to be uniform, we see that the estimate converges to the estimate with true prior within a few time steps. This indicates that the proposed method is robust to modeling uncertainties and lack of information in the initial state.

V-B Tracking ensembles over a network

In this example we consider the problem of tracking a number of indistinguishable agents over a network given measurements from sensors that are distributed around the network. This is inspired by [13], where an HMM is used to estimate the flow of a crowd in an urban environment based on observations generated when cell phones connect to Wi-Fi sensors. The environment is modeled as a network of nodes and arcs, where the arcs represent walking paths in the area and the nodes are the intersections between the paths.

For this application, the optimization problem (10) needs to be extended to allow for multiple measurements. To this end, let $\Phi_{st}$ be a set of observations from measurement unit $s$ at time point $t$, for $t=1,\dots,T$, and for $s=1,\dots,S$. We obtain the maximum likelihood solution as the optimal solution to111This is a convex optimization problem which in principle can be solved with off-the-shelf solvers. In this example we use an efficient algorithm in the spirit of Proposition 2, but due to lack of space we defer the exact algorithm to a forthcoming paper.

[TABLE]

Consider a hidden Markov model where the states $X=\{X_{1},\ldots,X_{n}\}$ are the edges in the directed graph $\mathcal{G}=(V,X)$, and where the edge $X_{i}=(V^{\rm in}_{i},V^{\rm out}_{i})$ goes from $V^{\rm in}_{i}\in V$ to $V^{\rm out}_{i}\in V$. In this example we will use the graph illustrated in Figure 3, consisting of 11 nodes and $n=28$ edges. For the true model, transition probabilities are defined according to weights in the graph that represent which walking paths are preferred by the pedestrians. For the model used in the estimation we assume that this information is not known and use uniform weights. More specifically, the transition probabilities are given by

[TABLE]

where $\{w_{ij}\}$ is a set of weights. For the transition matrix $A$ used in the estimation we assign uniform weights $w_{ij}=1$, for all $(i,j)$ with $V^{\rm in}_{j}=V_{i}^{\rm out}$. For the true transitions the weights are defined as

[TABLE]

for $(i,j)$ such that $V^{\rm in}_{j}=V_{i}^{\rm out}$, and where $\mathcal{W}$ is the set of edges highlighted in Figure 3. Note that the second case implies that agents do not transition to the reverse edge in the next time step.

The agents are observed by a set of $S=7$ sensors located at the positions indicated in Figure 3. The observation probability for an agent on edge $X_{i}$ to be detected by a given sensor $s$ is defined as $b^{s}_{i1}=\min(0.99,2e^{-5d})$, where $d$ denotes the Euclidean distance between the location of $s$ and the midpoint of $X_{i}$. Consequently the probability of not being detected is $b^{s}_{i2}=1-b^{s}_{i1}$.

Given an initial distribution of 100 agents on the edge $(1,3)$, the flow and the measurements for the true ensemble are computed for $T=20$ time steps using the true transition and observation probabilities. Then we estimate the flow by solving the optimization problem (11). The true and estimated particle distributions are compared for some time steps in Figure 4, where the width of each edge is proportional to the number of agents on it. As can be seen in the figure the proposed method provides a good estimate also for this example.

VI Conclusions and future directions

In this work, we propose a method for estimating the flow of an ensemble of particles on a hidden Markov chain. The estimation, which is formulated as a maximum likelihood problem, can be recast as a convex optimization problem, for which we provide an efficient algorithm.

There are several natural directions in which this work can be extended. One restriction in this paper is that the number of agents are fixed and known, thus extensions with a variable number of particles, such as birth/death processes, could be of interest (cf. [14]). Furthermore, the model may be extended to continuous-state dynamics. Another natural direction is to study the connections to reciprocal processes.

-A Proof of Proposition 1

Let $\mu_{0}^{(N)}$ and $M^{(N)}$ be as described in the statement of the proposition. Moreover, let $\mathcal{Z}^{(N)}=\{(i,j)\mid m^{(N)}_{ij}\neq 0\}$, which is non-empty, and let $\mathcal{Z}^{(N)}_{i}=\{j\mid(i,j)\in\mathcal{Z}^{(N)}\}$, which is non-empty if and only if $(\mu_{0}^{(N)})_{i}>0$. Furthermore, let $\mathcal{Y}^{(N)}=\{i\mid(\mu_{0}^{(N)})_{i}>0\}$. Then, using Stirling’s formula

[TABLE]

see, e.g., [30], for $i\in\mathcal{Y}^{(N)}$ the multinomial coefficient in (4) can be bounded from above by

[TABLE]

where $|\mathcal{Z}^{(N)}_{i}|$ denotes the cardinality of the set, and where the second inequality follows from the fact that $\sum_{j\in\mathcal{Z}^{(N)}_{i}}m_{ij}=\sum_{j=1}^{n}m_{ij}=(\mu_{0})_{i}$, and that $e^{-(|\mathcal{Z}^{(N)}_{i}|-1)}\leq 1$. Thus, the log-likelihood of the probability for a transfer plan $M^{(N)}$ can be upper-bounded as follows:

[TABLE]

where the last inequality comes from the second-to-last expression by i) identifying the first term as the KL divergence, ii) noting that the second term is nonpositive, and iii) overestimating the third term by taking $(\mu_{0}^{(N)})_{i}=N$ for $i=1,\ldots,n$.

Similarly, underestimating the multinomial coefficients gives that the log-likelihood can be bounded from below by

[TABLE]

By using the two inequalities, the result follows.

-B Proof of Proposition 2

We solve (10) by a block coordinate ascent in the dual. First note that since $\mu_{0}$, $A$, $B$, and $\Phi_{t}$ are all elementwise nonnegative, so will the optimal solution $M_{[1:T]}^{*},D_{[1:T]}^{*},\mu_{[1:T]}^{*}$ also be. We can therefore add the constraint $\mu_{[1:T]}\geq 0$ to (10) without changing the optimal solution. For this problem, we relax the constraints in (10) with corresponding dual variables $\lambda_{M_{t}},\nu_{M_{t}},\lambda_{D_{t}},\nu_{D_{t}}$. Let $\mathfrak{M}:=M_{[1:T]}$, $\mathfrak{D}:=D_{[1:T]}$, and define the Lagrangian

[TABLE]

Minimizing this with respect to the matrices $M_{t}$ and $D_{t}$ gives explicit expressions for the optimal solution in terms of $\mu_{[0:T]}$ and the dual variables, i.e.,

[TABLE]

where $u_{t}=\exp(\lambda_{M_{t}})$, $w_{t}=\exp(\nu_{M_{t}})$, $x_{t}=\exp(\lambda_{D_{t}})$ and $v_{t}=\exp(\nu_{D_{t}})$, for $t=1,\dots,T$. Plugging these into the Lagrangian, we get the modified Lagrangian

[TABLE]

Noting that since $\mu_{[1:T]}$ occurs linearly in $L$, for the modified dual functional $\inf_{\mu_{[1:T]}\geq 0}L(\mu_{[1:T]},\lambda_{\mathfrak{M}},\nu_{\mathfrak{M}},\lambda_{\mathfrak{D}},\nu_{\mathfrak{D}})$ to be bounded from below, the corresponding factors need to be elementwise nonnegative. In this case the corresponding terms will be zero when taking the infimum, and thus the dual problem is to maximize

[TABLE]

subject to

[TABLE]

for $t=1,\dots,T-1$, and

[TABLE]

Note that neither the objective function (12), nor the first $T-1$ constraints (13) depend on $x_{T}$. Thus the optimal choice of $x_{T}$ is the one creating the most slack in the last constraint, i.e., the one maximizing the first two terms of (14). This is achieved by $x_{T}=e{\bf 1}./(Bv_{T})$, and for this choice of $x_{T}$ the constraint (14) can be replaced by

[TABLE]

Similarly, the slack in (13) is maximized by selecting

[TABLE]

and thus the constraints (13) can be replaced by

[TABLE]

for $t=1,\ldots,T-1$. Next, note that if some of the constraints (15) and (16) are not fulfilled with equality, the objective function (12) can be improved by increasing the values of the corresponding $v_{t}$. Therefore, in an optimal point we must have that

[TABLE]

This gives an expression for $w_{1}$, that depends on $v_{t}$, $t=1,\ldots,T$. Inserting this into the objective (12) leads to an unconstrained problem that depends only on $u_{1}$ and $v_{t}$, $t=1,\dots,T$, which we solve using block coordinate ascent. To this end, the unconstrained objective is first maximized with respect to $u_{1}$, which gives

[TABLE]

Further, note that the gradient of the unconstrained objective with respect to $v_{t}$ is

[TABLE]

where $y_{t}$ is defined by the recursion

[TABLE]

Hence, maximization with respect to $v_{t}$ is achieved by

[TABLE]

As the unconstrained problem is convex and the objective function continuously differentiable, the block coordinate ascent method converges [3, Prop. 2.7.1]. In the limit point, the hidden marginals can be reconstructed from

[TABLE]