Initial density distribution provides a basis for understanding the complete evolution of cosmological density fluctuations. While reconstruction in our local Universe exploits the observations of galaxy surveys with large volumes, observations of high-redshift galaxies are performed with a small field of view and therefore can hardly be used for reconstruction. Here we propose to reconstruct the initial density field using the H I 21 cm and CO line intensity maps from the epoch of reionization. Observations of these two intensity maps provide complementary information of the density field -- the H I 21 cm field is a proxy of matter distributions in the neutral regions, while the CO line intensity maps are sensitive to the high-density, star-forming regions that host the sources for reionization. Technically, we employ the conjugate gradient method and develop the machinery for minimizing the cost function for the intensity mapping observations. Analytical expressions for the gradient of cost function are derived explicitly. We show that the resimulated intensity maps match the input maps of mock observations using semi-numerical simulations of reionization with an rms error ≲7% at all stages of reionization. This reconstruction is also robust at the same level of accuracy when a noise at the level of ≲1% of the standard deviation is applied to each map. Our proof-of-concept work demonstrates the robustness of the reconstruction method, thereby providing an effective technique for reconstructing the cosmological initial density distribution from high-redshift observations.
