Abstract: To accurately evaluate the probabilistic characteristics of seismic responses of structures, the temporal and spatial variability and the spectral compatibility need to be considered when simulating stochastic processes of ground motions. In classical simulations of multivariate fully non-stationary processes, the decomposition efficiency of the evolutionary power spectral density (EPSD) matrix is low due to the inseparability of frequency and temporal variables. To speed up the decomposition of the EPSD matrix, a novel Cholesky decomposition approach is proposed to simulate the non-stationary processes. The key of this approach is to separate the EPSD matrix into aphase and a modulus matrix. The modulus matrix will be further transformed to a lagged coherency matrix. The lagged coherence matrix is only related to the frequency, as opposed to the time-dependent EPSD matrix, which remarkably enhances the Cholesky decomposition efficiency. Furthermore, the lagged coherency matrix is more compatible with interpolation techniques. Finally, the novel Cholesky decomposition method and interpolation techniques are used in a stochastic method which is capable of generating spectrum-compatible ground motion samples. Results show that the novel Cholesky decomposition and interpolation techniques are valid for generating fully non-stationary and spectrum-compatible multivariate ground motion samples. Both the linear and cubic spline interpolations achieve a satisfactory level of resolution and accuracy, even with a small number of interpolation points.