The matrix product approximation for the dynamic cavity method

TL;DR

This paper introduces a matrix product edge message algorithm for the dynamic cavity method, enabling efficient and accurate analysis of stochastic dynamics on tree-like graphs with controllable precision.

Contribution

It presents a novel matrix product approximation technique for the dynamic cavity equations, improving computational efficiency and extending applicability to models with continuum time limits.

Findings

Computation costs grow linearly with time for Glauber-Ising dynamics.

The method outperforms Monte Carlo in error scaling, especially for low probability events.

Optimized truncation schemes enhance accuracy and efficiency.

Abstract

Stochastic dynamics of classical degrees of freedom, defined on vertices of locally tree-like graphs, can be studied in the framework of the dynamic cavity method which is exact for tree graphs. Such models correspond for example to spin-glass systems, Boolean networks, neural networks, and other technical, biological, and social networks. The central objects in the cavity method are edge messages -- conditional probabilities of two vertex variable trajectories. In this paper, we discuss a rather pedagogical derivation for the dynamic cavity method, give a detailed account of the novel matrix product edge message (MPEM) algorithm for the solution of the dynamic cavity equation as introduced in Phys. Rev. E 97, 010104(R) (2018), and present optimizations and extensions. Matrix product approximations of the edge messages are constructed recursively in an iteration over time. Computation costs and precision can be tuned by controlling the matrix dimensions of the MPEM in truncations. Without truncations, the dynamics is exact. Data for Glauber-Ising dynamics shows a linear growth of computation costs in time. In contrast to Monte Carlo simulations, the approach has a much better error scaling. Hence, it gives for example access to low probability events and decaying observables like temporal correlations. We discuss optimized truncation schemes and an extension that allows to capture models which have a continuum time limit.

Equations

$${\mathbf{\sigma }}^t:=({\sigma }_1^t,{\sigma }_2^t,\dots ).$$

$$P({\mathbf{\sigma }}^{t+1})=\sum _{{\mathbf{\sigma }}^t}W({\mathbf{\sigma }}^{t+1}|{\mathbf{\sigma }}^t)P({\mathbf{\sigma }}^t)\text{with}\sum _{\mathbf{\sigma }}W(\mathbf{\sigma }|{\mathbf{\sigma }}^{\prime })=1.$$

$$W({\mathbf{\sigma }}^{t+1}|{\mathbf{\sigma }}^t)=\prod _{i\in \mathcal{G}}{w}_i({\sigma }_i^{t+1}|{\mathbf{\sigma }}_{\partial i}^t)\text{with}\sum _{{\sigma }_i}{w}_i({\sigma }_i|{\mathbf{\sigma }}_{\partial i}^{\prime })=1,$$

$$P({\mathbf{\sigma }}^0{\mathbf{\sigma }}^1\dots {\mathbf{\sigma }}^{t-1}{\mathbf{\sigma }}^t)=\prod _{s=1}^{t}W({\mathbf{\sigma }}^s|{\mathbf{\sigma }}^{s-1})P({\mathbf{\sigma }}^0).$$

$$P({\mathbf{\sigma }}^0)=\prod _{i\in \mathcal{G}}{p}_i({\sigma }_i^0)\text{with}\sum _{{\sigma }_i}{p}_i({\sigma }_i)=1.$$

$${\bar{\sigma }}_i^{t+1}:=({\sigma }_i^0,{\sigma }_i^1,\dots ,{\sigma }_i^t,{\sigma }_i^{t+1}),$$

$$\mu_{i\to j}(\bar{\sigma}_{i}^{t+1}|\bar{\sigma}_{j}^{t}):=\sum_{\{\bar{\sigma}_{k}^{t+1}\}_{k\in\mathcal{G}_{i\to j}\setminus\{i\}}}\Big{[}\prod_{s=0}^{t}\prod_{q\in\mathcal{G}_{i\to j}}w_{q}(\sigma_{q}^{s+1}|{\boldsymbol{\sigma}}_{\partial q}^{s})\Big{]}\prod_{q\in\mathcal{G}_{i\to j}}p_{q}(\sigma_{q}^{0}).$$

$$\mu_{i\to j}(\bar{\sigma}_{i}^{t+1}|\bar{\sigma}_{j}^{t})=\sum_{\{\sigma_{k}^{s}\}_{(k,s)\in\bar{\mathcal{G}}^{t+1}_{i\to j}}}\Big{[}\prod_{s=0}^{t}\prod_{q\in\mathcal{G}^{t-s}_{i\to j}}w_{q}(\sigma_{q}^{s+1}|{\boldsymbol{\sigma}}_{\partial q}^{s})\Big{]}\prod_{q\in\mathcal{G}^{t+1}_{i\to j}}p_{q}(\sigma_{q}^{0}).$$

$$\mathcal{G}^{0}_{i\to j}:=\{i\},\quad\mathcal{G}^{s+1}_{i\to j}:=\mathcal{G}^{s}_{i\to j}\cup\Big{[}\bigcup_{k\in\mathcal{G}^{s}_{i\to j}}(\partial k\setminus\{j\})\Big{]}=\big{\{}k\in\mathcal{G}_{i\to j}\,|\,\operatorname{dist}(i,k)\leq s+1\big{\}}$$

$$\bar{\mathcal{G}}^{t}_{i\to j}:=\bigcup_{s=0}^{t}\Big{\{}(s,k)\,|\,k\in\mathcal{G}^{t-s}_{i\to j}\setminus\{i\}\Big{\}}.$$

$$P({\bar{\sigma }}_i^t,{\bar{\sigma }}_j^t)\equiv \sum _{\{{\bar{\sigma }}_k^t{\}}_{k\in \mathcal{G}\setminus \{i,j\}}}P({\mathbf{\sigma }}^0{\mathbf{\sigma }}^1\dots {\mathbf{\sigma }}^{t-1}{\mathbf{\sigma }}^t)={\mu }_{i\to j}({\bar{\sigma }}_i^t|{\bar{\sigma }}_j^{t-1}){\mu }_{j\to i}({\bar{\sigma }}_j^t|{\bar{\sigma }}_i^{t-1}).$$

$$\mu_{i\to j}(\bar{\sigma}_{i}^{t+1}|\bar{\sigma}_{j}^{t})=\sum_{\{\bar{\sigma}_{k}^{t}\}_{k\in\partial i\setminus\{i,j\}}}p_{i}(\sigma_{i}^{0})\Big{[}\prod_{s=0}^{t}w_{i}(\sigma_{i}^{s+1}|{\boldsymbol{\sigma}}_{\partial i}^{s})\Big{]}\Big{[}\prod_{k\in\partial i\setminus\{i,j\}}\mu_{k\to i}(\bar{\sigma}_{k}^{t}|\bar{\sigma}_{i}^{t-1})\Big{]}$$

$$|\psi ⟩=\sum _{n_1,\dots ,n_L}{\psi }_{n_1,\dots ,n_L}|\mathbf{n}⟩,\text{where}|\mathbf{n}⟩\equiv |n_1⟩\otimes |n_2⟩\otimes \dots \otimes |n_L⟩$$

$$|\psi ⟩\approx \sum _{\mathbf{n}}{A}_1^{n_1}{A}_2^{n_1}\cdots A_{L-1}^{n_{L-1}}{A}_L^{n_L}|\mathbf{n}⟩.$$

$$\mu_{i\to j}(\bar{\sigma}_{i}^{t}|\bar{\sigma}_{j}^{t-1})=A^{(0)}_{i\to j}(\sigma_{j}^{0})\Big{[}\prod_{s=1}^{t-1}A^{(s)}_{i\to j}(\sigma_{i}^{s-1}|\sigma_{j}^{s})\Big{]}A^{(t)}_{i\to j}(\sigma_{i}^{t-1})A^{(t+1)}_{i\to j}(\sigma_{i}^{t})$$

$$|\psi ⟩=\sum _{\mathbf{n}}{A}_1^{n_1}{A}_2^{n_1}\cdots A_{L-1}^{n_{L-1}}{A}_L^{n_L}|\mathbf{n}⟩,$$

$$|\psi ⟩=:\sum _{a,b}{\psi }_{a,b}|a⟩\otimes |b⟩.$$

$$|\psi ⟩=\sum _{a,b,k}{U}_{a,k}{\lambda }_k{V}_{b,k}^{\ast }|a⟩\otimes |b⟩=\sum _{k=1}^{M^{\prime }}{\lambda }_k|k_{\mathcal{L}}⟩\otimes |k_{\mathcal{R}}⟩,$$

$$|{\psi }_{\text{trunc}}⟩:=\sum _{k=1}^M{\lambda }_k|k_{\mathcal{L}}⟩\otimes |k_{\mathcal{R}}⟩\text{with error}{\Vert \psi -{\psi }_{\text{trunc}}\Vert }^2=\sum _{k=M+1}^{M^{\prime }}{\lambda }_k^2$$

$$|\psi ⟩=\sum _{\mathbf{n}}{Y}_1^{n_1}\cdots Y_{\ell }^{n_{\ell }}\cdot C\cdot Z_{\ell +1}^{n_{\ell +1}}\cdots Z_L^{n_L}|\mathbf{n}⟩\text{with}$$

$$\sum _n(Y_i^n{)}^{†}{Y}_i^n=\text{mathbbm}1\text{and}\sum _n{Z}_i^n(Z_i^n{)}^{†}=\text{mathbbm}1.$$

$$\displaystyle|a\rangle=\sum_{n_{1},\dotsc,n_{\ell}}\big{[}Y_{1}^{n_{1}}\dotsb Y_{\ell}^{n_{\ell}}\big{]}_{1,a}\,|n_{1}\dots n_{\ell}\rangle\quad\text{for}\ \mathcal{L}\ \text{and}$$

$$\displaystyle|b\rangle=\sum_{n_{\ell+1},\dotsc,n_{L}}\big{[}Z_{\ell+1}^{n_{\ell+1}}\dotsb Z_{L}^{n_{L}}\big{]}_{b,1}\,|n_{\ell+1}\dots n_{L}\rangle\quad\text{for}\ \mathcal{R}.$$

$$⟨a^{\prime }|a⟩$$

$$\displaystyle=\sum_{n_{2},\dotsc,n_{\ell}}\big{[}(Y_{\ell}^{n_{\ell}})^{\dagger}\dotsb(Y_{3}^{n_{3}})^{\dagger}(Y_{2}^{n_{2}})^{\dagger}Y_{2}^{n_{2}}Y_{3}^{n_{3}}\dotsb Y_{\ell}^{n_{\ell}}\big{]}_{a^{\prime},a}=\dots=[\mathbbm{1}]_{a^{\prime},a}=\delta_{a,a^{\prime}}$$

$$\mu_{i\to j}(\bar{\sigma}_{i}^{t+1}|\bar{\sigma}_{j}^{t})=B^{(0)}_{i\to j}(\sigma_{j}^{0})\Big{[}\prod_{s=1}^{t}B^{(s)}_{i\to j}(\sigma_{i}^{s-1}|\sigma_{j}^{s})\Big{]}B^{(t+1)}_{i\to j}(\sigma_{i}^{t})B^{(t+2)}_{i\to j}(\sigma_{i}^{t+1})$$

$$\mu_{i\to j}(\bar{\sigma}_{i}^{t+1}|\bar{\sigma}_{j}^{t})=C^{(0)}_{i\to j}(\sigma_{i}^{0})\Big{[}\prod_{s=1}^{t+1}C^{(s)}_{i\to j}(\sigma_{i}^{s}|\sigma_{j}^{s-1})\Big{]}.$$

$$C^{(s)}_{i\to j}(\sigma_{i}^{s}|\sigma_{j}^{s-1})=\sum_{\sigma_{k_{1}}^{s-1}\dotsc\sigma_{k_{z-1}}^{s-1}}w_{i}(\sigma_{i}^{s}|{\boldsymbol{\sigma}}_{\partial i}^{s-1})\Big{[}\bigotimes_{n=1}^{z-1}A^{(s)}_{k_{n}\to i}(\sigma_{k_{n}}^{s-1}|\sigma_{i}^{s})\Big{]}$$

$$\big{[}C^{(s)}_{i\to j}(\sigma_{i}^{s}|\sigma_{j}^{s-1})\big{]}_{(a_{1},a_{2}),(b_{1},b_{2})}=\sum_{\sigma_{k_{1}}^{s-1},\sigma_{k_{2}}^{s-1}}w_{i}(\sigma_{i}^{s}|\sigma_{j}^{s-1}\sigma_{k_{1}}^{s-1}\sigma_{k_{2}}^{s-1})\\[-3.00003pt] \times\big{[}A^{(s)}_{k_{1}\to i}(\sigma_{k_{1}}^{s-1}|\sigma_{i}^{s})\big{]}_{a_{1},b_{1}}\cdot\big{[}A^{(s)}_{k_{2}\to i}(\sigma_{k_{2}}^{s-1}|\sigma_{i}^{s})\big{]}_{a_{2},b_{2}}.$$

$$C_{i\to j}^{(s)}({\sigma }_i^s|{\sigma }_j^{s-1})U^{(s+1)}{\Lambda }^{(s+1)}\stackrel{\text{SVD}}{=:}{U}^{(s)}{\Lambda }^{(s)}{\tilde{C}}_{i\to j}^{(s)}({\sigma }_i^s|{\sigma }_j^{s-1}),$$

$$\mu_{i\to j}(\bar{\sigma}_{i}^{t+1}|\bar{\sigma}_{j}^{t})=\tilde{C}^{(0)}_{i\to j}(\sigma_{i}^{0})\Big{[}\prod_{s=1}^{t+1}\tilde{C}^{(s)}_{i\to j}(\sigma_{i}^{s}|\sigma_{j}^{s-1})\Big{]}\stackrel{{\scriptstyle\text{trunc}}}{{\approx}}\bar{C}^{(0)}_{i\to j}(\sigma_{i}^{0})\Big{[}\prod_{s=1}^{t+1}\bar{C}^{(s)}_{i\to j}(\sigma_{i}^{s}|\sigma_{j}^{s-1})\Big{]}.$$

$$\mu_{i\to j}(\bar{\sigma}_{i}^{t+1}|\bar{\sigma}_{j}^{t})\stackrel{{\scriptstyle\text{trunc}}}{{\approx}}\Big{[}\prod_{s=0}^{t}D^{(s)}_{i\to j}(\sigma_{i}^{s}|\sigma_{j}^{s})\Big{]}D^{(t+1)}_{i\to j}(\sigma_{i}^{t+1}).$$

Full text

1. I Introduction

2. II Stochastic dynamics on locally tree-like graphs

3. II.1 Global equation of motion

4. II.2 Probability for global trajectories

5. II.3 Motivation and definition of edge messages

6. II.4 Simplification of edge messages

7. II.5 Computing edge trajectory probabilities from edge messages

8. II.6 Dynamic cavity equation

9. III Matrix product edge messages (MPEM)

10. III.1 General idea and motivation

11. III.2 Canonical form

12. III.3 Controlled truncations

13. IV Algorithm for the time evolution of matrix product edge messages

14. IV.1 One exact MPEM evolution step

15. IV.2 Truncating MPEMs and recovering the canonical form

16. V Evaluation of observables

17. VI An improved truncation scheme

18. VII Computation costs in Glauber-Ising dynamics

19. VIII Models with vertex-state dependence

20. IX Discussion

I.   Introduction

Statistical physics provides a comprehensive framework for the study of many-body systems in equilibrium with their environment. It is a foundation of modern physics, tracing back the laws of thermodynamics to characteristics of the microscopic degrees of freedom. Beyond physics, its formalism is successfully applied in many different fields such as information theory, computer science, economy, and sociology. In contrast, our methodological toolbox for dynamics in probabilistic systems is rather limited, which restricts our understanding of stochastic dynamical phenomena.

Stochastic dynamics, i.e., dynamics that are governed by probabilistic instead of deterministic rules, are ubiquitous in nature as well as in social and technological systems VanKampen2007 ; Freund2010 ; Capasso2004 ; AitSahalia2014 . Here, we focus in particular on stochastic dynamics in networks Barrat2008 ; Newman2010 . Some challenging examples are the thermal dynamics in spin glasses, avalanche dynamics, the dynamics of infectious diseases, the evolution of opinions in social networks, the stability of functionality under perturbations in technological networks, or synchronization phenomena. The probabilistic features can be intrinsic or due to our ignorance of certain details that are not essential for the observed (macroscopic) phenomena. A simple method for the simulation of a stochastic system with states $\sigma_{i}$ on $N$ vertices of a network is to propagate the probability distribution $P(\sigma_{1}^{t},\dotsc,\sigma_{N}^{t})$ in time. Unfortunately, the corresponding computation costs are exponential in the system size $N$ and it is also difficult to asses temporal correlations in this way.

For the investigation of network systems, a strong simplification can be achieved when cycles in the interaction graph are rare or sufficiently long. This is the case for locally tree-like graphs [see Fig. 1(a)] such as random regular graphs, Erdős-Rényi graphs, and Gilbert graphs. For such random graphs with $N$ vertices, almost all cycles have length $\gtrsim\log N$ such that their effect is negligible in the thermodynamic limit Mezard2009 . For static problems, this feature is used in the cavity method Mezard1986-1 ; Mezard2001-20 , where conditional nearest-neighbor probabilities are computed iteratively within the Bethe-Peierls approximation. This efficient method has been applied very successfully to study, for example, equilibrium properties of spin glasses Mezard2001-20 , computationally hard satisfiability problems Mezard2002-297 ; Mezard2002-66 , and random matrix ensembles Rogers2008-78 .

Subsequently, the cavity method has been generalized to dynamical problems, resulting in the dynamic cavity method or dynamic belief propagation Neri2009-08 ; Karrer2010-82 . The central objects in this approach are not the time-dependent probabilities for global system states, but conditional probabilities $\mu_{i\to j}(\sigma^{0}_{i},\dotsc,\sigma^{t}_{i}|\sigma^{0}_{j},\dotsc,\sigma^{t-1}_{j})$ for state trajectories on neighboring sites $i$ and $j$ (edges $i\to j$). These so-called edge messages $\mu_{i\to j}$ are generated in an iteration over time. Consequently, the computational complexity is now linear instead of exponential in the system size $N$. Unfortunately, the number of possible trajectories and, hence, the computational complexity increase exponentially in time. Applications have thus been restricted to either very short times Neri2009-08 ; Kanoria2011-21 , oriented graphs Neri2009-08 , a no-backtracking approximation Shrestha2015-92 ; Castellano2018-98 , or unidirectional dynamics with local absorbing states Karrer2010-82 ; Lokhov2015-91 ; Lokhov2014-90 ; Altarelli2014-112 ; Altarelli2014-10 ; Shrestha2014-89 ; Lokhov2017-114 . In the latter case, one can exploit that vertex trajectories can be parametrized by a few switching times.

For general stochastic network dynamics, it is a challenging endeavor to find good approximative solutions to the dynamic cavity equations with polynomial computation costs. A drastic simplification is to neglect temporal correlations completely as in the one-time approximation (a.k.a. time factorization) Neri2009-08 ; Aurell2011-04 ; Aurell2012-85 ; Zhang2012-148 or to retain only short-time correlations as in the one-step Markov ansatz DelFerraro2015-92 . While this can be expected to work well for stationary states at high temperatures, such approximations are usually severe for short to intermediate times or low temperatures. Other approximative approaches are the cluster variational method Pelizzola2013-86 ; Vazquez2017-3 ; Pelizzola2017-7 (applicable for short-range spatio-temporal correlations) or perturbative schemes Roudi2011 ; Aurell2012-85 ; BachschmidRomano2016-49 , generating functional analysis Hatchett2004-37 ; Mimura2009-42 ; Mozeika2011-106 ; Mozeika2012-92 ; Coolen2012-92 , and the generalized mean field approximation Mahmoudi2014-7 .

In Ref. Barthel2018-97 , we have introduced a novel efficient algorithm for precise solutions of the parallel dynamic cavity equations for locally tree-like graphs and general bidirectional dynamics. The approach is based on matrix-product approximations for the edge messages $\mu_{i\to j}$. The computational complexity is decreased from exponential to polynomial in the duration of the dynamical process. As we demonstrated for non-equilibrium Glauber dynamics in the kinetic Ising model, one can obtain quasi-exact results with relatively small matrix dimensions. Computation costs and accuracy can be tuned through controlled truncations of the matrix dimensions. Matrix-product approximations are widely used in condensed matter theory to encode spatial correlations of quantum states Accardi1981 ; Fannes1992-144 ; Rommer1997 ; Schollwoeck2011-326 . Here, we employ matrix products to encode temporal correlations in the edge messages.

This paper gives a pedagogical introduction to the dynamic cavity method in Sec. II. Sections III-V give a detailed account of the novel matrix product edge message (MPEM) algorithm Barthel2018-97 for the solution of the dynamic cavity equation. Section VI provides an improved MPEM truncation scheme that substantially reduces computation costs. For Glauber-Ising dynamics, we demonstrate a linear growth of computation costs in Sec. VII. In Sec. VIII, the MPEM method is generalized to systems where the transition probability for the state of a vertex does not only depend on the state of its neighbors in the interaction graph, but also on the current state of the vertex itself.

II.   Stochastic dynamics on locally tree-like graphs

II II.1.   Global equation of motion

Let us consider a graph $\mathcal{G}$ and denote the state at vertex $i\in\mathcal{G}$ by $\sigma_{i}$. The state of the full system at time $t$ will be denoted by

[TABLE]

Given the state probabilities $P({\boldsymbol{\sigma}}^{t})$ for time $t$, the probabilities $P({\boldsymbol{\sigma}}^{t+1})$ for the subsequent time step $t+1$ are obtained by applying the transition matrix $W({\boldsymbol{\sigma}}^{t+1}|{\boldsymbol{\sigma}}^{t})$, where the notation of the arguments indicates that $W$ is a conditional probability.

[TABLE]

In correspondence with the structure of the graph $\mathcal{G}$, we assume that, in every time step, the probability for $\sigma_{i}^{t+1}$ only depends on the states $\sigma_{k}^{t}$ of nearest-neighbors $k$ at the previous time step such that the global transition matrix $W$ is a product of local transition matrices $w_{i}$,

[TABLE]

where the vicinity $\partial i:=\{k\in\mathcal{G}\,|\,\operatorname{dist}(i,k)\leq 1\}$ of vertex $i$, contains $i$ and its nearest neighbors. Equations (1) and (2) specify the parallel stochastic dynamics of the system.

II II.2.   Probability for global trajectories

With the definitions above, the probability for a trajectory ${\boldsymbol{\sigma}}^{0}\to{\boldsymbol{\sigma}}^{1}\to\dots\to{\boldsymbol{\sigma}}^{t}$ of states up to time $t$ takes the form

[TABLE]

For simplicity, let us assume initially uncorrelated states, i.e.,

[TABLE]

II II.3.   Motivation and definition of edge messages

An exact solution of the global equation of motion (1) is usually not possible. In principle, one can employ a simple iteration over time. However, the corresponding memory costs (to keep track of $P({\boldsymbol{\sigma}}^{t})$) and computation costs scale exponentially in the system size $|\mathcal{G}|$ such that this approach is limited to very small systems.

A very useful technique is the Markov chain Monte Carlo method. In this approach, one repeatedly generates states ${\boldsymbol{\sigma}}^{0}$ according to the $t=0$ probability $P({\boldsymbol{\sigma}}^{0})$, and propagates them in time by choosing a state ${\boldsymbol{\sigma}}^{t+1}$ with probability $W({\boldsymbol{\sigma}}^{t+1}|{\boldsymbol{\sigma}}^{t})$ for $t=0,1,2,\dotsc$. While the technique can be used to study many interesting problems, it has several disadvantages. First of all, one cannot address the thermodynamic limit directly and needs to perform a costly finite-size scaling analysis. Secondly, the number of required samples can be very large if temporal correlations are non-trivial. Also, the error scaling is not very favorable. When increasing the number of samples $N_{s}$, errors decrease slowly as $1/\sqrt{N_{s}}$. Especially if the absolute value of the desired observable is small, as is usually the case when studying temporal correlations, it is often impossible to achieve a sufficient precision.

To resolve these issues in our approach for locally tree-like graphs, we shift the focus from global-state probabilities $P({\boldsymbol{\sigma}}^{t})$ to so-called edge messages. These are conditional probabilities for the trajectory $\bar{\sigma}^{t+1}_{i}$ of the variable at vertex $i$,

[TABLE]

and the trajectory $\bar{\sigma}^{t}_{j}$ on a neighboring site $j$. If the graph is a tree, the exact evolution of the system can be formulated in terms of the edge messages and leads to the dynamic cavity equations Neri2009-08 ; Karrer2010-82 discussed in the following. If the graph contains some longer loops, i.e., is only locally tree-like, the dynamic cavity equations give an approximation to the exact dynamics with the accuracy depending on the number and lengths of loops.

When we remove an edge $(i,j)$ from a tree graph $\mathcal{G}$, it is decomposed into two parts $\mathcal{G}_{i\to j}$ and $\mathcal{G}_{j\to i}$ such that $\mathcal{G}_{i\to j}$ contains $i$ and the subgraph that is still connected to $i$ (and all subgraphs of $\mathcal{G}$ that are disconnected from $i$ and $j$) while $\mathcal{G}_{j\to i}$ contains $j$ and the subgraph that is still connected to $j$. See Fig. 1(a). Now, an edge message $\mu_{i\to j}(\bar{\sigma}_{i}^{t+1}|\bar{\sigma}_{j}^{t})$ for edge $(i,j)$ is the probability of the trajectory $\bar{\sigma}_{i}^{t+1}$ on vertex $i$ for the dynamics being restricted to the subgraph $\mathcal{G}_{i\to j}$ under the condition that we impose the trajectory $\bar{\sigma}_{j}^{t}$ on vertex $j$. Specifically, it is the product of all transition matrices $w_{k}(\sigma_{k}^{s+1}|{\boldsymbol{\sigma}}_{\partial k}^{s})$ and $p_{k}(\sigma_{k}^{0})$ for all $k\in\mathcal{G}_{i\to j}$ and $0\leq s\leq t$, summed over all $\bar{\sigma}_{k}^{t+1}$ with $k\in\mathcal{G}_{i\to j}\setminus\{i\}$, i.e.,

[TABLE]

A pictorial representation for a one-dimensional graph is given in Fig. 2.

II II.4.   Simplification of edge messages

Due to the normalization constraint $\sum_{\sigma_{i}}w_{i}(\sigma_{i}|{\boldsymbol{\sigma}}^{\prime}_{\partial i})=1$ of the transition matrices, the expression (6) for the edge message can be simplified substantially by removing all $w_{k}$ and $p_{k}$ outside the causal cone of edge $i\to j$ such that

[TABLE]

In this equation, we use the (recursively defined) subgraphs

[TABLE]

and the causal cone $\bar{\mathcal{G}}^{t+1}_{i\to j}$ of edge $i\to j$ which is a subset of $\mathbb{N}_{0}\times\mathcal{G}_{i\to j}$,

[TABLE]

See Figs. 3 and 4 for examples.

II II.5.   Computing edge trajectory probabilities from edge messages

Before showing that edge messages can be generated in an iteration (Sec. II.6), note that they allow us to evaluate easily many observables of interest. In particular, the product of edge messages $\mu_{i\to j}(\bar{\sigma}_{i}^{t}|\bar{\sigma}_{j}^{t-1})$ and $\mu_{j\to i}(\bar{\sigma}_{j}^{t}|\bar{\sigma}_{i}^{t-1})$ yields the joint probability for the occurrence of trajectories $\bar{\sigma}_{i}^{t}$ and $\bar{\sigma}_{j}^{t}$ on vertices $i$ and $j$,

[TABLE]

See Fig. 5.

By marginalizing Eq. (10) over certain subsets of the variables $\{\sigma_{i}^{s}\}_{0\leq s\leq t}$ and $\{\sigma_{j}^{s}\}_{0\leq s\leq t}$, one obtains the probabilities required for the evaluation of time-local observables or temporal correlation functions.

II II.6.   Dynamic cavity equation

Inspection of the simplified expression (7) for edge message $\mu_{i\to j}(\bar{\sigma}_{i}^{t+1}|\bar{\sigma}_{j}^{t})$ shows that it can be constructed from the time-$t$ edge messages $\mu_{k\to i}(\bar{\sigma}_{k}^{t}|\bar{\sigma}_{i}^{t-1})$ according to

[TABLE]

with $\mu_{i\to j}(\sigma_{i}^{0})\equiv p_{i}(\sigma_{i}^{0})$. This recursion relation for the edge messages, the dynamic cavity equation Neri2009-08 ; Karrer2010-82 , is very useful and central for the matrix product edge message time-evolution algorithm presented in Sec. IV. An illustration of Eq. (11) for the example of a Y-junction graph is given in Fig. 4.

III.   Matrix product edge messages (MPEM)

This section and Sec. IV give a detailed account of an evolution algorithm that employs precise approximations to edge messages (7) in matrix product form Barthel2018-97 . Extensions will be described in Secs. VI and VIII.

III III.1.   General idea and motivation

In the following let $d$ denote the number of different vertex states ($\sigma_{i}=1,\dotsc,d$). In Sec. II, we have seen that the stochastic dynamics on a tree graph $\mathcal{G}$ can not only be simulated by sampling trajectories $({\boldsymbol{\sigma}}^{0}{\boldsymbol{\sigma}}^{1}\dots\boldsymbol{\sigma}^{t-1}{\boldsymbol{\sigma}}^{t})$ of the full system according to their weights (1) using a Monte Carlo algorithm, but also by an iterative construction of edge messages $\mu_{i\to j}(\bar{\sigma}_{i}^{t}|\bar{\sigma}_{j}^{t-1})$ which give access to many observables of interest such as the evolution of time-local observables or temporal correlators. Importantly, the computation costs and memory requirements for evolving and storing an edge message ($\mathcal{O}(d^{2t-1})$) are independent of the system size. However, without approximations, the costs increase exponentially with time $t$. So, an exact treatment is limited to short times.

To resolve this problem, we can employ an idea from the density matrix renormalization group (DMRG) White1992-11 ; White1993-10 ; Schollwoeck2005 which is a numerical algorithm for the simulation of strongly-correlated quantum many-body systems (predominantly for one-dimensional lattices) which has been used very successfully for a wide range of applications in condensed matter physics, the physics of ultracold atomic gases, and quantum chemistry. In principle, the encoding of a many-body state for a lattice of $L$ sites,

[TABLE]

label orthonormal basis states ($\langle{\boldsymbol{n}}|{\boldsymbol{n}}^{\prime}\rangle=\delta_{{\boldsymbol{n}}{\boldsymbol{n}}^{\prime}}$), requires storage of the $d^{L}$ expansion coefficients $\psi_{n_{1},\dotsc,n_{L}}$. However, in ground states of typical quantum many-body systems, spatial correlations like $\langle\hat{S}^{+}_{x}\hat{S}^{-}_{x^{\prime}}\rangle$ decay quickly in the distance $|x-x^{\prime}|$, exponentially or algebraically. Due to this fact and corresponding entanglement properties, one can approximate $\psi_{n_{1},\dotsc,n_{L}}$ by much fewer effective degrees of freedom, given by the elements of $M_{x}\times M_{x+1}$ matrices $A_{x}^{n_{x}}$ in the approximation Hastings2007-76 ; Brandao2013-9 ; Verstraete2005-5 ; Barthel2017_08unused

[TABLE]

The $M_{x}$ are called bond dimensions. In order for the matrix product $A_{1}^{n_{1}}\dotsb A_{L}^{n_{L}}$ to yield a scalar, we require dimensions $M_{1}=M_{L+1}=1$ at the system boundaries. The elements of matrices $A_{x}^{n_{x}}$ are effective degrees of freedom that encode correlations around site $x$. The larger the bond dimensions $M_{x}$ are chosen, the more effective degrees of freedom are taken into account, the higher the computation costs, and the higher the precision of the approximation. In any case, the memory costs ($\mathcal{O}(dLM^{2})$) and computation costs ($\mathcal{O}(d^{k}LM^{3})$) are now linear in the system size $L$. The right-hand side of Eq. (13) is called a matrix product state (MPS) Accardi1981 ; Fannes1992-144 ; Rommer1997 ; Schollwoeck2011-326 . More recently, MPS are also discussed in the more mathematical literature, sometimes under the name tensor train decomposition Oseledets2011-33 . DMRG is a class of algorithms operating on MPS – most importantly, to compute ground state approximations, thermal equilibrium states, or to study non-equilibrium phenomena. Observables can often be computed to machine precision.

The idea is now to similarly exploit that (connected) temporal correlations in edge messages often decay quickly in time $t$ and/or in the time difference $|t-t^{\prime}|$ and to approximate edge messages $\mu_{i\to j}(\bar{\sigma}_{i}^{t}|\bar{\sigma}_{j}^{t-1})$ in the form of a matrix product.

III III.2.   Canonical form

Let us define the canonical form of a matrix product edge message (MPEM) as

[TABLE]

The graphical representation is shown in Fig. 6. The particular choice of assigning state labels $\{\sigma^{s}_{i}\}$ and $\{\sigma^{s}_{j}\}$ to the $M_{s}\times M_{s+1}$ matrices occurring in the matrix product (14) is advantageous for the implementation of the dynamic cavity equation (11) for MPEMs as will become clear in Sec. IV. In order for the matrix product (14) to yield a scalar, we set $M_{0}=M_{t+2}=1$.

III III.3.   Controlled truncations

When advancing an MPEM $\mu_{i\to j}$ one step in time according to the dynamic cavity equation (11), we need to contract 111With the contraction of two tensors, we refer to their product, summed over certain sets of common indices. several MPEMs, $\mu_{k\to i}(\bar{\sigma}_{i}^{t}|\bar{\sigma}_{j}^{t-1})$ with $k$ from the neighborhood $\partial i$ of vertex $i$, and local transition matrices $w_{i}$ and, subsequently, express the result again in MPEM form. In general, this generates an MPEM $\mu_{i\to j}(\bar{\sigma}_{i}^{t+1}|\bar{\sigma}_{j}^{t})$ with increased bond dimension $M^{\prime}\geq M$. In order to control the growth of the bond dimensions, and hence the computation costs, we would like to reduce the bond dimension of the new MPEM. This can indeed be done in a controlled fashion such that the resulting truncated MPEM is close to the original. Here, we discuss the most simple truncation scheme Barthel2018-97 and an optimized scheme is presented in Sec. VI. The iterative time evolution of MPEMs is described in Sec. IV.

Because the notation is a little simpler, let us demonstrate truncations using the example of an MPS for a quantum system of $L$ sites,

[TABLE]

Let us split the system into two parts $\mathcal{L}$, containing sites $[1,\ell]$, and $\mathcal{R}$, containing sites $[\ell+1,L]$. Let $\{|a\rangle\}$ and $\{|b\rangle\}$ be orthonormal bases for left and right parts, respectively, such that

[TABLE]

The objective is to find a good approximation $|\psi_{\text{trunc}}\rangle$ of $|\psi\rangle$ in a reduced vector space, where the 2-norm distance is used to quantify the accuracy. This can be done using a singular value decomposition (SVD) of $\psi_{a,b}$. In matrix form, it reads $\psi=U\Lambda V^{\dagger}$, where $U$ and $V$ are unitary matrices and $\Lambda$ is the diagonal matrix containing descendingly ordered singular values $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{M^{\prime}}$ such that

[TABLE]

The fact that this procedure yields the best rank-$M$ approximation of $|\psi\rangle$ for the given bipartition corresponds to the Eckart-Young theorem Golub1996 ; Stewart1993-35 .

So, to truncate in this fashion bond dimensions of an MPS $|\psi\rangle$ (and analogously for MPEM), we need to express it in suitable orthonormal bases for $\mathcal{L}$ and $\mathcal{R}$. This can be achieved by exploiting the freedom to replace two subsequent matrices $(A^{n_{x}}_{x},A^{n_{x+1}}_{x+1})$ in the matrix product (15) by $(A^{n_{x}}_{x}X^{-1},XA^{n_{x+1}}_{x+1})$. The inserted non-singular matrices $X$ and $X^{-1}$ clearly leave $|\psi\rangle$ invariant as they cancel in the matrix product. Using this gauge freedom, we can bring the matrix product (15) into the form

[TABLE]

The left and right orthonormality constraints (18b) are imposed by sequences of singular value decompositions, sweeping from site 1 to site $\ell$ and from site $L$ to site $\ell+1$ with details discussed in Sec. IV.2. The resulting matrix product (18) is in fact of the form (16), with $\psi_{a,b}\leftrightarrow C_{a,b}$, and orthonormal basis states

[TABLE]

The orthonormality of these states follows from Eq. (18b):

[TABLE]

and similarly for the states $\{|b\rangle\}$. With a singular value decomposition of the matrix $C$, we then obtain an optimally truncated state $|\psi_{\text{trunc}}\rangle$ [Eq. (17b)] in MPS form.

IV.   Algorithm for the time evolution of matrix product edge messages

Given MPEMs $\mu_{i\to j}(\bar{\sigma}_{i}^{t}|\bar{\sigma}_{j}^{t-1})$ [Eq. (14)] for all edges at time $t$, we want to do one time step according to the dynamic cavity equation (11) and obtain MPEM approximations

[TABLE]

for the edge messages at time $t+1$.

Let us assume for now that vertex $i$ is not member of its neighborhood $\partial i$, i.e., that the local transition matrix $w_{i}(\sigma_{i}^{s+1}|{\boldsymbol{\sigma}}_{\partial i}^{s})$ is independent of $\sigma_{i}^{s}$. The more general case can also be handled as well but requires a slightly more complicated algorithm. The corresponding extension is described in Sec. VIII.

IV IV.1.   One exact MPEM evolution step

First, we generate a non-canonical matrix product representation of the evolved edge message (11), in particular, choosing the non-canonical form

[TABLE]

To this purpose, the matrices $C^{(s)}_{i\to j}(\sigma_{i}^{s}|\sigma_{j}^{s-1})$ for the bulk, $0<s<t$, are obtained by contracting the transition matrix $w_{i}(\sigma_{i}^{s}|{\boldsymbol{\sigma}}_{\partial i}^{s-1})$ for vertex $i$ with the tensors $A^{(s)}_{k\to i}(\sigma_{k}^{s-1}|\sigma_{i}^{s})$ from the time-$t$ MPEMs for edges $k\to i$ with vertices $k\in\partial i\setminus\{j\}=:\{k_{1},\dotsc,k_{z-1}\}$, where $z$ is the degree of vertex $i$. In this contraction, we sum over the $z-1$ common indices $\sigma_{k}^{s-1}$ as illustrated in Fig. 7. The resulting matrices

[TABLE]

have left and right multi-indices of dimensions $\bar{M}_{s}=(M_{s})^{z-1}$ and $\bar{M}_{s+1}=(M_{s+1})^{z-1}$, corresponding to the direct products ${\boldsymbol{a}}:=(a_{k_{1}},\dotsc,a_{k_{z-1}})$ and ${\boldsymbol{b}}:=(b_{k_{1}},\dotsc,b_{k_{z-1}})$ of the left and right indices of the matrices $A^{(s)}_{k\to i}(\sigma_{k}^{s-1}|\sigma_{i}^{s})$. For notational simplicity, we have assumed here that the bond dimensions $M_{s}$ are the same for all time-$t$ MPEMs of the edges $k\to i$ with $k\in\{k_{1},\dotsc,k_{z-1}\}$.

Specifically, for the case $z=3$ with $\partial i=\{j,k_{1},k_{2}\}$, the contractions take for example the form

[TABLE]

At the extremal time slices $s=0$ and $s=t,t+1$, the contractions differ slightly. They are specified graphically in Fig. 8.

At this point, it has become clear why the specific choice for the assignment of state indices $\{\sigma_{i}^{s}\,|\,0\leq s\leq t\}$ and $\{\sigma_{j}^{s}\,|\,0\leq s\leq t-1\}$ to tensors $A_{i\to j}^{(s)}$ in the canonical MPEM form (14) is favorable. It allows for an entirely local construction of tensor $C_{i\to j}^{(s)}$ from tensors $A_{k\to i}^{(s)}$, i.e., only $A$-tensors of a single time-slice $s$ are involved.

IV IV.2.   Truncating MPEMs and recovering the canonical form

The exact progression by one time step yields the non-canonical MPEM (22) with bond dimensions being increased from $M_{s}$ to $\tilde{M}_{s}=(M_{s})^{z-1}$. If we would proceed without any approximation, the computation cost would increase exponentially in time and, hence, restrict the simulation to very short times. Therefore, we are faced with two objectives: (a) reducing bond dimensions by a controlled truncation of the MPEM as described in Sec. III.3 and motivated physically in Sec. III.1, and (b) rearranging the assignment of physical state indices $\{\sigma_{i}^{s}\}$ and $\{\sigma_{j}^{s}\}$ to attain the canonical form (21) of the evolved MPEM. This can be achieved in different ways: A relatively simple method Barthel2018-97 is described in the following and an optimized scheme is introduced in Sec. VI.

If we want to truncate bond $s$, i.e., reduce the bond dimension $\tilde{M}_{s}$ to something smaller, we need to take care of orthonormality as discussed in Sec. III.3. In particular, we need to express the edge message in orthonormal reduced bases for left block $[0,s]$ and right block $[s+1,t+1]$. This can be done by imposing left and right orthonormality constraints (18b) for the $C$-tensors. If we would not do so and truncate bond dimension $\tilde{M}_{s}$, say, through an SVD of tensor $C^{(s)}$ without having the other tensors in orthonormal form, the approximation error would be uncontrolled. Truncating the smallest singular values of $C^{(s)}$ would then not correspond to the best rank-$M_{s}$ approximation.

In a first sweep from right ($s=t+1$) to left ($s=0$), we can sequentially impose the right orthonormality constraint $\sum_{n}Z_{i}^{n}(Z_{i}^{n})^{\dagger}=\mathbbm{1}$ on all $C$-tensors. As the original $C$-tensors do not obey the orthonormality constraints, we can only truncate singular values $\lambda_{k}$ that are (up to machine precision) zero. At the right boundary, we start with the singular value decomposition $C^{(t+1)}_{i\to j}(\sigma_{i}^{t+1}|\sigma_{j}^{t})\stackrel{{\scriptstyle\text{SVD}}}{{=:}}U^{(t+1)}\Lambda^{(t+1)}\tilde{C}^{(t+1)}_{i\to j}(\sigma_{i}^{t+1}|\sigma_{j}^{t})$, continue for the bulk tensors $s=t,t-1,\dotsc,1$ with

[TABLE]

as shown in Fig. 9 and end at the left boundary with $\tilde{C}^{(0)}_{i\to j}(\sigma_{i}^{0}):=C^{(0)}_{i\to j}(\sigma_{i}^{0})U^{(1)}\Lambda^{(1)}$ such that all $\tilde{C}$-tensors except $\tilde{C}^{0}$ now obey the right orthonormality constraint. The computation cost for each such SVD is proportional to $\tilde{M}^{3}$.

In a subsequent sweep from left to right, again based on singular value decompositions, we can now do the actual truncations to decrease bond dimensions ($\tilde{C}\to\bar{C}$), e.g., by setting a threshold $\varepsilon$ and truncating all singular values $\lambda_{k}$ with $\lambda_{k}^{2}<\varepsilon\cdot\sum_{q}\lambda_{q}^{2}$. All $\bar{C}$-tensors except $\bar{C}^{(t+1)}$ now obey the left orthonormality constraint and

[TABLE]

We now need to reorder the indices $\{\sigma_{i}^{s}\}$ and $\{\sigma_{j}^{s}\}$ of the vertex states. In a sweep from right to left, we go from the index assignment $(\sigma_{i}^{0})(\sigma_{i}^{1}|\sigma_{j}^{0})(\sigma_{i}^{2}|\sigma_{j}^{1})\dots(\sigma_{i}^{t+1}|\sigma_{j}^{t})$ in the truncated and orthonormalized version of the matrix product in Eq. (25) to the assignment $(\sigma_{i}^{0}\sigma_{j}^{0})(\sigma_{i}^{1}|\sigma_{j}^{1})\dots(\sigma_{i}^{t}|\sigma_{j}^{t})(\sigma_{i}^{t+1})$ in the matrix product

[TABLE]

At the right boundary, we start with a singular value decomposition and controlled truncation (SVT) $\bar{C}^{(t+1)}_{i\to j}(\sigma_{i}^{t+1}|\sigma_{j}^{t})\stackrel{{\scriptstyle\text{SVT}}}{{\approx:}}U^{(t+1)}(\sigma_{j}^{t})\Lambda^{(t+1)}D^{(t+1)}_{i\to j}(\sigma_{i}^{t+1})$ as described in footnote 222The resulting matrices $\mathcal{U}^{\dagger}$ and $\mathcal{D}$, defined as $\mathcal{U}_{(a,\sigma),k}:=[U^{(t+1)}(\sigma)]_{a,k}$ and $\mathcal{D}_{k,(\sigma,b)}:=[D^{(t+1)}_{i\to j}(\sigma)]_{k,b}$ are isometric in the sense that $\mathcal{U}^{\dagger}\mathcal{U}=\mathbbm{1}=\mathcal{D}\mathcal{D}^{\dagger}$. These matrices are obtained from the singular value decomposition $\bar{\mathcal{C}}\stackrel{{\scriptstyle\text{SVT}}}{{\approx}}\mathcal{U}\Lambda\mathcal{D}$ and subsequent truncation of small singular values $\lambda_{k}$, where $\bar{\mathcal{C}}_{(a,\sigma),(b,\sigma^{\prime})}:=[\bar{C}^{(t+1)}_{i\to j}(\sigma|\sigma^{\prime})]_{a,b}$. The isometric property of $\mathcal{D}$, corresponds to the right orthonormality constraint as defined in the second part of Eq. (18b)., continue for the bulk tensors $s=t,t-1,\dotsc,1$ with

[TABLE]

and end at the left boundary with $D^{(0)}_{i\to j}(\sigma_{i}^{0}|\sigma_{j}^{0}):=\bar{C}^{(0)}_{i\to j}(\sigma_{i}^{0})U^{(1)}(\sigma_{j}^{0})\Lambda^{(1)}$.

To finally obtain the evolved edge message in canonical form (21), in a sweep from left to right, we go from the index assignment $(\sigma_{i}^{0}\sigma_{j}^{0})(\sigma_{i}^{1}|\sigma_{j}^{1})\dots(\sigma_{i}^{t}|\sigma_{j}^{t})(\sigma_{i}^{t+1})$ in Eq. (26) to the canonical assignment $(\sigma_{j}^{0})(\sigma_{i}^{0}|\sigma_{j}^{1})(\sigma_{i}^{1}|\sigma_{j}^{2})\dots(\sigma_{i}^{t-1}|\sigma_{j}^{t})(\sigma_{i}^{t})(\sigma_{i}^{t+1})$. At the left boundary, we start with $D^{(0)}_{i\to j}(\sigma_{i}^{0}|\sigma_{j}^{0})\stackrel{{\scriptstyle\text{SVT}}}{{\approx:}}B^{(0)}_{i\to j}(\sigma_{j}^{0})\Lambda^{(0)}V^{(0)}(\sigma_{i}^{0})$, continue for the bulk tensors $s=1,2,\dotsc,t$ with

[TABLE]

and end at the left boundary with $\Lambda^{(t)}V^{(t)}(\sigma_{i}^{t})D^{(t+1)}_{i\to j}(\sigma_{i}^{t+1})\stackrel{{\scriptstyle\text{SVT}}}{{\approx:}}B^{(t+1)}_{i\to j}(\sigma_{i}^{t})\Lambda^{(t+1)}V^{(t+1)}(\sigma_{i}^{t+1})$ and $B^{(t+2)}_{i\to j}(\sigma_{i}^{t+1}):=\Lambda^{(t+1)}V^{(t+1)}(\sigma_{i}^{t+1})$. The matrix products for the different stages in the evolution and truncation of the edge message $\mu_{i\to j}$ are illustrated in Fig. 10.

V.   Evaluation of observables

As described in Sec. II.5, the joint probability of trajectories $\bar{\sigma}_{i}^{t}$ and $\bar{\sigma}_{j}^{t-1}$ for the vertices of an edge $(i,j)$ is given by the product of the two corresponding edge messages, $P(\bar{\sigma}_{i}^{t},\bar{\sigma}_{j}^{t})=\mu_{i\to j}(\bar{\sigma}_{i}^{t}|\bar{\sigma}_{j}^{t-1})\mu_{j\to i}(\bar{\sigma}_{j}^{t}|\bar{\sigma}_{i}^{t-1})$. Given matrix product representations (approximations) of these edge messages in canonical form (14), time-local observables and temporal correlation functions can be evaluated efficiently. In order to evaluate, for example, the probability $P(\sigma_{i}^{t},\sigma_{j}^{t})$ of the edge being in state $(\sigma_{i}^{t},\sigma_{j}^{t})$ at time $t$, we simply contract all indices that occur in both MPEMs as depicted in Fig. 11. The contractions can be started at the left boundary ($s=0$) in such a way that the total computation cost scales as $\mathcal{O}(tM^{3})$.

To this purpose we start at the left boundary, i.e., at time slice $s=0$, with the multiplication E^{(0)}(\sigma_{i}^{0},\sigma_{j}^{0}):=A^{(0)}_{i\to j}(\sigma_{j}^{0})\big{[}A^{(0)}_{j\to i}(\sigma_{i}^{0})\big{]}^{\intercal}, continue for the bulk tensors $s=0,1,\dotsc,t-2$ with

[TABLE]

and finish with E^{(t)}:=\sum_{\sigma_{i}^{t-1}}A^{(t)}_{i\to j}(\sigma_{i}^{t-1})\big{(}\sum_{\sigma_{j}^{t-1}}E^{(t-1)}(\sigma_{i}^{t-1},\sigma_{j}^{t-1})\big{[}A^{(t)}_{j\to i}(\sigma_{j}^{t-1})\big{]}^{\intercal}\big{)} and the final step E^{(t+1)}(\sigma_{i}^{t},\sigma_{j}^{t}):=A^{(t+1)}_{i\to j}(\sigma_{i}^{t})\big{(}E^{(t)}\big{[}A^{(t+1)}_{j\to i}(\sigma_{j}^{t})\big{]}^{\intercal}\big{)}. Now, $E^{(t+1)}(\sigma,\sigma^{\prime})$ are scalars ($1\times 1$ matrices because $M_{t+2}=1$) that yield the desired joint probability

[TABLE]

During the time-evolution of the MPEMs (Sec. IV) it is natural to normalize the edge messages according to the 2-norm. In that case, it can be necessary to add a renormalization in Eq. (29) in order to avoid the generation of very large matrix elements, e.g., by replacing in each step $E^{(s)}$ by $E^{(s)}/\left\|E^{(s)}\right\|$.

VI.   An improved truncation scheme

In order to recast the evolved state (22) into canonical form (21) and truncate bond dimensions, the first step in Sec. IV.2 was to reorthonormalize the tensors ($C\to\tilde{C}$) in a sweep from right to left before doing truncations in a subsequent sweep from left to right ($\tilde{C}\to\bar{C}$). The preparatory first sweep can be avoided as described in the following. This also reduces computation costs substantially from \mathcal{O}\big{(}M^{3z-3}\big{)} to $\mathcal{O}(M^{2z-1})$, where $M$ denotes MPEM bond dimensions and $z$ denotes vertex degrees.

In extension of Sec. III.3, let us first discuss how controlled truncations (of Schmidt coefficients) can be done when the basis $\{|a\rangle\}$ of the left part $\mathcal{L}$, containing sites $[1,\ell]$, is orthonormal while the basis $\{|b\rangle\}$ for the right part $\mathcal{R}$, containing sites $[\ell+1,L]$, is not orthonormal. Recall that the purpose of the preparatory first sweep in Sec. IV.2 was to orthonormalize both bases. For notational simplicity, let us again use the example of a quantum many-body system with $L$ lattice sites. We are given a state

[TABLE]

with orthonormal states $\langle a|a^{\prime}\rangle=\delta_{a,a^{\prime}}$ for the left part and arbitrary states $|b\rangle$ for the right part. Now one can find a good approximation $|\psi_{\text{trunc}}\rangle$ of $|\psi\rangle$ (with respect to the 2-norm distance) by diagonalization of the reduced density matrix $\hat{\varrho}_{\mathcal{L}}$ for the left part. This density matrix,

[TABLE]

is obtained by a partial trace over $\mathcal{R}$. With the overlap matrix $F_{b,b^{\prime}}:=\langle b^{\prime}|b\rangle$ we have $[\varrho_{\mathcal{L}}]_{a,a^{\prime}}:=\langle a|\hat{\varrho}_{\mathcal{L}}|a^{\prime}\rangle=[\psi F\psi^{\dagger}]_{a,a^{\prime}}$ which can be diagonalized according to

[TABLE]

where, as in Sec. III.3, $U$ is a unitary matrix and $\Lambda$ is the diagonal matrix containing descendingly ordered Schmidt coefficients $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{M^{\prime}}$ (square roots of density matrix eigenvalues).

[TABLE]

The approximation after truncating all Schmidt coefficients $\lambda_{k}$ with $k>M$ is

[TABLE]

The orthogonality and norms of the states $\{|\tilde{k}_{\mathcal{R}}\rangle\}$ follow as

[TABLE]

The truncation (34b) is suitable for an MPS $|\psi\rangle$ (and analogously for MPEM) given in the form

[TABLE]

This is the form (31) with orthonormal basis states |a\rangle=\sum_{n_{1},\dotsc,n_{\ell}}\big{[}Y_{1}^{n_{1}}\dotsb Y_{\ell}^{n_{\ell}}\big{]}_{1,a}\,|n_{1}\dots n_{\ell}\rangle for $\mathcal{L}$ [cf. Eq. (20)], non-orthogonal states |b\rangle=\sum_{n_{\ell+1},\dotsc,n_{L}}\big{[}A_{\ell+1}^{n_{\ell+1}}\dotsb A_{L}^{n_{L}}\big{]}_{b,1}\,|n_{\ell+1}\dots n_{L}\rangle for $\mathcal{R}$, and $\psi=\mathbbm{1}$. The overlap matrix $F$ can be computed in an iteration with $F_{L}:=\mathbbm{1}$ and $F_{j}:=\sum_{n}A^{n}_{j+1}F_{j+1}(A^{n}_{j+1})^{\dagger}$ for $j=L-1,\dotsc,\ell$. With a diagonalization of $F\equiv F_{\ell}$, as in Eq. (33), we then obtain an optimally truncated state $|\psi_{\text{trunc}}\rangle$ in MPS form.

With this type of truncation, the preparatory first orthonormalization sweep for the evolved MPEM (right to left, $C\to\tilde{C}$), described in Sec. IV.2, can be avoided and computation costs can be reduced substantially. Going from right to left, one can first compute all overlap matrices $F^{(s)}$ for the evolved MPEM (22) with

[TABLE]

as illustrated in Fig. 12(a). With these, we can directly truncate the evolved MPEM, in a sweep from left to right. At the left boundary, we start with

[TABLE]

where we diagonalize and truncate (diagt) with the diagonalization as in Eq. (33) and the truncation as in Eq. (34b). The new tensor $\bar{C}^{(0)}_{i\to j}$ obeys the left orthonormality constraint [cf. Eq. (36)] and we need to include the matrix $X^{(0)}$ right of $\bar{C}^{(0)}_{i\to j}$ to keep the matrix product invariant such that $C^{(0)}_{i\to j}(\sigma_{i}^{0})\stackrel{{\scriptstyle\text{trunc}}}{{\approx}}\bar{C}^{(0)}_{i\to j}(\sigma_{i}^{0})X^{(0)}$ in Eq. (22). For the bulk tensors $s=1,2,\dotsc,t$, we continue with

[TABLE]

See Fig. 12(b). We end at the right boundary with $\bar{C}^{(t+1)}_{i\to j}(\sigma_{i}|\sigma_{j}):=X^{(t)}C^{(t+1)}_{i\to j}(\sigma_{i}|\sigma_{j})$ such that all $\bar{C}$-tensors except $\bar{C}^{(t+1)}$ now obey the left orthonormality constraint.

Truncating in this way, we arrive at Eq. (25) of the conceptually simpler truncation scheme of Sec. IV.2 and can continue in the same way as described there, rearranging vertex variables in the evolved MPEM to bring it into the canonical form (21). The major advantage of the more elaborate density matrix truncation scheme is the reduced computation cost. With bond dimensions $M$ of the MPEM (14) before the evolution step, the $C$-tensors of the exactly evolved MPEM (22) have increased bond dimensions $M^{z-1}$, where $z$ is the degree of vertex $i$. The computationally most expensive step in the simple truncation scheme of Sec. IV.2 is the singular value decomposition (24) for the orthonormalization of the $C$-tensors; the cost scales as \mathcal{O}\big{(}M^{3z-3}\big{)}. The most expensive steps in the more efficient density matrix truncation scheme are the computations of overlap matrices (37) and the reduced density matrices (39a). Exploiting the structure of the $C$-tensors (23), these operations can be accomplished with a cost $\mathcal{O}(M^{2z-1})$ if we assume that the bond dimensions of tensors $\bar{C}$ of the evolved MPEM after truncation are similar to that of the $A$-tensors in the original MPEM (14), i.e., approximately $M$. This is generally the case. When we fix a truncation threshold $\varepsilon$ for the discarded 2-norm weight as discussed in Sec. IV.2, bond dimensions evolve smoothly in time.

VII.   Computation costs in Glauber-Ising dynamics

As the number of possible trajectories on a vertex increases exponentially in time, the computation costs for exact dynamic belief propagation (11) also grow exponentially. As discussed in the introduction, this has, so far, substantially limited the applicability of the dynamic cavity method. The MPEM approach allows for a controlled reduction of the computational complexity, exploiting the decay of temporal correlations to truncate unimportant components of the edge messages.

Concerning computation costs, a decisive question is now, what bond dimensions $M_{s}$ in MPEMs are required to achieve a certain approximation accuracy in comparison to the exact evolution. Also, how does the required bond dimension evolve with time? Generally, it is to be expected that the maximum bond dimension, $\max_{s}M_{s}$, will converge to a constant as a function of time if connected temporal correlations decay exponentially. In simulations, it is often favorable not to fix the truncation dimension, but to fix instead a threshold $\varepsilon$ and to truncate all singular values $\lambda_{k}$ with $\lambda_{k}^{2}<\varepsilon\cdot\sum_{q}\lambda_{q}^{2}$. The threshold $\varepsilon$ controls the 2-norm loss in each truncation. Allowing bond dimensions to evolve accordingly, avoids wasting computation time and gives a measure for the information-theoretic complexity of the edge messages.

As an example, we consider Glauber dynamics of the kinetic Ising model Glauber1963-4 on random regular graphs with vertex degree $z=3$ in the thermodynamic limit. Specifically, Ising spins interact ferromagnetically with local transition matrices

[TABLE]

At time $t=0$, all spins have magnetization $\langle\sigma^{0}_{i}\rangle=1/2$, i.e., $p_{i}(\uparrow{}\!)=3/4$.

Figure 13 shows the evolution of maximum bond dimensions with time. They increase with decreasing truncation threshold $\varepsilon$. As expected, they become constant at longer times. The number of tensors $A^{(0)},\dotsc,A^{(t+1)}$ in an MPEM increases linearly with $t$, and we need to sweep a few times through the matrix product in each iteration. With converging bond dimensions, this implies that the computation cost per iteration grow only linearly in time, instead of exponentially. For fixed $\varepsilon$, the largest bond dimensions are needed for inverse temperature $\beta=1/2$. This is so, because it is the $\beta$ closest to the phase transition in the Ising model. For this temperature, the bond dimensions do actually not yet show saturation on the time interval displayed in Fig. 13. They will converge at larger $t$.

VIII.   Models with vertex-state dependence

Except for the algorithms described in Secs. IV and VI our description was generic in the sense that the local transition matrices $w_{i}=w_{i}(\sigma_{i}^{t+1}|{\boldsymbol{\sigma}}_{\partial i}^{t})$ were allowed to depend on the time-$t$ state $\sigma_{i}^{t}$ of vertex $i$ itself, in addition to the states on neighboring vertices. However, in Secs. IV and VI we considered the case where $\partial i$ does not contain $i$ itself, i.e., that $w_{i}$ is independent of $\sigma_{i}^{t}$. Here, we generalize to MPEM algorithms that allow for the dependence on $\sigma_{i}^{t}$. This is important for many applications like stochastic models for the dynamics of infectious diseases Murray1989 ; Dangerfield2009-6 ; Karrer2010-82 or of opinions is social networks Castellano2009-81 . Naturally, transition matrices with dependence on $\sigma_{i}^{t}$ also arise in time-discretized versions of stochastic continuous-time dynamics, i.e., all models that have a well-defined continuum-time limit as, for example, the Glauber dynamics of Ising spin systems Glauber1963-4 . In particular, when decreasing the time step $\Delta t$, one should have

[TABLE]

Vertex $i$ is now contained in its neighborhood $\partial i=\{i,j,k_{1},\dotsc,k_{z-1}\}$, where $z$ denotes the vertex degree. In this more general scenario, the evolved MPEM (22) cannot be constructed in an entirely local fashion anymore. Instead, the $\tilde{C}$-tensors are constructed in a sweep from right ($s=t+1$) to left ($s=0$), imposing at the same time the right orthonormality constraint, \sum_{\sigma,\sigma^{\prime}}\tilde{C}^{(s)}_{i\to j}(\sigma|\sigma^{\prime})\big{[}\tilde{C}^{(s)}_{i\to j}(\sigma|\sigma^{\prime})\big{]}^{\dagger}=\mathbbm{1} [cf. Eq. (18b)].

We start at the right boundary ($s=t+1$) by doing a singular value decomposition (SVD) of the tensor

[TABLE]

in order to obtain the tensor $\tilde{C}^{(t+1)}(\sigma_{i}^{t+1}|\sigma_{j}^{t})$, where $\Lambda^{(t+1)}$ is the diagonal matrix of singular values, $\tilde{C}^{(t+1)}_{i\to j}$ obeys the right orthonormality constraint, and, similarly, \sum_{\sigma}\big{[}U^{(t+1)}(\sigma)\big{]}^{\dagger}U^{(t+1)}(\sigma)=\mathbbm{1}. For time $s=t$, the process continues with

[TABLE]

For the bulk tensors, $s=t,t-1,\dotsc,1$, the corresponding equations read very similarly

[TABLE]

as illustrated in Fig. 14. The process ends at the left boundary with the assignment

[TABLE]

The remaining truncation and vertex variable reordering sweeps, done in order to transform the non-canonical form (22) into the canonical form (21), can be executed exactly as discussed in the somewhat simpler situation covered in Sec. IV.2.

The more efficient density matrix truncation scheme discussed in Sec. VI, similarly, can be adapted to the situation where transition matrix $w_{i}=w_{i}(\sigma_{i}^{t+1}|{\boldsymbol{\sigma}}_{\partial i}^{t})$ depends on $\sigma_{i}^{t}$. Like the unitaries $U^{(s)}$ in Eqs. (42)-(44), the $X^{(s)}$-matrices in Eqs. (38) and (39) of the density matrix truncation scheme will then depend on $\sigma_{i}^{s-1}$.

Continuous-time stochastic dynamics on locally tree-like networks can be simulated with the described approach after discretization of the time axis using a small time step $\Delta t$. It should also be possible to, alternatively, work with continuous matrix products Verstraete2010-104 ; Osborne2010-105 as introduced in the context of condensed matter physics . This however, entails some technical difficulties and one can expect that using discrete-time MPEMs and a small time step should be the best option. Similarly, quantum many-body systems in continuous real-space have been treated efficiently using MPS with a sufficiently fine space discretization Stoudenmire2012-109 ; Dolfi2012-109 . Other approximative schemes to simulate continuous-time stochastic dynamics are the dynamical replica analysis Mozeika2008-41 that captures macroscopic observables and the cavity master equation method Aurell2017-95 that operates on local conditional probabilities.

IX.   Discussion

The described MPEM algorithm for the simulation of stochastic dynamics is based on the cavity method, applicable for locally tree-like interaction graphs, and on the matrix product approximation for edge messages which exploits the decay of (connected) temporal correlations. The MPEM method lifts restrictions of earlier approaches for the solution of dynamic cavity equations, mentioned in the introduction. It can also be used to simulate in the thermodynamic limit. Compared to Monte Carlo simulations, errors decrease much faster as a function of the invested computation time. This has been demonstrated for Glauber dynamics of the kinetic Ising model in Ref. Barthel2018-97 . Hence, the MPEM method is particularly suited for the precise analysis of temporal correlations, decay processes, and low-probability events. For the Glauber-Ising dynamics and fixed truncation thresholds, we have seen here that MPEM bond dimensions converge to a constant as a function of time. This should generally apply when connected temporal correlations decay exponentially. It implies that required computation costs per iteration grow only linearly instead of exponentially in time.

There are several ways in which the MPEM method can be developed further. Here, we have discussed a more efficient truncation scheme that reduces computation costs considerably from $\mathcal{O}(M^{3z-3})$ to $\mathcal{O}(M^{2z-1})$, and we have generalized the method to models where the transition probabilities for a vertex depend explicitly on the vertex state at the previous time step, in addition to the states for its nearest neighbors.

Efficient codes for the MPEM algorithms, including the optimized truncation scheme of Sec. VI, are available from the author and under www.manyparticle.org/$\sim$barthel/mpem. They were, for example, employed for the simulations in Ref. Barthel2018-97 . We are also happy to collaborate on specific applications.

I gratefully acknowledge discussions with Silvio Franz and Caterina De Bacco that initiated this work and support through US Department of Energy grant DE-SC0019449.