The theory for skeletonized least squares wave equation migration (LSM) is presented. The key idea is, for an assumed velocity model, the source-side Green's function and the geophone-side Green's function are computed by a numerical solution of the wave equation. Only the early-arrivals of these Green's functions are saved and skeletonized to form the migration Green's function (MGF) by convolution. Then the migration image is obtained by a dot product between the recorded shot gathers and the MGF for every trial image point. The key to an efficient implementation of iterative LSM is that at each conjugate gradient iteration, the MGF is reused and no new finitedifference (FD) simulations are needed to get the updated migration image. It is believed that this procedure combined with phase-encoded multi-source technology will allow for the efficient computation of wave equation LSM images in less time than that of conventional reverse time migration (RTM).