  A theoretically interesting and practically important question in cosmology
is the reconstruction of the initial density distribution provided a late-time
density field. This is a long-standing question with a revived interest
recently, especially in the context of optimally extracting the baryonic
acoustic oscillation (BAO) signals from observed galaxy distributions. We
present a new efficient method to carry out this reconstruction, which is based
on numerical solutions to the nonlinear partial differential equation that
governs the mapping between the initial Lagrangian and final Eulerian
coordinates of particles in evolved density fields. This is motivated by
numerical simulations of the quartic Galileon gravity model, which has similar
equations that can be solved effectively by multigrid Gauss-Seidel relaxation.
The method is based on mass conservation, and does not assume any specific
cosmological model. Our test shows that it has a performance comparable to that
of state-of-the-art algorithms which were very recently put forward in the
literature, with the reconstructed density field over $\sim80\%$ ($50\%$)
correlated with the initial condition at $k\lesssim0.6h/{\rm Mpc}$ ($1.0h/{\rm
Mpc}$). With an example, we demonstrate that this method can significantly
improve the accuracy of BAO reconstruction.
