Estimating functional brain networks from fMRI data has been the focus of much research in recent years. Low sample sizes (time-points) and high dimensionality of fMRI has restricted estimation to a temporally averaged connectivity matrix per subject, due to which the dynamics of functional connectivity is largely unknown. In this paper, we propose a novel method based on constrained matrix factorization that addresses two major issues. Firstly, it finds a set of basis networks that are the semantic parts of the time-varying whole-brain functional networks. The whole-brain network at any point in time, for any subject, is a non-negative combination of these basis networks. Secondly, significant dimensionality reduction is achieved by projecting the data onto this basis, facilitating subsequent analysis of temporal dynamics. Results on simulated fMRI data show that our method can effectively recover underlying basis networks. We apply this method on a normative dataset of resting state fMRI scans. Results indicate that the functional connectivity of a subject at any point during the scan is composed of combinations of overlapping task-positive/negative pairs of sub-networks.