Home > manopt > manifolds > fixedranktensors > fixedrankfactory_tucker_preconditioned.m

fixedrankfactory_tucker_preconditioned

PURPOSE ^

Manifold of fixed multilinear rank tensors in Tucker format.

SYNOPSIS ^

function M = fixedrankfactory_tucker_preconditioned(tensor_size, tensor_rank)

DESCRIPTION ^

 Manifold of fixed multilinear rank tensors in Tucker format.

 function M = fixedrankfactory_tucker_preconditioned(tensor_size, tensor_rank)

 n1 = tensor_size(1);
 n2 = tensor_size(2);
 n3 = tensor_size(3);
 r1 = tensor_rank(1);
 r2 = tensor_rank(2);
 r3 = tensor_rank(3);

 A point X on the manifold is represented as a structure with four
 fields: U1, U2, U3 and G. The matrices U1 (n1-by-r1), U2 (n2-by-r2),
 and U3 (n3-by-r3) are orthogonal matrices. G (r1-by-r2-by-r3) is a 
 multidimensional array.

 Tangent vectors are represented as a structure with four fields: 
 U1, U2, U3, and G.

 We exploit the quotient nature of Tucker decompositions to impose a
 scaled inner product on the manifold. This suits least-squares problems.
 For details, refer to the technical report:
 "{R}iemannian preconditioning for tensor completion",
 H. Kasai and B. Mishra, Arxiv preprint arXiv:1506.02159, 2015.

 Paper link: http://arxiv.org/abs/1506.02159.

 Please cite the Manopt paper as well as the research paper:
     @TechReport{kasai2015precon,
       Title   = {{R}iemannian preconditioning for tensor completion},
       Author  = {Kasai, H. and Mishra, B.},
       Journal = {Arxiv preprint arXiv:1506.02159},
       Year    = {2015}
     }

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = fixedrankfactory_tucker_preconditioned(tensor_size, tensor_rank)
0002 % Manifold of fixed multilinear rank tensors in Tucker format.
0003 %
0004 % function M = fixedrankfactory_tucker_preconditioned(tensor_size, tensor_rank)
0005 %
0006 % n1 = tensor_size(1);
0007 % n2 = tensor_size(2);
0008 % n3 = tensor_size(3);
0009 % r1 = tensor_rank(1);
0010 % r2 = tensor_rank(2);
0011 % r3 = tensor_rank(3);
0012 %
0013 % A point X on the manifold is represented as a structure with four
0014 % fields: U1, U2, U3 and G. The matrices U1 (n1-by-r1), U2 (n2-by-r2),
0015 % and U3 (n3-by-r3) are orthogonal matrices. G (r1-by-r2-by-r3) is a
0016 % multidimensional array.
0017 %
0018 % Tangent vectors are represented as a structure with four fields:
0019 % U1, U2, U3, and G.
0020 %
0021 % We exploit the quotient nature of Tucker decompositions to impose a
0022 % scaled inner product on the manifold. This suits least-squares problems.
0023 % For details, refer to the technical report:
0024 % "{R}iemannian preconditioning for tensor completion",
0025 % H. Kasai and B. Mishra, Arxiv preprint arXiv:1506.02159, 2015.
0026 %
0027 % Paper link: http://arxiv.org/abs/1506.02159.
0028 %
0029 % Please cite the Manopt paper as well as the research paper:
0030 %     @TechReport{kasai2015precon,
0031 %       Title   = {{R}iemannian preconditioning for tensor completion},
0032 %       Author  = {Kasai, H. and Mishra, B.},
0033 %       Journal = {Arxiv preprint arXiv:1506.02159},
0034 %       Year    = {2015}
0035 %     }
0036 
0037 % Original authors: Hiroyuki Kasai and Bamdev Mishra, June 5, 2015.
0038 % Contributors:
0039 % Change log:
0040 %
0041 %    Apr. 17, 2018 (NB):
0042 %        Removed dependency on lyap.
0043 %
0044 %    Sep.  6, 2018 (NB):
0045 %        Removed M.exp() as it was not implemented.
0046 
0047     if length(tensor_rank) > 3
0048         error('Bad usage of fixedrankfactory_tucker_preconditioned. Currently, only handles 3-order tensors.');
0049     end
0050     
0051     % Tensor size
0052     n1 = tensor_size(1);
0053     n2 = tensor_size(2);
0054     n3 = tensor_size(3);
0055     
0056     % Core size or multilinear rank
0057     r1 = tensor_rank(1);
0058     r2 = tensor_rank(2);
0059     r3 = tensor_rank(3);
0060     
0061     
0062     speyer1 = speye(r1); % Sparse version of identity that is used in M.proj
0063     speyer2 = speye(r2);
0064     speyer3 = speye(r3);
0065     
0066 
0067     M.name = @() sprintf('G x U1 x U2 x U3 quotient Tucker manifold of %d-by-%d-by-%d tensor of rank %d-by-%d-by-%d.', n1, n2, n3, r1, r2, r3);
0068     
0069     M.dim = @() n1*r1-r1^2 + n2*r2-r2^2 + n3*r3-r3^2 + r1*r2*r3;
0070     
0071     % Some precomputations at point X to be used in the inner product (and
0072     % pretty much everywhere else)
0073     function X = prepare(X)
0074         if ~all(isfield(X,{'G1G1t','G1',...
0075                 'G2G2t','G2', ...
0076                 'G3G3t','G3'}) == 1)
0077             
0078             X.G1 =  reshape(X.G, r1, r2*r3);
0079             X.G1G1t = X.G1*X.G1'; % Positive definite
0080             
0081             
0082             X.G2 = reshape(permute(X.G, [2 1 3]), r2, r1*r3); 
0083             X.G2G2t = X.G2*X.G2'; % Positive definite
0084             
0085             
0086             X.G3 = reshape(permute(X.G, [3 1 2]), r3, r1*r2);  
0087             X.G3G3t = X.G3*X.G3'; % Positive definite
0088         end
0089         
0090         
0091     end
0092     
0093     % Choice of metric is motivated by symmetry and tuned to least-squares
0094     % cost function
0095     M.inner = @iproduct;
0096     function ip = iproduct(X, eta, zeta)
0097         X = prepare(X);
0098         ip =  trace(X.G1G1t*(eta.U1'*zeta.U1)) ...
0099             + trace(X.G2G2t*(eta.U2'*zeta.U2)) ...
0100             + trace(X.G3G3t*(eta.U3'*zeta.U3)) ...
0101             + (eta.G(:)'*zeta.G(:));
0102     end
0103     M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));
0104     
0105     M.dist = @(x, y) error('fixedrankfactory_tucker_preconditioned.dist not implemented yet.');
0106     
0107     M.typicaldist = @() 10*n1*r1; % BM: To do
0108     
0109     skew = @(X) .5*(X-X');
0110     symm = @(X) .5*(X+X');
0111     
0112     M.egrad2rgrad = @egrad2rgrad;
0113     function rgrad = egrad2rgrad(X, egrad)
0114         X = prepare(X); % Reuse already computed terms
0115         
0116         SSU1 = X.G1G1t;
0117         ASU1 = 2*symm(SSU1*(X.U1' * egrad.U1));
0118         
0119         SSU2 = X.G2G2t;
0120         ASU2 = 2*symm(SSU2*(X.U2' * egrad.U2));
0121         
0122         SSU3 = X.G3G3t;
0123         ASU3 = 2*symm(SSU3*(X.U3' * egrad.U3));
0124         
0125         
0126         BU1 = lyapunov_symmetric(SSU1, ASU1);
0127         BU2 = lyapunov_symmetric(SSU2, ASU2);
0128         BU3 = lyapunov_symmetric(SSU3, ASU3);
0129         
0130         % The sylvester solutions ensure that the Riemannian gradient rgrad
0131         % is now on the tangent space. From the Riemannian submersion
0132         % theory, it also belongs to the horizontal space. Therefore,
0133         % no need to further project it on the horizontal space.
0134         
0135         rgrad.U1 = (egrad.U1 - X.U1*BU1)/X.G1G1t;
0136         rgrad.U2 = (egrad.U2 - X.U2*BU2)/X.G2G2t;
0137         rgrad.U3 = (egrad.U3 - X.U3*BU3)/X.G3G3t;
0138         rgrad.G = egrad.G;
0139 
0140         
0141     end
0142     
0143     
0144     
0145     M.ehess2rhess = @ehess2rhess;
0146     function Hess = ehess2rhess(X, egrad, ehess, eta) 
0147         X = prepare(X); % Reuse already computed terms
0148         
0149         % Riemannian gradient
0150         SSU1 = X.G1G1t;
0151         ASU1 = 2*symm(SSU1*(X.U1' * egrad.U1));
0152         SSU2 = X.G2G2t;
0153         ASU2 = 2*symm(SSU2*(X.U2' * egrad.U2));
0154         SSU3 = X.G3G3t;
0155         ASU3 = 2*symm(SSU3*(X.U3' * egrad.U3));
0156         
0157         BU1 = lyapunov_symmetric(SSU1, ASU1);
0158         BU2 = lyapunov_symmetric(SSU2, ASU2);
0159         BU3 = lyapunov_symmetric(SSU3, ASU3);
0160         
0161         rgrad.U1 = (egrad.U1 - X.U1*BU1)/X.G1G1t;
0162         rgrad.U2 = (egrad.U2 - X.U2*BU2)/X.G2G2t;
0163         rgrad.U3 = (egrad.U3 - X.U3*BU3)/X.G3G3t;
0164         rgrad.G = egrad.G;
0165         
0166         % Directional derivative of Riemannian gradient
0167         
0168         eta_G1 = reshape(eta.G, r1, r2*r3); % double(tenmat(eta.G,1));
0169         eta_G2 = reshape(permute(eta.G, [2 1 3]), r2, r1*r3); % double(tenmat(eta.G,2));
0170         eta_G3 = reshape(permute(eta.G, [3 1 2]), r3, r1*r2); % double(tenmat(eta.G,3));
0171         egrad_G1 = reshape(egrad.G, r1, r2*r3); % double(tenmat(egrad.G,1));
0172         egrad_G2 = reshape(permute(egrad.G, [2 1 3]), r2, r1*r3); % double(tenmat(egrad.G,2));
0173         egrad_G3 = reshape(permute(egrad.G, [3 1 2]), r3, r1*r2); % double(tenmat(egrad.G,3));
0174         ehess_G1 = reshape(ehess.G, r1, r2*r3); % double(tenmat(ehess.G,1));
0175         ehess_G2 = reshape(permute(ehess.G, [2 1 3]), r2, r1*r3); % double(tenmat(ehess.G,2));
0176         ehess_G3 = reshape(permute(ehess.G, [3 1 2]), r3, r1*r2); % double(tenmat(ehess.G,3));
0177         rgrad_G1 = reshape(rgrad.G, r1, r2*r3); % double(tenmat(rgrad.G,1));
0178         rgrad_G2 = reshape(permute(rgrad.G, [2 1 3]), r2, r1*r3); % double(tenmat(rgrad.G,2));
0179         rgrad_G3 = reshape(permute(rgrad.G, [3 1 2]), r3, r1*r2); % double(tenmat(rgrad.G,3));
0180         
0181         ASU1dot = 2*symm((2*symm(X.G1*eta_G1')*(egrad_G1*X.G1')) + X.G1G1t*(ehess_G1*X.G1' + egrad_G1*eta_G1')) - 4*symm(symm(eta_G1*X.G1')*BU1);
0182         ASU2dot = 2*symm((2*symm(X.G2*eta_G2')*(egrad_G2*X.G2')) + X.G2G2t*(ehess_G2*X.G2' + egrad_G2*eta_G2')) - 4*symm(symm(eta_G2*X.G2')*BU2);
0183         ASU3dot = 2*symm((2*symm(X.G3*eta_G3')*(egrad_G3*X.G3')) + X.G3G3t*(ehess_G3*X.G3' + egrad_G3*eta_G3')) - 4*symm(symm(eta_G3*X.G3')*BU3);
0184         
0185         
0186         SSU1dot = X.G1G1t;
0187         SSU2dot = X.G2G2t;
0188         SSU3dot = X.G3G3t;
0189         BU1dot = lyapunov_symmetric(SSU1dot, ASU1dot);
0190         BU2dot = lyapunov_symmetric(SSU2dot, ASU2dot);
0191         BU3dot = lyapunov_symmetric(SSU3dot, ASU3dot);
0192         
0193         
0194         Hess.U1 = (ehess.U1 - eta.U1*BU1 - X.U1*BU1dot - 2*rgrad.U1*symm(eta_G1*X.G1'))/X.G1G1t;
0195         Hess.U2 = (ehess.U2 - eta.U2*BU2 - X.U2*BU2dot - 2*rgrad.U2*symm(eta_G2*X.G2'))/X.G2G2t;
0196         Hess.U3 = (ehess.U3 - eta.U3*BU3 - X.U3*BU3dot - 2*rgrad.U3*symm(eta_G3*X.G3'))/X.G3G3t;
0197         Hess.G = ehess.G;
0198         
0199         
0200         
0201         % BM: we need a correction factor for the non-constant metric
0202         % The correction factor owes itself to the Koszul formula.
0203         % This is the Riemannian connection in the Euclidean space with the
0204         % scaled metric.
0205         
0206         
0207         Hess.U1 = Hess.U1 + (eta.U1*symm(rgrad_G1*X.G1') + rgrad.U1*symm(eta_G1*X.G1'))/X.G1G1t;
0208         Hess.U2 = Hess.U2 + (eta.U2*symm(rgrad_G2*X.G2') + rgrad.U2*symm(eta_G2*X.G2'))/X.G2G2t;
0209         Hess.U3 = Hess.U3 + (eta.U3*symm(rgrad_G3*X.G3') + rgrad.U3*symm(eta_G3*X.G3'))/X.G3G3t;
0210         Hess.G = Hess.G  - permute(reshape(symm(rgrad.U1'*eta.U1)*X.G1,r1,r2,r3), [1 2 3]) ...
0211             - permute(reshape(symm(rgrad.U2'*eta.U2)*X.G2,r2,r1,r3), [2 1 3]) ...
0212             - permute(reshape(symm(rgrad.U3'*eta.U3)*X.G3,r3,r1,r2), [2 3 1]);
0213         
0214         % The Riemannian connection on the quotient space is the
0215         % projection on the tangent space of the total space and then onto the horizontal
0216         % space. This is accomplished with the following operation.
0217         
0218         Hess = M.proj(X, Hess);
0219         
0220         
0221     end
0222     
0223     
0224     
0225     
0226     M.proj = @projection;
0227     function etaproj = projection(X, eta)
0228         X = prepare(X); % Reuse already computed terms
0229         
0230         % First, projection onto tangent space of total space
0231         SSU1 = X.G1G1t;
0232         ASU1 = 2*symm(X.G1G1t*(X.U1'*eta.U1)*X.G1G1t);
0233         BU1 = lyapunov_symmetric(SSU1, ASU1);
0234         eta.U1 = eta.U1 - X.U1*(BU1/X.G1G1t);
0235         
0236         SSU2 = X.G2G2t;
0237         ASU2 = 2*symm(X.G2G2t*(X.U2'*eta.U2)*X.G2G2t);
0238         BU2 = lyapunov_symmetric(SSU2, ASU2);
0239         eta.U2 = eta.U2 - X.U2*(BU2/X.G2G2t);
0240         
0241         SSU3 = X.G3G3t;
0242         ASU3 = 2*symm(X.G3G3t*(X.U3'*eta.U3)*X.G3G3t);
0243         BU3 = lyapunov_symmetric(SSU3, ASU3);
0244         eta.U3 = eta.U3 - X.U3*(BU3/X.G3G3t);
0245         
0246 
0247         eta_G1 = reshape(eta.G, r1, r2*r3); 
0248         eta_G2 = reshape(permute(eta.G, [2 1 3]), r2, r1*r3); 
0249         eta_G3 = reshape(permute(eta.G, [3 1 2]), r3, r1*r2);
0250         
0251         
0252         % Project onto the horizontal space.
0253         PU1 = skew((X.U1'*eta.U1)*X.G1G1t) + skew(X.G1*eta_G1');
0254         PU2 = skew((X.U2'*eta.U2)*X.G2G2t) + skew(X.G2*eta_G2');
0255         PU3 = skew((X.U3'*eta.U3)*X.G3G3t) + skew(X.G3*eta_G3');
0256         
0257         % Calculate Omega1, Omega2, Omega3 that are required in finding the
0258         % horizontal component.
0259         % We use the Matlab's pcg function to solve the system efficiently.
0260         % We exploit the structure by designing a good preconditioner as well.
0261         % The preconditioner takes the block positive definite part of the
0262         % linear system.
0263         
0264         % Options for PCG
0265         tol_omegax_pcg = 1e-6; % BM: standard tolerance as suggested in PCG.
0266         max_iterations_pcg = 15;% BM: fix this to 15 for simulations. In practice, it requires 7 to 10 iteraions.
0267         
0268         % Preconditioner for PCG
0269         M1 = kron(speyer1,SSU1) + kron(SSU1, speyer1);
0270         M2 = kron(speyer2,SSU2) + kron(SSU2, speyer2);
0271         M3 = kron(speyer3,SSU3) + kron(SSU3, speyer3);
0272         
0273         Mprecon_pcg = sparse(zeros(r1^2 + r2^2 + r3^2));
0274         Mprecon_pcg(1 : r1^2, 1 : r1^2 ) = M1;
0275         Mprecon_pcg(1 + r1^2 : r1^2 + r2^2, 1 + r1^2 : r1^2 + r2^2) = M2;
0276         Mprecon_pcg(1 + r1^2 + r2^2 : end, 1 + r1^2 + r2^2 : end) = M3;
0277         
0278         % Call PCG
0279         [Omegaxsol, unused] = pcg(@compute_residual, [PU1(:); PU2(:); PU3(:)],  tol_omegax_pcg, max_iterations_pcg, Mprecon_pcg); %#ok<ASGLU>
0280         
0281         Omega1 = reshape(Omegaxsol(1:r1^2), r1, r1);
0282         Omega2 = reshape(Omegaxsol(1 + r1^2 : r1^2 + r2^2), r2, r2);
0283         Omega3 = reshape(Omegaxsol(1 + r1^2 + r2^2 : end), r3, r3);
0284             
0285         function AOmegax = compute_residual(Omegax)
0286             Omegax1 = reshape(Omegax(1:r1^2), r1, r1);
0287             Omegax2 = reshape(Omegax(1 + r1^2 : r1^2 + r2^2), r2, r2);
0288             Omegax3 = reshape(Omegax(1 + r1^2 + r2^2 : end), r3, r3);
0289             
0290             OffsetU1 = X.G1*((kron(speyer3,Omegax2) + kron(Omegax3, speyer2))*X.G1');
0291             OffsetU2 = X.G2*((kron(speyer3,Omegax1) + kron(Omegax3, speyer1))*X.G2');
0292             OffsetU3 = X.G3*((kron(speyer2,Omegax1) + kron(Omegax2, speyer1))*X.G3');
0293             
0294             residual1 = Omegax1*SSU1 + SSU1*Omegax1 - OffsetU1;
0295             residual2 = Omegax2*SSU2 + SSU2*Omegax2 - OffsetU2;
0296             residual3 = Omegax3*SSU3 + SSU3*Omegax3 - OffsetU3;
0297             
0298             AOmegax = [residual1(:); residual2(:); residual3(:)];
0299         end
0300         
0301         
0302         % Calculate projection along U1, U2, and U3
0303         etaproj.U1 = eta.U1 - (X.U1*Omega1);
0304         etaproj.U2 = eta.U2 - (X.U2*Omega2);
0305         etaproj.U3 = eta.U3 - (X.U3*Omega3);
0306         
0307         % Calculate projection algong G
0308         GOmega1 = reshape(Omega1*X.G1, r1, r2, r3);
0309         GOmega2 = permute(reshape(Omega2*X.G2, r2, r1, r3), [2 1 3]);
0310         GOmega3 = permute(reshape(Omega3*X.G3, r3, r1, r2), [2 3 1]); 
0311         etaproj.G = eta.G -(-(GOmega1+GOmega2+GOmega3));
0312         
0313     end
0314     
0315     
0316     
0317     M.tangent = M.proj;
0318     M.tangent2ambient = @(X, eta) eta;
0319     
0320     M.retr = @retraction;
0321     function Y = retraction(X, eta, t)
0322         if nargin < 3
0323             t = 1.0;
0324         end
0325         
0326         Y.G = (X.G + t*eta.G);
0327         Y.U1 = uf((X.U1 + t*eta.U1)); % U factor of Polar factorization
0328         Y.U2 = uf((X.U2 + t*eta.U2));
0329         Y.U3 = uf((X.U3 + t*eta.U3));
0330         
0331         Y = prepare(Y);
0332     end
0333     
0334     
0335     M.hash = @(X) ['z' hashmd5([sum(X.U1(:)) ; sum(X.U2(:)); sum(X.U3(:)); sum(X.G(:)) ])]; % Efficient, suggested by Bart Vandereycken.
0336     % M.hash = @(X) ['z' hashmd5([X.U1(:); X.U2(:); X.U3(:); X.G(:)])];
0337     
0338     M.rand = @random;
0339     function X = random()
0340         %         % Random generator on the total space
0341         %         % Factors U1, U2, and U3 are on Stiefel manifolds, hence we reuse
0342         %         % their random generator.
0343         %         stiefell = stiefelfactory(n1, r1);
0344         %         stiefelm = stiefelfactory(n2, r2);
0345         %         stiefeln = stiefelfactory(n3, r3);
0346         %
0347         %         X.U1 = stiefell.rand();
0348         %         X.U2 = stiefelm.rand();
0349         %         X.U3 = stiefeln.rand();
0350         %
0351         %         % Random initialization: generalization of randn(r1, r1 = r2) in the
0352         %         % matrix case.
0353         %         X.G = randn(r1,r2,r3);
0354         
0355         
0356         %  Random generator on the fixed-rank space from a uniform distribution on [0, 1].
0357         [U1, R1] = qr(rand(n1, r1), 0);
0358         [U2, R2] = qr(rand(n2, r2), 0);
0359         [U3, R3] = qr(rand(n3, r3), 0);
0360         C  = rand(r1, r2, r3);
0361         
0362         C1 = reshape(C, r1, r2*r3);
0363         CR1 = reshape(R1*C1, r1, r2, r3); % Multplication by R1
0364         
0365         C2 = reshape(permute(CR1, [2 1 3]), r2, r1*r3);
0366         CR1R2 = permute(reshape(R2*C2, r2, r1, r3), [2 1 3]); % Multplication by R2
0367         
0368         C3 = reshape(permute(CR1R2, [3 1 2]), r3, r1*r2);
0369         CR1R2R3 = permute(reshape(R3*C3, r3, r1, r2), [2 3 1]); % Multplication by R3
0370         
0371         X.U1 = U1;
0372         X.U2 = U2;
0373         X.U3 = U3;
0374         X.G = CR1R2R3;
0375     
0376         
0377         % Compute some terms that are used subsequently.
0378         X = prepare(X);
0379         
0380     end
0381     
0382     M.randvec = @randomvec;
0383     function eta = randomvec(X)
0384         % A random vector on the horizontal space
0385         eta.U1 = randn(n1, r1);
0386         eta.U2 = randn(n2, r2);
0387         eta.U3 = randn(n3, r3);
0388         eta.G = randn(r1, r2, r3);
0389         eta = projection(X, eta);
0390         nrm = M.norm(X, eta);
0391         eta.U1 = eta.U1 / nrm;
0392         eta.U2 = eta.U2 / nrm;
0393         eta.U3 = eta.U3 / nrm;
0394         eta.G = eta.G / nrm;
0395     end
0396     
0397     M.lincomb = @lincomb;
0398     
0399     M.zerovec = @(X) struct('U1', zeros(n1, r1), 'U2', zeros(n2, r2), ...
0400         'U3', zeros(n3, r3), 'G', zeros(r1, r2, r3));
0401     
0402     M.transp = @(x1, x2, d) projection(x2, d);
0403     
0404     % vec and mat are not isometries, because of the scaled metric.
0405     M.vec = @(X, U1) [U1.U1(:); U1.U2(:); U1.U3(:); U1.G(:)];
0406     M.mat = @(X, u) struct ...
0407         ('U1', reshape(u(1  : n1*r1), n1, r1), ...
0408         'U2', reshape(u(n1*r1 + 1 : n1*r1 + n2*r2), n2, r2), ...
0409         'U3', reshape(u(n1*r1 + n2*r2 + 1 : n1*r1 + n2*r2 + n3*r3), n3, r3), ...
0410         'G', reshape(u(n1*r1 + n2*r2 + n3*r3 + 1 : end), r1, r2, r3));
0411     M.vecmatareisometries = @() false;
0412     
0413 end
0414 
0415 % Linear combination of tangent vectors
0416 function d = lincomb(X, a1, d1, a2, d2) %#ok<INUSL>
0417     
0418     if nargin == 3
0419         d.U1 = a1*d1.U1;
0420         d.U2 = a1*d1.U2;
0421         d.U3 = a1*d1.U3;
0422         d.G = a1*d1.G;
0423     elseif nargin == 5
0424         d.U1 = a1*d1.U1 + a2*d2.U1;
0425         d.U2 = a1*d1.U2 + a2*d2.U2;
0426         d.U3 = a1*d1.U3 + a2*d2.U3;
0427         d.G = a1*d1.G + a2*d2.G;
0428     else
0429         error('Bad use of fixedrankfactory_tucker_preconditioned.lincomb.');
0430     end
0431     
0432 end
0433 
0434 function U = uf(A) % U factor of Polar factorization of a tall matrix A.
0435     [L, unused, R] = svd(A, 0); %#ok
0436     U = L*R';
0437 end
0438 
0439 
0440

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