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

- 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 = retraction_orthographic(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'. 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