A numerically stable and efficient generalized least squares (GLS) algorithm is presented, based on the minimization of an alternative functional. The method relies on orthogonal rotations, is highly parallel and equivalent to a deterministic algorithm that does not `square' the covariance square root to compute the Kalman gain. Computation of eigenvalue and singular value decompositions is not required. We present numerical results for the Lorenz (2005, J. Atmos. Sci.) two-scale dynamics, using perfect and imperfect models. Compared to the Ensemble Transform Kalman Filter (ETKF), Monte-Carlo filters, and other square-root formulations, the GLS-EnKF algorithm results in lower mean square errors and is less susceptible to filter divergence. Although the vector form of the GLS algorithm and the two-step serial Ensemble Adjustment Kalman Filter (EAKF) described by Anderson (2003, Monthly Weather Review) are fundamentally different, they both lead to the same error levels.