Manifold of m-by-n matrices of rank k with polar quotient geometry. function M = fixedrankfactory_3factors(m, n, k) The first-order geometry follows the balanced quotient geometry described in the paper, "Linear regression under fixed-rank constraints: a Riemannian approach", G. Meyer, S. Bonnabel and R. Sepulchre, ICML 2011. Paper link: http://www.icml-2011.org/papers/350_icmlpaper.pdf. The second-order geometry follows from the paper "Fixed-rank matrix factorizations and Riemannian low-rank optimization", B. Mishra, R. Meyer, S. Bonnabel and R. Sepulchre, Computational Statistics, 29(3 - 4), pp. 591 - 621, 2014. A point X on the manifold is represented as a structure with three fields: L, S and R. The matrices L (mxk) and R (nxk) are orthonormal, while the matrix S (kxk) is a symmetric positive definite full rank matrix. Tangent vectors are represented as a structure with three fields: L, S and R. For first-order geometry, please cite the Manopt paper as well as the research paper: @InProceedings{meyer2011linear, Title = {Linear regression under fixed-rank constraints: a {R}iemannian approach}, Author = {Meyer, G. and Bonnabel, S. and Sepulchre, R.}, Booktitle = {{28th International Conference on Machine Learning}}, Year = {2011}, Organization = {{ICML}} } For second-order geometry, 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} } See also fixedrankembeddedfactory fixedrankfactory_2factors fixedrankfactory_3factors_preconditioned
0001 function M = fixedrankfactory_3factors(m, n, k) 0002 % Manifold of m-by-n matrices of rank k with polar quotient geometry. 0003 % 0004 % function M = fixedrankfactory_3factors(m, n, k) 0005 % 0006 % The first-order geometry follows the balanced quotient geometry described 0007 % in the paper, 0008 % "Linear regression under fixed-rank constraints: a Riemannian approach", 0009 % G. Meyer, S. Bonnabel and R. Sepulchre, ICML 2011. 0010 % 0011 % Paper link: http://www.icml-2011.org/papers/350_icmlpaper.pdf. 0012 % 0013 % The second-order geometry follows from the paper 0014 % "Fixed-rank matrix factorizations and Riemannian low-rank optimization", 0015 % B. Mishra, R. Meyer, S. Bonnabel and R. Sepulchre, 0016 % Computational Statistics, 29(3 - 4), pp. 591 - 621, 2014. 0017 % 0018 % A point X on the manifold is represented as a structure with three 0019 % fields: L, S and R. The matrices L (mxk) and R (nxk) are orthonormal, 0020 % while the matrix S (kxk) is a symmetric positive definite full rank 0021 % matrix. 0022 % 0023 % Tangent vectors are represented as a structure with three fields: L, S 0024 % and R. 0025 % 0026 % 0027 % For first-order geometry, please cite the Manopt paper as well as the research paper: 0028 % @InProceedings{meyer2011linear, 0029 % Title = {Linear regression under fixed-rank constraints: a {R}iemannian approach}, 0030 % Author = {Meyer, G. and Bonnabel, S. and Sepulchre, R.}, 0031 % Booktitle = {{28th International Conference on Machine Learning}}, 0032 % Year = {2011}, 0033 % Organization = {{ICML}} 0034 % } 0035 % For second-order geometry, please cite the Manopt paper as well as the research paper: 0036 % @Article{mishra2014fixedrank, 0037 % Title = {Fixed-rank matrix factorizations and {Riemannian} low-rank optimization}, 0038 % Author = {Mishra, B. and Meyer, G. and Bonnabel, S. and Sepulchre, R.}, 0039 % Journal = {Computational Statistics}, 0040 % Year = {2014}, 0041 % Number = {3-4}, 0042 % Pages = {591--621}, 0043 % Volume = {29}, 0044 % Doi = {10.1007/s00180-013-0464-z} 0045 % } 0046 % 0047 % 0048 % See also fixedrankembeddedfactory fixedrankfactory_2factors fixedrankfactory_3factors_preconditioned 0049 0050 % This file is part of Manopt: www.manopt.org. 0051 % Original author: Bamdev Mishra, Dec. 30, 2012. 0052 % Contributors: 0053 % Change log: 0054 % 0055 % Apr. 18, 2018 (NB): 0056 % Removed dependency on lyap. 0057 % 0058 % Sep. 6, 2018 (NB): 0059 % Removed M.exp() as it was not implemented. 0060 0061 M.name = @() sprintf('LSR'' quotient manifold of %dx%d matrices of rank %d', m, n, k); 0062 0063 M.dim = @() (m+n-k)*k; 0064 0065 % Choice of the metric on the orthnormal space is motivated by the symmetry present in the 0066 % space. The metric on the positive definite space is its natural metric. 0067 M.inner = @(X, eta, zeta) eta.L(:).'*zeta.L(:) + eta.R(:).'*zeta.R(:) ... 0068 + trace( (X.S\eta.S) * (X.S\zeta.S) ); 0069 0070 M.norm = @(X, eta) sqrt(M.inner(X, eta, eta)); 0071 0072 M.dist = @(x, y) error('fixedrankfactory_3factors.dist not implemented yet.'); 0073 0074 M.typicaldist = @() 10*k; 0075 0076 skew = @(X) .5*(X-X'); 0077 symm = @(X) .5*(X+X'); 0078 stiefel_proj = @(L, H) H - L*symm(L'*H); 0079 0080 M.egrad2rgrad = @egrad2rgrad; 0081 function rgrad = egrad2rgrad(X, egrad) 0082 rgrad.L = stiefel_proj(X.L, egrad.L); 0083 rgrad.S = X.S*symm(egrad.S)*X.S; 0084 rgrad.R = stiefel_proj(X.R, egrad.R); 0085 end 0086 0087 0088 M.ehess2rhess = @ehess2rhess; 0089 function Hess = ehess2rhess(X, egrad, ehess, eta) 0090 0091 % Riemannian gradient for the factor S. 0092 rgrad.S = X.S*symm(egrad.S)*X.S; 0093 0094 % Directional derivatives of the Riemannian gradient. 0095 Hess.L = ehess.L - eta.L*symm(X.L'*egrad.L); 0096 Hess.L = stiefel_proj(X.L, Hess.L); 0097 0098 Hess.R = ehess.R - eta.R*symm(X.R'*egrad.R); 0099 Hess.R = stiefel_proj(X.R, Hess.R); 0100 0101 Hess.S = X.S*symm(ehess.S)*X.S + 2*symm(eta.S*symm(egrad.S)*X.S); 0102 0103 % Correction factor for the non-constant metric on the factor S. 0104 Hess.S = Hess.S - symm(eta.S*(X.S\rgrad.S)); 0105 0106 % Projection onto the horizontal space. 0107 Hess = M.proj(X, Hess); 0108 end 0109 0110 0111 M.proj = @projection; 0112 function etaproj = projection(X, eta) 0113 % First, projection onto the tangent space of the total space. 0114 eta.L = stiefel_proj(X.L, eta.L); 0115 eta.R = stiefel_proj(X.R, eta.R); 0116 eta.S = symm(eta.S); 0117 0118 % Then, projection onto the horizontal space. 0119 SS = X.S*X.S; 0120 AS = X.S*(skew(X.L'*eta.L) + skew(X.R'*eta.R) - 2*skew(X.S\eta.S))*X.S; 0121 omega = lyapunov_symmetric(SS, AS); 0122 0123 etaproj.L = eta.L - X.L*omega; 0124 etaproj.S = eta.S - (X.S*omega - omega*X.S); 0125 etaproj.R = eta.R - X.R*omega; 0126 end 0127 0128 M.tangent = M.proj; 0129 M.tangent2ambient = @(X, eta) eta; 0130 0131 M.retr = @retraction; 0132 function Y = retraction(X, eta, t) 0133 if nargin < 3 0134 t = 1.0; 0135 end 0136 0137 L = chol(X.S); 0138 Y.S = L'*expm(L'\(t*eta.S)/L)*L; 0139 Y.L = uf(X.L + t*eta.L); 0140 Y.R = uf(X.R + t*eta.R); 0141 end 0142 0143 0144 M.hash = @(X) ['z' hashmd5([X.L(:) ; X.S(:) ; X.R(:)])]; 0145 0146 M.rand = @random; 0147 % Factors L and R are on Stiefel manifolds, hence we reuse 0148 % their random generators. 0149 stiefelm = stiefelfactory(m, k); 0150 stiefeln = stiefelfactory(n, k); 0151 function X = random() 0152 X.L = stiefelm.rand(); 0153 X.R = stiefeln.rand(); 0154 X.S = diag(1+rand(k, 1)); 0155 end 0156 0157 M.randvec = @randomvec; 0158 function eta = randomvec(X) 0159 % A random vector on the horizontal space. 0160 eta.L = randn(m, k); 0161 eta.R = randn(n, k); 0162 eta.S = randn(k, k); 0163 eta = projection(X, eta); 0164 nrm = M.norm(X, eta); 0165 eta.L = eta.L / nrm; 0166 eta.R = eta.R / nrm; 0167 eta.S = eta.S / nrm; 0168 end 0169 0170 M.lincomb = @lincomb; 0171 0172 M.zerovec = @(X) struct('L', zeros(m, k), 'S', zeros(k, k), ... 0173 'R', zeros(n, k)); 0174 0175 M.transp = @(x1, x2, d) projection(x2, d); 0176 0177 % vec and mat are not isometries, because of the scaled inner metric. 0178 M.vec = @(X, U) [U.L(:) ; U.S(:); U.R(:)]; 0179 M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ... 0180 'S', reshape(u((m*k+1): m*k + k*k), k, k), ... 0181 'R', reshape(u((m*k+ k*k + 1):end), n, k)); 0182 M.vecmatareisometries = @() false; 0183 0184 end 0185 0186 % Linear combination of tangent vectors. 0187 function d = lincomb(x, a1, d1, a2, d2) %#ok<INLSL> 0188 0189 if nargin == 3 0190 d.L = a1*d1.L; 0191 d.R = a1*d1.R; 0192 d.S = a1*d1.S; 0193 elseif nargin == 5 0194 d.L = a1*d1.L + a2*d2.L; 0195 d.R = a1*d1.R + a2*d2.R; 0196 d.S = a1*d1.S + a2*d2.S; 0197 else 0198 error('Bad use of fixedrankfactory_3factors.lincomb.'); 0199 end 0200 0201 end 0202 0203 function A = uf(A) 0204 [L, unused, R] = svd(A, 0); %#ok 0205 A = L*R'; 0206 end