Home > manopt > manifolds > fixedrank > fixedrankfactory_3factors_preconditioned.m

fixedrankfactory_3factors_preconditioned

PURPOSE ^

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

SYNOPSIS ^

function M = fixedrankfactory_3factors_preconditioned(m, n, k)

DESCRIPTION ^

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

 function M = fixedrankfactory_3factors_preconditioned(m, n, k)

 This geometry is tuned to least squares problems such as low-rank matrix
 completion with ell-2 loss.

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

 Tangent vectors are represented as a structure with three fields: L, S
 and R.

 Please cite the Manopt paper as well as the research paper:
     @InProceedings{mishra2014r3mc,
       Title        = {{R3MC}: A {R}iemannian three-factor algorithm for low-rank matrix completion},
       Author       = {Mishra, B. and Sepulchre, R.},
       Booktitle    = {{53rd IEEE Conference on Decision and Control}},
       Year         = {2014},
       Organization = {{IEEE CDC}}
     }


 See also: fixedrankfactory_3factors fixedrankfactory_2factors_preconditioned

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = fixedrankfactory_3factors_preconditioned(m, n, k)
0002 % Manifold of m-by-n matrices of rank k with three factor quotient geometry.
0003 %
0004 % function M = fixedrankfactory_3factors_preconditioned(m, n, k)
0005 %
0006 % This geometry is tuned to least squares problems such as low-rank matrix
0007 % completion with ell-2 loss.
0008 %
0009 % A point X on the manifold is represented as a structure with three
0010 % fields: L, S and R. The matrices L (mxk) and R (nxk) are orthonormal,
0011 % while the matrix S (kxk) is a full rank matrix such that X = L*S*R'.
0012 %
0013 % Tangent vectors are represented as a structure with three fields: L, S
0014 % and R.
0015 %
0016 % Please cite the Manopt paper as well as the research paper:
0017 %     @InProceedings{mishra2014r3mc,
0018 %       Title        = {{R3MC}: A {R}iemannian three-factor algorithm for low-rank matrix completion},
0019 %       Author       = {Mishra, B. and Sepulchre, R.},
0020 %       Booktitle    = {{53rd IEEE Conference on Decision and Control}},
0021 %       Year         = {2014},
0022 %       Organization = {{IEEE CDC}}
0023 %     }
0024 %
0025 %
0026 % See also: fixedrankfactory_3factors fixedrankfactory_2factors_preconditioned
0027 
0028 % This file is part of Manopt: www.manopt.org.
0029 % Original author: Bamdev Mishra, Dec. 30, 2012.
0030 % Contributors:
0031 % Change log:
0032 %
0033 %    April 04, 2015 (BM):
0034 %       Cosmetic changes including avoiding storing the inverse of a kxk matrix.
0035 
0036     
0037     M.name = @() sprintf('LSR'' (tuned for least square problems) quotient manifold of %dx%d matrices of rank %d', m, n, k);
0038     
0039     M.dim = @() (m+n-k)*k;
0040     
0041     % Some precomputations at the point X that are to be used in the inner product (and
0042     % pretty much everywhere else).
0043     function X = prepare(X)
0044         if ~all(isfield(X,{'StS','SSt'}) == 1)
0045             X.SSt = X.S*X.S';
0046             X.StS = X.S'*X.S;
0047         end
0048     end
0049     
0050     % The choice of metric is motivated by symmetry and tuned to least square
0051     % objective function.
0052     M.inner = @iproduct;
0053     function ip = iproduct(X, eta, zeta)
0054         X = prepare(X);
0055         
0056         ip = trace(X.SSt*(eta.L'*zeta.L)) + trace(X.StS*(eta.R'*zeta.R)) ...
0057             + trace(eta.S'*zeta.S);
0058     end
0059     
0060     M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));
0061     
0062     M.dist = @(x, y) error('fixedrankfactory_3factors_preconditioned.dist not implemented yet.');
0063     
0064     M.typicaldist = @() 10*k;
0065     
0066     skew = @(X) .5*(X-X');
0067     symm = @(X) .5*(X+X');
0068     
0069     M.egrad2rgrad = @egrad2rgrad;
0070     function rgrad = egrad2rgrad(X, egrad)
0071         X = prepare(X);
0072         
0073         SSL = X.SSt;
0074         ASL = 2*symm(SSL*(egrad.S*X.S'));
0075         
0076         SSR = X.StS;
0077         ASR = 2*symm(SSR*(egrad.S'*X.S));
0078         
0079         [BL, BR] = tangent_space_lyap(X.S, ASL, ASR); % It computes the solution without calling Matlab's Lyap.
0080         
0081         rgrad.L = (egrad.L - X.L*BL)/X.SSt;
0082         rgrad.R = (egrad.R - X.R*BR)/X.StS;
0083         rgrad.S = egrad.S;
0084         
0085         % Debug
0086         %         BL1 = lyap(SSL, -ASL); % Alternate way
0087         %         BR1 = lyap(SSR, -ASR);
0088         %         norm(skew(X.SSt*(rgrad.L'*X.L) + rgrad.S*X.S'), 'fro')
0089         %         norm(skew(X.StS*(rgrad.R'*X.R) - X.S'*rgrad.S), 'fro')
0090         
0091     end
0092     
0093     
0094     
0095     M.ehess2rhess = @ehess2rhess;
0096     function Hess = ehess2rhess(X, egrad, ehess, eta)
0097         X = prepare(X);
0098         
0099         % Riemannian gradient.
0100         SSL = X.SSt;
0101         ASL = 2*symm(SSL*(egrad.S*X.S'));
0102         SSR = X.StS;
0103         ASR = 2*symm(SSR*(egrad.S'*X.S));
0104         [BL, BR] = tangent_space_lyap(X.S, ASL, ASR);
0105         
0106         rgrad.L = (egrad.L - X.L*BL)/X.SSt;
0107         rgrad.R = (egrad.R - X.R*BR)/X.StS;
0108         rgrad.S = egrad.S;
0109         
0110         % Directional derivative of the Riemannian gradient.
0111         ASLdot = 2*symm((2*symm(X.S*eta.S')*(egrad.S*X.S')) + X.SSt*(ehess.S*X.S' + egrad.S*eta.S')) - 4*symm(symm(eta.S*X.S')*BL);
0112         ASRdot = 2*symm((2*symm(X.S'*eta.S)*(egrad.S'*X.S)) + X.StS*(ehess.S'*X.S + egrad.S'*eta.S)) - 4*symm(symm(eta.S'*X.S)*BR);
0113         
0114         %         SSLdot = X.SSt;
0115         %         SSRdot = X.StS;
0116         %         BLdot = lyap(SSLdot, -ASLdot);
0117         %         BRdot = lyap(SSRdot, -ASRdot);
0118         
0119         [BLdot, BRdot] = tangent_space_lyap(X.S, ASLdot, ASRdot);
0120         
0121         Hess.L = (ehess.L - eta.L*BL - X.L*BLdot - 2*rgrad.L*symm(eta.S*X.S'))/X.SSt;
0122         Hess.R = (ehess.R - eta.R*BR - X.R*BRdot - 2*rgrad.R*symm(eta.S'*X.S))/X.StS;
0123         Hess.S = ehess.S;
0124         
0125         
0126         
0127         % BM: Till this, everything seems correct.
0128         % We still need a correction factor for the non-constant metric
0129         % that is imposed.
0130         % The computation of the correction factor owes itself to the Koszul formula.
0131         % This corresponds to the Riemannian connection in the Euclidean space with the
0132         % scaled metric.
0133         Hess.L = Hess.L + (eta.L*symm(rgrad.S*X.S') + rgrad.L*symm(eta.S*X.S'))/X.SSt;
0134         Hess.R = Hess.R + (eta.R*symm(rgrad.S'*X.S) + rgrad.R*symm(eta.S'*X.S))/X.StS;
0135         Hess.S = Hess.S - symm(rgrad.L'*eta.L)*X.S - X.S*symm(rgrad.R'*eta.R);
0136         
0137         % The Riemannian connection on the quotient space is the
0138         % projection of the Riemmian connection in the ambient space onto the tangent space of the total space and
0139         % then onto the horizontal space.
0140         % This is accomplished by the following operation.
0141         Hess = M.proj(X, Hess);
0142         
0143         % Debug
0144         %         norm(skew(X.SSt*(Hess.L'*X.L) + Hess.S*X.S'))
0145         %         norm(skew(X.StS*(Hess.R'*X.R) - X.S'*Hess.S))
0146         
0147     end
0148     
0149     
0150     
0151     
0152     M.proj = @projection;
0153     function etaproj = projection(X, eta)
0154         X = prepare(X);
0155         
0156         % First, projection onto the tangent space of the total space.
0157         SSL = X.SSt;
0158         ASL = 2*symm(X.SSt*(X.L'*eta.L)*X.SSt);
0159         BL = lyap(SSL, -ASL);
0160         eta.L = eta.L - X.L*(BL/X.SSt);
0161         
0162         SSR = X.StS;
0163         ASR = 2*symm(X.StS*(X.R'*eta.R)*X.StS);
0164         BR = lyap(SSR, -ASR);
0165         eta.R = eta.R - X.R*(BR/X.StS);
0166         
0167         % Project onto the horizontal space
0168         PU = skew((X.L'*eta.L)*X.SSt) + skew(X.S*eta.S');
0169         PV = skew((X.R'*eta.R)*X.StS)  + skew(X.S'*eta.S);
0170         [Omega1, Omega2] = coupled_lyap(X.S, PU, PV);
0171         %         norm(2*skew(Omega1*X.SSt) - PU -(X.S*Omega2*X.S'),'fro' )
0172         %         norm(2*skew(Omega2*X.StS) - PV -(X.S'*Omega1*X.S),'fro' )
0173         %
0174         
0175         etaproj.L = eta.L - (X.L*Omega1);
0176         etaproj.S = eta.S - (X.S*Omega2 - Omega1*X.S) ;
0177         etaproj.R = eta.R - (X.R*Omega2);
0178         
0179         
0180         % Debug
0181         %         norm(skew(X.SSt*(etaproj.L'*X.L) + etaproj.S*X.S'))
0182         %         norm(skew(X.StS*(etaproj.R'*X.R) - X.S'*etaproj.S))
0183         %
0184         %         norm(skew(X.SSt*(etaproj.L'*X.L) - X.S*etaproj.S'))
0185         %         norm(skew(X.StS*(etaproj.R'*X.R) + etaproj.S'*X.S))
0186         
0187     end
0188     
0189     
0190     M.tangent = M.proj;
0191     M.tangent2ambient = @(X, eta) eta;
0192     
0193     M.retr = @retraction;
0194     function Y = retraction(X, eta, t)
0195         if nargin < 3
0196             t = 1.0;
0197         end
0198         
0199         Y.S = (X.S + t*eta.S);
0200         Y.L = uf((X.L + t*eta.L));
0201         Y.R = uf((X.R + t*eta.R));
0202         
0203         Y = prepare(Y);
0204     end
0205     
0206     M.exp = @exponential;
0207     function Y = exponential(X, eta, t)
0208         if nargin < 3
0209             t = 1.0;
0210         end
0211         Y = retraction(X, eta, t);
0212         warning('manopt:fixedrankfactory_3factors_preconditioned:exp', ...
0213             ['Exponential for fixed rank ' ...
0214             'manifold not implemented yet. Used retraction instead.']);
0215     end
0216     
0217     M.hash = @(X) ['z' hashmd5([X.L(:) ; X.S(:) ; X.R(:)])];
0218     
0219     M.rand = @random;
0220     % Factors L and R live on Stiefel manifolds, hence we will reuse
0221     % their random generator.
0222     stiefelm = stiefelfactory(m, k);
0223     stiefeln = stiefelfactory(n, k);
0224     function X = random()
0225         X.L = stiefelm.rand();
0226         X.R = stiefeln.rand();
0227         X.S = diag(1+rand(k, 1));
0228         
0229         X = prepare(X);
0230     end
0231     
0232     M.randvec = @randomvec;
0233     function eta = randomvec(X)
0234         % A random vector on the horizontal space
0235         eta.L = randn(m, k);
0236         eta.R = randn(n, k);
0237         eta.S = randn(k, k);
0238         eta = projection(X, eta);
0239         nrm = M.norm(X, eta);
0240         eta.L = eta.L / nrm;
0241         eta.R = eta.R / nrm;
0242         eta.S = eta.S / nrm;
0243     end
0244     
0245     M.lincomb = @lincomb;
0246     
0247     M.zerovec = @(X) struct('L', zeros(m, k), 'S', zeros(k, k), ...
0248         'R', zeros(n, k));
0249     
0250     M.transp = @(x1, x2, d) projection(x2, d);
0251     
0252     % vec and mat are not isometries, because of the unusual inner metric.
0253     M.vec = @(X, U) [U.L(:) ; U.S(:); U.R(:)];
0254     M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ...
0255         'S', reshape(u((m*k+1): m*k + k*k), k, k), ...
0256         'R', reshape(u((m*k+ k*k + 1):end), n, k));
0257     M.vecmatareisometries = @() false;
0258     
0259 end
0260 
0261 % Linear combination of tangent vectors
0262 function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>
0263     
0264     if nargin == 3
0265         d.L = a1*d1.L;
0266         d.R = a1*d1.R;
0267         d.S = a1*d1.S;
0268     elseif nargin == 5
0269         d.L = a1*d1.L + a2*d2.L;
0270         d.R = a1*d1.R + a2*d2.R;
0271         d.S = a1*d1.S + a2*d2.S;
0272     else
0273         error('Bad use of fixedrankfactory_3factors_preconditioned.lincomb.');
0274     end
0275     
0276 end
0277 
0278 function A = uf(A)
0279     [L, unused, R] = svd(A, 0); %#ok
0280     A = L*R';
0281 end
0282 
0283 function[BU, BV] = tangent_space_lyap(R, E, F)
0284     % We intent to solve a linear system    RR^T  BU + BU RR^T  = E
0285     %                                       R^T R BV + BV R^T R = F
0286     % for BU and BV.
0287     %
0288     % This can be solved using two calls to the Matlab's lyap.
0289     % However, we can still have a more efficient implementation
0290     % that does not require the full functionaliyt of Matlab's lyap.
0291     
0292     [U, Sigma, V] = svd(R);
0293     E_mod = U'*E*U;
0294     F_mod = V'*F*V;
0295     b1 = E_mod(:);
0296     b2 = F_mod(:);
0297     
0298     r = size(Sigma, 1);
0299     sig = diag(Sigma); % all the singular values in a vector
0300     sig1 = sig*ones(1, r); % columns repeat
0301     sig1t = sig1'; % rows repeat
0302     s1 = sig1(:);
0303     s2 = sig1t(:);
0304     
0305     % The block elements
0306     a =  s1.^2 + s2.^2; % a column vector
0307     
0308     % Solve the linear system of equations
0309     cu = b1./a; %a.\b1;
0310     cv = b2./a; %a.\b2;
0311     
0312     % Matricize
0313     CU = reshape(cu, r, r);
0314     CV = reshape(cv, r, r);
0315     
0316     % Do the similarity transforms
0317     BU = U*CU*U';
0318     BV = V*CV*V';
0319     
0320     % %% Debug
0321     %
0322     % norm(R*R'*BU + BU*R*R' - E, 'fro');
0323     % norm((Sigma.^2)*CU + CU*(Sigma.^2) - E_mod, 'fro');
0324     % norm(a.*cu - b1, 'fro');
0325     %
0326     % norm(R'*R*BV + BV*R'*R - F, 'fro');
0327     %
0328     % BU1 = lyap(R*R', - E);
0329     % norm(R*R'*BU1 + BU1*R*R' - E, 'fro');
0330     %
0331     % BV1 = lyap(R'*R, - F);
0332     % norm(R'*R*BV1 + BV1*R'*R - F, 'fro');
0333     %
0334     % % as accurate as the lyap
0335     % norm(BU - BU1, 'fro')
0336     % norm(BV - BV1, 'fro')
0337 end
0338 
0339 
0340 
0341 function[Omega1, Omega2] = coupled_lyap(R, E, F)
0342     % We intent to solve the coupled system of Lyapunov equations
0343     %
0344     % RR^T Omega1 + Omega1 RR^T  - R Omega2 R^T = E
0345     % R^T R Omega2 + Omega1 R^T R  - R^T Omega2 R = F,
0346     %
0347     % for Omega1 and Omega2, both are skew symmetric matrices.
0348     %
0349     % Below is an efficient implementation
0350     
0351     [U, Sigma, V] = svd(R);
0352     E_mod = U'*E*U;
0353     F_mod = V'*F*V;
0354     b1 = E_mod(:);
0355     b2 = F_mod(:);
0356     
0357     r = size(Sigma, 1);
0358     sig = diag(Sigma); % All the singular values in a vector
0359     sig1 = sig*ones(1, r); % Columns repeat
0360     sig1t = sig1'; % Rows repeat
0361     s1 = sig1(:);
0362     s2 = sig1t(:);
0363     
0364     % The block elements
0365     a =  s1.^2 + s2.^2; % A column vector
0366     c = s1.*s2;
0367     
0368     % Solve directly using the formula
0369     % A = diag(a);
0370     % C = diag(c);
0371     % Y1_sol = (A*(C\A) - C) \ (b2 + A*(C\b1));
0372     % Y2_sol = A\(b2 + C*Y1_sol);
0373     
0374     Y1_sol = (b2 + (a./c).*b1) ./ ((a.^2)./c - c);
0375     Y2_sol = (b2 + c.*Y1_sol)./a;
0376     
0377     % Matricize
0378     Omega1 = reshape(Y1_sol, r, r);
0379     Omega2 = reshape(Y2_sol, r, r);
0380     
0381     % Do the similarity transforms
0382     Omega1 = U*Omega1*U';
0383     Omega2 = V*Omega2*V';
0384     
0385     % %% Debug: whether we have the right solution.
0386     % norm(R*R'*Omega1 + Omega1*R*R'  - R*Omega2*R' - E, 'fro')
0387     % norm(R'*R*Omega2 + Omega2*R'*R  - R'*Omega1*R - F, 'fro')
0388 end

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