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:
• stiefelfactory Returns a manifold structure to optimize over orthonormal matrices.
• hashmd5 Computes the MD5 hash of input data.
• lincomb Computes a linear combination of tangent vectors in the Manopt framework.
This function is called by:

## 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     M.name = @() sprintf('Manifold of %dx%d matrices of rank %d', m, n, k);
0082
0083     M.dim = @() (m+n-k)*k;
0084
0085     M.inner = @(x, d1, d2) d1.M(:).'*d2.M(:) + d1.Up(:).'*d2.Up(:) ...
0086                                              + d1.Vp(:).'*d2.Vp(:);
0087
0088     M.norm = @(x, d) sqrt(M.inner(x, d, d));
0089
0090     M.dist = @(x, y) error('fixedrankembeddedfactory.dist not implemented yet.');
0091
0092     M.typicaldist = @() M.dim();
0093
0094     % Given Z in tangent vector format, projects the components Up and Vp
0095     % such that they satisfy the tangent space constraints up to numerical
0096     % errors. If Z was indeed a tangent vector at X, this should barely
0097     % affect Z (it would not at all if we had infinite numerical accuracy).
0098     M.tangent = @tangent;
0099     function Z = tangent(X, Z)
0100         Z.Up = Z.Up - X.U*(X.U'*Z.Up);
0101         Z.Vp = Z.Vp - X.V*(X.V'*Z.Vp);
0102     end
0103
0104     % For a given ambient vector Z, applies it to a matrix W. If Z is given
0105     % as a matrix, this is straightfoward. If Z is given as a structure
0106     % with fields U, S, V such that Z = U*S*V', the product is executed
0107     % efficiently.
0108     function ZW = apply_ambient(Z, W)
0109         if ~isstruct(Z)
0110             ZW = Z*W;
0111         else
0112             ZW = Z.U*(Z.S*(Z.V'*W));
0113         end
0114     end
0115
0116     % Same as apply_ambient, but applies Z' to W.
0117     function ZtW = apply_ambient_transpose(Z, W)
0118         if ~isstruct(Z)
0119             ZtW = Z'*W;
0120         else
0121             ZtW = Z.V*(Z.S'*(Z.U'*W));
0122         end
0123     end
0124
0125     % Orthogonal projection of an ambient vector Z represented as an mxn
0126     % matrix or as a structure with fields U, S, V to the tangent space at
0127     % X, in a tangent vector structure format.
0128     M.proj = @projection;
0129     function Zproj = projection(X, Z)
0130
0131         ZV = apply_ambient(Z, X.V);
0132         UtZV = X.U'*ZV;
0133         ZtU = apply_ambient_transpose(Z, X.U);
0134
0135         Zproj.M = UtZV;
0136         Zproj.Up = ZV  - X.U*UtZV;
0137         Zproj.Vp = ZtU - X.V*UtZV';
0138
0139     end
0140
0141     M.egrad2rgrad = @projection;
0142
0143     % Code supplied by Bart.
0144     % Given the Euclidean gradient at X and the Euclidean Hessian at X
0145     % along H, where egrad and ehess are vectors in the ambient space and H
0146     % is a tangent vector at X, returns the Riemannian Hessian at X along
0147     % H, which is a tangent vector.
0148     M.ehess2rhess = @ehess2rhess;
0149     function rhess = ehess2rhess(X, egrad, ehess, H)
0150
0151         % Euclidean part
0152         rhess = projection(X, ehess);
0153
0154         % Curvature part
0155         T = apply_ambient(egrad, H.Vp)/X.S;
0156         rhess.Up = rhess.Up + (T - X.U*(X.U'*T));
0157         T = apply_ambient_transpose(egrad, H.Up)/X.S;
0158         rhess.Vp = rhess.Vp + (T - X.V*(X.V'*T));
0159
0160     end
0161
0162     % Transforms a tangent vector Z represented as a structure (Up, M, Vp)
0163     % into a structure with fields (U, S, V) that represents that same
0164     % tangent vector in the ambient space of mxn matrices, as U*S*V'.
0165     % This matrix is equal to X.U*Z.M*X.V' + Z.Up*X.V' + X.U*Z.Vp'. The
0166     % latter is an mxn matrix, which could be too large to build
0167     % explicitly, and this is why we return a low-rank representation
0168     % instead. Note that there are no guarantees on U, S and V other than
0169     % that USV' is the desired matrix. In particular, U and V are not (in
0170     % general) orthonormal and S is not (in general) diagonal.
0171     % (In this implementation, S is identity, but this might change.)
0172     M.tangent2ambient = @tangent2ambient;
0173     function Zambient = tangent2ambient(X, Z)
0174         Zambient.U = [X.U*Z.M + Z.Up, X.U];
0175         Zambient.S = eye(2*k);
0176         Zambient.V = [X.V, Z.Vp];
0177     end
0178
0179     % This retraction is second order, following general results from
0180     % Absil, Malick, "Projection-like retractions on matrix manifolds",
0181     % SIAM J. Optim., 22 (2012), pp. 135-158.
0182     M.retr = @retraction;
0183     function Y = retraction(X, Z, t)
0184         if nargin < 3
0185             t = 1.0;
0186         end
0187
0188         % See personal notes June 28, 2012 (NB)
0189         [Qu, Ru] = qr(Z.Up, 0);
0190         [Qv, Rv] = qr(Z.Vp, 0);
0191
0192         % Calling svds or svd should yield the same result, but BV
0193         % advocated svd is more robust, and it doesn't change the
0194         % asymptotic complexity to call svd then trim rather than call
0195         % svds. Also, apparently Matlab calls ARPACK in a suboptimal way
0196         % for svds in this scenario.
0197         % [Ut St Vt] = svds([X.S+t*Z.M , t*Rv' ; t*Ru , zeros(k)], k);
0198         [Ut, St, Vt] = svd([X.S+t*Z.M , t*Rv' ; t*Ru , zeros(k)]);
0199
0200         Y.U = [X.U Qu]*Ut(:, 1:k);
0201         Y.V = [X.V Qv]*Vt(:, 1:k);
0202         Y.S = St(1:k, 1:k) + eps*eye(k);
0203
0204         % equivalent but very slow code
0205         % [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);
0206         % Y.U = U; Y.V = V; Y.S = S;
0207
0208     end
0209
0210
0211     % Orthographic retraction provided by Teng Zhang. One interst of the
0212     % orthographic retraction is that if matrices are represented in full
0213     % size, it can be computed without any SVDs. If for an application it
0214     % makes sense to represent the matrices in full size, this may be a
0215     % good idea, but it won't shine in the present implementation of the
0216     % manifold.
0217     M.retr_ortho = @retraction_orthographic;
0218     function Y = retraction_orthographic(X, Z, t)
0219         if nargin < 3
0220             t = 1.0;
0221         end
0222
0223         % First, write Y (the output) as U1*S0*V1', where U1 and V1 are
0224         % orthogonal matrices and S0 is of size r by r.
0225         [U1, ~] = qr(t*(X.U*Z.M  + Z.Up) + X.U*X.S, 0);
0226         [V1, ~] = qr(t*(X.V*Z.M' + Z.Vp) + X.V*X.S, 0);
0227         S0 = (U1'*X.U)*(X.S + t*Z.M)*(X.V'*V1) + ...
0228              t*((U1'*Z.Up)*(X.V'*V1) + (U1'*X.U)*(Z.Vp'*V1));
0229
0230         % Then, obtain the singular value decomposition of Y.
0231         [U2, S2, V2] = svd(S0);
0232         Y.U = U1*U2;
0233         Y.S = S2;
0234         Y.V = V1*V2;
0235
0236     end
0237
0238
0239     M.exp = @exponential;
0240     function Y = exponential(X, Z, t)
0241         if nargin < 3
0242             t = 1.0;
0243         end
0244         Y = retraction(X, Z, t);
0245         warning('manopt:fixedrankembeddedfactory:exp', ...
0246                ['Exponential for fixed rank ' ...
0247                 'manifold not implemented yet. Used retraction instead.']);
0248     end
0249
0250     % Less safe but much faster checksum, June 24, 2014.
0251     % Older version right below.
0252     M.hash = @(X) ['z' hashmd5([sum(X.U(:)) ; sum(X.S(:)); sum(X.V(:)) ])];
0253     %M.hash = @(X) ['z' hashmd5([X.U(:) ; X.S(:) ; X.V(:)])];
0254
0255     M.rand = @random;
0256     % Factors U and V live on Stiefel manifolds, hence we will reuse
0257     % their random generator.
0258     stiefelm = stiefelfactory(m, k);
0259     stiefeln = stiefelfactory(n, k);
0260     function X = random()
0261         X.U = stiefelm.rand();
0262         X.V = stiefeln.rand();
0263         X.S = diag(sort(rand(k, 1), 1, 'descend'));
0264     end
0265
0266     % Generate a random tangent vector at X.
0267     % TODO: consider a possible imbalance between the three components Up,
0268     % Vp and M, when m, n and k are widely different (which is typical).
0269     M.randvec = @randomvec;
0270     function Z = randomvec(X)
0271         Z.Up = randn(m, k);
0272         Z.Vp = randn(n, k);
0273         Z.M  = randn(k);
0274         Z = tangent(X, Z);
0275         nrm = M.norm(X, Z);
0276         Z.Up = Z.Up / nrm;
0277         Z.Vp = Z.Vp / nrm;
0278         Z.M  = Z.M  / nrm;
0279     end
0280
0281     M.lincomb = @lincomb;
0282
0283     M.zerovec = @(X) struct('Up', zeros(m, k), 'M', zeros(k, k), ...
0284                                                         'Vp', zeros(n, k));
0285
0286     % New vector transport on June 24, 2014 (as indicated by Bart)
0287     % Reference: Absil, Mahony, Sepulchre 2008 section 8.1.3:
0288     % For Riemannian submanifolds of a Euclidean space, it is acceptable to
0289     % transport simply by orthogonal projection of the tangent vector
0290     % translated in the ambient space.
0291     M.transp = @project_tangent;
0292     function Z2 = project_tangent(X1, X2, Z1)
0293         Z2 = projection(X2, tangent2ambient(X1, Z1));
0294     end
0295
0296
0297     M.vec = @vec;
0298     function Zvec = vec(X, Z)
0299         Zamb = tangent2ambient(X, Z);
0300         Zamb_mat = Zamb.U*Zamb.S*Zamb.V';
0301         Zvec = Zamb_mat(:);
0302     end
0303     M.mat = @(X, Zvec) projection(X, reshape(Zvec, [m, n]));
0304     M.vecmatareisometries = @() true;
0305
0306 end
0307
0308 % Linear combination of tangent vectors
0309 function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>
0310
0311     if nargin == 3
0312         d.Up = a1*d1.Up;
0313         d.Vp = a1*d1.Vp;
0314         d.M  = a1*d1.M;
0315     elseif nargin == 5
0316         d.Up = a1*d1.Up + a2*d2.Up;
0317         d.Vp = a1*d1.Vp + a2*d2.Vp;
0318         d.M  = a1*d1.M  + a2*d2.M;
0319     else
0320         error('fixedrank.lincomb takes either 3 or 5 inputs.');
0321     end
0322
0323 end

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005