Home > manopt > manifolds > fixedrank > fixedrankfactory_2factors_subspace_projection.m

fixedrankfactory_2factors_subspace_projection

PURPOSE ^

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

SYNOPSIS ^

function M = fixedrankfactory_2factors_subspace_projection(m, n, k)

DESCRIPTION ^

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

 function M = fixedrankfactory_2factors_subspace_projection(m, n, k)

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

 Tangent vectors are represented as a structure with two fields: L, R.

 Note: L is orthonormal, i.e., columns are orthogonal to each other.
 Such a geometry might be of interest where the left factor has a
 subspace interpretation. A motivation is in Sections 3.3 and 6.4 of the
 paper below.

 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: fixedrankfactory_2factors fixedrankembeddedfactory fixedrankfactory_2factors_preconditioned

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = fixedrankfactory_2factors_subspace_projection(m, n, k)
0002 % Manifold of m-by-n matrices of rank k with two factor quotient geometry.
0003 %
0004 % function M = fixedrankfactory_2factors_subspace_projection(m, n, k)
0005 %
0006 % A point X on the manifold is represented as a structure with two
0007 % fields: L and R. The matrix L (mxk) is orthonormal,
0008 % while the matrix R (nxk) is a full column-rank
0009 % matrix such that X = L*R'.
0010 %
0011 % Tangent vectors are represented as a structure with two fields: L, R.
0012 %
0013 % Note: L is orthonormal, i.e., columns are orthogonal to each other.
0014 % Such a geometry might be of interest where the left factor has a
0015 % subspace interpretation. A motivation is in Sections 3.3 and 6.4 of the
0016 % paper below.
0017 %
0018 % Please cite the Manopt paper as well as the research paper:
0019 %     @Article{mishra2014fixedrank,
0020 %       Title   = {Fixed-rank matrix factorizations and {Riemannian} low-rank optimization},
0021 %       Author  = {Mishra, B. and Meyer, G. and Bonnabel, S. and Sepulchre, R.},
0022 %       Journal = {Computational Statistics},
0023 %       Year    = {2014},
0024 %       Number  = {3-4},
0025 %       Pages   = {591--621},
0026 %       Volume  = {29},
0027 %       Doi     = {10.1007/s00180-013-0464-z}
0028 %     }
0029 %
0030 % See also: fixedrankfactory_2factors fixedrankembeddedfactory fixedrankfactory_2factors_preconditioned
0031 
0032 
0033 
0034 % This file is part of Manopt: www.manopt.org.
0035 % Original author: Bamdev Mishra, Dec. 30, 2012.
0036 % Contributors:
0037 % Change log:
0038     
0039     M.name = @() sprintf('LR'' quotient manifold of %dx%d matrices of rank %d', m, n, k);
0040     
0041     M.dim = @() (m+n-k)*k;
0042     
0043     % Some precomputations at the point X to be used in the inner product (and
0044     % pretty much everywhere else).
0045     function X = prepare(X)
0046         if ~all(isfield(X,{'RtR'}) == 1)
0047             X.RtR = X.R'*X.R;
0048         end
0049     end
0050     
0051     % The choice of the metric is motivated by symmetry and scale
0052     % invariance in the total space.
0053     M.inner = @iproduct;
0054     function ip = iproduct(X, eta, zeta)
0055         X = prepare(X);
0056         
0057         ip = eta.L(:).'*zeta.L(:)  + trace(X.RtR\(eta.R'*zeta.R));
0058     end
0059     
0060     M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));
0061     
0062     M.dist = @(x, y) error('fixedrankfactory_2factors_subspace_projection.dist not implemented yet.');
0063     
0064     M.typicaldist = @() 10*k;
0065     
0066     skew = @(X) .5*(X-X');
0067     symm = @(X) .5*(X+X');
0068     stiefel_proj = @(L, H) H - L*symm(L'*H);
0069     
0070     M.egrad2rgrad = @egrad2rgrad;
0071     function rgrad = egrad2rgrad(X, egrad)
0072         X = prepare(X);
0073         
0074         rgrad.L = stiefel_proj(X.L, egrad.L);
0075         rgrad.R = egrad.R*X.RtR;
0076     end
0077     
0078     
0079     M.ehess2rhess = @ehess2rhess;
0080     function Hess = ehess2rhess(X, egrad, ehess, eta)
0081         X = prepare(X);
0082         
0083         % Riemannian gradient.
0084         rgrad = egrad2rgrad(X, egrad);
0085         
0086         % Directional derivative of the Riemannian gradient.
0087         Hess.L = ehess.L - eta.L*symm(X.L'*egrad.L);
0088         Hess.L = stiefel_proj(X.L, Hess.L);
0089         
0090         Hess.R = ehess.R*X.RtR + 2*egrad.R*symm(eta.R'*X.R);
0091         
0092         % Correction factor for the non-constant metric on the factor R.
0093         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)));
0094         
0095         % Projection onto the horizontal space.
0096         Hess = M.proj(X, Hess);
0097     end
0098     
0099     
0100     M.proj = @projection;
0101     function etaproj = projection(X, eta)
0102         X = prepare(X);
0103         
0104         eta.L = stiefel_proj(X.L, eta.L); % On the tangent space.
0105         SS = X.RtR;
0106         AS1 = 2*X.RtR*skew(X.L'*eta.L)*X.RtR;
0107         AS2 = 2*skew(X.RtR*(X.R'*eta.R));
0108         AS  = skew(AS1 + AS2);
0109         
0110         Omega = nested_sylvester(SS,AS);
0111         etaproj.L = eta.L - X.L*Omega;
0112         etaproj.R = eta.R - X.R*Omega;
0113     end
0114     
0115     M.tangent = M.proj;
0116     M.tangent2ambient = @(X, eta) eta;
0117     
0118     M.retr = @retraction;
0119     function Y = retraction(X, eta, t)
0120         if nargin < 3
0121             t = 1.0;
0122         end
0123         Y.L = uf(X.L + t*eta.L);
0124         Y.R = X.R + t*eta.R;
0125         
0126         % These are reused in the computation of the gradient and Hessian.
0127         Y = prepare(Y);
0128     end
0129     
0130     M.exp = @exponential;
0131     function R = exponential(X, eta, t)
0132         if nargin < 3
0133             t = 1.0;
0134         end
0135         
0136         R = retraction(X, eta, t);
0137         warning('manopt:fixedrankfactory_2factors_subspace_projection:exp', ...
0138             ['Exponential for fixed rank ' ...
0139             'manifold not implemented yet. Lsed retraction instead.']);
0140     end
0141     
0142     M.hash = @(X) ['z' hashmd5([X.L(:) ; X.R(:)])];
0143     
0144     M.rand = @random;
0145     % Factors L lives on Stiefel manifold, hence we will reuse
0146     % its random generator.
0147     stiefelm = stiefelfactory(m, k);
0148     function X = random()
0149         X.L = stiefelm.rand();
0150         X.R = randn(n, k);
0151     end
0152     
0153     M.randvec = @randomvec;
0154     function eta = randomvec(X)
0155         eta.L = randn(m, k);
0156         eta.R = randn(n, k);
0157         eta = projection(X, eta);
0158         nrm = M.norm(X, eta);
0159         eta.L = eta.L / nrm;
0160         eta.R = eta.R / nrm;
0161     end
0162     
0163     M.lincomb = @lincomb;
0164     
0165     M.zerovec = @(X) struct('L', zeros(m, k),...
0166         'R', zeros(n, k));
0167     
0168     M.transp = @(x1, x2, d) projection(x2, d);
0169     
0170     % vec and mat are not isometries, because of the scaled inner metric.
0171     M.vec = @(X, U) [U.L(:) ; U.R(:)];
0172     M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ...
0173         'R', reshape(u((m*k+1):end), n, k));
0174     M.vecmatareisometries = @() false;
0175     
0176     
0177 end
0178 
0179 % Linear combination of tangent vectors.
0180 function d = lincomb(x, a1, d1, a2, d2) %#ok<INLSL>
0181     
0182     if nargin == 3
0183         d.L = a1*d1.L;
0184         d.R = a1*d1.R;
0185     elseif nargin == 5
0186         d.L = a1*d1.L + a2*d2.L;
0187         d.R = a1*d1.R + a2*d2.R;
0188     else
0189         error('Bad use of fixedrankfactory_2factors_subspace_projection.lincomb.');
0190     end
0191     
0192 end
0193 
0194 function A = uf(A)
0195     [L, unused, R] = svd(A, 0); %#ok
0196     A = L*R';
0197 end
0198 
0199 function omega = nested_sylvester(sym_mat, asym_mat)
0200     % omega=nested_sylvester(sym_mat,asym_mat)
0201     % This function solves the system of nested Sylvester equations:
0202     %
0203     %     X*sym_mat + sym_mat*X = asym_mat
0204     %     Omega*sym_mat+sym_mat*Omega = X
0205     % Mishra, Meyer, Bonnabel and Sepulchre, 'Fixed-rank matrix factorizations and Riemannian low-rank optimization'
0206     
0207     % Uses built-in lyap function, but does not exploit the fact that it's
0208     % twice the same sym_mat matrix that comes into play.
0209     
0210     X = lyap(sym_mat, -asym_mat);
0211     omega = lyap(sym_mat, -X);
0212     
0213 end
0214 
0215 
0216

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005