Abstract: Nonnegative Matrix Factorization (NMF) is a recently developed technique for nonlinearly finding purely additive, parts-based, linear, and low-dimension representations of nonnegative multivariate data to consequently reveal the latent structure, feature or pattern in the data. Although NMF has been successfully applied to several research fields, it is confronted with two main problems (unsatisfactory accuracy and bad generality) while the processed is a matrix-set, because the processed objects of NMF are intrinsically vectors and the necessary vectorization for every matrix in the processed matrix-set often make corresponding NMF learning to be typical small-sample learning. In this paper, Nonnegative Matrix-Set Factorization (NMSF) is conceived to overcome the problems and to retain NMFs good properties. As opposed to NMF, NMSF directly processes original data matrices rather than vectorization results of them. Theoretical analysis shows that while processing a data matrix-set, NMSF should be more accurate and has better generality than NMF. To show how to implement NMSF, and to validate NMSFs properties by experiments, Bilinear Form-Based NMSF (BFBNMSF) algorithm, as an implementation mode of NMSF, is formulated. Results of comparison experiments between BFBNMSF and NMF stably support the theoretical analysis. It is worth noting that higher accuracy and better generality actually means that NMSF is better at extracting essential features of data matrices than NMF.