We denote the set of rank-$k$ matrices of size $m\times n$ as $\mathbb{R}_k ^{m \times n}$. We parameterize it as $X = LR^T$ where $L \in \mathbb{R}^{m\times k}_k$ and $R \in \mathbb{R}^{n\times k}_k$. Notice that this factorization is not unique. Indeed, for any invertible matrix $M$ of size $k$, $(L, R)$ and $(LM^{-1}, RM^T)$ are mapped to the same $X$. This observation leads to the identification of $\mathbb{R}_k ^{m \times n}$ with the quotient space $(\mathbb{R}_k ^{m \times k} \times \mathbb{R}_k ^{n \times k}) / \mathbb{R}_k ^{k \times k}$. Points on this quotient space are equivalence classes $[(L, R)] = \{(LM^{-1}, RM^T) : M \in \mathbb{R}_k^{k\times k}\}$. Conceptually, we move on the quotient space but matrix representations are done in the total space, $\mathbb{R}_k ^{m \times k} \times \mathbb{R}_k ^{n \times k} $.

A function $\phi :\mathbb{R}_k ^{m \times n} \rightarrow \mathbb{R}: X \mapsto \phi(X) $ induces a function $f: \mathbb{R}_k ^{m \times k} \times \mathbb{R}_k ^{n \times k} \rightarrow \mathbb{R}: (L, R) \mapsto f(L, R)$ where $f(L, R) = \phi(LR^T)$ (though in some developments belows we use the notation $f(X)$ for simplicity). The total space $\mathbb{R}_k ^{m \times k} \times \mathbb{R}_k ^{n \times k} $ is endowed with a scaled metric that is different from the standard Euclidean metric. This choice leads, in turn, to a Riemannian quotient geometry for $\mathbb{R}_k^{m \times n}$.

The following material is referenced from the papers [AAM12] and [MMBS12] (see below).

Factory call: M = fixedrankfactory_2factors(m,n,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.

This is an example of low-rank matrix approximation. Let $A\in\mathbb{R}^{m\times n}$ be given. The best rank-$k$ approximation of $A$ w.r.t. the Frobenius norm can be computed via SVD of $A$. Here, instead, we will compute it using optimization on the manifold we just described. There is no claim that this is a good way to proceed for this particular problem: this is merely an example to illustrate the use of the manifold.

We wish to minimize the following cost function:

$$\phi(X) = \frac{1}{2}\| X - A \|_F ^2,$$

such that $X \in \mathbb{R}_k ^{m \times n}$. $X$ is factorized as $LR ^T$ with $L \in \mathbb{R}_k ^{m \times k}$ and $R \in \mathbb{R}_k ^{n \times k}$. The function $f:\mathbb{R}_k ^{m \times k} \times \mathbb{R}_k ^{n \times k} \rightarrow \mathbb{R}$ is defined as $f(L, R) = (1/2) \cdot \| LR^T - A \|_F ^2 $. Compute the Euclidean gradient and Hessian of $f$ w.r.t. $L$ and $R$:

$$\nabla _L f(X) = (LR^T - A)R,$$

$$\nabla _R f(X) = (LR^T - A)^T L,$$

$$(\nabla^2 f(X)[U])_L = (U_L R^T + LU_R^T )R + (LR^T - A)U_R, $$

$$(\nabla^2 f(X)[U])_R = (U_L R^T + LU_R^T )^T L + (LR^T - A)^T U_L.$$

The Riemannian gradient and Hessian are obtained automatically by applying the M.egrad2rgrad and M.ehess2rhess operators. After the prototyping stage, it may be helpful though to check what the expressions for the Riemannian gradient and Hessian turn out to be, as it is often possible to compute them more efficiently then with the generic conversion tools.

function fixedrank_test()

% Generate the problem data. 

% Problem
m = 500;
n = 1000;
A = randn(m, n);

% We aim for a rank k approximation of A
k = 5;

% Create the problem structure.
manifold = fixedrankfactory_2factors(m,n,k);
problem.M = manifold;

% Define the problem cost function and its gradient.

problem.cost = @cost;
    function f = cost(X)
        f = .5*norm(X.L*X.R' - A, 'fro')^2;
    end

problem.grad = @(X) problem.M.egrad2rgrad(X, egrad(X)); 
    function g = egrad(X)
        P = (X.L*X.R' - A);
        g.L = P*X.R;        
        g.R = P'*X.L;
    end

problem.hess = @(X, U) problem.M.ehess2rhess(X, egrad(X), ehess(X, U), U); 
    function Hess = ehess(X, eta)
        P = (X.L*X.R' - A);
        Pdot = (eta.L*X.R' + X.L*eta.R');
        Hess.L = Pdot*X.R + P*eta.R;
        Hess.R = Pdot'*X.L + P'*eta.L;
    end


% Numerically check the differentials.
checkgradient(problem); pause;
checkhessian(problem); pause;

X = trustregions(problem);

end

Notice that there are a lot of redundant computations in this implementation. For example, computing the gradient at $(L, R)$ will recompute the residue $LR^T - A$ which was probably already computed when the cost function was evaluated at the same point. See the tutorial for ways to avoid redundancy. This can result in significant speed ups.

Notice that the checkhessian test fails: the slope is not right. This is because the retraction is not second-order compatible with the Riemannian exponential on this manifold, making the checkhessian tool unusable. The Hessian is correct though. The warning about the exponential map can be disabled with this command: warning('off', 'manopt:fixedrankfactory_2factors:exp').