The efficient computation of the matrix exponential \(e^{A t}\) is one fundamental building block for the simulation of time-dependent science and engineering problems. Motivated by applications in quantum thermodynamics and reaction–diffusion biological systems, we assume in this paper that \(A = T + V\) where both T and V  have compressible features, e.g., diagonal or low-rank structures in the physical or frequency domains. We compare two different numerical approaches that utilize the operator splitting technique for efficiently evaluating \(e^{A t }\) . In the first approach, the Lie–Trotter–Suzuki fractal decompositions are applied when \(e^{A t }\) is efficiently approximated by the products of \(e^{\alpha T}\) and \(e^{\beta V}\) terms for different choices of \(\alpha \) and \(\beta \) values. In the second approach, we introduce a polynomial-based pseudo-spectral Gauss Runge–Kutta (GRK) collocation formulation for \(Y'(t)=A Y(t)\) with initial \(Y(0)\) being either a vector containing the initial conditions or an Identity matrix. The spectral deferred correction method is then applied to efficiently solve the resulting GRK formulation in each time marching step. For both approaches, matrix–vector or matrix–matrix multiplications are accelerated using the compressible features of T, V , \(e^{\alpha T}\) , and \(e^{\beta V}\) . Analytical and numerical results are presented to demonstrate the accuracy and efficiency properties of the two approaches and to provide guidelines for choosing an appropriate method for a given system. By combining randomized rank-revealing and low-rank decomposition algorithms, we demonstrate how the developed tools can be used to study the interactions of the low-rank structure with diagonal structure when time t evolves in the matrix exponential.