Matrix Product Edge Message (MPEM) Algorithm

• Matrix Product Edge Message (MPEM) algorithm is a numerical method that approximates high-dimensional edge messages with matrix-product states for efficient simulation of stochastic dynamics on tree-like graphs.
• It converts the exponential scaling of dynamic cavity equations into a tractable process using local tensor contractions and SVD-based truncation techniques.
• The method outperforms Monte Carlo simulations by accurately capturing rare events and long-range temporal correlations in models such as the kinetic Ising model and neural networks.

The Matrix Product Edge Message (MPEM) algorithm is a numerical method designed to efficiently solve the dynamic cavity equations describing the stochastic dynamics of classical degrees of freedom on graphs, particularly for locally tree-like topologies. Its core innovation is the approximation of high-dimensional edge messages—conditional probability distributions of entire vertex trajectories—by matrix-product states, enabling long-time, high-accuracy simulations of processes such as Glauber dynamics, neural network evolution, or Boolean network dynamics. The approach draws upon concepts from quantum many-body theory to enable controlled truncation and efficient updating of these trajectory distributions, providing improved scaling and accuracy over traditional Monte Carlo techniques, especially for low-probability events and long-range temporal observables (Barthel, 2019, Barthel et al., 2015).

1. Edge Messages and the Dynamic Cavity Equation

In the dynamic cavity framework for stochastic processes on graphs, each vertex $i$ assumes a state $\sigma_i^t \in \{1,\dots,d\}$ at discrete time $t$. The full trajectory of a vertex up to time $t$ is denoted $\sigma_i^{0\ldots t}$. The transition dynamics are specified by local transition probabilities $w_i(\sigma_i^{t+1} \mid \sigma_{\partial i}^t)$, where $\partial i$ are the neighbors of $i$, and initial distributions $p_i(\sigma_i^0)$. When a directed edge $(i \to j)$ is removed, the graph splits into complementary subgraphs, and the edge message $$\mu_{i\to j}(\sigma_i^{0\ldots t+1} \mid \sigma_j^{0\ldots t})$$ encodes the conditional probability of a full trajectory of $i$, given that of $j$, on the subtree devoid of $j$'s influence.

The exact update for edge messages under the Bethe (locally tree-like) approximation is defined recursively: $$\mu_{i\to j}(\sigma_i^{0\ldots t+1}\mid\sigma_j^{0\ldots t}) =\sum_{\{\sigma_k^{0\ldots t}\}_{k\in\partial i\setminus\{j\}}} p_i(\sigma_i^0) \prod_{s=0}^t w_i(\sigma_i^{s+1} \mid \sigma_{\partial i}^s) \prod_{k\in\partial i\setminus\{j\}} \mu_{k\to i}(\sigma_k^{0\ldots t} \mid \sigma_i^{0\ldots t-1})$$ This relation is exact for trees and constitutes the dynamic cavity or Bethe–Peierls equation for general (locally tree-like) graphs. However, the necessary storage and computational cost for these objects scales exponentially as $\mathcal{O}(d^{2t+1})$ with time $t$.

2. Matrix-Product Edge Message Ansatz

To mitigate exponential scaling, the MPEM algorithm approximates edge messages as matrix-product states over the history of the system. The canonical MPEM decomposition is: $$\mu_{i\to j}(\sigma_i^{0\ldots t} \mid \sigma_j^{0\ldots t-1}) \approx A^{(0)}(\sigma_j^0) \left[\prod_{s=1}^{t-1}A^{(s)}(\sigma_i^{s-1}\mid\sigma_j^s)\right] A^{(t)}(\sigma_i^{t-1})A^{(t+1)}(\sigma_i^t)$$ Each $A^{(s)}$ is a three-index tensor of size $M_s \times d \times d \times M_{s+1}$, where the $M_s$ are the bond dimensions controlling the approximation accuracy. This ansatz enables representing the exponentially large conditional probability as a product of lower-dimensional tensors, allowing controlled truncation of the information retained at each time "bond" by adjusting the bond dimensions (Barthel, 2019, Barthel et al., 2015).

3. Recursive Update and Truncation Procedures

The MPEM update proceeds by local tensor contractions and SVD-based truncations:

1. Evolution step: The new, non-canonical tensors $C^{(s)}$ for time $t+1$ are constructed by contracting the local transition kernel $w_i$ with all neighboring edge MPEMs $A_{k\to i}^{(s)}$:

$$C^{(s)}_{i\to j}(\sigma_i^{s}\mid\sigma_j^{s-1}) = \sum_{\{\sigma_k^{s-1}\}} w_i(\sigma_i^{s} \mid \sigma_j^{s-1},\{\sigma_k^{s-1}\}) \bigotimes_{k\in\partial i\setminus\{j\}} A_{k\to i}^{(s)}(\sigma_k^{s-1}\mid\sigma_i^s)$$

1. Canonicalization and truncation: As contraction increases the bond dimension exponentially (from $M_s$ to $M_s^{z-1}$ for degree $z$), SVD-based truncations restore manageable bond dimension and canonical form. Two schemes are used:
   • Simple SVD sweep: Alternating right-to-left and left-to-right sweeps orthonormalize and truncate singular values below a fixed threshold $\epsilon$.
   • Optimized density-matrix truncation: Employs the computation of overlap matrices to select an optimal subspace, significantly reducing computational overhead to $\mathcal{O}(M^{2z-1})$ per edge (from $\mathcal{O}(M^{3z-3})$).

The $2$-norm truncation error per step is bounded by the sum of discarded squared singular values $\sum_{k>M_s}\lambda_k^2 \le \epsilon$. One typically fixes $\epsilon$ (e.g., $10^{-8}$ to $10^{-12}$) and allows $M_s$ to adapt dynamically.

4. Computational Complexity and Error Scaling

The computational cost per directed edge per time step is $\mathcal{O}(d^2 M^{2z-1})$, where $d$ is the local state space dimension and $z$ is the average degree. In sparse, locally tree-like graphs, the total cost scales linearly with time: $\mathcal{O}(|E| d^2 M^{2z-1} T)$, where $|E|$ is the number of edges and $T$ the number of time steps.

Empirical results demonstrate that for systems such as the Glauber–Ising model, the bond dimension $M_s$ saturates over time for fixed $\epsilon$, resulting in linear overall scaling with $T$ rather than exponential growth. MPEM truncation yields error scaling as $\mathcal{O}(\epsilon T)$ in the $2$-norm, contrasting favorably with Monte Carlo's statistical error $\mathcal{O}(1/\sqrt{N_s})$ for $N_s$ samples, especially in resolving long-time correlations and rare events (Barthel, 2019, Barthel et al., 2015).

5. Extraction of Observables

With MPEMs for all directed edges, local and multi-site observables are computed by tensor contractions:

• Single-site marginals: For vertex $i$,

$$P(\sigma_i^t) \propto \sum_{\substack{\sigma_i^{0:t-1} \ \{\sigma_j^{0:t-1}\}_{j\in\partial i}}} p_i(\sigma_i^0) \prod_{s=0}^{t-1} w_i(\sigma_i^{s+1} \mid \{\sigma_j^s\}) \prod_{j\in\partial i} \mu_{j\to i}(\sigma_j^{0:t-1} \mid \sigma_i^{0:t-2})$$

The magnetization is $m_i(t) = \langle \sigma_i^t \rangle$.

• Two-point spatial correlations: For edge $(i,j)$,

$$P(\sigma_i^t, \sigma_j^t) = \sum_{\sigma_i^{0:t-1},\sigma_j^{0:t-1}} \mu_{i\to j}(\sigma_i^{0:t} \mid \sigma_j^{0:t-1}) \mu_{j\to i}(\sigma_j^{0:t} \mid \sigma_i^{0:t-1})$$

Computation proceeds by sequential contraction of the MPEM tensors at all time steps.

• Temporal correlations: By grooming the MPEM, quantities such as $$C(t,s) = \langle \sigma_i^t \sigma_i^s \rangle - \langle \sigma_i^t \rangle \langle \sigma_i^s \rangle$$ are accessed.

Evaluation of these observables is deterministic, non-stochastic, and scales efficiently in time and bond dimension (Barthel, 2019, Barthel et al., 2015).

6. Applications, Optimizations, and Extensions

MPEM has been applied to prototypical spin systems such as the kinetic Ising model with Glauber dynamics on random regular graphs. The transition kernel for the Ising application is: $$w_i(\sigma_i^{t+1}\mid\sigma_{\partial i}^t) = \frac{1}{Z} \exp\left[\beta \sigma_i^{t+1} \sum_{j\in\partial i} \sigma_j^t\right]$$ with initial conditions $p_i(\sigma_i^0)=\frac12(1+\sigma_i^0 m_0)$.

Algorithmic optimizations include the adoption of density-matrix truncation to lower cost, the extension to continuous-time dynamics via right-to-left SVD sweeps including the self-state index, and easy parallelization over edges or vertices.

Empirical results indicate that for a threshold $\epsilon=10^{-10}$ and system size $T=200$, each edge update requires only milliseconds. For comparable computational effort, MPEM achieves much lower errors ($|\delta m(t)| \sim 10^{-6}$ with $M\sim 20$) than Monte Carlo, especially at late times and for small signals, which would require $10^6$ trajectories in MC (Barthel, 2019).

7. Significance and Impact

The MPEM algorithm provides a controlled, efficient, and noise-free framework for simulating the stochastic dynamics of classical systems on graphs with locally tree-like structure. Unlike traditional Monte Carlo, it enables direct and accurate computation of low-probability observables and is particularly well-suited to long-time studies where MC error is prohibitive. Its matrix-product formulation is closely related to quantum many-body methods for open system dynamics (e.g., matrix product states and operators), making it a point of intersection between statistical mechanics and tensor network techniques. The approach's public implementation further facilitates adoption for a wide range of network dynamics models (Barthel, 2019, Barthel et al., 2015).