We consider the computational problem arising in magnetoencephalography (MEG), where the goal is to estimate the electric activity within the brain noninvasively from extracranial measurements of the magnetic field components. The problem is severely ill-posed due to the intrinsic nonuniqueness of the solution, and suffers further from the challenges of a weak data signal, its high dimensionality, and complexity of the noise, part of which is due to the brain itself. In this work, we suggest a new algorithm that is based on a truncated conjugate gradient algorithm for least squares with statistically inspired left and right preconditioners. We demonstrate that by carefully accounting for the spatiotemporal statistical structure of the brain noise and by adopting a suitable prior within the Bayesian framework, we can design a robust and efficient method for the numerical solution of the MEG inverse problem which can improve the spatial and temporal resolution of events of short duration. The effectiveness of the proposed method is demonstrated on a synthetic example of localization of spiking simulating the focal onset of epileptic seizures.