We present a multilevel algorithm for the solution of a source identification problem in which the forward problem is an elliptic partial differential equation on the 2D unit box. The Hessian corresponds to a Tikhonov-regularized first-kind Fredholm equation. Our method uses an approximate Hessian operator for which, first, the spectral decomposition is known, and second, there exists a fast algorithm that can perform the spectral transform. Based on this decomposition we propose a conjugate gradients solver which we precondition with a multilevel subspace projection scheme. The coarse-level preconditioner is an exact solve and the finer-levels preconditioner is one step of the scaled Richardson iteration. As a model problem, we consider the 2D-Neumann Poisson problem with variable coefficients and partial observations. The approximate Hessian for this case is the Hessian related to a problem with constant coefficients and full observations. We can use a fast cosine transform to compute the spectral transforms. We examine the effect of using Galerkin or level-discretized Hessian operators and we provide results from numerical experiments that indicate the effectiveness of the method for full and partial observations.