Home > manopt > manifolds > fixedrank > fixedrankfactory_2factors.m

fixedrankfactory_2factors

PURPOSE ^

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

SYNOPSIS ^

function M = fixedrankfactory_2factors(m, n, k)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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