Home > manopt > manifolds > ttfixedrank > TTeMPS_1.1 > algorithms > eigenvalue > block_eigenvalue.m

# block_eigenvalue

## PURPOSE

BLOCK_EIGENVALUE Calculate p smallest eigenvalues of a TTeMPS operator

## SYNOPSIS

function [X, C, evalue, residuums, micro_res, objective, elapsed_time] = block_eigenvalue(A, p, rr, opts)

## DESCRIPTION

```BLOCK_EIGENVALUE Calculate p smallest eigenvalues of a TTeMPS operator

[X, C, evalue, residuums, micro_res, objective, elapsed_time] = block_eigenvalue(A, P, RR, OPTS)
performs a block-eigenvalue optimization scheme to compute the p smallest eigenvalues of A
using the algorithm described in [1].

RR defines the starting rank and should usually be set to ones(1,d) where d is the dimensionality.
If p == 1, the algorithm equals a standard ALS-procedure for the eigenvalue problem, which is NOT
rank adaptive. Hence, in this case RR should be taken to be the expected rank of the solution or,
if unknown, the highest affordable rank.

Specify the wanted options using the OPTS argument. All options have
default values (denoted in brackets). Hence, you only have to specify those you want
to change.
The OPTS struct is organized as follows:
OPTS.maxiter        Maximum number of full sweeps to perform        [3]
OPTS.maxrank        Maximum rank during the iteration               [40]
OPTS.tol            Tolerance for the shift from one core
to the next                                     [1e-8]
OPTS.tolLOBPCG      Tolerance for inner LOBPCG solver               [1e-6]
OPTS.maxiterLOBPCG  Max. number of iterations for LOBPCG            [2000]
OPTS.verbose        Show iteration information (true/false)         [true]
OPTS.precInner      Precondition the LOBPCG (STRONGLY RECOMMENDED!) [true]

NOTE: To run block_eigenvalue, Knyazev's implementation of LOBPCG is required. The corresponding
http://www.mathworks.com/matlabcentral/fileexchange/48-lobpcg-m

References:
[1] S.V. Dolgov, B.N. Khoromskij, I.V. Oseledets, D.V. Savostyanov
Computation of extreme eigenvalues in higher dimensions using block tensor train format
Computer Physics Communications 185 (2014) 1207-1216.```

## CROSS-REFERENCE INFORMATION

This function calls:
• disp DISP Display TT/MPS tensor.
• innerprod INNERPROD Inner product between two TT/MPS tensors.
• norm NORM Norm of a TT/MPS tensor.
• orthogonalize ORTHOGONALIZE Orthogonalize tensor.
• disp DISP Display TT/MPS block-mu tensor.
• innerprod INNERPROD Inner product between two TT/MPS tensors.
• norm NORM Norm of a TT/MPS block-mu tensor.
• orthogonalize ORTHOGONALIZE Orthogonalize TT/MPS Block-mu tensor.
• apply APPLY Application of TT/MPS operator to a TT/MPS tensor
• disp DISP Display TT/MPS operator.
• apply APPLY Application of TT/MPS Laplace-like operator to a TT/MPS tensor
• constr_precond_inner TTeMPS Toolbox.
• disp DISP Display TT/MPS operator.
• TTeMPS_rand TTEMPS_RAND Create random TTeMPS tensor
• tensorprod_ttemps TENSORPROD_TTEMPS Tensor-times-Matrix product.
• trunc_singular REL_TRUNC_SINGULAR Helper routine to truncate singular values
• orthogonalize Orthonormalizes a basis of tangent vectors in the Manopt framework.
This function is called by:

## SOURCE CODE

```0001 function [X, C, evalue, residuums, micro_res, objective, elapsed_time] = block_eigenvalue(A, p, rr, opts)
0002     %BLOCK_EIGENVALUE Calculate p smallest eigenvalues of a TTeMPS operator
0003     %
0004     %   [X, C, evalue, residuums, micro_res, objective, elapsed_time] = block_eigenvalue(A, P, RR, OPTS)
0005     %       performs a block-eigenvalue optimization scheme to compute the p smallest eigenvalues of A
0006     %       using the algorithm described in [1].
0007     %
0008     %   RR defines the starting rank and should usually be set to ones(1,d) where d is the dimensionality.
0009     %   If p == 1, the algorithm equals a standard ALS-procedure for the eigenvalue problem, which is NOT
0010     %   rank adaptive. Hence, in this case RR should be taken to be the expected rank of the solution or,
0011     %   if unknown, the highest affordable rank.
0012     %
0013     %   Specify the wanted options using the OPTS argument. All options have
0014     %   default values (denoted in brackets). Hence, you only have to specify those you want
0015     %   to change.
0016     %   The OPTS struct is organized as follows:
0017     %       OPTS.maxiter        Maximum number of full sweeps to perform        [3]
0018     %       OPTS.maxrank        Maximum rank during the iteration               [40]
0019     %       OPTS.tol            Tolerance for the shift from one core
0020     %                           to the next                                     [1e-8]
0021     %       OPTS.tolLOBPCG      Tolerance for inner LOBPCG solver               [1e-6]
0022     %       OPTS.maxiterLOBPCG  Max. number of iterations for LOBPCG            [2000]
0023     %       OPTS.verbose        Show iteration information (true/false)         [true]
0024     %       OPTS.precInner      Precondition the LOBPCG (STRONGLY RECOMMENDED!) [true]
0025     %
0026     %
0027     %   NOTE: To run block_eigenvalue, Knyazev's implementation of LOBPCG is required. The corresponding
0029     %               http://www.mathworks.com/matlabcentral/fileexchange/48-lobpcg-m
0030     %
0031     %    References:
0032     %    [1] S.V. Dolgov, B.N. Khoromskij, I.V. Oseledets, D.V. Savostyanov
0033     %        Computation of extreme eigenvalues in higher dimensions using block tensor train format
0034     %        Computer Physics Communications 185 (2014) 1207-1216.
0035
0036     %   TTeMPS Toolbox.
0037     %   Michael Steinlechner, 2013-2014
0038     %   Questions and contact: michael.steinlechner@epfl.ch
0040
0041
0042 nn = A.size_col;
0043 d = A.order;
0044
0045 if ~isfield( opts, 'maxiter');       opts.maxiter = 3;           end
0046 if ~isfield( opts, 'maxrank');       opts.maxrank = 40;          end
0047 if ~isfield( opts, 'tol');           opts.tol = 1e-8;            end
0048 if ~isfield( opts, 'tolLOBPCG');     opts.tolLOBPCG = 1e-6;      end
0049 if ~isfield( opts, 'maxiterLOBPCG'); opts.maxiterLOBPCG = 2000;  end
0050 if ~isfield( opts, 'verbose');       opts.verbose = 1;           end
0051 if ~isfield( opts, 'precInner');     opts.precInner = true;      end
0052 if ~isfield( opts, 'X0');            useX0 = 0; else useX0 = 1;  end
0053
0054 if p == 1
0055     tolLOBPCGmod = opts.tolLOBPCG;
0056 else
0057     tolLOBPCGmod = 0.00001;
0058 end
0059
0060
0061 if ~useX0
0062     X = TTeMPS_rand( rr, nn );
0063     % Left-orthogonalize the tensor:
0064     X = orthogonalize( X, X.order );
0065
0066     X.U{d} = rand( X.rank(d), X.size(d), X.rank(d+1), p );
0067     C = cell( 1, p );
0068     C{1} = X.U{d}(:,:,:,1);
0069     for i = 2:p
0070         C{i} = X.U{d}(:,:,:,i);
0071     end
0072
0073 else
0074     X = opts.X0;
0075     X.U{d} = rand( X.rank(d), X.size(d), X.rank(d+1), p );
0076     C = cell( 1, p );
0077     C{1} = X.U{d}(:,:,:,1);
0078     for i = 2:p
0079         C{i} = X.U{d}(:,:,:,i);
0080     end;
0081 end
0082
0083 evalue = [];
0084 residuums = zeros(p,2*opts.maxiter);
0085 micro_res = [];
0086 resi_norm = zeros(p,1);
0087 objective = [];
0088 tic;
0089 elapsed_time = [];
0090
0091
0092 for i = 1:opts.maxiter
0093
0094     disp(sprintf('Iteration %i:', i));
0095     % right-to-left sweep
0096     fprintf(1, 'RIGHT-TO-LEFT SWEEP. ---------------\n')
0097
0098     for mu = d:-1:2
0099         sz = [X.rank(mu), X.size(mu), X.rank(mu+1)];
0100         % shows itertaion information
0101         if opts.verbose
0102             fprintf(1,'Current core: %i. Current iterate (first eigenvalue): \n', mu);
0103             disp(X);
0104             fprintf(1,'Running LOBPCG: system size %i, tol %g ... ', prod(sz), max(opts.tolLOBPCG, min( tolLOBPCGmod, sqrt(sum(residuums(:,2*(i-1)+1).^2))/sqrt(p)*tolLOBPCGmod )))
0105         else
0106             fprintf(1,'%i  ',mu);
0107         end
0108
0109         [left, right] = Afun_prepare( A, X, mu );
0110
0111         if opts.precInner
0112             if prod(sz) < 5*p;
0113                 Amat = zeros( prod(sz) );
0114                 for i_mat = 1:prod(sz)
0115                     Amat(:,i_mat) = Afun_block_optim( A, [zeros(i_mat-1,1); 1; zeros(prod(sz)-i_mat,1)], sz, left, right, mu);
0116                 end
0117                 [vmat,emat] = eig(Amat);
0118                 [emat,sortind] = sort(diag(emat),'ascend');
0119                 vmat = vmat(:, sortind);
0120                 L = emat(1:p);
0121                 V = vmat(:, 1:p);
0122
0123                 disp('Reduced system too small for LOBPCG, using eig instead');
0124                 disp(['Current Eigenvalue approx: ', num2str( L(1:p)')]);
0125
0126
0127             else
0128                 expB = constr_precond_inner( A, X, mu );
0129                 [V,L,failureFlag,Lhist ] = lobpcg( rand( prod(sz), p), ...
0130                                                    @(y) Afun_block_optim( A, y, sz, left, right, mu), [], ...
0131                                                    @(y) apply_local_precond( A, y, sz, expB), ...
0132                                                    max(opts.tolLOBPCG, min( tolLOBPCGmod, sqrt(sum(residuums(:,2*(i-1)+1).^2))/sqrt(p)*tolLOBPCGmod )), ...
0133                                                    opts.maxiterLOBPCG, 0);
0134                 if opts.verbose
0135                     if failureFlag
0136                         fprintf(1,'NOT CONVERGED within %i steps!\n', opts.maxiterLOBPCG)
0137                     else
0138                         fprintf(1,'converged after %i steps!\n', size(Lhist,2));
0139                     end
0140                     disp(['Current Eigenvalue approx: ', num2str( L(1:p)')]);
0141                 end
0142             end
0143         else
0144             if prod(sz) < 25;
0145                 Amat = zeros( prod(sz) );
0146                 for i_mat = 1:prod(sz)
0147                     Amat(:,i_mat) = Afun_block_optim( A, [zeros(i_mat-1,1); 1; zeros(prod(sz)-i_mat,1)], sz, left, right, mu);
0148                 end
0149                 [vmat,emat] = eig(Amat);
0150                 [emat,sortind] = sort(diag(emat),'ascend');
0151                 vmat = vmat(:, sortind);
0152                 L = emat(1:p);
0153                 V = vmat(:, 1:p);
0154
0155                 disp('Condition number!')
0156                 cond(Amat)
0157
0158                 disp('Reduced system too small for LOBPCG, using eig instead');
0159                 disp(['Current Eigenvalue approx: ', num2str( L(1:p)')]);
0160
0161             else
0162                 X0 = rand( prod(sz), p);
0163                 X0 = orth(X0);
0164                 [V,L,failureFlag,Lhist ] = lobpcg( X0, ...
0165                                                @(y) Afun_block_optim( A, y, sz, left, right, mu), ...
0166                                                opts.tolLOBPCG, opts.maxiterLOBPCG, 0);
0167             end
0168         end
0169
0170         evalue = [evalue, L(1:p)];
0171         objective = [objective, sum(evalue(:,end))];
0172         elapsed_time = [elapsed_time, toc];
0173
0174         X.U{mu} = reshape( V, [sz, p] );
0175         lamX = X;
0176         lamX.U{mu} = repmat( reshape(L, [1 1 1 p]), [X.rank(mu), X.size(mu), X.rank(mu+1), 1]).*lamX.U{mu};
0177         res_new = apply(A, X) - lamX
0178
0179         for j = 1:p
0180             X.U{mu} = reshape( V(:,j), sz );
0181             res = apply(A, X) - L(j)*X;
0182             resi_norm(j) = norm(res);
0183             residuums(j,2*(i-1)+1) = resi_norm(j);
0184
0185         end
0186
0187
0188         micro_res = [micro_res, resi_norm];
0189
0190         % split new core
0191         V = reshape( V, [sz, p] );
0192         V = permute( V, [1, 4, 2, 3] );
0193         V = reshape( V, [sz(1)*p, sz(2)*sz(3)] );
0194
0195         [U,S,V] = svd( V, 'econ' );
0196
0197         if p == 1
0198             s = length(diag(S));
0199         else
0200             s = trunc_singular( diag(S), opts.tol, true, opts.maxrank );
0201         end
0202
0203         V = V(:,1:s)';
0204         X.U{mu} = reshape( V, [s, sz(2), sz(3)] );
0205
0206         W = U(:,1:s)*S(1:s,1:s);
0207         W = reshape( W, [sz(1), p, s]);
0208         W = permute( W, [1, 3, 2]);
0209         for k = 1:p
0210             C{k} = tensorprod_ttemps( X.U{mu-1}, W(:,:,k)', 3);
0211         end
0212
0213         X.U{mu-1} = C{1};
0214
0215         if opts.verbose
0216             disp( ['Augmented system of size (', num2str( [sz(1)*p, sz(2)*sz(3)]), '). Cut-off tol: ', num2str(opts.tol) ])
0217             disp( sprintf( 'Number of SVs: %i. Truncated to: %i => %g %%', length(diag(S)), s, s/length(diag(S))*100))
0218             disp(' ')
0219         end
0220     end
0221
0222     % calculate current residuum
0223
0224     fprintf(1, '---------------    finshed sweep.    ---------------\n')
0225     disp(['Current residuum: ', num2str(residuums(:,2*(i-1)+1).')])
0226     disp(' ')
0227     disp(' ')
0228     fprintf(1, '--------------- LEFT-TO-RIGHT SWEEP. ---------------\n')
0229     % left-to-right sweep
0230     for mu = 1:d-1
0231         sz = [X.rank(mu), X.size(mu), X.rank(mu+1)];
0232
0233         if opts.verbose
0234             fprintf(1,'Current core: %i. Current iterate (first eigenvalue): \n', mu);
0235             disp(X);
0236             fprintf(1,'Running LOBPCG: system size %i, tol %g ... ', prod(sz), max(opts.tolLOBPCG, min( tolLOBPCGmod, sqrt(sum(residuums(:,2*(i-1)+2).^2))/sqrt(p)*tolLOBPCGmod )))
0237         else
0238             fprintf(1,'%i  ',mu);
0239         end
0240
0241         [left, right] = Afun_prepare( A, X, mu );
0242
0243         if opts.precInner
0244             expB = constr_precond_inner( A, X, mu );
0245             X0 = rand( prod(sz), p);
0246             X0 = orth(X0);
0247             [U,L,failureFlag,Lhist ] = lobpcg( X0 , ...
0248                                                @(y) Afun_block_optim( A, y, sz, left, right, mu), [], ...
0249                                                @(y) apply_local_precond( A, y, sz, expB), ...
0250                                                max(opts.tolLOBPCG, min( tolLOBPCGmod, sqrt(sum(residuums(:,2*(i-1)+2).^2))/sqrt(p)*tolLOBPCGmod )), ...
0251                                                opts.maxiterLOBPCG, 0);
0252         else
0253             if prod(sz) < 25;
0254                 Amat = zeros( prod(sz) );
0255                 for i_mat = 1:prod(sz)
0256                     Amat(:,i_mat) = Afun_block_optim( A, [zeros(i_mat-1,1); 1; zeros(prod(sz)-i_mat,1)], sz, left, right, mu);
0257                 end
0258                 [vmat,emat] = eig(Amat);
0259                 [emat,sortind] = sort(diag(emat),'ascend');
0260                 vmat = vmat(:, sortind);
0261                 L = emat(1:p);
0262                 U = vmat(:, 1:p);
0263
0264                 disp('Reduced system too small for LOBPCG, using eig instead');
0265                 disp(['Current Eigenvalue approx: ', num2str( L(1:p)')]);
0266             else
0267                 X0 = rand( prod(sz), p);
0268                 X0 = orth(X0);
0269                 [U,L,failureFlag,Lhist ] = lobpcg( X0 , ...
0270                                                 @(y) Afun_block_optim( A, y, sz, left, right, mu), ...
0271                                                 opts.tolLOBPCG, opts.maxiterLOBPCG, 0);
0272             end
0273         end
0274
0275         if opts.verbose
0276             if failureFlag
0277                 fprintf(1,'NOT CONVERGED within %i steps!\n', opts.maxiterLOBPCG)
0278             else
0279                 fprintf(1,'converged after %i steps!\n', size(Lhist,2));
0280             end
0281             disp(['Current Eigenvalue approx: ', num2str( L(1:p)')]);
0282         end
0283         evalue = [evalue, L(1:p)];
0284         objective = [objective, sum(evalue(:,end))];
0285         elapsed_time = [elapsed_time, toc];
0286
0287         for j = 1:p
0288             X.U{mu} = reshape( U(:,j), sz );
0289             res = apply(A, X) - L(j)*X;
0290             resi_norm(j) = norm(res);
0291             residuums(j,2*(i-1)+2) = resi_norm(j);
0292         end
0293         micro_res = [micro_res, resi_norm];
0294
0295         % split new core
0296         U = reshape( U, [sz, p] );
0297         U = permute( U, [1, 2, 4, 3] );
0298         U = reshape( U, [sz(1)*sz(2), p*sz(3)] );
0299
0300         [U,S,V] = svd( U, 'econ' );
0301         if p == 1
0302             s = length(diag(S));
0303         else
0304             s = trunc_singular( diag(S), opts.tol, true, opts.maxrank );
0305         end
0306
0307         U = U(:,1:s);
0308         X.U{mu} = reshape( U, [sz(1), sz(2), s] );
0309         W = S(1:s,1:s)*V(:,1:s)';
0310         W = reshape( W, [s, p, sz(3)]);
0311         W = permute( W, [1, 3, 2]);
0312
0313         for k = 1:p
0314             C{k} = tensorprod_ttemps( X.U{mu+1}, W(:,:,k), 1);
0315         end
0316
0317         X.U{mu+1} = C{1};
0318
0319         if opts.verbose
0320             disp( ['Augmented system of size (', num2str( [sz(1)*sz(2), p*sz(3)]), '). Cut-off tol: ', num2str(opts.tol) ])
0321             disp( sprintf( 'Number of SVs: %i. Truncated to: %i => %g %%', length(diag(S)), s, s/length(diag(S))*100))
0322             disp(' ')
0323         end
0324     end
0325
0326     fprintf(1, '---------------    finshed sweep.    ---------------\n')
0327     disp(['Current residuum: ', num2str(residuums(:,2*(i-1)+2).')])
0328     disp(' ')
0329     disp(' ')
0330 end
0331
0332 evs = zeros(p,1);
0333 % Evaluate rayleigh quotient for the p TT/MPS tensors
0334 for i=1:p
0335     evec = X;
0336     evec.U{d} = C{i};
0337     evs(i) = innerprod( evec, apply(A, evec));
0338 end
0339 evalue = [evalue, evs];
0340
0341 end
0342
0343 function [left, right] = Afun_prepare( A, x, idx )
0344     y = A.apply(x);
0345     if idx == 1
0346         right = innerprod( x, y, 'RL', idx+1 );
0347         left = [];
0348     elseif idx == x.order
0349         left = innerprod( x, y, 'LR', idx-1 );
0350         right = [];
0351     else
0352         left = innerprod( x, y, 'LR', idx-1 );
0353         right = innerprod( x, y, 'RL', idx+1 );
0354     end
0355 end
0356
0357 function res = Afun_block_optim( A, U, sz, left, right, mu )
0358
0359     p = size(U, 2);
0360
0361     V = reshape( U, [sz, p] );
0362     V = A.apply( V, mu );
0363
0364     if mu == 1
0365         tmp = reshape( permute( V, [3 1 2 4] ), [size(V, 3), sz(1)*sz(2)*p]);
0366         tmp = right * tmp;
0367         tmp = reshape( tmp, [size(right, 1), sz(1), sz(2), p]);
0368         tmp = ipermute( tmp, [3 1 2 4] );
0369     elseif mu == A.order
0370         tmp = reshape( V, [size(V,1), sz(2)*sz(3)*p]);
0371         tmp = left * tmp;
0372         tmp = reshape( tmp, [size(left, 1), sz(2), sz(3), p]);
0373     else
0374         tmp = reshape( permute( V, [3 1 2 4] ), [size(V, 3), size(V, 1)*sz(2)*p]);
0375         tmp = right * tmp;
0376         tmp = reshape( tmp, [size(right, 1), size(V, 1), sz(2), p]);
0377         tmp = ipermute( tmp, [3 1 2 4] );
0378
0379         tmp = reshape( tmp, [size(V, 1), sz(2)*sz(3)*p]);
0380         tmp = left * tmp;
0381         tmp = reshape( tmp, [size(left, 1), sz(2), sz(3), p]);
0382
0383     end
0384
0385     res = reshape( tmp, [prod(sz), p] );
0386
0387
0388 end
0389
0390 function res = apply_local_precond( A, U, sz, expB)
0391
0392     p = size(U, 2);
0393
0394     x = reshape( U, [sz, p] );
0395     res = zeros( [sz, p] );
0396
0397     for i = 1:size( expB, 1)
0398         tmp = reshape( x, [sz(1), sz(2)*sz(3)*p] );
0399         tmp = reshape( expB{1,i}*tmp, [sz(1), sz(2), sz(3), p] );
0400
0401         tmp = reshape( permute( tmp, [2 1 3 4] ), [sz(2), sz(1)*sz(3)*p] );
0402         tmp = ipermute( reshape( expB{2,i}*tmp, [sz(2), sz(1), sz(3), p] ), [2 1 3 4] );
0403
0404         tmp = reshape( permute( tmp, [3 1 2 4] ), [sz(3), sz(1)*sz(2)*p] );
0405         tmp = ipermute( reshape( expB{3,i}*tmp, [sz(3), sz(1), sz(2), p] ), [3 1 2 4] );
0406
0407         res = res + tmp;
0408     end
0409     res = reshape( res, [prod(sz), p] );
0410
0411 end
0412```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005