Home > manopt > manifolds > fixedrank > fixedrankMNquotientfactory.m

fixedrankMNquotientfactory

PURPOSE ^

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

SYNOPSIS ^

function M = fixedrankMNquotientfactory(m, n, k)

DESCRIPTION ^

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

 function M = fixedrankMNquotientfactory(m, n, k)

 This follows the quotient geometry described in the following paper:
 P.-A. Absil, L. Amodei and G. Meyer,
 "Two Newton methods on the manifold of fixed-rank matrices endowed
  with Riemannian quotient geometries", arXiv, 2012.

 Paper link: http://arxiv.org/abs/1209.0068

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

 Tangent vectors are represented as a structure with two fields (M, N).

 Please cite the Manopt paper as well as the research paper:
     @Article{absil2014fixedrank,
       Title   = {Two Newton methods on the manifold of fixed-rank matrices endowed with Riemannian quotient geometries},
       Author  = {Absil, P.-A. and Amodei, L. and Meyer, G.},
       Journal = {Computational Statistics},
       Year    = {2014},
       Number  = {3-4},
       Pages   = {569--590},
       Volume  = {29},
       Doi     = {10.1007/s00180-013-0441-6}
     }

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = fixedrankMNquotientfactory(m, n, k)
0002 % Manifold of m-by-n matrices of rank k with two factor quotient geometry.
0003 %
0004 % function M = fixedrankMNquotientfactory(m, n, k)
0005 %
0006 % This follows the quotient geometry described in the following paper:
0007 % P.-A. Absil, L. Amodei and G. Meyer,
0008 % "Two Newton methods on the manifold of fixed-rank matrices endowed
0009 %  with Riemannian quotient geometries", arXiv, 2012.
0010 %
0011 % Paper link: http://arxiv.org/abs/1209.0068
0012 %
0013 % A point X on the manifold is represented as a structure with two
0014 % fields: M and N. The matrix M (mxk) is orthonormal, while the matrix N
0015 % (nxk) is full-rank such that X = M*N';
0016 %
0017 % Tangent vectors are represented as a structure with two fields (M, N).
0018 %
0019 % Please cite the Manopt paper as well as the research paper:
0020 %     @Article{absil2014fixedrank,
0021 %       Title   = {Two Newton methods on the manifold of fixed-rank matrices endowed with Riemannian quotient geometries},
0022 %       Author  = {Absil, P.-A. and Amodei, L. and Meyer, G.},
0023 %       Journal = {Computational Statistics},
0024 %       Year    = {2014},
0025 %       Number  = {3-4},
0026 %       Pages   = {569--590},
0027 %       Volume  = {29},
0028 %       Doi     = {10.1007/s00180-013-0441-6}
0029 %     }
0030 
0031 % This file is part of Manopt: www.manopt.org.
0032 % Original author: Nicolas Boumal, Dec. 30, 2012.
0033 % Contributors:
0034 % Change log:
0035     
0036     
0037     M.name = @() sprintf('MN'' quotient manifold of %dx%d matrices of rank %d', m, n, k);
0038     
0039     M.dim = @() (m+n-k)*k;
0040     
0041     % Choice of the metric is motivated by the symmetry present in the
0042     % space.
0043     M.inner = @(X, eta, zeta) eta.M(:).'*zeta.M(:) + eta.N(:).'*zeta.N(:);
0044     
0045     M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));
0046     
0047     M.dist = @(x, y) error('fixedrankMNquotientfactory.dist not implemented yet.');
0048     
0049     M.typicaldist = @() 10*k;
0050     
0051     symm = @(X) .5*(X+X');
0052     stiefel_proj = @(M, H) H - M*symm(M'*H);
0053     
0054     M.egrad2rgrad = @egrad2rgrad;
0055     function eta = egrad2rgrad(X, eta)
0056         eta.M = stiefel_proj(X.M, eta.M);
0057     end
0058     
0059     M.ehess2rhess = @ehess2rhess;
0060     function Hess = ehess2rhess(X, egrad, ehess, eta)
0061         
0062         % Directional derivative of the Riemannian gradient.
0063         Hess.M = ehess.M - eta.M*symm(X.M'*egrad.M);
0064         Hess.M = stiefel_proj(X.M, Hess.M);
0065         
0066         Hess.N = ehess.N;
0067         
0068         % Projection onto the horizontal space.
0069         Hess = M.proj(X, Hess);
0070     end
0071     
0072     
0073     M.proj = @projection;
0074     function etaproj = projection(X, eta)
0075         
0076         % Start by projecting the vector from Rmp x Rnp to the tangent
0077         % space to the total space, that is, eta.M should be in the
0078         % tangent space to Stiefel at X.M and eta.N is arbitrary.
0079         eta.M = stiefel_proj(X.M, eta.M);
0080         
0081         % Now project from the tangent space to the horizontal space, that
0082         % is, take care of the quotient.
0083         
0084         % First solve a Sylvester equation (A symm., B skew-symm.)
0085         A = X.N'*X.N + eye(k);
0086         B = eta.M'*X.M + eta.N'*X.N;
0087         B = B-B';
0088         omega = lyap(A, -B);
0089         
0090         % And project along the vertical space to the horizontal space.
0091         etaproj.M = eta.M + X.M*omega;
0092         etaproj.N = eta.N + X.N*omega;
0093         
0094     end
0095     
0096     M.exp = @exponential;
0097     function Y = exponential(X, eta, t)
0098         if nargin < 3
0099             t = 1.0;
0100         end
0101         
0102         A = t*X.M'*eta.M;
0103         S = t^2*eta.M'*eta.M;
0104         Y.M = [X.M t*eta.M]*expm([A -S ; eye(k) A])*eye(2*k, k)*expm(-A);
0105         
0106         % re-orthonormalize (seems necessary from time to time).
0107         [Q R] = qr(Y.M, 0);
0108         Y.M = Q * diag(sign(diag(R)));
0109         
0110         Y.N = X.N + t*eta.N;
0111         
0112     end
0113     
0114     % Factor M lives on the Stiefel manifold, hence we will reuse its
0115     % random generator.
0116     stiefelm = stiefelfactory(m, k);
0117     
0118     M.retr = @retraction;
0119     function Y = retraction(X, eta, t)
0120         if nargin < 3
0121             t = 1.0;
0122         end
0123         
0124         Y.M = uf(X.M + t*eta.M); % This is a valid retraction
0125         Y.N = X.N + t*eta.N;
0126     end
0127     
0128     M.hash = @(X) ['z' hashmd5([X.M(:) ; X.N(:)])];
0129     
0130     M.rand = @random;
0131     function X = random()
0132         X.M = stiefelm.rand();
0133         X.N = randn(n, k);
0134     end
0135     
0136     M.randvec = @randomvec;
0137     function eta = randomvec(X)
0138         eta.M = randn(m, k);
0139         eta.N = randn(n, k);
0140         eta = projection(X, eta);
0141         nrm = M.norm(X, eta);
0142         eta.M = eta.M / nrm;
0143         eta.N = eta.N / nrm;
0144     end
0145     
0146     M.lincomb = @lincomb;
0147     
0148     M.zerovec = @(X) struct('M', zeros(m, k), 'N', zeros(n, k));
0149     
0150     M.transp = @(x1, x2, d) projection(x2, d);
0151     
0152 end
0153 
0154 
0155 % Linear combination of tangent vectors
0156 function d = lincomb(x, a1, d1, a2, d2) %#ok<INMSL>
0157     
0158     if nargin == 3
0159         d.M = a1*d1.M;
0160         d.N = a1*d1.N;
0161     elseif nargin == 5
0162         d.M = a1*d1.M + a2*d2.M;
0163         d.N = a1*d1.N + a2*d2.N;
0164     else
0165         error('Bad use of fixedrankMNquotientfactory.lincomb.');
0166     end
0167     
0168 end
0169 
0170 
0171 function A = uf(A)
0172     [L, unused, R] = svd(A, 0);
0173     A = L*R';
0174 end

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005