As mentioned in the introduction, it is not possible to compute the entire eigenspectrum of huge sparse matrices using standard dense solvers based on Givens rotation or Householder reflections for memory limitations. Instead, sparse and iterative eigenvalue solvers [1]}, [2]}, [3]} are to be implemented to extract a limited portion of the eigenspectrum and to approximate the evolution operator by its projection on the reduced eigenbasis [4]}. In practice, the basis encompasses the \(k\)  eigenvectors with smallest eigenvalues. By doing so, a reduction of the model is performed and the symmetric transition rate matrix \(\mathbf {A}^a\)  in Eq.(REF ) is approximated by
\(\mathbf {A}^{(k)} = \mathbf {SR} \left[\sum _{h=0}^k \mathbf {\varphi }_h\mathbf {\varphi }^T_h\lambda _h \right] \mathbf {SR}^{-1},\) 
