A great interest has been given to the Nonnegative Matrix Factorization (NMF) due to its ability of extracting highly-interpretable parts from data sets. Nonetheless, its usage is hindered by the computational complexity when processing large matrices. In this paper, we present three implementations of NMF to analyze the benefits that parallelism can provide to this method on data sets of different size. Our first implementation is based on Message Passing Interface (MPI). Input data is distributed among multiple processors. The second version uses CUDA in Graphics Processing Units (GPUs). Large data sets are blockwise transferred and processed. Finally, we combine both paradigms. Compared to the single-GPU version, it achieves super linear speedups when data portions assigned to each GPU fit into their memory and can be transferred only once.