Resting-state functional magnetic resonance imaging (rs-fMRI) provides a valuable tool to study spontaneous brain activity. Using rs-fMRI, researchers have extensively studied the organization of the brain functional network and found several consistent communities consisting of functionally connected but spatially separated brain regions across subjects. However, increasing evidence in many disciplines has shown that most realistic complex networks have overlapping community structure. Only recently has the overlapping community structure drawn increasing interest in the domain of brain network studies. Another issue is that the inter-subject variability is often not directly reflected in the process of community detection at the group level. In this paper, we propose a novel method called collective sparse symmetric non-negative matrix factorization (cssNMF) to address these issues. The cssNMF approach identifies the group-level overlapping communities across subjects and in the meantime preserves the information of individual variation in brain functional network organization. To comprehensively validate cssNMF, a simulated fMRI dataset with ground-truth, a real rs-fMRI dataset with two repeated sessions and another different real rs-fMRI dataset have been used for performance comparison in the experiment. Experimental results show that the proposed cssNMF method accurately and stably identifies group-level overlapping communities across subjects as well as individual differences in network organization with neurophysiologically meaningful interpretations. This research extends our understanding of the common underlying community structures and individual differences in community strengths in brain functional network organization.