This is Part II of a two-part work on the estimation for a multi-layer generalized linear model (ML-GLM) in large system limits. In Part I, we had analyzed the asymptotic performance of an exact MMSE estimator, and obtained a set of coupled equations that could characterize its MSE performance. To work around the implementation difficulty of the exact estimator, this paper continues to propose an approximate solution, ML-GAMP, which could be derived by blending a moment-matching projection into the Gaussian approximated loopy belief propagation. The ML-GAMP estimator is then shown to enjoy a great simplicity in its implementation, where its per-iteration complexity is as low as GAMP. Further analysis on its asymptotic performance also reveals that, in large system limits, its dynamical MSE behavior is fully characterized by a set of simple one-dimensional iterating equations, termed state evolution (SE). Interestingly, this SE of ML-GAMP share exactly the same fixed points with an exact MMSE estimator whose fixed points were obtained in Part I via a replica analysis. Given the Bayes-optimality of the exact implementation, this proposed estimator (if converged) is optimal in the MSE sense.