We establish a scalable manifold learning method and theory, motivated by the problem of estimating fMRI activation manifolds in the Human Connectome Project (HCP). Our primary contribution is the development of an efficient estimation technique for heat kernel Gaussian processes in the exponential family model. This approach handles large sample sizes $n$, preserves the intrinsic geometry of data, and significantly reduces computational complexity from $\mathcal{O}(n^3)$ to $\mathcal{O}(n)$ via a novel reduced-rank approximation of the graph Laplacian's transition matrix and a Truncated Singular Value Decomposition for the eigenpair computation. The numerical experiments demonstrate the scalability and improved accuracy of our method for manifold learning tasks involving complex large-scale data.