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

 See also fixedranktensorembeddedfactory

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

Generated on Tue 19-May-2020 18:46:12 by m2html © 2005