Despite recent progress in linear scaling (LS) density function theory (DFT), the computational cost of the existing LS methods remains too high for a widespread adoption at present. In this work, we exploit nonorthogonal localized molecular orbitals to develop a series of LS methods for molecular systems with a low computational overhead. High efficiency of the proposed methods is achieved with a new robust two-stage variational procedure or by replacing the optimization altogether with an accurate nonself-consistent approach. We demonstrate that, even for challenging condensed-phase systems, the implemented LS methods are capable of extending the range of accurate DFT simulations to molecular systems that are an order of magnitude larger than those previously treated.