Cancer progression can be described by continuous-time Markov chains whose state space grows exponentially in the number of somatic mutations. The age of a tumor at diagnosis is typically unknown. Therefore, the quantity of interest is the time-marginal distribution over all possible genotypes of tumors, defined as the transient distribution integrated over an exponentially distributed observation time. It can be obtained as the solution of a large linear system. However, the sheer size of this system renders classical solvers infeasible. We consider Markov chains whose transition rates are separable functions, allowing for an efficient low-rank tensor representation of the linear system’s operator. Thus we can reduce the computational complexity from exponential to linear. We derive a convergent iterative method using low-rank formats whose result satisfies the normalization constraint of a distribution. We also perform numerical experiments illustrating that the marginal distribution is well approximated with low rank.

The solution of parameter-dependent linear systems, by classical methods, leads to an arithmetic effort that grows exponentially in the number of parameters. This renders the multigrid method, which has a well understood convergence theory, infeasible. A parameter-dependent representation, e.g., a low-rank tensor format, can avoid this exponential dependence, but in these it is unknown how to calculate the inverse directly within the representation. The combination of these representations with the multigrid method requires a parameter-dependent version of the classical multigrid theory and a parameter-dependent representation of the linear system, the smoother, the prolongation and the restriction. A derived parameter-dependent version of the smoothing property, fulfilled by parameter-dependent versions of the Richardson and Jacobi methods, together with the approximation property prove the convergence of the multigrid method for arbitrary parameter-dependent representations. For a model problem low-rank tensor formats represent the parameter-dependent linear system, prolongation and restriction. The smoother, a damped Jacobi method, is directly approximated in the low-rank tensor format by using exponential sums. Proving the smoothing property for this approximation guarantees the convergence of the parameter-dependent method. Numerical experiments for the parameter-dependent model problem, with bounded parameter value range, indicate a grid size independent convergence rate.

Continuous-time Markov chains describing interacting processes exhibit a state space that grows exponentially in the number of processes. This state-space explosion renders the computation or storage of the time-marginal distribution, which is defined as the solution of a certain linear system, infeasible using classical methods. We consider Markov chains whose transition rates are separable functions, which allows for an efficient low-rank tensor representation of the operator of this linear system. Typically, the right-hand side also has low-rank structure, and thus we can reduce the cost for computation and storage from exponential to linear. Previously known iterative methods also allow for low-rank approximations of the solution but are unable to guarantee that its entries sum up to one as required for a probability distribution. We derive a convergent iterative method using low-rank formats satisfying this condition. We also perform numerical experiments illustrating that the marginal distribution is well approximated with low rank.

Motivation: We consider continuous-time Markov chains that describe the stochastic evolution of a dynamical system by a transition-rate matrix Q which depends on a parameter θ. Computing the probability distribution over states at time t requires the matrix exponential exp(tQ), and inferring θ from data requires its derivative ∂ exp(tQ)/∂θ. Both are challenging to compute when the state space and hence the size of Q is huge. This can happen when the state space consists of all combinations of the values of several interacting discrete variables. Often it is even impossible to store Q. However, when Q can be written as a sum of tensor products, computing exp(tQ) becomes feasible by the uniformization method, which does not require explicit storage of Q.Results: Here we provide an analogous algorithm for computing ∂ exp(tQ)/∂θ, the differentiated uniformization method. We demonstrate our algorithm for the stochastic SIR model of epidemic spread, for which we show that Q can be written as a sum of tensor products. We estimate monthly infection and recovery rates during the first wave of the COVID-19 pandemic in Austria and quantify their uncertainty in a full Bayesian analysis. Availability: Implementation and data are available at https://github.com/spang-lab/TenSIR.