[Paper Review] Matrix Product Belief Propagation for reweighted stochastic dynamics over graphs

Stefano Crotti, Alfredo Braunstein|arXiv (Cornell University)|Mar 30, 2023

Complex Network Analysis Techniques4 citations

TL;DR

This paper introduces Matrix Product Belief Propagation (MPBP), a scalable method for computing rare-event statistics in Markov processes on graphs by extending the matrix product cavity method to handle arbitrary reweighting factors. It reduces the computational cost of node updates from exponential to polynomial in node degree, enabling efficient inference and large deviation analysis in models like SIRS epidemics and kinetic Ising systems.

ABSTRACT

Stochastic processes on graphs can describe a great variety of phenomena ranging from neural activity to epidemic spreading. While many existing methods can accurately describe typical realizations of such processes, computing properties of extremely rare events is a hard task. Particularly so in the case of recurrent models, in which variables may return to a previously visited state. Here, we build on the matrix product cavity method, extending it fundamentally in two directions: first, we show how it can be applied to Markov processes biased by arbitrary reweighting factors that concentrate most of the probability mass on rare events. Second, we introduce an efficient scheme to reduce the computational cost of a single node update from exponential to polynomial in the node degree. Two applications are considered: inference of infection probabilities from sparse observations within the SIRS epidemic model, and the computation of both typical observables and large deviations of several kinetic Ising models.

Motivation & Objective

• Address the challenge of computing observables for rare events in recurrent Markov processes on large networks, where standard Monte Carlo methods fail due to exponential sampling costs.
• Overcome the computational intractability of dynamic cavity methods on general graphs by exploiting matrix product state (MPS) representations to compress trajectory distributions.
• Enable the study of large deviations and atypical trajectories by incorporating arbitrary non-negative reweighting factors into the belief propagation framework.
• Develop a computationally efficient scheme that reduces the cost of a single node update from exponential to polynomial in node degree for a broad class of models.
• Provide a scalable framework applicable to both epidemic models (e.g., SIRS) and statistical physics models (e.g., Glauber dynamics) with practical inference and observables computation.

Proposed method

• Formulate the joint distribution of stochastic trajectories on graphs using a factorized form that includes transition probabilities and reweighting factors for rare-event biasing.
• Adopt the Matrix Product State (MPS) representation to parametrize the trajectory distributions of individual nodes, enabling compact encoding of high-dimensional trajectory spaces.
• Develop a belief propagation (BP) formalism—Matrix Product Belief Propagation (MPBP)—that operates on messages encoded as matrices, with fixed-point iteration to compute marginals.
• Introduce a recursive message-passing scheme that leverages intermediate aggregate variables (e.g., number of infected neighbors or local field) to reduce computational complexity.
• Apply tensor network contraction techniques to efficiently compute single-variable and pair marginals, time autocorrelations, and Bethe free energy for model learning.
• Use population dynamics in the thermodynamic limit to extend the method to heterogeneous or random graphs with variable degrees and parameters.

Figure 1: Marginal probabilities of free dynamics under the SIS model, comparison with models mentioned in the text. The main panels correspond to marginals for a single node of the graph, insets show the average absolute error over all nodes with respect to Monte Carlo simulations. Panels (a-c) com

Experimental results

Research questions

• RQ1 Can the matrix product cavity method be extended to handle arbitrary reweighting factors that concentrate probability on rare trajectories in recurrent stochastic processes?
• RQ2 What is the computational cost of node updates in existing MPS-based dynamic cavity methods, and can it be reduced from exponential to polynomial in node degree?
• RQ3 For which classes of models (e.g., epidemic or Ising-type) does the intermediate variable assumption hold, enabling efficient message computation?
• RQ4 How does the proposed MPBP method compare to standard belief propagation and Monte Carlo in computing both typical and large deviation observables?
• RQ5 Can the Bethe free energy computed via MPBP be used for maximum-likelihood parameter estimation in stochastic dynamics on graphs?

Key findings

• The proposed MPBP method reduces the computational cost of a single node update from 𝒪(M^{3z-3}) to 𝒪(M^6) for a wide class of models, where M is the matrix dimension and z is the node degree.
• The method enables accurate computation of both typical observables and large deviations in SIRS epidemic models with sparse observations, outperforming standard Monte Carlo in rare-event regimes.
• For kinetic Ising models, MPBP successfully computes time autocorrelations and magnetization dynamics under various reweighting schemes, including those favoring atypical spin configurations.
• The Bethe free energy approximation computed via MPBP allows for maximum-likelihood learning of model parameters such as infection rates or inverse temperature.
• The method achieves polynomial scaling in node degree by exploiting conditional independence structures, such as the local field or infection count as sufficient statistics.
• The Julia implementation of MPBP is publicly available, enabling reproducibility and extension to other stochastic processes on complex networks.

Figure 2: MPBP (solid line) with bond dimension $3$ correctly computes marginals of an SIS model defined on an Erdos-Renyi graph with $23$ nodes and average connectivity $4$ , $\lambda=0.15,\rho=0.12,\gamma=0.13$ . The state of a random half of the variables was observed at final time $T=10$ and use