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.

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:
• hashmd5 Computes the MD5 hash of input data.
• lincomb Computes a linear combination of tangent vectors in the Manopt framework.
• lyapunov_symmetric Solves AX + XA = C when A = A', as a pseudo-inverse.
This function is called by:

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

Generated on Sun 05-Sep-2021 17:57:00 by m2html © 2005