In electronic structure methods based on the correction of approximate density-functional theory (DFT) for systematic inaccuracies, Hubbard $U$ parameters may be used to quantify and amend the self-interaction errors ascribed to selected subspaces. Here, in order to enable the accurate, computationally convenient calculation of $U$ by means of DFT algorithms that locate the ground-state by direct total-energy minimization, we introduce a reformulation of the successful linear-response method for $U$ in terms of the fully-relaxed constrained ground-state density. Defining $U$ as an implicit functional of the ground-state density implies the comparability of DFT + Hubbard $U$ (DFT+$U$) total-energies, and related properties, as external parameters such as ionic positions are varied together with their corresponding first-principles $U$ values. Our approach provides a framework in which to address the partially unresolved question of self-consistency over $U$, for which plausible schemes have been proposed, and to precisely define the energy associated with subspace many-body self-interaction error. We demonstrate that DFT+$U$ precisely corrects the total energy for self-interaction error under ideal conditions, but only if a simple self-consistency condition is applied. Such parameters also promote to first-principles a recently proposed DFT+$U$ based method for enforcing Koopmans theorem.