Nonnegative matrix factorization (NMF) is an advanced method for nonnegative feature extraction, with widespread applications. However, the NMF solution often entails to solve a global optimization problem with a nonconvex objective function and nonnegativity constraints. This paper presents a collective neurodynamic optimization (CNO) approach to this challenging problem. The proposed collective neurodynamic system consists of a population of recurrent neural networks (RNNs) at the lower level and a particle swarm optimization (PSO) algorithm with wavelet mutation at the upper level. The RNNs act as search agents carrying out precise local searches according to their neurodynamics and initial conditions. The PSO algorithm coordinates and guides the RNNs with updated initial states toward global optimal solution(s). A wavelet mutation operator is added to enhance PSO exploration diversity. Through iterative interaction and improvement of the locally best solutions of RNNs and global best positions of the whole population, the population-based neurodynamic systems are almost sure able to achieve the global optimality for the NMF problem. It is proved that the convergence of the group-best state to the global optimal solution with probability one. The experimental results substantiate the efficacy and superiority of the CNO approach to bound-constrained global optimization with several benchmark nonconvex functions and NMF-based clustering with benchmark data sets in comparison with the state-of-the-art algorithms.