Hyperspectral unmixing is an important technique for remote sensing image analysis. Among various unmixing techniques, nonnegative matrix factorization (NMF) shows unique advantage in providing a unified solution with well physical interpretation. In order to explore the geometric information of the hyperspectral data, graph regularization is often used to improve the NMF unmixing performance. It groups neighboring pixels, uses groups as graph vertices, and then assigns weights to connected vertices. The construction of neighborhood and the weights are normally determined by k-nearest neighbors or heat kernel in a deterministic process, which do not fully reveal the structural relationships among data. In this article, we introduce an adaptive graph to regularize a multilayer NMF (AGMLNMF) model for hyperspectral unmixing. In AGMLNMF, a graph is constructed based on the probabilities between neighbors. This enables the optimal neighborhood be automatically determined. Moreover, the weights of the graph are assigned based on the relationships among neighbors, which reflects the intrinsic structure of the complex data. Experiments on both synthetic and real datasets show that this method has outperformed several state-of-the-art unmixing approaches.