block_eigenvalue

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
         m-file can be downloaded from 
               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.

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
0028     %         m-file can be downloaded from
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
0039     %   BSD 2-clause license, see LICENSE.txt
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