Home > manopt > manifolds > fixedrank > fixedrankfactory_2factors_subspace_projection.m

# fixedrankfactory_2factors_subspace_projection

## PURPOSE

Manifold of m-by-n matrices of rank k with two factor quotient geometry.

## SYNOPSIS

function M = fixedrankfactory_2factors_subspace_projection(m, n, k)

## DESCRIPTION

``` Manifold of m-by-n matrices of rank k with two factor quotient geometry.

function M = fixedrankfactory_2factors_subspace_projection(m, n, k)

A point X on the manifold is represented as a structure with two
fields: L and R. The matrix L (mxk) is orthonormal,
while the matrix R (nxk) is a full column-rank
matrix such that X = L*R'.

Tangent vectors are represented as a structure with two fields: L, R.

Note: L is orthonormal, i.e., columns are orthogonal to each other.
Such a geometry might be of interest where the left factor has a
subspace interpretation. A motivation is in Sections 3.3 and 6.4 of the
paper below.

Please cite the Manopt paper as well as the research paper:
@Article{mishra2014fixedrank,
Title   = {Fixed-rank matrix factorizations and {Riemannian} low-rank optimization},
Author  = {Mishra, B. and Meyer, G. and Bonnabel, S. and Sepulchre, R.},
Journal = {Computational Statistics},
Year    = {2014},
Number  = {3-4},
Pages   = {591--621},
Volume  = {29},
Doi     = {10.1007/s00180-013-0464-z}
}

## CROSS-REFERENCE INFORMATION

This function calls:
• stiefelfactory Returns a manifold structure to optimize over orthonormal matrices.
• hashmd5 Computes the MD5 hash of input data.
• lincomb Computes a linear combination of tangent vectors in the Manopt framework.
• lyapunov_symmetric_eig Solves AX + XA = C when A = A', as a pseudo-inverse, given eig(A).
This function is called by:

## SOURCE CODE

```0001 function M = fixedrankfactory_2factors_subspace_projection(m, n, k)
0002 % Manifold of m-by-n matrices of rank k with two factor quotient geometry.
0003 %
0004 % function M = fixedrankfactory_2factors_subspace_projection(m, n, k)
0005 %
0006 % A point X on the manifold is represented as a structure with two
0007 % fields: L and R. The matrix L (mxk) is orthonormal,
0008 % while the matrix R (nxk) is a full column-rank
0009 % matrix such that X = L*R'.
0010 %
0011 % Tangent vectors are represented as a structure with two fields: L, R.
0012 %
0013 % Note: L is orthonormal, i.e., columns are orthogonal to each other.
0014 % Such a geometry might be of interest where the left factor has a
0015 % subspace interpretation. A motivation is in Sections 3.3 and 6.4 of the
0016 % paper below.
0017 %
0018 % Please cite the Manopt paper as well as the research paper:
0019 %     @Article{mishra2014fixedrank,
0020 %       Title   = {Fixed-rank matrix factorizations and {Riemannian} low-rank optimization},
0021 %       Author  = {Mishra, B. and Meyer, G. and Bonnabel, S. and Sepulchre, R.},
0022 %       Journal = {Computational Statistics},
0023 %       Year    = {2014},
0024 %       Number  = {3-4},
0025 %       Pages   = {591--621},
0026 %       Volume  = {29},
0027 %       Doi     = {10.1007/s00180-013-0464-z}
0028 %     }
0029 %
0031
0032
0033
0034 % This file is part of Manopt: www.manopt.org.
0035 % Original author: Bamdev Mishra, Dec. 30, 2012.
0036 % Contributors:
0037 % Change log:
0038 %
0039 %    Apr. 18, 2018 (NB):
0040 %        Removed lyap dependency.
0041 %    Aug. 31, 2018 (NB):
0042 %        Improved efficiency of nested_sylvester using lyapunov_symmetric_eig.
0043 %    Sep.  6, 2018 (NB):
0044 %        Removed M.exp() as it was not implemented.
0045
0046     M.name = @() sprintf('LR'' quotient manifold of %dx%d matrices of rank %d', m, n, k);
0047
0048     M.dim = @() (m+n-k)*k;
0049
0050     % Some precomputations at the point X to be used in the inner product (and
0051     % pretty much everywhere else).
0052     function X = prepare(X)
0053         if ~all(isfield(X,{'RtR'}) == 1)
0054             X.RtR = X.R'*X.R;
0055         end
0056     end
0057
0058     % The choice of the metric is motivated by symmetry and scale
0059     % invariance in the total space.
0060     M.inner = @iproduct;
0061     function ip = iproduct(X, eta, zeta)
0062         X = prepare(X);
0063
0064         ip = eta.L(:).'*zeta.L(:)  + trace(X.RtR\(eta.R'*zeta.R));
0065     end
0066
0067     M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));
0068
0069     M.dist = @(x, y) error('fixedrankfactory_2factors_subspace_projection.dist not implemented yet.');
0070
0071     M.typicaldist = @() 10*k;
0072
0073     skew = @(X) .5*(X-X');
0074     symm = @(X) .5*(X+X');
0075     stiefel_proj = @(L, H) H - L*symm(L'*H);
0076
0079         X = prepare(X);
0080
0083     end
0084
0085
0086     M.ehess2rhess = @ehess2rhess;
0087     function Hess = ehess2rhess(X, egrad, ehess, eta)
0088         X = prepare(X);
0089
0092
0093         % Directional derivative of the Riemannian gradient.
0094         Hess.L = ehess.L - eta.L*symm(X.L'*egrad.L);
0095         Hess.L = stiefel_proj(X.L, Hess.L);
0096
0097         Hess.R = ehess.R*X.RtR + 2*egrad.R*symm(eta.R'*X.R);
0098
0099         % Correction factor for the non-constant metric on the factor R.
0101
0102         % Projection onto the horizontal space.
0103         Hess = M.proj(X, Hess);
0104     end
0105
0106
0107     M.proj = @projection;
0108     function etaproj = projection(X, eta)
0109         X = prepare(X);
0110
0111         eta.L = stiefel_proj(X.L, eta.L); % On the tangent space.
0112         SS = X.RtR;
0113         AS1 = 2*X.RtR*skew(X.L'*eta.L)*X.RtR;
0114         AS2 = 2*skew(X.RtR*(X.R'*eta.R));
0115         AS  = skew(AS1 + AS2);
0116
0117         Omega = nested_sylvester(SS, AS);
0118         etaproj.L = eta.L - X.L*Omega;
0119         etaproj.R = eta.R - X.R*Omega;
0120     end
0121
0122     M.tangent = M.proj;
0123     M.tangent2ambient = @(X, eta) eta;
0124
0125     M.retr = @retraction;
0126     function Y = retraction(X, eta, t)
0127         if nargin < 3
0128             t = 1.0;
0129         end
0130         Y.L = uf(X.L + t*eta.L);
0131         Y.R = X.R + t*eta.R;
0132
0133         % These are reused in the computation of the gradient and Hessian.
0134         Y = prepare(Y);
0135     end
0136
0137
0138     M.hash = @(X) ['z' hashmd5([X.L(:) ; X.R(:)])];
0139
0140     M.rand = @random;
0141     % Factors L lives on Stiefel manifold, hence we will reuse
0142     % its random generator.
0143     stiefelm = stiefelfactory(m, k);
0144     function X = random()
0145         X.L = stiefelm.rand();
0146         X.R = randn(n, k);
0147     end
0148
0149     M.randvec = @randomvec;
0150     function eta = randomvec(X)
0151         eta.L = randn(m, k);
0152         eta.R = randn(n, k);
0153         eta = projection(X, eta);
0154         nrm = M.norm(X, eta);
0155         eta.L = eta.L / nrm;
0156         eta.R = eta.R / nrm;
0157     end
0158
0159     M.lincomb = @lincomb;
0160
0161     M.zerovec = @(X) struct('L', zeros(m, k),...
0162         'R', zeros(n, k));
0163
0164     M.transp = @(x1, x2, d) projection(x2, d);
0165
0166     % vec and mat are not isometries, because of the scaled inner metric.
0167     M.vec = @(X, U) [U.L(:) ; U.R(:)];
0168     M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ...
0169         'R', reshape(u((m*k+1):end), n, k));
0170     M.vecmatareisometries = @() false;
0171
0172
0173 end
0174
0175 % Linear combination of tangent vectors.
0176 function d = lincomb(x, a1, d1, a2, d2) %#ok<INLSL>
0177
0178     if nargin == 3
0179         d.L = a1*d1.L;
0180         d.R = a1*d1.R;
0181     elseif nargin == 5
0182         d.L = a1*d1.L + a2*d2.L;
0183         d.R = a1*d1.R + a2*d2.R;
0184     else
0186     end
0187
0188 end
0189
0190 function A = uf(A)
0191     [L, unused, R] = svd(A, 0); %#ok
0192     A = L*R';
0193 end
0194
0195 function omega = nested_sylvester(sym_mat, asym_mat)
0196     % omega = nested_sylvester(sym_mat, asym_mat)
0197     % This function solves the system of nested Sylvester equations:
0198     %
0199     %     X*sym_mat + sym_mat*X = asym_mat
0200     %     omega*sym_mat+sym_mat*omega = X
0201     %
0202     % Mishra, Meyer, Bonnabel and Sepulchre,
0203     % 'Fixed-rank matrix factorizations and Riemannian low-rank optimization'
0204
0205     % Solve each Lyapunov equation efficiently, exploiting the fact
0206     % that it is twice the same sym_mat matrix that comes into play.
0207     [V, lambda] = eig(sym_mat, 'vector');
0208     X = lyapunov_symmetric_eig(V, lambda, asym_mat);
0209     omega = lyapunov_symmetric_eig(V, lambda, X);
0210
0211 end```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005