In this work, three versions of self-consistent field/Kohn-Sham density functional theory (SCF/KS-DFT) orbital optimization are described and benchmarked. The methods are a modified version of the geometry version of the direct inversion in the iterative subspace approach (which we call r-GDIIS), the modified restricted step rational function optimization method (RS-RFO), and the novel subspace gradient-enhanced Kriging method combined with restricted variance optimization (S-GEK/RVO). The modifications introduced are aimed at improving the robustness and computational scaling of the procedures. In particular, the subspace approach in S-GEK/RVO allows the application to SCF/KS-DFT optimization of a machine learning technique that has proven to be successful in geometry optimizations. The performance of the three methods is benchmarked for a large number of small- to medium-sized organic molecules, at equilibrium structures and close to a transition state, and a second set of molecules containing closed- and open-shell transition metals. The results indicate the importance of the resetting technique in boosting the performance of the r-GDIIS procedure. Moreover, it is demonstrated that already at the inception of the subspace version of GEK to optimize SCF wave functions, it displays superior and robust convergence properties as compared to those of the standard state-of-the-art SCF/KS-DFT optimization methods.