  We leverage powerful mathematical tools stemming from optimal transport
theory and transform them into an efficient algorithm to reconstruct the
fluctuations of the primordial density field, built on solving the
Monge-Amp\`ere-Kantorovich equation. Our algorithm computes the optimal
transport between an initial uniform continuous density field, partitioned into
Laguerre cells, and a final input set of discrete point masses, linking the
early to the late Universe. While existing early universe reconstruction
algorithms based on fully discrete combinatorial methods are limited to a few
hundred thousand points, our algorithm scales up well beyond this limit, since
it takes the form of a well-posed smooth convex optimization problem, solved
using a Newton method. We run our algorithm on cosmological $N$-body
simulations, from the AbacusCosmos suite, and reconstruct the initial positions
of $\mathcal{O}(10^7)$ particles within a few hours with an off-the-shelf
personal computer. We show that our method allows a unique, fast and precise
recovery of subtle features of the initial power spectrum, such as the baryonic
acoustic oscillations.
