The set $S_+(n, k)$ of rank-$k$ symmetric positive semidefinite matrices $X$ of size $n\times n$ is parametrized as $X = YY^T$ where $Y \in \mathbb{R}^{n\times k}_*$ is a full column rank matrix of size $n\times k$. Notice that for all orthogonal matrices $O$ of size $k$, the transformation $Y \rightarrow YO$ leaves $X$ unchanged. As a result, we may consider all matrices of the form $YO$ as being equivalent for our purpose. Quotienting this equivalence relation out results in a quotient space, $\mathbb{R}^{n \times k}_* / {\rm O}(k)$, with ${\rm O}(k)$ the group of orthogonal matrices. A point on the quotient space is (conceptually) an equivalence class $[Y] = \{YO : O \in {\rm O}(k)\}$. Numerically, we represent that equivalence class with (any) matrix $Y$ in $[Y]$.
A function $\phi : S_+(n,k) \to \mathbb{R}: X \mapsto \phi(X) $ induces a function $f: \mathbb{R}^{n \times k}_* \rightarrow \mathbb{R}: Y \mapsto f(Y)$ defined as $f(Y) = \phi(YY^T)$. Of course, $f$ has the appropriate invariance $f(Y) = f(YO)$ for all orthogonal $O$. Thus, $f$ is defined on the total space $\mathbb{R}^{n\times k}_*$ and descends as a well-defined function on the quotient space $\mathbb{R}^{n \times k}_* / {\rm O}(k)$. Endowing the total space with the usual Riemannian structure of a Euclidean space (the inner product $\langle H_1, H_2 \rangle = \operatorname{trace}(H_1^T H_2)$), a Riemannian structure for the quotient space follows.
See the paper [JBAS10] referenced below for details on this geometry, and [AMS08] for general background.
Factory call: M = symfixedrankYYfactory(n, k).
| Name | Formula | Numerical representation |
| Set | $S_+(n,k) = \{X \in \mathbb{R}^{n\times n} : X=X^T\succeq 0, {\rm rank}(X) = k\} \\ = \{ YY^T: Y \in \mathbb{R}^{n\times k}, {\rm rank}(Y) = k\}$ | $X$ is represented by means of $Y$ as a matrix Y of size $n\times k$. |
| Vertical space at $Y$ | $\{ Y\Omega : \Omega \in \mathbb{R}^{k\times k}, \Omega + \Omega^T = 0 \}$ | A vertical vector $U$ at $Y$ is represented as a matrix U of size $n\times k$: it is tangent to the equivalence classes. |
| Horizontal space at $Y$ | $\{ U \in \mathbb{R}^{n\times k}: U^T Y = Y^T U \}$ | A horizontal vector $U$ at $Y$ is represented as a matrix U of size $n\times k$: it is orthogonal to all vertical vectors. |
| Ambient space | $\mathbb{R}^{n\times k}$ | Points and vectors in the ambient space are, naturally, represented as matrices of size $n\times k$. |
The following table shows some of the nontrivial available functions in the structure M. The tutorial page gives more details about the functionality implemented by each function.
| Name | Field usage | Formula |
| Dimension | M.dim() |
$\operatorname{dim}\mathcal{M} = kn - k(k-1)/2$ |
| Metric | M.inner(Y, U, V) |
$\langle U, V\rangle_Y = \operatorname{trace}(U^T V)$ |
| Norm | M.norm(Y, U) |
$\|U\|_Y = \sqrt{\langle U, U \rangle_Y}$ |
| Horizontal space projector | M.proj(Y, H) |
$P_Y(H) = H - Y\Omega$, where H represents a vector in the ambient space. $\Omega$ is a skew-symmetric matrix obtained as the unique solution to the Sylvester equation $\Omega (Y^T Y) + (Y^T Y)\Omega = Y^T H - H^T Y $ |
| Euclidean to Riemannian gradient | M.egrad2rgrad(Y, egrad) |
$\operatorname{grad} f(Y) = \nabla f(Y)$, where egrad represents the Euclidean gradient $\nabla f(Y)$, which is a vector in the ambient space. There is indeed nothing to do, since the Euclidean gradient of $f$ is already a horizontal vector thanks to the fact that $f(Y) = f(YO)$ for all orthogonal $O$. |
| Euclidean to Riemannian Hessian | M.ehess2rhess(Y, egrad, ehess, U) |
$\operatorname{Hess} f(Y)[U] = P_Y(\nabla^2 f(Y)[U])$, where egrad represents the Euclidean gradient $\nabla f(Y)$ (which happens to not be necessary in this case) and ehess represents the Euclidean Hessian $\nabla^2 f(Y)[U]$, both being vectors in the ambient space. |
| Retraction | M.retr(Y, U, t) |
$\operatorname{Retr}_Y(tU) = Y + tU$ |
| Random point | M.rand() |
Returns a point at random as follows: generate $Y$ with i.i.d. normal entries; return $Y$. With high probability, $Y$ will indeed be full-rank. |
| Random vector | M.randvec(Y) |
Returns a horizontal vector at $Y$ with uniformly random direction, obtained as follows: generate $H$ with i.i.d. normal entries; return $P_Y(H)/\|P_Y(H)\|_Y$. |
| Vector transport | M.transp(Y, Z, U) |
$\operatorname{Transp}_{Z\leftarrow Y}(U) = P_Z(U)$, where $U$ is a horizontal vector at $Y$ that is transported to the horizontal space at $Z$. |
This is an example of low-rank matrix approximation. Let $A\in\mathbb{R}^{n\times n}$ be a symmetric matrix. We minimize the following cost function:
$$\phi(X) =\| X - A \|_F ^2,$$
such that $X \in S_+(n, k)$. $X$ is factorized as $YY^T$ with $Y \in \mathbb{R}^{n \times k}_*$. The function $f:\mathbb{R}^{n \times k}_* \rightarrow \mathbb{R}$ is defined as $f(Y) = \| YY^T - A \|_F ^2 $. Compute the Euclidean gradient and Hessian of $f$:
$$\nabla f(Y) = 4(YY^T - A)Y,$$
$$\nabla^2 f(Y)[U] = 4(YU^T + UY^T)Y + 4(YY^T - A)U.$$
The Riemannian gradient and Hessian are obtained by applying the M.egrad2rgrad and M.ehess2rhess operators:
$$\operatorname{grad} f(x) = 4(YY^T - A)Y,$$
$$\operatorname{Hess} f(Y)[U] = \nabla^2 f(Y)[U] - Y\Omega$$ where $\Omega$ is obtained by solving the Sylvester equation, $$ \Omega (Y^T Y) + (Y^T Y )\Omega = Y^T \nabla^2 f(Y)[U] - \nabla^2 f(Y)[U] ^T Y. $$
The Sylvester equation can be solved efficiently using lyap of MATLAB. Notice that there is no need to compute these explicitly: it suffices to write code for the Euclidean quantities and to apply the conversion tools on them to obtain the Riemannian quantities, as in the following code:
% Generate the problem data. n = 1000; k = 5; Y_org = randn(n,k); A = Y_org*Y_org'; % Create the problem structure. manifold = symfixedrankYYfactory(n,k); problem.M = manifold; % Define the problem cost function and its gradient.
problem.cost = @(Y) norm( Y*Y' - A, 'fro')^2;
egrad = @(Y) 4*(Y*Y' - A)*Y;
ehess = @(Y, U) 4*(Y*U' + U*Y')*Y + 4*(Y*Y' - A)*U; problem.grad = @(Y) manifold.egrad2rgrad(Y, egrad(Y));
problem.hess = @(Y, U) manifold.ehess2rhess(Y, egrad(Y), ehess(Y, U), U);
% Numerically check the differentials.
checkgradient(problem); pause;
checkhessian(problem); pause;
% Execute the optimization
Y = trustregions(problem);
We are still working on building a collection of examples for Manopt. The relevant ones will be referenced here when the time comes.