An approach to handling colored observation noise in large least-squares (LS) problems is presented. The handling of colored noise is reduced to the problem of solving a Toeplitz system of linear equations. The colored noise is represented as an auto regressive moving-average (ARMA) process. Stability and invertability of the ARMA model allows the solution of the Toeplitz system to be reduced to two successive filtering operations using the inverse transfer function of the ARMA model. The numerical complexity of the algorithm scales proportionally to the order of the ARMA model times the number of observations. This makes the algorithm particularly suited for LS problems with millions of observations. It can be used in combination with direct and iterative algorithms for the solution of the normal equations. The performance of the algorithm is demonstrated for the computation of a model of the Earth’s gravity field from simulated satellite-derived gravity gradients up to spherical harmonic degree 300.