When applying the Bayesian manifold regularization method to function estimation problem with manifold constraints, the direct implementation has computational complexity O(N3), where N is the number of input-output data measurements. This becomes particularly costly when N is large. In this paper, we propose a more efficient implementation based on the Kalman filter and smoother using a state-space model realization of the underlying Gaussian process. Moreover, we explore the sequentially semi-separable structure of the Laplacian matrix and the posterior covariance matrix. Our proposed implementation has computational complexity O(N) and thus can be applied to large data problems. We exemplify the effectiveness of our proposed implementation through numerical simulations.