Described here are systems and methods for quantitative susceptibility mapping (“QSM”) using magnetic resonance imaging (“MRI”). Susceptibility maps are reconstructed from phase images using an automatic regularization technique based in part on variable splitting. Two different regularization parameters are used, one, λ, that controls the smoothness of the final susceptibility map and one, μ, that controls the convergence speed of the reconstruction. For instance, the regularization parameters can be determined using an L-curve heuristic to find the parameters that yield the maximum curvature on the L-curve. The μ parameter can be determined based on an l2-regularization and the λ parameter can be determined based on the iterative l1-regularization used to reconstruct the susceptibility map.