Home > manopt > manifolds > fixedrank > fixedrankembeddedfactory.m

fixedrankembeddedfactory

PURPOSE ^

Manifold struct to optimize fixed-rank matrices w/ an embedded geometry.

SYNOPSIS ^

function M = fixedrankembeddedfactory(m, n, k)

DESCRIPTION ^

 Manifold struct to optimize fixed-rank matrices w/ an embedded geometry.

 function M = fixedrankembeddedfactory(m, n, k)

 Manifold of m-by-n real matrices of fixed rank k. This follows the
 embedded geometry described in Bart Vandereycken's 2013 paper:
 "Low-rank matrix completion by Riemannian optimization".
 
 Paper link: http://arxiv.org/pdf/1209.3834.pdf

 A point X on the manifold is represented as a structure with three
 fields: U, S and V. The matrices U (mxk) and V (nxk) are orthonormal,
 while the matrix S (kxk) is any /diagonal/, full rank matrix.
 Following the SVD formalism, X = U*S*V'. Note that the diagonal entries
 of S are not constrained to be nonnegative.

 Tangent vectors are represented as a structure with three fields: Up, M
 and Vp. The matrices Up (mxk) and Vp (mxk) obey Up'*U = 0 and Vp'*V = 0.
 The matrix M (kxk) is arbitrary. Such a structure corresponds to the
 following tangent vector in the ambient space of mxn matrices:
   Z = U*M*V' + Up*V' + U*Vp'
 where (U, S, V) is the current point and (Up, M, Vp) is the tangent
 vector at that point.

 Vectors in the ambient space are best represented as mxn matrices. If
 these are low-rank, they may also be represented as structures with
 U, S, V fields, such that Z = U*S*V'. There are no resitrictions on what
 U, S and V are, as long as their product as indicated yields a real, mxn
 matrix.

 The chosen geometry yields a Riemannian submanifold of the embedding
 space R^(mxn) equipped with the usual trace (Frobenius) inner product.


 Please cite the Manopt paper as well as the research paper:
     @Article{vandereycken2013lowrank,
       Title   = {Low-rank matrix completion by {Riemannian} optimization},
       Author  = {Vandereycken, B.},
       Journal = {SIAM Journal on Optimization},
       Year    = {2013},
       Number  = {2},
       Pages   = {1214--1236},
       Volume  = {23},
       Doi     = {10.1137/110845768}
     }

 See also: fixedrankfactory_2factors fixedrankfactory_3factors

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = fixedrankembeddedfactory(m, n, k)
0002 % Manifold struct to optimize fixed-rank matrices w/ an embedded geometry.
0003 %
0004 % function M = fixedrankembeddedfactory(m, n, k)
0005 %
0006 % Manifold of m-by-n real matrices of fixed rank k. This follows the
0007 % embedded geometry described in Bart Vandereycken's 2013 paper:
0008 % "Low-rank matrix completion by Riemannian optimization".
0009 %
0010 % Paper link: http://arxiv.org/pdf/1209.3834.pdf
0011 %
0012 % A point X on the manifold is represented as a structure with three
0013 % fields: U, S and V. The matrices U (mxk) and V (nxk) are orthonormal,
0014 % while the matrix S (kxk) is any /diagonal/, full rank matrix.
0015 % Following the SVD formalism, X = U*S*V'. Note that the diagonal entries
0016 % of S are not constrained to be nonnegative.
0017 %
0018 % Tangent vectors are represented as a structure with three fields: Up, M
0019 % and Vp. The matrices Up (mxk) and Vp (mxk) obey Up'*U = 0 and Vp'*V = 0.
0020 % The matrix M (kxk) is arbitrary. Such a structure corresponds to the
0021 % following tangent vector in the ambient space of mxn matrices:
0022 %   Z = U*M*V' + Up*V' + U*Vp'
0023 % where (U, S, V) is the current point and (Up, M, Vp) is the tangent
0024 % vector at that point.
0025 %
0026 % Vectors in the ambient space are best represented as mxn matrices. If
0027 % these are low-rank, they may also be represented as structures with
0028 % U, S, V fields, such that Z = U*S*V'. There are no resitrictions on what
0029 % U, S and V are, as long as their product as indicated yields a real, mxn
0030 % matrix.
0031 %
0032 % The chosen geometry yields a Riemannian submanifold of the embedding
0033 % space R^(mxn) equipped with the usual trace (Frobenius) inner product.
0034 %
0035 %
0036 % Please cite the Manopt paper as well as the research paper:
0037 %     @Article{vandereycken2013lowrank,
0038 %       Title   = {Low-rank matrix completion by {Riemannian} optimization},
0039 %       Author  = {Vandereycken, B.},
0040 %       Journal = {SIAM Journal on Optimization},
0041 %       Year    = {2013},
0042 %       Number  = {2},
0043 %       Pages   = {1214--1236},
0044 %       Volume  = {23},
0045 %       Doi     = {10.1137/110845768}
0046 %     }
0047 %
0048 % See also: fixedrankfactory_2factors fixedrankfactory_3factors
0049 
0050 % This file is part of Manopt: www.manopt.org.
0051 % Original author: Nicolas Boumal, Dec. 30, 2012.
0052 % Contributors:
0053 % Change log:
0054 %
0055 %    Feb. 20, 2014 (NB):
0056 %       Added function tangent to work with checkgradient.
0057 %
0058 %   June 24, 2014 (NB):
0059 %       A couple modifications following
0060 %       Bart Vandereycken's feedback:
0061 %       - The checksum (hash) was replaced for a faster alternative: it's a
0062 %         bit less "safe" in that collisions could arise with higher
0063 %         probability, but they're still very unlikely.
0064 %       - The vector transport was changed.
0065 %       The typical distance was also modified, hopefully giving the
0066 %       trustregions method a better initial guess for the trust region
0067 %       radius, but that should be tested for different cost functions too.
0068 %
0069 %    July 11, 2014 (NB):
0070 %       Added ehess2rhess and tangent2ambient, supplied by Bart.
0071 %
0072 %    July 14, 2014 (NB):
0073 %       Added vec, mat and vecmatareisometries so that hessianspectrum now
0074 %       works with this geometry. Implemented the tangent function.
0075 %       Made it clearer in the code and in the documentation in what format
0076 %       ambient vectors may be supplied, and generalized some functions so
0077 %       that they should now work with both accepted formats.
0078 %       It is now clearly stated that for a point X represented as a
0079 %       triplet (U, S, V), the matrix S needs to be diagonal.
0080 %
0081 %    Sep.  6, 2018 (NB):
0082 %       Removed M.exp() as it was not implemented.
0083 
0084     M.name = @() sprintf('Manifold of %dx%d matrices of rank %d', m, n, k);
0085     
0086     M.dim = @() (m+n-k)*k;
0087     
0088     M.inner = @(x, d1, d2) d1.M(:).'*d2.M(:) + d1.Up(:).'*d2.Up(:) ...
0089                                              + d1.Vp(:).'*d2.Vp(:);
0090     
0091     M.norm = @(x, d) sqrt(M.inner(x, d, d));
0092     
0093     M.dist = @(x, y) error('fixedrankembeddedfactory.dist not implemented yet.');
0094     
0095     M.typicaldist = @() M.dim();
0096     
0097     % Given Z in tangent vector format, projects the components Up and Vp
0098     % such that they satisfy the tangent space constraints up to numerical
0099     % errors. If Z was indeed a tangent vector at X, this should barely
0100     % affect Z (it would not at all if we had infinite numerical accuracy).
0101     M.tangent = @tangent;
0102     function Z = tangent(X, Z)
0103         Z.Up = Z.Up - X.U*(X.U'*Z.Up);
0104         Z.Vp = Z.Vp - X.V*(X.V'*Z.Vp);
0105     end
0106 
0107     % For a given ambient vector Z, applies it to a matrix W. If Z is given
0108     % as a matrix, this is straightfoward. If Z is given as a structure
0109     % with fields U, S, V such that Z = U*S*V', the product is executed
0110     % efficiently.
0111     function ZW = apply_ambient(Z, W)
0112         if ~isstruct(Z)
0113             ZW = Z*W;
0114         else
0115             ZW = Z.U*(Z.S*(Z.V'*W));
0116         end
0117     end
0118 
0119     % Same as apply_ambient, but applies Z' to W.
0120     function ZtW = apply_ambient_transpose(Z, W)
0121         if ~isstruct(Z)
0122             ZtW = Z'*W;
0123         else
0124             ZtW = Z.V*(Z.S'*(Z.U'*W));
0125         end
0126     end
0127     
0128     % Orthogonal projection of an ambient vector Z represented as an mxn
0129     % matrix or as a structure with fields U, S, V to the tangent space at
0130     % X, in a tangent vector structure format.
0131     M.proj = @projection;
0132     function Zproj = projection(X, Z)
0133             
0134         ZV = apply_ambient(Z, X.V);
0135         UtZV = X.U'*ZV;
0136         ZtU = apply_ambient_transpose(Z, X.U);
0137 
0138         Zproj.M = UtZV;
0139         Zproj.Up = ZV  - X.U*UtZV;
0140         Zproj.Vp = ZtU - X.V*UtZV';
0141 
0142     end
0143 
0144     M.egrad2rgrad = @projection;
0145     
0146     % Code supplied by Bart.
0147     % Given the Euclidean gradient at X and the Euclidean Hessian at X
0148     % along H, where egrad and ehess are vectors in the ambient space and H
0149     % is a tangent vector at X, returns the Riemannian Hessian at X along
0150     % H, which is a tangent vector.
0151     M.ehess2rhess = @ehess2rhess;
0152     function rhess = ehess2rhess(X, egrad, ehess, H)
0153         
0154         % Euclidean part
0155         rhess = projection(X, ehess);
0156         
0157         % Curvature part
0158         T = apply_ambient(egrad, H.Vp)/X.S;
0159         rhess.Up = rhess.Up + (T - X.U*(X.U'*T));
0160         T = apply_ambient_transpose(egrad, H.Up)/X.S;
0161         rhess.Vp = rhess.Vp + (T - X.V*(X.V'*T));
0162         
0163     end
0164 
0165     % Transforms a tangent vector Z represented as a structure (Up, M, Vp)
0166     % into a structure with fields (U, S, V) that represents that same
0167     % tangent vector in the ambient space of mxn matrices, as U*S*V'.
0168     % This matrix is equal to X.U*Z.M*X.V' + Z.Up*X.V' + X.U*Z.Vp'. The
0169     % latter is an mxn matrix, which could be too large to build
0170     % explicitly, and this is why we return a low-rank representation
0171     % instead. Note that there are no guarantees on U, S and V other than
0172     % that USV' is the desired matrix. In particular, U and V are not (in
0173     % general) orthonormal and S is not (in general) diagonal.
0174     % (In this implementation, S is identity, but this might change.)
0175     M.tangent2ambient_is_identity = false;
0176     M.tangent2ambient = @tangent2ambient;
0177     function Zambient = tangent2ambient(X, Z)
0178         Zambient.U = [X.U*Z.M + Z.Up, X.U];
0179         Zambient.S = eye(2*k);
0180         Zambient.V = [X.V, Z.Vp];
0181     end
0182     
0183     % This retraction is second order, following general results from
0184     % Absil, Malick, "Projection-like retractions on matrix manifolds",
0185     % SIAM J. Optim., 22 (2012), pp. 135-158.
0186     %
0187     % Notice that this retraction is only locally smooth. Indeed, the
0188     % following code exhibits a discontinuity when retracting from
0189     % X = [1 0 ; 0 0] along V = [0 0 ; 0 1]:
0190     %
0191     % M = fixedrankembeddedfactory(2, 2, 1);
0192     % X = struct('U', [1;0], 'V', [1;0], 'S', 1);
0193     % V = struct('Up', [0;1], 'Vp', [0;1], 'M', 1);
0194     % entry = @(M) M(1, 1);
0195     % mat = @(X) X.U*X.S*X.V';
0196     % g = @(t) entry(mat(M.retr(X, V, t)));
0197     % ezplot(g, [-2, 2]);
0198     M.retr = @retraction;
0199     function Y = retraction(X, Z, t)
0200         if nargin < 3
0201             t = 1.0;
0202         end
0203 
0204         % See personal notes June 28, 2012 (NB)
0205         [Qu, Ru] = qr(Z.Up, 0);
0206         [Qv, Rv] = qr(Z.Vp, 0);
0207         
0208         % Calling svds or svd should yield the same result, but BV
0209         % advocated svd is more robust, and it doesn't change the
0210         % asymptotic complexity to call svd then trim rather than call
0211         % svds. Also, apparently Matlab calls ARPACK in a suboptimal way
0212         % for svds in this scenario.
0213         % [Ut St Vt] = svds([X.S+t*Z.M , t*Rv' ; t*Ru , zeros(k)], k);
0214         [Ut, St, Vt] = svd([X.S+t*Z.M , t*Rv' ; t*Ru , zeros(k)]);
0215         
0216         Y.U = [X.U Qu]*Ut(:, 1:k);
0217         Y.V = [X.V Qv]*Vt(:, 1:k);
0218         Y.S = St(1:k, 1:k) + eps*eye(k);
0219         
0220         % equivalent but very slow code
0221         % [U S V] = svds(X.U*X.S*X.V' + t*(X.U*Z.M*X.V' + Z.Up*X.V' + X.U*Z.Vp'), k);
0222         % Y.U = U; Y.V = V; Y.S = S;
0223         
0224     end
0225 
0226 
0227     % Orthographic retraction provided by Teng Zhang. One interst of the
0228     % orthographic retraction is that if matrices are represented in full
0229     % size, it can be computed without any SVDs. If for an application it
0230     % makes sense to represent the matrices in full size, this may be a
0231     % good idea, but it won't shine in the present implementation of the
0232     % manifold.
0233     M.retr_ortho = @retraction_orthographic;
0234     function Y = retraction_orthographic(X, Z, t)
0235         if nargin < 3
0236             t = 1.0;
0237         end
0238         
0239         % First, write Y (the output) as U1*S0*V1', where U1 and V1 are
0240         % orthogonal matrices and S0 is of size r by r.
0241         [U1, ~] = qr(t*(X.U*Z.M  + Z.Up) + X.U*X.S, 0);
0242         [V1, ~] = qr(t*(X.V*Z.M' + Z.Vp) + X.V*X.S, 0);
0243         S0 = (U1'*X.U)*(X.S + t*Z.M)*(X.V'*V1) + ...
0244              t*((U1'*Z.Up)*(X.V'*V1) + (U1'*X.U)*(Z.Vp'*V1));
0245         
0246         % Then, obtain the singular value decomposition of Y.
0247         [U2, S2, V2] = svd(S0);
0248         Y.U = U1*U2;
0249         Y.S = S2;
0250         Y.V = V1*V2;
0251         
0252     end
0253 
0254 
0255     % Less safe but much faster checksum, June 24, 2014.
0256     % Older version right below.
0257     M.hash = @(X) ['z' hashmd5([sum(X.U(:)) ; sum(X.S(:)); sum(X.V(:)) ])];
0258     %M.hash = @(X) ['z' hashmd5([X.U(:) ; X.S(:) ; X.V(:)])];
0259     
0260     M.rand = @random;
0261     % Factors U and V live on Stiefel manifolds, hence we will reuse
0262     % their random generator.
0263     stiefelm = stiefelfactory(m, k);
0264     stiefeln = stiefelfactory(n, k);
0265     function X = random()
0266         X.U = stiefelm.rand();
0267         X.V = stiefeln.rand();
0268         X.S = diag(sort(rand(k, 1), 1, 'descend'));
0269     end
0270     
0271     % Generate a random tangent vector at X.
0272     % TODO: consider a possible imbalance between the three components Up,
0273     % Vp and M, when m, n and k are widely different (which is typical).
0274     M.randvec = @randomvec;
0275     function Z = randomvec(X)
0276         Z.Up = randn(m, k);
0277         Z.Vp = randn(n, k);
0278         Z.M  = randn(k);
0279         Z = tangent(X, Z);
0280         nrm = M.norm(X, Z);
0281         Z.Up = Z.Up / nrm;
0282         Z.Vp = Z.Vp / nrm;
0283         Z.M  = Z.M  / nrm;
0284     end
0285     
0286     M.lincomb = @lincomb;
0287     
0288     M.zerovec = @(X) struct('Up', zeros(m, k), 'M', zeros(k, k), ...
0289                                                         'Vp', zeros(n, k));
0290     
0291     % New vector transport on June 24, 2014 (as indicated by Bart)
0292     % Reference: Absil, Mahony, Sepulchre 2008 section 8.1.3:
0293     % For Riemannian submanifolds of a Euclidean space, it is acceptable to
0294     % transport simply by orthogonal projection of the tangent vector
0295     % translated in the ambient space.
0296     M.transp = @project_tangent;
0297     function Z2 = project_tangent(X1, X2, Z1)
0298         Z2 = projection(X2, tangent2ambient(X1, Z1));
0299     end
0300 
0301 
0302     M.vec = @vec;
0303     function Zvec = vec(X, Z)
0304         Zamb = tangent2ambient(X, Z);
0305         Zamb_mat = Zamb.U*Zamb.S*Zamb.V';
0306         Zvec = Zamb_mat(:);
0307     end
0308     M.mat = @(X, Zvec) projection(X, reshape(Zvec, [m, n]));
0309     M.vecmatareisometries = @() true;
0310 
0311 end
0312 
0313 % Linear combination of tangent vectors
0314 function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>
0315 
0316     if nargin == 3
0317         d.Up = a1*d1.Up;
0318         d.Vp = a1*d1.Vp;
0319         d.M  = a1*d1.M;
0320     elseif nargin == 5
0321         d.Up = a1*d1.Up + a2*d2.Up;
0322         d.Vp = a1*d1.Vp + a2*d2.Vp;
0323         d.M  = a1*d1.M  + a2*d2.M;
0324     else
0325         error('fixedrank.lincomb takes either 3 or 5 inputs.');
0326     end
0327 
0328 end

Generated on Mon 10-Sep-2018 11:48:06 by m2html © 2005