Hyperspectral images contain mixed pixels due to low spatial resolution of hyperspectral sensors. Spectral unmixing problem refers to decomposing mixed pixels into a set of endmembers and abundance fractions. Due to nonnegativity constraint on abundance fractions, nonnegative matrix factorization (NMF) methods have been widely used for solving spectral unmixing problem. In this letter we proposed using multilayer NMF (MLNMF) for the purpose of hyperspectral unmixing. In this approach, spectral signature matrix can be modeled as a product of sparse matrices. In fact MLNMF decomposes the observation matrix iteratively in a number of layers. In each layer, we applied sparseness constraint on spectral signature matrix as well as on abundance fractions matrix. In this way signatures matrix can be sparsely decomposed despite the fact that it is not generally a sparse matrix. The proposed algorithm is applied on synthetic and real datasets. Synthetic data is generated based on endmembers from USGS spectral library. AVIRIS Cuprite dataset has been used as a real dataset for evaluation of proposed method. Results of experiments are quantified based on SAD and AAD measures. Results in comparison with previously proposed methods show that the multilayer approach can unmix data more effectively.