Not yet reviewed by Pith; the record is open.

Our research objective in this paper is to reconstruct an initial linear density field, which follows the multivariate Gaussian distribution with variances given by the linear power spectrum of the current CDM model and evolves through gravitational instability to the present-day density field in the local Universe. For this purpose, we develop a Hamiltonian Markov Chain Monte Carlo method to obtain the linear density field from a posterior probability function that consists of two components: a prior of a Gaussian density field with a given linear spectrum, and a likelihood term that is given by the current density field. The present-day density field can be reconstructed from galaxy groups using the method developed in Wang et al. (2009a). Using a realistic mock SDSS DR7, obtained by populating dark matter haloes in the Millennium simulation with galaxies, we show that our method can effectively and accurately recover both the amplitudes and phases of the initial, linear density field. To examine the accuracy of our method, we use $N$-body simulations to evolve these reconstructed initial conditions to the present day. The resimulated density field thus obtained accurately matches the original density field of the Millennium simulation in the density range 0.3 <= rho/rho_mean <= 20 without any significant bias. Especially, the Fourier phases of the resimulated density fields are tightly correlated with those of the original simulation down to a scale corresponding to a wavenumber of ~ 1 h/Mpc, much smaller than the translinear scale, which corresponds to a wavenumber of ~ 0.15 h\Mpc.