Kohn-Sham density-functional calculations are used in many branches of science to obtain information about the electronic structure of molecular systems and materials. Conventional algorithms for minimization of the Kohn-Sham energy have certain deficiencies, however, that may cause divergence or, worse, convergence to unphysical saddle points. We here present a three-level hierarchical minimization strategy which is both more efficient and robust than the conventional algorithms and which does not suffer from the flaws of these algorithms. Using the three-level minimization strategy, the molecular density is built up in a hierarchical fashion in accordance with chemical insight: First, the molecular density is composed by a superposition of atomic densities; next, bonds are formed by performing a simple valence-shell optimization; finally, the molecular description is refined by an optimization in the full molecular basis. Importantly, the density matrix generated at each of the two lower levels in this hierarchy is transferred to the next without loss of information. Examples demonstrate the efficacy and robustness of the proposed scheme.