This study presents a new method of centroid moment tensor (CMT) data inversion to estimate time-dependent regional stress fields. A Gaussian process (GP) is applied to resolve a computational difficulty of the existing basis function expansion method in analyzing high-dimensional data including time dependence. A critical step in the formulation is an analytical derivation of the relationship of the covariance function, which is a key ingredient of GP, between CMT data and a model stress field based on an observation equation. Applications to CMT data in and around Japan after the 2011 Tohoku earthquake show the efficiency and validity of the method, which clarifies that the stress field has small-scale heterogeneity in space and long-term stability in time for most regions. Additionally, significant temporal variations are observed around the margin of the focal region of the 2011 event, the sense of which is opposite in the landward side and the oceanward side. GP would be particularly effective in geophysical inversions of high-dimensional data distributed in a broad region.