The problem of estimating a low-dimensional subspace from a collection of experimentally measured subspaces arises in many applications of statistical signal processing. In this paper we address this problem, and give a solution for the average subspace that minimizes an extrinsic mean-squared error, defined by the squared Frobenius norm between projection matrices. The solution automatically returns the dimension of the optimal average subspace, which is the novel result of the paper. The proposed order fitting rule is based on thresholding the eigenvalues of the average projection matrix, and thus it is free of penalty terms or other tuning parameters commonly used by other rank estimation techniques. Several numerical examples demonstrate the usefulness and applicability of the proposed criterion, showing how the dimension of the average subspace captures the variability of the measured subspaces.