10 pages, 6 figures ; Density-functional theory has been one of the most successful approaches ever to address the electronic-structure problem; nevertheless, since its implementations are by necessity approximate, they can suffer from a number of fundamental qualitative shortcomings, often rooted in the remnant electronic self-interaction present in the approximate energy functionals adopted. Functionals that strive to correct for such self-interaction errors, such as those obtained by imposing the Perdew- Zunger self-interaction correction [Phys. Rev. B 23, 5048 (1981)] or the generalized Koopmans' condition [Phys. Rev. B 82, 115121 (2010)], become orbital dependent or orbital-density depen- dent, and provide a very promising avenue to go beyond density-functional theory, especially when studying electronic, optical and dielectric properties, charge-transfer excitations, and molecular dis- sociations. Unlike conventional density functionals, these functionals are not invariant under unitary transformations of occupied electronic states, which leave the total charge density intact, and this added complexity has greatly inhibited both their development and their practical applicability. Here, we first recast the minimization problem for non-unitary invariant energy functionals into the language of ensemble density-functional theory [Phys. Rev. Lett. 79, 1337 (1997)], decoupling the variational search into an inner loop of unitary transformations that minimize the energy at fixed orbital subspace, and an outer-loop evolution of the orbitals in the space orthogonal to the occupied manifold. Then, we show that the potential energy surface in the inner loop is far from convex parabolic in the early stages of the minimization and hence minimization schemes based on these assumptions are unstable, and present an approach to overcome such difficulty. The overall for- mulation allows for a stable, robust, and efficient variational minimization of non-unitary-invariant functionals, essential to study complex materials and molecules, and to investigate the bulk thermo- dynamic limit, where orbitals converge typically to localized Wannier functions. In particular, using maximally localized Wannier functions as an initial guess can greatly reduce the computational costs needed to reach the energy minimum while not affecting or improving the convergence efficiency.