Symmetric Nonnegative Matrix Factorization (symNMF) is a special case of the standard Nonnegative Matrix Factorization (NMF) method which is the most popular linear dimensionality reduction technique for analyzing nonnegative data. Examples of symmetric matrices that arise in real-life applications include covariance matrices in finance, adjacency matrices associated with undirected graphs, and distance matrices arising from video and other media summarization technologies. The advantages of having sparse factors in symNMF include saving a great deal of storage space and enhancing the extraction of localized basis features that better represent the latent features of the original data. In order to obtain sparser basis factors, we formulate a new sparse symNMF model by imposing an l1-norm based sparsity constraint on the symNMF problem. We use the concept of rank-one nonnegative matrix approximation to decouple this non-convex optimization problem into convex nonnegative least squares sub-problems which are easier to solve. We develop a new and efficient coordinate-descent based algorithm and test it on four well-known databases of facial images. The proposed algorithm is shown to have fast convergence, well-suited for solving the sparse symNMF problem, and extract localized basis features of image datasets that facilitate interpretation.