How to impose orthonormality constraints by method of Lagrange multipliers I want to find the matrix $\Phi: \mathbb{R}^n \to \mathbb{R}^m$, $m<n$ that minimizes
$$V={\rm tr}(\Phi R \Phi^T)$$
subject to the orthonormality constraint
$$\Phi\Phi^T=I$$
where $R: \mathbb{R}^n \to \mathbb{R}^n$ is a given symmetric positive definite matrix.
How do I apply the Lagrange multiplier method to this constrained optimization problem?
I tried:
$$\tilde{V}={\rm tr}(\Phi R \Phi^T) - {\rm tr}\left(\Lambda(\Phi\Phi^T-I)\right)$$
where $\Lambda:\mathbb{R}^m \to \mathbb{R}^m$ is a general matrix of Lagrange multipliers, and this yields the necessary conditions for minimality:
$$\frac{\partial \tilde{V}}{\partial \Lambda}=\Phi\Phi^T-I=0$$
$$\frac{\partial \tilde{V}}{\partial \Phi}=R\Phi^T-\Phi^T\Lambda=0$$
The first equation is the orthonormality condition, so far so good. But how is the second equation supposed to find minimum candidates/critical points? Intuitively, I think I know that the global minimum argument $\Phi$ must be the basis of the eigenspace of the lowest $m$ eigenvalues of $R$. But $\Lambda$ is a general $m\times m$ matrix and nothing seems to constrain it to a diagonal matrix with $m$ eigenvalues (let alone the lowest) on the diagonal.
Or am I applying the Lagrange multiplier method the wrong way?
By the way: it seems easy for $m=1$. Then $\Phi^T$ simply becomes a vector and $\Lambda=\lambda$ becomes a scalar. This yields the eigenvalue problem for one eigenvector, which is 'diagonal' in the trivial sense:
$$R\Phi^T-\lambda\Phi^T=0$$
 A: The Lagrange multiplier method does put a special constraint on the structure of $\Lambda$, but that isn't what you expected.
From $\Phi\Phi^T=I$ and $R\Phi^T=\Phi^T\Lambda$, we obtain $\Lambda=\Phi\Phi^T\Lambda=\Phi R\Phi^T$. The constraint placed on the structure of $\Lambda$ is not that $\Lambda$ must be diagonal, but that $\Lambda$ must be symmetric.
Since $\Lambda$ is symmetric, it can be orthogonally diagonalised as $QDQ^T$. Therefore, $R\Phi^T=\Phi^T\Lambda$ implies that $R(\Phi^TQ)=(\Phi^TQ)D$. The eigenvectors of $R$ are the columns of $\Phi^TQ$ rather than the columns of $\Phi^T$.
This makes sense if you look at the original objective function. Since $\Phi R\Phi^T$ has the same trace as $(Q^T\Phi) R(\Phi^TQ)$ for every $Q\in SO(m)$, there is no reason why the columns of any $\Phi^T$ that minimises $\operatorname{tr}(\Phi R\Phi^T)$ must be eigenvectors of $R$.
A: It is impossible to achieve the diagonalization of $\Lambda$ as you wish to do. 
Simply, the relation $\Phi\Phi^T=1$ imposes in total $\frac{m(m+1)}{2}$ constraints. Therefore one needs to use exactly the same amount of independent Lagrange multipliers to impose all these constraints on the minimization functional. The minimal choice of a multiplier matrix is the one presented above, where $\Lambda$ is a $m\times m$ symmetric matrix. Notice that this matrix contains exactly $\frac{m(m+1)}{2}$ independent Lagrange multipliers. If one was permitted to diagonalize this matrix, that would imply that there would only be $m$ independent constraints to impose, and that the expression for $\Lambda$ would be redundant, and hence subject to transformation to one with fewer degrees of freedom.
$\textbf{Example:}$
Let's attack the simplest case that is non-trivial, $m=2$. Then express $\Phi$ as follows:
$$\Phi=\begin{pmatrix}
\mathbf{q}_1\\\mathbf{q}_2
\end{pmatrix}~~,~~\Lambda=\begin{pmatrix}
\lambda_1&\lambda_3\\ \lambda_3&\lambda_2
\end{pmatrix}$$
where $q$'s are horizontal line vectors 
$$\mathbf{q}_a=(q_{a1}, q_{a2},..., q_{an})$$
With this notation in place, and dropping the boldface notation for vectors henceforth, we can show that the minimization functional reduces for $m=2$ to the following:
$$V=q_1^TRq_1+q_2^TRq_2-\lambda_1(q_1^Tq_1-1)-\lambda_2(q_2^Tq_2-1)-2\lambda_3 q_1^Tq_2$$
and minimizing with respect to the vectors $\mathbf{q}$ we obtain the equations:
$$(R-\lambda_1I)q_1=\lambda_3 q_2\\
(R-\lambda_2I)q_2=\lambda_3 q_1\\$$
Eliminating $q_2$ from the equations we find that $q_1$ satisfies the following nonlinear eigenvalue problem:
$$[R^2-(\lambda_1+\lambda_2)R+(\lambda_1\lambda_2-\lambda_3^2)]q_1=0$$
We get lucky however, because this problem is factorizable into the form
$$(R-r_1(\lambda))(R-r_2(\lambda))q_1=0~\\r_1+r_2=\lambda_1+\lambda_2~,~r_1r_2=\lambda_1\lambda_2-\lambda_3^2$$
which basically an eigenvalue problem. We denote the set of eigenvalues of $R$ as $E_R=\{\omega_1,..., \omega_n\}$ and the set of normalized eigenvectors $V_R=\{e_{1R},...,e_{nR}\}$. It is easy to show that $r_1, r_2 \in E_R$, and therefore the set of solutions to the above equation is discrete: $$(r_1,r_2)=(\omega_i,\omega_j), j\geq i$$
Note that there exist exactly $\frac{m(m+1)}{2}$ different solutions and not $m^2$ due to the symmetry $r_1\to r_2$. Now the question is, can we pin down the values of the Lagrange multipliers using these equations? The answer is positive, at least in principle. 
Choose one of the solutions by fixing $r_1=\omega_i, r_2=\omega_j$ for some, appropriately chosen $i,j$. The general solution to the factorized eigenvalue problem is:
$$q_1=Ae_{iR}+Be_{jR}\\
q_2=\frac{A(\omega_i-\lambda_1)}{\lambda_3}e_{iR}+\frac{B(\omega_j-\lambda_1)}{\lambda_3}e_{jR}$$
where $A,B$ are arbitrary real parameters. We wish to solve for the matrices $\Lambda, \Phi$ but we only have partial information on them so far. It turns out we can write down exactly as many equations as necessary to determine all of their elements. We just need to apply the 3 constraints and the fact that the $\Lambda$ matrix is related to eigenvalues as follows:
$$\begin{align}
&\lambda_1+\lambda_2=\omega_{i}+\omega_{j}\\
&\lambda_1\lambda_2-\lambda^2_3=\omega_i\omega_j
\\& q_1^Tq_1=A^2+B^2=1
 \\&q_2^Tq_2=A^2(\frac{\omega_i-\lambda_1}{\lambda_3})^2+B^2(\frac{\omega_j-\lambda_1}{\lambda_3})^2=1
\\&q_2^Tq_1=A^2(\frac{\omega_i-\lambda_1}{\lambda_3})+B^2(\frac{\omega_j-\lambda_1}{\lambda_3})=0
\end{align}$$ 
These are five equations for five unknowns $A,B,\lambda_1,\lambda_2, \lambda_3$ and therefore they can in principle be solved. 
$\textbf{Generalization:}$
These equations are highly nonlinear for general $m$  and it is not clear for general even whether a solution exists or not, because the vectors are supposed to be real. In a complex vector space where orthonormal matrices are replaced by Hermitian matrices, it's easier to guarantee that the problem has a solution. However, for general values of $m$ the problem can still be reduced to an eigenvalue problem. Consider the Cayley-Hamilton polynomial of the matrix $\Lambda$, defined by $P_{\Lambda}(x)=\det(\Lambda-xI)$. It may be written in the form:
$$P_{\Lambda}(x)=x^m-\text{tr}(\Lambda)x^{m-1}+\frac{\text{tr}(\Lambda^2)-\text{tr}(\Lambda)^2}{2}x^{m-2}+...+(-1)^m \det(\Lambda)$$ 
It may be shown that the vectors $\mathbf{q}_1,...,\mathbf{q}_m$ all satisfy the equation
$$P_{\Lambda}(R)\mathbf{q_i}=0~~,~~ i={1,...,m}$$
which can be subsequently factorized and solved for $q_1$ at least:
$$q_1=\sum_{a=1}^{m-1} A_a e_{i_a R}~,~ i_a\geq i_{a+1}~,~ i_a\in \{1,..n\}$$
The other vectors can be determined from $q_1$. One can also see easily that the factorization provides $m-1$ constraints and the total constraints are $\frac{m(m+3)}{2}-1$. These are exactly enough of them to determine the coefficients of each eigenvector attached to each possible eigenvalue combination and the Lagrange multipliers, as illustrated in the $m=2$ example.  
