Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It only takes a minute to sign up.
Sign up to join this community
Given symmetric matrix $A \in \Bbb R^{n \times n}$, how does one solve the following quadratic matrix equation in matrix $X \in \Bbb R^{m \times n}$?
$$ X^T X = A $$
I know this must be simple, but I have been wrestling with it for some time now. Thanks in advance for any help.
Since $A$ is symmetric $n \times n$, we can write it as $A=Q \Lambda Q^T$, where $Q$ is an orthogonal matrix whose columns are the eigenvectors of $A$, and $\Lambda$ is a diagonal matrix whose entries are the eigenvalues of $A$.
We want to decompose $A=X^TX$, where $X$ is $m \times n$.
Let $\mathbf{rank}(A)=r<\min(m, n)$, therefore, $A$ is not full rank and have $n-r$ eigenvalues equal to $0$.
Thus, we can form and $n \times r$ matrix $\tilde{Q}$ by deleting the columns of $Q$ corresponding to the eigenvectors associated to the eigenvalues equal to $0$, and also form $\tilde{\Lambda}$, an $r \times r$ diagonal matrix whose entries are the non-zero eigenvalues of $A$. Thus, we have $$A=\tilde{Q} \tilde{\Lambda} \tilde{Q}^T = \tilde{Q} \sqrt{\tilde{\Lambda}} \sqrt{\tilde{\Lambda}} \tilde{Q}^T = (\sqrt{\tilde{\Lambda}} \tilde{Q}^T)^T \sqrt{\tilde{\Lambda}} \tilde{Q}^T = X^TX,$$
where $X := \sqrt{\tilde{\Lambda}} \tilde{Q}^T$, with dimensions $r \times n$.
The only problem is that $\mathbf{rank}(A)=r=\min(m, n)$, for $m=r$. To solve this issue, when forming $\tilde{Q}$, delete all but one of the columns of $Q$ corresponding eigenvectors associated to the eigenvalues equal to $0$. Doing this, $X$ will end up with dimensions $(r+1) \times n$, satisfying the inequality $\mathbf{rank}(A)=r<\min(m, n)$, for $m=r+1$.
The solution can only exist for symmetric $A$. We can then write $A=O^T DO$ with $D$ diagonal (and with non-negative eigenvalues) and $O$ orthogonal. We can then define a square root of $D$ in the obvious way, obtaining $X=\sqrt{D}O$ as one solution. For orthogonal $O'$, a more general solution is $X=O'\sqrt{D}O$.
Assume $m=n$. The rows of $X$ are some $n$ vectors in $\mathbb{R}^n$, say, $v_1,\dots, v_n$. The $(i,j)$ element of $X^tX$ is the inner product $v_i\cdot v_j$. So there are at most ${n\choose 2}+n$ different entries in $X^tX$, but $n^2$ possible entries in a general $n\times n$ matrix $A$. As a result, one cannot expect to solve the equation $X^tX=A$ for every $A$.
I propose a solution of the case where potentially $m \gg n$.
Assume that $\mathbf A$ is positive semidefinite. By eigendecomposition, $\mathbf A = \mathbf Q\mathbf S\mathbf Q^\top$ where $\mathbf Q$ is an orthonormal matrix of shape $n \times n$. Consider another decomposition of $\mathbf A$ such that $\mathbf A = \mathbf P|\mathbf M|^2\mathbf P^\top$ where $|\mathbf M|^2$ is diagonal and $\mathbf P$ is an $n \times m$ matrix ($m > n$) defined as $\mathbf P \triangleq \mathbf Q\boldsymbol\Pi^\top$. Thus, if $\mathbf X^\top \mathbf X = \mathbf A$, then $\mathbf X = |\mathbf M|\boldsymbol\Pi\mathbf Q^\top$. Since $\mathbf Q$ is known given $\mathbf A$, we need only to find $|\mathbf M|$ and $\boldsymbol\Pi$.
$$ \mathbf Q\mathbf S\mathbf Q^\top = \mathbf P|\mathbf M|^2\mathbf P^\top = \mathbf Q\boldsymbol\Pi^\top|\mathbf M|^2\boldsymbol\Pi\mathbf Q^\top $$
Naturally, our objective is to minimize $L(|\mathbf M|,\boldsymbol\Pi)=\|\boldsymbol\Pi^\top|\mathbf M|^2\boldsymbol\Pi-\mathbf S\|_F^2$, where $\|\cdot\|_F$ is the Frobenius norm. We need also apply an orthogonal constraint $\boldsymbol\Pi^\top\boldsymbol\Pi=\mathbf I$ to ensure that the column space of $\mathbf P$ is necessarily isometric to that of $\mathbf Q$. To apply orthogonal constraint in gradient method, see Section 1.1 Constraint-Preserving Update of this paper (too long to put it here).
The solution to such optimization could be nonunique (or must be nonunique?), but we may find a particular solution using coordinate gradient descent, namely first find a feasible $\boldsymbol\Pi_0$, fix $\boldsymbol\Pi_0$ and find a better $|\mathbf M_1|$, fix $|\mathbf M_1|$ and find a better $\boldsymbol\Pi_1$ subject to the constraint, fix $\boldsymbol\Pi_1$ and find a better $|\mathbf M_2|$, fix $|\mathbf M_2|$ and find a better $\boldsymbol\Pi_2$, etc. We may use one-step gradient descent (given a small learning rate) for each optimization step. To find the first feasible $\boldsymbol\Pi_0$, one may look into the doc of scipy.stats.ortho_group.
For the implementation, we may resort to PyTorch, which does the gradient things for us.
For example, given
$$ \mathbf A = \begin{pmatrix} 0.38672121 & -0.72476649\\ -0.72476649 & 2.58435428\\ \end{pmatrix}\,, $$
here is what such algorithm gives for $m=10$:
$$ \mathbf X = \begin{pmatrix} 3.02239738\times 10^{-1}& 8.15764615\times 10^{-2}\\ 1.55106782\times 10^{-1}& 1.21846205\times 10^{-1}\\ 4.31693560\times 10^{-2}&-3.22591722\times 10^{-1}\\ 9.10323361\times 10^{-2}&-2.12030684\times 10^{-1}\\ -1.48292760\times 10^{-2}&-7.20156447\times 10^{-2}\\ -1.26908969\times 10^{-1}&-1.17936127\times 10^{-1}\\ 1.47863458\times 10^{-1}&-2.93733723\times 10^{-1}\\ -5.72824350\times 10^{-2}&-6.42074711\times 10^{-2}\\ 4.68723305\times 10^{-1}&-1.51760316\times 10^{+0}\\ -9.26707466\times 10^{-4}& 3.47160090\times 10^{-2}\\ \end{pmatrix}\,. $$
One may verify that $\mathbf X^\top \mathbf X$ is numerically closed to $\mathbf A$.
