Charge and energy transfer in biological and synthetic organic materials are strongly influenced by the coupling of electronic states to high-frequency underdamped vibrations under dephasing noise. Non-perturbative simulations of these systems require a substantial computational effort and current methods can only be applied to large systems with severely coarse-grained environmental structures. In this work, we introduce a dissipation-assisted matrix product factorization (DAMPF) method based on a memory-efficient matrix product operator (MPO) representation of the vibronic state at finite temperature. In this approach, the correlations between environmental vibrational modes can be controlled by the MPO bond dimension, allowing for systematic interpolation between approximate and numerically exact dynamics. Crucially, by subjecting the vibrational modes to damping, we show that one can significantly reduce the bond dimension required to achieve a desired accuracy, and also consider a continuous, highly structured spectral density in a non-perturbative manner. We demonstrate that our method can simulate large vibronic systems consisting of 10-50 sites coupled with 100-1000 underdamped modes in total and for a wide range of parameter regimes. An analytical error bound is provided which allows one to monitor the accuracy of the numerical results. This formalism will facilitate the investigation of spatially extended systems with applications to quantum biology, organic photovoltaics and quantum thermodynamics. arXiv:1903.05443v2 [quant-ph]