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'. Their 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

- 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.

- low_rank_matrix_completion Given partial observation of a low rank matrix, attempts to complete it.

- function Z = tangent(X, Z)
- function ZW = apply_ambient(Z, W)
- function ZtW = apply_ambient_transpose(Z, W)
- function Zproj = projection(X, Z)
- function rhess = ehess2rhess(X, egrad, ehess, H)
- function Zambient = tangent2ambient(X, Z)
- function Y = retraction(X, Z, t)
- function Y = exponential(X, Z, t)
- function X = random()
- function Z = randomvec(X)
- function Z2 = project_tangent(X1, X2, Z1)
- function Zvec = vec(X, Z)
- function d = lincomb(x, a1, d1, a2, d2)

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'. Their 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 M.exp = @exponential; 0211 function Y = exponential(X, Z, t) 0212 if nargin < 3 0213 t = 1.0; 0214 end 0215 Y = retraction(X, Z, t); 0216 warning('manopt:fixedrankembeddedfactory:exp', ... 0217 ['Exponential for fixed rank ' ... 0218 'manifold not implemented yet. Used retraction instead.']); 0219 end 0220 0221 % Less safe but much faster checksum, June 24, 2014. 0222 % Older version right below. 0223 M.hash = @(X) ['z' hashmd5([sum(X.U(:)) ; sum(X.S(:)); sum(X.V(:)) ])]; 0224 %M.hash = @(X) ['z' hashmd5([X.U(:) ; X.S(:) ; X.V(:)])]; 0225 0226 M.rand = @random; 0227 % Factors U and V live on Stiefel manifolds, hence we will reuse 0228 % their random generator. 0229 stiefelm = stiefelfactory(m, k); 0230 stiefeln = stiefelfactory(n, k); 0231 function X = random() 0232 X.U = stiefelm.rand(); 0233 X.V = stiefeln.rand(); 0234 X.S = diag(sort(rand(k, 1), 1, 'descend')); 0235 end 0236 0237 % Generate a random tangent vector at X. 0238 % TODO: consider a possible imbalance between the three components Up, 0239 % Vp and M, when m, n and k are widely different (which is typical). 0240 M.randvec = @randomvec; 0241 function Z = randomvec(X) 0242 Z.Up = randn(m, k); 0243 Z.Vp = randn(n, k); 0244 Z.M = randn(k); 0245 Z = tangent(X, Z); 0246 nrm = M.norm(X, Z); 0247 Z.Up = Z.Up / nrm; 0248 Z.Vp = Z.Vp / nrm; 0249 Z.M = Z.M / nrm; 0250 end 0251 0252 M.lincomb = @lincomb; 0253 0254 M.zerovec = @(X) struct('Up', zeros(m, k), 'M', zeros(k, k), ... 0255 'Vp', zeros(n, k)); 0256 0257 % New vector transport on June 24, 2014 (as indicated by Bart) 0258 % Reference: Absil, Mahony, Sepulchre 2008 section 8.1.3: 0259 % For Riemannian submanifolds of a Euclidean space, it is acceptable to 0260 % transport simply by orthogonal projection of the tangent vector 0261 % translated in the ambient space. 0262 M.transp = @project_tangent; 0263 function Z2 = project_tangent(X1, X2, Z1) 0264 Z2 = projection(X2, tangent2ambient(X1, Z1)); 0265 end 0266 0267 0268 M.vec = @vec; 0269 function Zvec = vec(X, Z) 0270 Zamb = tangent2ambient(X, Z); 0271 Zamb_mat = Zamb.U*Zamb.S*Zamb.V'; 0272 Zvec = Zamb_mat(:); 0273 end 0274 M.mat = @(X, Zvec) projection(X, reshape(Zvec, [m, n])); 0275 M.vecmatareisometries = @() true; 0276 0277 end 0278 0279 % Linear combination of tangent vectors 0280 function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL> 0281 0282 if nargin == 3 0283 d.Up = a1*d1.Up; 0284 d.Vp = a1*d1.Vp; 0285 d.M = a1*d1.M; 0286 elseif nargin == 5 0287 d.Up = a1*d1.Up + a2*d2.Up; 0288 d.Vp = a1*d1.Vp + a2*d2.Vp; 0289 d.M = a1*d1.M + a2*d2.M; 0290 else 0291 error('fixedrank.lincomb takes either 3 or 5 inputs.'); 0292 end 0293 0294 end

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