Nonnegative matrix factorization (NMF) is a prominent technique for data dimensionality reduction that has been widely used for text mining, computer vision, pattern discovery, and bioinformatics. In this paper, a framework called ARkNLS (Alternating Rank-k Nonnegativity constrained Least Squares) is proposed for computing NMF. First, a recursive formula for the solution of the rank-k nonnegativity-constrained least squares (NLS) is established. This recursive formula can be used to derive the closed-form solution for the Rank-k NLS problem for any integer k ≥ 1. As a result, each subproblem for an alternating rank-k nonnegative least squares framework can be obtained based on this closed form solution. Assuming that all matrices involved in rank-k NLS in the context of NMF computation are of full rank, two of the currently best NMF algorithms HALS (hierarchical alternating least squares) and ANLS-BPP (Alternating NLS based on Block Principal Pivoting) can be considered as special cases of ARkNLS with k = 1 and k = r for rank r NMF, respectively. This paper is then focused on the framework with k = 3, which leads to a new algorithm for NMF via the closed-form solution of the rank-3 NLS problem. Furthermore, a new strategy that efficiently overcomes the potential singularity problem in rank-3 NLS within the context of NMF computation is also presented. Extensive numerical comparisons using real and synthetic data sets demonstrate that the proposed algorithm provides state-of-the-art performance in terms of computational accuracy and cpu time.