Manifold of m-by-n matrices of rank k with balanced quotient geometry. function M = fixedrankfactory_2factors(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 two fields: L and R. The matrices L (mxk) and R (nxk) are full column-rank matrices such that X = L*R'. Tangent vectors are represented as a structure with two fields: L, 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_3factors fixedrankfactory_2factors_preconditioned
0001 function M = fixedrankfactory_2factors(m, n, k) 0002 % Manifold of m-by-n matrices of rank k with balanced quotient geometry. 0003 % 0004 % function M = fixedrankfactory_2factors(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 two 0019 % fields: L and R. The matrices L (mxk) and R (nxk) are full column-rank 0020 % matrices such that X = L*R'. 0021 % 0022 % Tangent vectors are represented as a structure with two fields: L, R. 0023 % 0024 % For first-order geometry, please cite the Manopt paper as well as the research paper: 0025 % @InProceedings{meyer2011linear, 0026 % Title = {Linear regression under fixed-rank constraints: a {R}iemannian approach}, 0027 % Author = {Meyer, G. and Bonnabel, S. and Sepulchre, R.}, 0028 % Booktitle = {{28th International Conference on Machine Learning}}, 0029 % Year = {2011}, 0030 % Organization = {{ICML}} 0031 % } 0032 % 0033 % For second-order geometry, please cite the Manopt paper as well as the research paper: 0034 % @Article{mishra2014fixedrank, 0035 % Title = {Fixed-rank matrix factorizations and {Riemannian} low-rank optimization}, 0036 % Author = {Mishra, B. and Meyer, G. and Bonnabel, S. and Sepulchre, R.}, 0037 % Journal = {Computational Statistics}, 0038 % Year = {2014}, 0039 % Number = {3-4}, 0040 % Pages = {591--621}, 0041 % Volume = {29}, 0042 % Doi = {10.1007/s00180-013-0464-z} 0043 % } 0044 % 0045 % 0046 % See also fixedrankembeddedfactory fixedrankfactory_3factors fixedrankfactory_2factors_preconditioned 0047 0048 % This file is part of Manopt: www.manopt.org. 0049 % Original author: Bamdev Mishra, Dec. 30, 2012. 0050 % Contributors: 0051 % Change log: 0052 % 0053 % July 10, 2013 (NB): 0054 % Added vec, mat, tangent, tangent2ambient. 0055 % 0056 % July 03, 2015 (BM): 0057 % Cosmetic changes including avoiding storing the inverse of a 0058 % k-by-k matrix. 0059 % 0060 % Apr. 17, 2018 (NB): 0061 % Using built-in sylvester instead of lyap, which requires a toolbox. 0062 % 0063 % Sep. 6, 2018 (NB): 0064 % Removed M.exp() as it was not implemented. 0065 0066 0067 M.name = @() sprintf('LR'' quotient manifold of %dx%d matrices of rank %d', m, n, k); 0068 0069 M.dim = @() (m+n-k)*k; 0070 0071 % Some precomputations at the point X to be used in the inner product 0072 % (and pretty much everywhere else). 0073 function X = prepare(X) 0074 if ~all(isfield(X,{'LtL','RtR'})) 0075 L = X.L; 0076 R = X.R; 0077 X.LtL = L'*L; 0078 X.RtR = R'*R; 0079 end 0080 end 0081 0082 % Choice of the metric is motivated by the symmetry present in the 0083 % space. The metric is the natural Grassmannian metric on L and R. 0084 M.inner = @iproduct; 0085 function ip = iproduct(X, eta, zeta) 0086 X = prepare(X); 0087 ip = trace(X.LtL\(eta.L'*zeta.L)) + trace( X.RtR\(eta.R'*zeta.R)); 0088 end 0089 0090 M.norm = @(X, eta) sqrt(M.inner(X, eta, eta)); 0091 0092 M.dist = @(x, y) error('fixedrankfactory_2factors.dist not implemented yet.'); 0093 0094 M.typicaldist = @() 10*k; 0095 0096 symm = @(M) .5*(M+M'); 0097 0098 M.egrad2rgrad = @egrad2rgrad; 0099 function rgrad = egrad2rgrad(X, egrad) 0100 X = prepare(X); 0101 rgrad.L = egrad.L*X.LtL; 0102 rgrad.R = egrad.R*X.RtR; 0103 end 0104 0105 M.ehess2rhess = @ehess2rhess; 0106 function Hess = ehess2rhess(X, egrad, ehess, eta) 0107 X = prepare(X); 0108 0109 % Riemannian gradient computation. 0110 rgrad = egrad2rgrad(X, egrad); 0111 0112 % Directional derivative of the Riemannian gradient. 0113 Hess.L = ehess.L*X.LtL + 2*egrad.L*symm(eta.L'*X.L); 0114 Hess.R = ehess.R*X.RtR + 2*egrad.R*symm(eta.R'*X.R); 0115 0116 % We need a correction term for the non-constant metric. 0117 Hess.L = Hess.L - rgrad.L*(X.LtL\(symm(X.L'*eta.L))) - eta.L*(X.LtL\(symm(X.L'*rgrad.L))) + X.L*(X.LtL\(symm(eta.L'*rgrad.L))); 0118 Hess.R = Hess.R - rgrad.R*(X.RtR\(symm(X.R'*eta.R))) - eta.R*(X.RtR\(symm(X.R'*rgrad.R))) + X.R*(X.RtR\(symm(eta.R'*rgrad.R))); 0119 0120 % Projection onto the horizontal space. 0121 Hess = M.proj(X, Hess); 0122 end 0123 0124 M.proj = @projection; 0125 % Projection of the vector eta in the ambient space onto the horizontal space. 0126 function etaproj = projection(X, eta) 0127 X = prepare(X); 0128 0129 SS = (X.LtL)*(X.RtR); 0130 AS = (X.LtL)*(X.R'*eta.R) - (eta.L'*X.L)*(X.RtR); 0131 Omega = sylvester_nochecks(SS, SS, AS); % could be sped up since SS appears twice 0132 etaproj.L = eta.L + X.L*Omega'; 0133 etaproj.R = eta.R - X.R*Omega; 0134 end 0135 0136 M.tangent = M.proj; 0137 M.tangent2ambient = @(X, eta) eta; 0138 0139 M.retr = @retraction; 0140 function Y = retraction(X, eta, t) 0141 if nargin < 3 0142 t = 1.0; 0143 end 0144 0145 Y.L = X.L + t*eta.L; 0146 Y.R = X.R + t*eta.R; 0147 0148 % Numerical conditioning step: A simpler version. 0149 % We need to ensure that L and R do not have very relative 0150 % skewed norms. 0151 0152 scaling = norm(X.L, 'fro')/norm(X.R, 'fro'); 0153 scaling = sqrt(scaling); 0154 Y.L = Y.L / scaling; 0155 Y.R = Y.R * scaling; 0156 0157 % These are reused in the computation of the gradient and Hessian. 0158 Y = prepare(Y); 0159 end 0160 0161 0162 M.hash = @(X) ['z' hashmd5([X.L(:) ; X.R(:)])]; 0163 0164 M.rand = @random; 0165 function X = random() 0166 % A random point on the total space. 0167 X.L = randn(m, k); 0168 X.R = randn(n, k); 0169 X = prepare(X); 0170 end 0171 0172 M.randvec = @randomvec; 0173 function eta = randomvec(X) 0174 % A random vector in the horizontal space. 0175 eta.L = randn(m, k); 0176 eta.R = randn(n, k); 0177 eta = projection(X, eta); 0178 nrm = M.norm(X, eta); 0179 eta.L = eta.L / nrm; 0180 eta.R = eta.R / nrm; 0181 end 0182 0183 M.lincomb = @lincomb; 0184 0185 M.zerovec = @(X) struct('L', zeros(m, k),'R', zeros(n, k)); 0186 0187 M.transp = @(x1, x2, d) projection(x2, d); 0188 0189 % vec and mat are not isometries, because of the unusual inner metric. 0190 M.vec = @(X, U) [U.L(:) ; U.R(:)]; 0191 M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ... 0192 'R', reshape(u((m*k+1):end), n, k)); 0193 M.vecmatareisometries = @() false; 0194 0195 end 0196 0197 % Linear combination of tangent vectors. 0198 function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL> 0199 0200 if nargin == 3 0201 d.L = a1*d1.L; 0202 d.R = a1*d1.R; 0203 elseif nargin == 5 0204 d.L = a1*d1.L + a2*d2.L; 0205 d.R = a1*d1.R + a2*d2.R; 0206 else 0207 error('Bad use of fixedrankfactory_2factors.lincomb.'); 0208 end 0209 0210 end 0211 0212 0213 0214 0215