# How do I solve the matrix equality constrained optimization problem using Lagrangian multipliers?

Solve the following minimization problem in $$\mathbf{X} \in \mathbb{R}^{m \times n}$$

$$\begin{array}{ll} \text{minimize} & \frac 12 \| \mathbf{X}\mathbf{X}^T -\mathbf{A} \|^2_\mathcal{F}\\ \text{subject to} & \mathbf{X}^T \mathbf{X} = \mathbf{I}_n\end{array}$$

where $$\mathbf{I}_n$$ is the $$n \times n$$ identity matrix, $$\mathbf{A} = \mathbf{B}\mathbf{B}^T$$ where $$\mathbf{B} \in \mathbb{R}^{m\times n}$$.

I was not sure how to write the Lagrangian and differentiate. Most probably, the final solution should be an eigenvalue problem but I am not sure how to arrive at that. Here is my attempt - $$\mathcal{L}(\mathbf{X},\mathbf{Z})=\frac 12 \| \mathbf{X}\mathbf{X}^T -\mathbf{A} \|^2_\mathcal{F} + <\mathbf{Z},\mathbf{X}^T \mathbf{X}-\mathbf{I}>_\mathcal{F}$$ where $$\mathbf{Z}$$ is the Lagrangian multiplier. Now differentiation with respect to $$\mathbf{Z}$$ yields $$\mathbf{X}^T\mathbf{X} = \mathbf{I}$$. Now to take differentiation with respect to $$\mathbf{X}$$ we expand the Lagrangian as follows - $$\mathcal{L}=\frac{1}{2}\text{Tr}((\mathbf{X}\mathbf{X}^T-\mathbf{A})^T(\mathbf{X}\mathbf{X}^T-\mathbf{A}))+ \text{Tr}(\mathbf{Z}^T(\mathbf{X}^T \mathbf{X}-\mathbf{I}))$$ $$\text{or} \quad \mathcal{L}=\frac{1}{2}\text{Tr}(\mathbf{X}\mathbf{X}^T\mathbf{X}\mathbf{X}^T-\mathbf{X}\mathbf{X}^T\mathbf{A}-\mathbf{A}^T\mathbf{X}\mathbf{X}^T+\mathbf{Z}^T\mathbf{X}^T \mathbf{X})+ \mathcal{L}_0$$ where $$\mathcal{L}_0$$ are terms not containing $$\mathbf{X}$$. Then I differentiate $$\mathcal{L}$$ with respect to $$\mathbf{X}$$ in order to obtain the second set of equations. So, I was wondering if there are alternative ways to solve the problem.

Edit : Using properties of trace operation and the constraint, I converted the problem into maximizing $$\text{Tr}(\mathbf{X}^T\mathbf{A}\mathbf{X})$$ with the same constraints.

• Are there any assumptions on what A is? Like is it positive semi-definite, positive definite, or anything like that? – spektr Apr 12 at 21:04
• A is of the form B times B transpose – Buna Apr 12 at 21:08
• Are you solving for X? If so, the problem seems poorly formed since any the constraint on X means the problem can't affect the choice of value. – Richard Apr 12 at 23:47
• Yes I am solving for X – Buna Apr 13 at 3:35
• Hint: the Frobenius norm is invariant under orthogonal transformations. You can diagonalize $A$ with no more than $m$ nonzero eigenvalues. – Brian Borchers Apr 13 at 5:15

I believe that there is a good alternative that benefits from the geometry of the parameter space and completely eliminates the need for constrained optimization. If you explicitly wanted to make use of Lagrangians, I will definitely not be answering the question, but I thought it might be worthwhile to consider the perspective I will describe. In particular, the approach I will be presenting uses Riemannian manifold optimization tools.

Note that the parameter of interest $$\mathbf{X}\in \mathbb{R}^{n \times k}$$ is simply an element of the Stiefel Manifold $$\mathcal{M}\equiv V_k(\mathbb{R}^n)$$ (the set of all orthonormal $$k$$-frames in $$\mathbb{R}^n$$): \begin{align} \mathbf{X} &\in V_k(\mathbb{R}^n) \subset \mathbb{R}^{n\times k}\\ V_k(\mathbb{R}^n) = \{\mathbf{X}&\in\mathbb{R}^{n\times k} : \mathbf{X}^\top\mathbf{X}=\mathbf{I}\}. \end{align}

With this definition, we can benefit from the Riemannian structure of $$V_k(\mathbb{R}^n)$$ and optimize the energy with any unconstrained optimization algorithm, such as the Riemannian-gradient descent:

\begin{align} &\text{ While } \mathbf{X}_{k} \text{ does not sufficiently minimize } f \\ &\text{\quad- Pick a gradient related descent direction }\boldsymbol{\eta}_k\in T_{\mathbf{X}_k}\mathcal{M}\\ &\text{\quad- Choose a retraction } R_{\mathbf{X}_k}:T_{\mathbf{X}_k}\mathcal{M}\rightarrow \mathcal{M}.\\ &\text{\quad- Choose a step length } \tau_k\in \mathbb{R}.\\ &\text{\quad- Set } \mathbf{X}_{k+1}\gets R_{\mathbf{X}_k}(\tau_k\boldsymbol{\eta}_k).\\ &\text{\quad- } k\gets k+1 \end{align}

Here $$f$$ denotes the function we optimize: $$f = \frac{1}{2}\|\mathbf{X}\mathbf{X}^\top - \mathbf{A} \|_\mathcal{F}^2 = \frac{1}{2}\text{tr}\big((\mathbf{X}\mathbf{X}^\top-\mathbf{A})(\mathbf{X}\mathbf{X}^\top-\mathbf{A})^\top\big)$$ whose gradient reads: $$\nabla_{\mathbf{X}}f = (\mathbf{X} \mathbf{X}^\top -\mathbf{A})\mathbf{X}+(\mathbf{X}\mathbf{X}^\top +(-\mathbf{A})^\top )\mathbf{X}.$$

We generally pick the gradient related direction $$\boldsymbol{\eta}$$ as the projection of the Euclidean gradient onto the tangent space of the manifold (or for a broad class of manifolds including the Stiefel manifold, we can instead use the logarithmic-map):

$$\boldsymbol{\eta}\triangleq \text{grad }_{\mathbf{X}} f = \Pi_{\mathbf{X}} \Big( -\nabla_{\mathbf{X}} f \Big)$$

For now let us assume that $$\tau_k$$ is fixed, e.g. $$\tau_k=0.1$$. The retraction operator $$R(\cdot)$$ of the Stiefel manifold is analytically defined or in other words the true exponential map is available. More details on that are given in this math-se post and of course in the seminal paper of Edelman et. al:

Edelman, Alan, Tomás A. Arias, and Steven T. Smith. "The geometry of algorithms with orthogonality constraints." SIAM journal on Matrix Analysis and Applications 20.2 (1998): 303-353.

All the Riemannian operators for Stiefel manifold are included in the toolboxes such as Manopt or ROPTLIB. Finally, we can update the obtain the new value $$\mathbf{X}_{k+1}$$ by: \begin{align} \mathbf{X}_{k+1} = R_{\mathbf{X}_k}(\tau_k \boldsymbol{\eta}_k) \end{align} Note that $$\mathbf{X}_{k+1}$$ will be an orthonormal $$k$$-frame obeying the aforementioned constrained naturally. Riemennian gradient descent is the simplest choice of Riemannian optimization and there are many others such as Riemannian-LBFGS or Riemannian-Trust Region. Many of those choices vary in how they compute the step size $$\tau_k$$, usually by a form of line-search.

More recently, Hu et al. considered a very similar problem to the one of this questions where the minimization with orthonormality constraints (again, Stiefel manifold) is made efficient. This method too, uses the Riemannian structure of the problem.