We consider the inverse medium problem for the time-harmonic wave equation with broadband and multi-point illumination in the low frequency regime. Such a problem finds many applications in geosciences (e.g. ground penetrating radar), non-destructive evaluation (acoustics), and medicine (optical tomography). We use an integral-equation (Lippmann-Schwinger) formulation, which we discretize using a quadrature method. We consider only small perturbations of the background medium(Born approximation). To solve this inverse problem, we use a least squares formulation that is regularized with the truncated Singular Value Decomposition (SVD). If Nfr is the number of excitation frequencies, Ns the number of incoming waves, Nd the number of detectors, and N the parameterization for the scatterer, a dense singular value decomposition for the overall input-output map will have [min(Ns Nfr Nd, N)]2 \times max(Ns Nfr Nd, N) cost. We have developed a fast SVD approach that brings the cost down to O( N Nfr Nd +N Nfr Ns) thus, providing orders of magnitude improvements over a black-box dense SVD. We provide numerical results that demonstrate the scalability of the method.