We consider the problem of image-domain least-squares migration (LSM) based on efficiently constructing the Hessian matrix with sparse beam data. Specifically, we use the ultra-wide-band phase space beam summation method, in which beams are used as local basis functions to represent scattered data collected at the surface. The beam domain data are sparse. One can identify seismic events with significant contributions so that only beams with nonnegligible amplitudes need to be used to image the subsurface. In addition, due to the beams’ spectral localization, only beams that pass near an imaging point need to be taken into account. These two properties reduce the computational complexity of computing the Hessian matrix — an essential ingredient for LSM. As a result, we can efficiently construct the Hessian matrix based on analyzing the sparse beam domain data.