Non-line-of-sight (NLoS) imaging reveals a hidden scene using indirect diffuse reflections. A common choice for analyzing the time-of-flight (ToF) data from a non-confocal system is an ellipsoid model whose operator is high-dimensional, leading to a computationally arduous task. In this Letter, the product-convolution expansions method is utilized to formulate the operator and its adjoint based on the observation of a shift-variant point spread function (PSF) in the ToF data. The operator and its adjoint are locally approximated as a convolution, which allows the forward and backward procedure to be computed efficiently through fast Fourier transform (FFT). Moreover, the low-rank approximation of the operator is obtained by matrix decompositions, further improving the computational efficiency. The proposed method is validated using publicly accessible datasets.