An Algorithm for Finding an Optimal Projection of a Symmetric Matrix onto a Diagonal Matrix

Abstract

This work is motivated by two two-sided optimization problems that arise in atomic chemistry when one is looking for a minimal set of localized orbitals with prescribed occupation numbers, respectively, a prescribed total number of electrons. We first propose an optimal analytic solution of the first problem, and then show that an optimal solution of the second problem can be found by solving a convex quadratic programming problem with box constraints and $p$ unknowns. We prove that the latter problem is solved by the active-set method in at most $2p$ iterations. These investigations reveal that the optimal solutions of both problems are projections that are in $\mathcal{C}=\{Y\in\mathbb{R}^{n\times p}:Y^TY=I_p, Y^TN Y=\Delta\}$ for $N$ symmetric and $\Delta$ diagonal and that these solutions are generally nonunique. The question of how one can then select a particular solution out of the set $\mathcal{C}$ leads to the main problem of this paper: the optimization of an arbitrary smooth function whose minimizer describes the solution of interest over $\mathcal{C}$. To solve this problem, we first find that a slight modification of $\mathcal{C}$ is a Riemannian manifold for which geometric objects can be derived that are required to allow optimization over this manifold. Using these geometric tools we propose an augmented Lagrangian-based algorithm that guarantees global convergence to a stationary point of our main problem. The latter is shown by investigating when the LICQ is satisfied. Finally we compare this algorithm numerically with a similar algorithm that, however, does not apply these geometric tools and that is, to our knowledge, not guaranteed to converge. Our results show that our algorithm yields a significantly better performance.