  We present a method to reconstruct the initial conditions of the universe
using observed galaxy positions and luminosities under the assumption that the
luminosities can be calibrated with weak lensing to give the mean halo mass.
Our method relies on following the gradients of forward model and since the
standard way to identify halos is non-differentiable and results in a discrete
sample of objects, we propose a framework to model the halo position and mass
field starting from the non-linear matter field using Neural Networks. We
evaluate the performance of our model with multiple metrics. Our model is more
than $95\%$ correlated with the halo-mass fields up to $k\sim 0.7 {\rm h/Mpc}$
and significantly reduces the stochasticity over the Poisson shot noise. We
develop a data likelihood model that takes our modeling error and intrinsic
scatter in the halo mass-light relation into account and show that a displaced
log-normal model is a good approximation to it. We optimize over the
corresponding loss function to reconstruct the initial density field and
develop an annealing procedure to speed up and improve the convergence. We
apply the method to halo number densities of $\bar{n} = 2.5\times 10^{-4}
-10^{-3}({\rm h/Mpc})^3$, typical of current and future redshift surveys, and
recover a Gaussian initial density field, mapping all the higher order
information in the data into the power spectrum. We show that our
reconstruction improves over the standard reconstruction. For baryonic acoustic
oscillations (BAO) the gains are relatively modest because BAO is dominated by
large scales where standard reconstruction suffices. We improve upon it by
$\sim 15-20\%$ in terms of error on BAO peak as estimated by Fisher analysis at
$z=0$. We expect larger gains will be achieved when applying this method to the
broadband linear power spectrum reconstruction on smaller scales.
