Nearest Reversible Markov Chains with Sparsity Constraints: An Optimization Approach

Technical Summary: Nearest Reversible Markov Chains with Sparsity Constraints

Problem Formulation
Markov chains are fundamental in probability, statistics, and computational science, particularly within Markov Chain Monte Carlo (MCMC) methods like Metropolis–Hastings, where reversibility is a critical property ensuring convergence to a target distribution. However, many practical applications—such as those in numerical linear algebra, computational chemistry, and data-driven modeling—yield non-reversible chains due to noisy observations or specific structural constraints.

The central problem addressed is how to modify a given non-reversible, sparse transition matrix $P$ as minimally as possible to obtain a reversible transition matrix $X$. This is formulated as a matrix nearness problem:
$$\min_{X} \frac{1}{2} \|X - P\|_F^2$$
subject to $X$ being a stochastic matrix ($X \ge 0, X\mathbf{1} = \mathbf{1}$), satisfying the detailed balance equations with respect to a stationary distribution $\pi$ ($D_\pi X = X^\top D_\pi$), and preserving a prescribed sparsity pattern $\Omega_M$. The authors specifically focus on the sparse case, where the number of non-zero entries is proportional to $n$ rather than $n^2$, which is typical in high-dimensional applications.

Methodology
The authors propose a quadratic programming (QP) approach to solve this optimization problem. The methodology proceeds through several key steps:

1. Variable Transformation: To handle the detailed balance constraint efficiently, the problem is reformulated using a change of variables $Y = D_{\hat{\pi}} X D_{\hat{\pi}}^{-1}$, where $\hat{\pi} = \pi^{1/2}$. This transforms the reversibility condition into a symmetry constraint ($Y = Y^\top$) and the stochasticity condition into an eigenvector constraint ($Y\hat{\pi} = \hat{\pi}$).
2. Sparsity Handling: The problem is restricted to matrices sharing the sparsity pattern of a given symmetric, irreducible binary matrix $M$. The authors define a vectorization strategy that reduces the number of decision variables by considering only the unique entries in the upper triangular part of the sparse matrix.
3. Convexity and Feasibility: The paper proves that the feasible set of matrices satisfying the sparsity, symmetry, and stochasticity constraints is convex and non-empty (Propositions 1 and 2). Furthermore, the objective function is shown to be strongly convex (Proposition 3), guaranteeing a unique global solution.
4. Handling Reducibility: For reducible Markov chains (those with multiple ergodic classes or transient states), the algorithm first decomposes the chain into ergodic classes. The optimization is then solved independently for each ergodic class, as reversibility is well-defined within these blocks. Transient states are discarded as they do not affect the detailed balance condition on the support of the stationary distribution.
5. Algorithmic Implementation: The reformulated problem is cast into a standard QP form and solved using off-the-shelf solvers, specifically MATLAB's quadprog and the commercial Gurobi solver. The authors provide a MATLAB implementation (nearest_sparse_reversible.m) that automates the decomposition, optimization, and reassembly of the reversible chain.

Key Results
Numerical experiments were conducted on both synthetic sparse matrices and real-world molecular dynamics data (butane torsional motion and Fs-peptide folding).

• Solver Performance: Comparisons between quadprog and Gurobi revealed that Gurobi is significantly more efficient, achieving solve times approximately 3–4 times faster than quadprog (median times of ~0.007s vs. ~0.023s). Both solvers satisfied reversibility and stationarity constraints to machine precision.
• Sparsity vs. Reversibility: The results confirm that enforcing reversibility on a sparse matrix generally increases the number of non-zero entries. The optimization preserves the original sparsity pattern but redistributes weights to satisfy detailed balance, often requiring the activation of symmetric counterparts to existing transitions.
• Comparison with Heuristics: The proposed optimization approach yields a perturbation $\Delta = X - P$ with a significantly smaller Frobenius norm compared to standard reversibilization heuristics like Metropolis–Hastings or Barker's algorithm. In the molecular dynamics examples, the optimization reduced the perturbation norm by factors ranging from 2 to 5 compared to these heuristics.
• Scalability: The computational cost exhibits subquadratic scaling with respect to the problem dimension, making the approach viable for moderate-scale problems (up to hundreds of states) and suitable for real-time or near-real-time applications.
• Comparison with Dense Methods: When compared to a Riemannian optimization approach designed for dense matrices (without sparsity constraints), the proposed sparse method achieved a comparable perturbation norm (5.54 vs. 5.46) but was orders of magnitude faster (0.07s vs. 37.8s) due to the reduced number of variables and constraints.

Significance and Claims
The paper claims to provide a "principled way" to enforce reversibility and sparsity patterns in Markov chains. Its primary contributions are:

• Theoretical Rigor: Establishing the strong convexity of the problem, which ensures the uniqueness of the solution and independence from the specific optimization algorithm used.
• Practical Efficiency: Demonstrating that a QP-based approach is computationally superior to both heuristic reversibilization methods (which do not minimize perturbation) and dense Riemannian optimization methods (which ignore sparsity).
• Applicability: Validating the framework on data-driven modeling scenarios where transition matrices are estimated from noisy simulation data, showing that the method effectively recovers reversible chains while maintaining the computational efficiency required for sparse structures.

The authors conclude that the trade-off between sparsity and reversibility is an intrinsic structural property of Markov chains, and their framework offers a robust, scalable tool for managing this trade-off in applications ranging from MCMC to computational chemistry. Future work is suggested to explore regularized interior point methods and the simultaneous optimization of the transition probabilities and the underlying graph topology.