  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.
