We present an approach to the DFT+U method (Density Functional Theory + Hubbard model) within which the computational effort for calculation of ground state energies and forces scales linearly with system size. We employ a formulation of the Hubbard model using nonorthogonal projector functions to define the localized subspaces, and apply it to a local-orbital DFT method including in situ orbital optimization. The resulting approach thus combines linear-scaling and systematic variational convergence. We demonstrate the scaling of the method by applying it to nickel oxide nano-clusters with sizes exceeding 7,000 atoms.