We develop new approximations of the Mori-Zwanzig equation based on operator series expansions of the orthogonal dynamics propagator. In particular, we study the Faber polynomial series, which yields asymptotically optimal approximations converging at least R-superlinearly with the polynomial order. Numerical applications are presented for random wave propagation and harmonic chains of oscillators interacting on the Bethe lattice and on graphs with arbitrary topology. Most of the analysis and numerical results presented in this paper are for linear dynamical systems evolving from random initial states. The extension to nonlinear systems is discussed as well the challenges in the implementation of the proposed approximation methods.