Fast and Numerically Stable Implementation of Rate Constant Matrix Contraction Method

TL;DR

This work tackles the challenge of fast and numerically stable simulation of chemical kinetics on large reaction-path networks by accelerating the greedy first step of the reformulated rate constant matrix contraction (RCMC) method. It reframes Step 1 as a greedy MAP inference for determinantal point processes and introduces StableLazyFastGreedy, a numerically stable modification that avoids subtracting like-signed numbers, aided by segment-tree based compression to enable efficient updates. To boost performance, the authors also develop RelaxedStableLazyFastGreedy, which permits controlled like-sign subtractions under an error bound, substantially speeding up the greedy phase without compromising stability. Across six large-scale chemical networks, the proposed methods demonstrate clear numerical stability and significant speedups over the original RCMC, making the approach practical for networks with $n$ in the tens of thousands and massive condition numbers.

Abstract

The rate constant matrix contraction (RCMC) method, proposed by Sumiya et al. (2015, 2017), enables fast and numerically stable simulations of chemical kinetics on large-scale reaction path networks. Later, Iwata et al. (2023) mathematically reformulated the RCMC method as a numerical algorithm to solve master equations whose coefficient matrices, known as rate constant matrices, possess the detailed balance property. This paper aims to accelerate the RCMC method. The bottleneck in the RCMC method lies in the greedy selection of steady states, which is actually equivalent to the greedy algorithm for the MAP inference in DPPs under cardinality constraints. Hemmi et al. (2022) introduced a fast implementation of the greedy DPP MAP inference, called LazyFastGreedy, by combining the greedy algorithm of Chen et al. (2018) with the lazy greedy algorithm by Minoux (1978), a practically efficient greedy algorithm that exploits the submodularity of the objective function. However, for instances arising from chemical kinetics, the straightforward application of LazyFastGreedy suffers from catastrophic cancellations due to the wide range of reaction time scales. To address this numerical instability, we propose a modification to LazyFastGreedy that avoids the subtraction of like-sign numbers by leveraging the properties of rate constant matrices and the connection of the DPP MAP inference to Cholesky decomposition. For faster implementation, we utilize a segment tree, a data structure that manages one-dimensional arrays of elements in a semigroup. We also analyze the increase in relative errors caused by like-sign subtractions and permit such subtractions when they do not lead to catastrophic cancellations, aiming to further accelerate the process. Using real instances, we confirm that the proposed algorithm is both numerically stable and significantly faster than the original RCMC method.

Figures (3)

• Figure 1: The iteration counts $k$ over the matrices size $n$ versus the running times of \ref{['step:1']}. The cross marks with the horizontal coordinate $k$ indicate the smallest $k$ where the algorithm selected an incorrect element as $s^{(k)}$ due to catastrophic cancellation. The vertical black lines indicate the value of $k$ corresponding to $t_\mathrm{max} = 86400sec$.
• Figure 3: The iteration counts $k$ over the matrices size $n$ versus $M_\mathrm{offdiag}$ (defined in \ref{['sec:lazy-greedy-runtime']}). The gray areas indicate the ranges of possible values of $M_\mathrm{offdiag}$ for every $k$ given by \ref{['eq:mlazy-bound']}. The cross marks and vertical lines indicate the same as \ref{['fig:runtime']}.
• Figure : pivoted Gaussian elimination

Theorems & Definitions (30)

• Proposition 2.1: Iwata2023
• Theorem 2.2
• proof
• Proposition 2.3: Iwata2023
• Theorem 3.1
• Lemma 3.2
• proof
• proof : Proof of \ref{['thm:equivalence-cholesky-crout']}
• Proposition 3.3
• Proposition 3.4