Nonlinear multi-scale problems are ubiquitous in materials science and biology. Complicated interactions between nonlinearities and (nonseparable) multiple scales pose a major challenge for analysis and simulation. In this paper, we study the numerical homogenization for multi-scale elliptic PDEs with monotone nonlinearity, in particular the Leray-Lions problem (a prototypical example is the p-Laplacian equation), where the nonlinearity cannot be parameterized with low dimensional parameters, and the linearization error is non-negligible. We develop the iterated numerical homogenization scheme by combining numerical homogenization methods for linear equations, and the so-called quasi-norm based iterative approach for monotone nonlinear equation. We propose a residual regularized nonlinear iterative method, and in addition, develop the sparse updating method for the efficient update of coarse spaces. A number of numerical results are presented to complement the analysis and valid the numerical method.