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 we 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 nanoclusters with sizes exceeding 7000 atoms.