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} }
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 % NB, April 17, 2018: added M.tangent 0036 % NB, April 18, 2018: removed lyap dependency 0037 0038 0039 M.name = @() sprintf('MN'' quotient manifold of %dx%d matrices of rank %d', m, n, k); 0040 0041 M.dim = @() (m+n-k)*k; 0042 0043 % Choice of the metric is motivated by the symmetry present in the 0044 % space. 0045 M.inner = @(X, eta, zeta) eta.M(:).'*zeta.M(:) + eta.N(:).'*zeta.N(:); 0046 0047 M.norm = @(X, eta) sqrt(M.inner(X, eta, eta)); 0048 0049 M.dist = @(x, y) error('fixedrankMNquotientfactory.dist not implemented yet.'); 0050 0051 M.typicaldist = @() 10*k; 0052 0053 symm = @(X) .5*(X+X'); 0054 stiefel_proj = @(M, H) H - M*symm(M'*H); 0055 0056 M.egrad2rgrad = @egrad2rgrad; 0057 function eta = egrad2rgrad(X, eta) 0058 eta.M = stiefel_proj(X.M, eta.M); 0059 end 0060 0061 M.ehess2rhess = @ehess2rhess; 0062 function Hess = ehess2rhess(X, egrad, ehess, eta) 0063 0064 % Directional derivative of the Riemannian gradient. 0065 Hess.M = ehess.M - eta.M*symm(X.M'*egrad.M); 0066 Hess.M = stiefel_proj(X.M, Hess.M); 0067 0068 Hess.N = ehess.N; 0069 0070 % Projection onto the horizontal space. 0071 Hess = M.proj(X, Hess); 0072 end 0073 0074 0075 M.proj = @projection; 0076 function etaproj = projection(X, eta) 0077 0078 % Start by projecting the vector from Rmp x Rnp to the tangent 0079 % space to the total space, that is, eta.M should be in the 0080 % tangent space to Stiefel at X.M and eta.N is arbitrary. 0081 eta.M = stiefel_proj(X.M, eta.M); 0082 0083 % Now project from the tangent space to the horizontal space, that 0084 % is, take care of the quotient. 0085 0086 % First solve a Sylvester equation (A symm., B skew-symm.) 0087 A = X.N'*X.N + eye(k); 0088 B = eta.M'*X.M + eta.N'*X.N; 0089 B = B-B'; 0090 omega = lyapunov_symmetric(A, B); 0091 0092 % And project along the vertical space to the horizontal space. 0093 etaproj.M = eta.M + X.M*omega; 0094 etaproj.N = eta.N + X.N*omega; 0095 0096 end 0097 0098 M.tangent = M.proj; 0099 0100 M.exp = @exponential; 0101 function Y = exponential(X, eta, t) 0102 if nargin < 3 0103 t = 1.0; 0104 end 0105 0106 A = t*X.M'*eta.M; 0107 S = t^2*eta.M'*eta.M; 0108 Y.M = [X.M t*eta.M]*expm([A -S ; eye(k) A])*eye(2*k, k)*expm(-A); 0109 0110 % re-orthonormalize (seems necessary from time to time). 0111 [Q R] = qr(Y.M, 0); 0112 Y.M = Q * diag(sign(diag(R))); 0113 0114 Y.N = X.N + t*eta.N; 0115 0116 end 0117 0118 % Factor M lives on the Stiefel manifold, hence we will reuse its 0119 % random generator. 0120 stiefelm = stiefelfactory(m, k); 0121 0122 M.retr = @retraction; 0123 function Y = retraction(X, eta, t) 0124 if nargin < 3 0125 t = 1.0; 0126 end 0127 0128 Y.M = uf(X.M + t*eta.M); % This is a valid retraction 0129 Y.N = X.N + t*eta.N; 0130 end 0131 0132 M.hash = @(X) ['z' hashmd5([X.M(:) ; X.N(:)])]; 0133 0134 M.rand = @random; 0135 function X = random() 0136 X.M = stiefelm.rand(); 0137 X.N = randn(n, k); 0138 end 0139 0140 M.randvec = @randomvec; 0141 function eta = randomvec(X) 0142 eta.M = randn(m, k); 0143 eta.N = randn(n, k); 0144 eta = projection(X, eta); 0145 nrm = M.norm(X, eta); 0146 eta.M = eta.M / nrm; 0147 eta.N = eta.N / nrm; 0148 end 0149 0150 M.lincomb = @lincomb; 0151 0152 M.zerovec = @(X) struct('M', zeros(m, k), 'N', zeros(n, k)); 0153 0154 M.transp = @(x1, x2, d) projection(x2, d); 0155 0156 end 0157 0158 0159 % Linear combination of tangent vectors 0160 function d = lincomb(x, a1, d1, a2, d2) %#ok<INMSL> 0161 0162 if nargin == 3 0163 d.M = a1*d1.M; 0164 d.N = a1*d1.N; 0165 elseif nargin == 5 0166 d.M = a1*d1.M + a2*d2.M; 0167 d.N = a1*d1.N + a2*d2.N; 0168 else 0169 error('Bad use of fixedrankMNquotientfactory.lincomb.'); 0170 end 0171 0172 end 0173 0174 0175 function A = uf(A) 0176 [L, unused, R] = svd(A, 0); 0177 A = L*R'; 0178 end