<p>Covariance operators form the foundation of functional principal component analysis (FPCA) and related methods in functional data analysis. For large-scale or distributed datasets, however, empirical covariance operator estimation becomes computationally infeasible. We introduce a Poisson subsampling framework for covariance operator estimation that combines unbiasedness, statistical efficiency, and computational scalability. By minimizing Hilbert–Schmidt error, we derive optimal inclusion probabilities and establish non-asymptotic operator-norm concentration bounds. Extensions to FPCA provide guarantees for eigenvalues, eigenfunctions, and projection operators. To address decentralized environments, we also propose a gossip-based federated algorithm that aggregates locally subsampled estimates while respecting communication constraints. Simulation studies demonstrate robustness to heavy-tailed data and improved estimation of covariance operator and non-dominant eigenfunctions. Applications to stellar spectra and financial time series highlight the method’s effectiveness in large-scale functional data analysis.</p>