amen_eigenvalue

AMEN_EIGENVALUE Calculate p smallest eigenvalues of a TTeMPS operator

   [X, C, evalue, residuums, micro_res, objective, elapsed_time] = amen_eigenvalue(A, PREC, P, RR, OPTS)
       performs a rank-adaptive AMEn-type optimization scheme to compute the p smallest eigenvalues of A
       using the algorithm described in [1].

   A preconditioner can be given by passing the PREC argument:
       []      no preconditioner used
       PR      with PR of type TTeMPS_op: PR is applied as a preconditioner
       true    exponential sum preconditioner described in [1] (recommended setting)    

   RR defines the starting rank and should usually be set to ones(1,d) where d is the dimensionality.
   
   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.maxrankRes     Maximum rank of residual (0 = No limit)         [0]
       OPTS.tol            Tolerance for the shift from one core 
                           to the next                                     [1e-8]
       OPTS.tolOP          Tolerance for truncation of residual            [1e-6]
       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 amen_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] D. Kressner, M. Steinlechner, A. Uschmajew:
        Low-rank tensor methods with subspace correction for symmetric eigenvalue problems
        SIAM J. Sci. Comput., 36(5):A2346-A2368, 2014.

0001 function [X, C, evalue, residuums, micro_res, objective, elapsed_time] = amen_eigenvalue(A, prec, p, rr, opts)
0002     %AMEN_EIGENVALUE Calculate p smallest eigenvalues of a TTeMPS operator
0003     %
0004     %   [X, C, evalue, residuums, micro_res, objective, elapsed_time] = amen_eigenvalue(A, PREC, P, RR, OPTS)
0005     %       performs a rank-adaptive AMEn-type optimization scheme to compute the p smallest eigenvalues of A
0006     %       using the algorithm described in [1].
0007     %
0008     %   A preconditioner can be given by passing the PREC argument:
0009     %       []      no preconditioner used
0010     %       PR      with PR of type TTeMPS_op: PR is applied as a preconditioner
0011     %       true    exponential sum preconditioner described in [1] (recommended setting)
0012     %
0013     %   RR defines the starting rank and should usually be set to ones(1,d) where d is the dimensionality.
0014     %
0015     %   Specify the wanted options using the OPTS argument. All options have
0016     %   default values (denoted in brackets). Hence, you only have to specify those you want
0017     %   to change.
0018     %   The OPTS struct is organized as follows:
0019     %       OPTS.maxiter        Maximum number of full sweeps to perform        [3]
0020     %       OPTS.maxrank        Maximum rank during the iteration               [40]
0021     %       OPTS.maxrankRes     Maximum rank of residual (0 = No limit)         [0]
0022     %       OPTS.tol            Tolerance for the shift from one core
0023     %                           to the next                                     [1e-8]
0024     %       OPTS.tolOP          Tolerance for truncation of residual            [1e-6]
0025     %       OPTS.tolLOBPCG      Tolerance for inner LOBPCG solver               [1e-6]
0026     %       OPTS.maxiterLOBPCG  Max. number of iterations for LOBPCG            [2000]
0027     %       OPTS.verbose        Show iteration information (true/false)         [true]
0028     %       OPTS.precInner      Precondition the LOBPCG (STRONGLY RECOMMENDED!) [true]
0029     %
0030     %
0031     %   NOTE: To run amen_eigenvalue, Knyazev's implementation of LOBPCG is required. The corresponding
0032     %         m-file can be downloaded from
0033     %               http://www.mathworks.com/matlabcentral/fileexchange/48-lobpcg-m
0034     %
0035     %    References:
0036     %    [1] D. Kressner, M. Steinlechner, A. Uschmajew:
0037     %        Low-rank tensor methods with subspace correction for symmetric eigenvalue problems
0038     %        SIAM J. Sci. Comput., 36(5):A2346-A2368, 2014.
0039 
0040     %   TTeMPS Toolbox.
0041     %   Michael Steinlechner, 2013-2016
0042     %   Questions and contact: michael.steinlechner@epfl.ch
0043     %   BSD 2-clause license, see LICENSE.txt
0044 
0045 
0046 nn = A.size_col;
0047 d = A.order;
0048 
0049 if ~isempty(prec)
0050     if isa(prec,'TTeMPS_op')
0051         precondition = 1;
0052     else
0053         precondition = 2;
0054     end
0055 else
0056     precondition = 0;
0057 end
0058 
0059 if ~isfield( opts, 'maxiter');       opts.maxiter = 3;           end
0060 if ~isfield( opts, 'maxrank');       opts.maxrank = 40;          end
0061 if ~isfield( opts, 'maxrankRes');    opts.maxrankRes = 0;        end
0062 if ~isfield( opts, 'tol');           opts.tol = 1e-8;            end
0063 if ~isfield( opts, 'tolOP');         opts.tolOP = 1e-6;          end
0064 if ~isfield( opts, 'tolLOBPCG');     opts.tolLOBPCG = 1e-6;      end
0065 if ~isfield( opts, 'maxiterLOBPCG'); opts.maxiterLOBPCG = 2000;  end
0066 if ~isfield( opts, 'verbose');       opts.verbose = true;        end
0067 if ~isfield( opts, 'precInner');     opts.precInner = true;     end
0068 
0069 if p == 1
0070     tolLOBPCGmod = opts.tolLOBPCG;
0071 else
0072     tolLOBPCGmod = 0.01;
0073 end
0074 
0075 X = TTeMPS_rand( rr, nn ); 
0076 % Left-orthogonalize the tensor:
0077 X = orthogonalize( X, X.order );
0078 
0079 C = cell( 1, p );
0080 C{1} = X.U{d};
0081 for i = 2:p
0082     C{i} = rand(size(X.U{d}));
0083 end
0084 
0085 X.U{d} = rand( X.rank(d), X.size(d), X.rank(d+1), p );
0086 
0087 evalue = [];
0088 residuums = zeros(p,2*opts.maxiter);
0089 micro_res = [];
0090 resi_norm = zeros(p,1);
0091 objective = [];
0092 tic;
0093 elapsed_time = [];
0094 
0095 for i = 1:opts.maxiter
0096     
0097     disp(sprintf('Iteration %i:', i));
0098     % right-to-left sweep
0099     fprintf(1, 'RIGHT-TO-LEFT SWEEP. ---------------\n')
0100 
0101     for mu = d:-1:2
0102         
0103         sz = [X.rank(mu), X.size(mu), X.rank(mu+1)];
0104     
0105         if opts.verbose
0106             fprintf(1,'Current core: %i. Current iterate (first eigenvalue): \n', mu);
0107             disp(X);
0108             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 ))) %opts.tolLOBPCG)
0109         else
0110             fprintf(1,'%i  ',mu);
0111         end
0112        
0113         [left, right] = Afun_prepare( A, X, mu );
0114         
0115         if opts.precInner
0116             if prod(sz) < 5*p
0117                 Amat = zeros( prod(sz) );
0118                 for i_mat = 1:prod(sz)
0119                     Amat(:,i_mat) = Afun_block_optim( A, [zeros(i_mat-1,1); 1; zeros(prod(sz)-i_mat,1)], sz, left, right, mu);
0120                 end
0121                 [vmat,emat] = eig(Amat);
0122                 [emat,sortind] = sort(diag(emat),'ascend');
0123                 vmat = vmat(:, sortind);
0124                 L = emat(1:p);
0125                 V = vmat(:, 1:p);
0126                 
0127                 disp('Reduced system too small for LOBPCG, using eig instead');
0128                 disp(['Current Eigenvalue approx: ', num2str( L(1:p)')]);
0129 
0130             else
0131                 expB = constr_precond_inner( A, X, mu );
0132                 [V,L,failureFlag,Lhist ] = lobpcg( rand( prod(sz), p), ...
0133                                                    @(y) Afun_block_optim( A, y, sz, left, right, mu), [], ...
0134                                                    @(y) apply_local_precond( y, sz, expB), ...
0135                                                    max(opts.tolLOBPCG, min( tolLOBPCGmod, sqrt(sum(residuums(:,2*(i-1)+1).^2))/sqrt(p)*tolLOBPCGmod )), opts.maxiterLOBPCG, 0);
0136                 if opts.verbose
0137                     if failureFlag
0138                         fprintf(1,'NOT CONVERGED within %i steps!\n', opts.maxiterLOBPCG)
0139                     else
0140                         fprintf(1,'converged after %i steps!\n', size(Lhist,2));
0141                     end
0142                     disp(['Current Eigenvalue approx: ', num2str( L(1:p)')]);
0143                 end
0144             end
0145         else
0146             [V,L,failureFlag,Lhist ] = lobpcg( rand( prod(sz), p), ...
0147                                                @(y) Afun_block_optim( A, y, sz, left, right, mu), ...
0148                                                opts.tolLOBPCG, opts.maxiterLOBPCG, 0);
0149         end
0150         
0151         evalue = [evalue, L(1:p)];
0152         objective = [objective, sum(evalue(:,end))];
0153         elapsed_time = [elapsed_time, toc];
0154 
0155         X.U{mu} = reshape( V, [sz, p] );
0156         lamX = X;
0157         lamX.U{mu} = repmat( reshape(L, [1 1 1 p]), [X.rank(mu), X.size(mu), X.rank(mu+1), 1]).*lamX.U{mu};
0158         
0159         res = apply(A, X) - lamX;
0160         
0161         if opts.maxrankRes ~= 0
0162             res_sz = [size(res.U{mu},1), size(res.U{mu},2), size(res.U{mu},3), size(res.U{mu},4)];
0163             res.U{mu} = reshape( permute( res.U{mu}, [1 2 4 3]), [ res_sz(1), res_sz(2)*p, res_sz(3)]);
0164             res = truncate( res, [1 opts.maxrankRes*ones(1,d-1) 1] );
0165             res.U{mu} = ipermute( reshape( res.U{mu}, [res.rank(mu), res_sz(2), p, res.rank(mu+1)] ), [1 2 4 3]);
0166         end
0167 
0168         % A not so nice way to calculate the residual norm
0169         for j = 1:p
0170             tmp = res;
0171             tmp.U{mu} = tmp.U{mu}(:,:,:,j);
0172             resi_norm(j) = norm(tmp);
0173             residuums(j,2*(i-1)+1) = resi_norm(j);
0174         end
0175         micro_res = [micro_res, resi_norm];
0176         
0177         if precondition == 1
0178             res = apply(prec, res);
0179         end
0180         
0181         res = contract( X, res, [mu-1, mu]);
0182         res_combined = unfold(res{1},'left') * reshape(res{2}, [size(res{2},1), size(res{2},2)*size(res{2},3)*p]);
0183 
0184         if precondition == 2
0185             expBglobal = constr_precond_outer( A, X, mu-1, mu );
0186             res_size = [size(res{1},1), size(res{1},2), size(res{2},2), size(res{2},3), p];
0187             res_combined = reshape( res_combined, res_size );
0188             res_combined = apply_global_precond( res_combined, res_size, expBglobal );
0189         end
0190 
0191         res_combined = reshape( res_combined, [size(res{1},1)*size(res{1},2), size(res{2},2)*size(res{2},3)*p]); 
0192 
0193         [uu,ss,vv] = svd( res_combined, 'econ');
0194         s = trunc_singular( diag(ss), opts.tolOP, true, opts.maxrankRes );
0195 
0196         res{1} = reshape( uu(:,1:s), [size(res{1},1), size(res{1},2), s]);
0197         res{2} = reshape( ss(1:s,1:s)*vv(:,1:s)', [s, size(res{2},2), size(res{2},3), p]);
0198 
0199         left = cat(3, X.U{mu-1}, res{1} );
0200         V = cat(1, X.U{mu}, res{2});
0201 
0202         tmp = permute( V, [1, 4, 2, 3] );
0203         V = reshape( tmp, [size(tmp,1)*p, size(tmp,3)*size(tmp,4)] );
0204 
0205         [U,S,V] = svd( V, 'econ' );
0206         s = trunc_singular( diag(S), opts.tol, true, opts.maxrank );
0207 
0208         V = V(:,1:s)';
0209         X.U{mu} = reshape( V, [s, sz(2), sz(3)] );
0210 
0211         W = U(:,1:s)*S(1:s,1:s);
0212         W = reshape( W, [size(tmp,1), p, s]);
0213         W = permute( W, [1, 3, 2]);
0214         for k = 1:p
0215             C{k} = tensorprod_ttemps( left, W(:,:,k)', 3);
0216         end
0217         X.U{mu-1} = C{1};
0218 
0219         if opts.verbose
0220             disp( ['Augmented system of size (', num2str( [sz(1)*p, sz(2)*sz(3)]), '). Cut-off tol: ', num2str(opts.tol) ])
0221             disp( sprintf( 'Number of SVs: %i. Truncated to: %i => %g %%', length(diag(S)), s, s/length(diag(S))*100))
0222             disp(' ')
0223         end
0224 
0225     end
0226 
0227     % calculate current residuum
0228 
0229 
0230     fprintf(1, '---------------    finshed sweep.    ---------------\n') 
0231     disp(['Current residuum: ', num2str(residuums(:,2*(i-1)+1).')])
0232     disp(' ')
0233     disp(' ')
0234     fprintf(1, '--------------- LEFT-TO-RIGHT SWEEP. ---------------\n') 
0235     % left-to-right sweep
0236     for mu = 1:d-1
0237         sz = [X.rank(mu), X.size(mu), X.rank(mu+1)];
0238 
0239         if opts.verbose
0240             fprintf(1,'Current core: %i. Current iterate (first eigenvalue): \n', mu);
0241             disp(X);
0242             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 )))
0243         else
0244             fprintf(1,'%i  ',mu);
0245         end
0246 
0247         [left, right] = Afun_prepare( A, X, mu );
0248         if opts.precInner
0249             expB = constr_precond_inner( A, X, mu );
0250             [U,L,failureFlag,Lhist ] = lobpcg( rand( prod(sz), p), ...
0251                                                @(y) Afun_block_optim( A, y, sz, left, right, mu), [], ...
0252                                                @(y) apply_local_precond( y, sz, expB), ...
0253                                                max(opts.tolLOBPCG, min( tolLOBPCGmod, sqrt(sum(residuums(:,2*(i-1)+2).^2))/sqrt(p)*tolLOBPCGmod )), ...
0254                                                opts.maxiterLOBPCG, 0);
0255         else
0256             [U,L,failureFlag,Lhist ] = lobpcg( rand( prod(sz), p), ... 
0257                                                @(y) Afun_block_optim( A, y, sz, left, right, mu), ...
0258                                                opts.tolLOBPCG, opts.maxiterLOBPCG, 0);
0259         end
0260        
0261         if opts.verbose
0262             if failureFlag
0263                 fprintf(1,'NOT CONVERGED within %i steps!\n', opts.maxiterLOBPCG)
0264             else
0265                 fprintf(1,'converged after %i steps!\n', size(Lhist,2));
0266             end
0267             disp(['Current Eigenvalue approx: ', num2str( L(1:p)')]);
0268         end
0269 
0270         evalue = [evalue, L(1:p)];
0271         objective = [objective, sum(evalue(:,end))];
0272         elapsed_time = [elapsed_time, toc];
0273 
0274         X.U{mu} = reshape( U, [sz, p] );
0275         lamX = X;
0276         lamX.U{mu} = repmat( reshape(L, [1 1 1 p]), [X.rank(mu), X.size(mu), X.rank(mu+1), 1]).*lamX.U{mu};
0277         res = apply(A, X) - lamX;
0278         
0279         if opts.maxrankRes ~= 0
0280             res_sz = [size(res.U{mu},1), size(res.U{mu},2), size(res.U{mu},3), size(res.U{mu},4)];
0281             res.U{mu} = reshape( permute( res.U{mu}, [1 2 4 3]), [ res_sz(1), res_sz(2)*p, res_sz(3)]);
0282             res = truncate( res, [1 opts.maxrankRes*ones(1,d-1) 1] );
0283             res.U{mu} = ipermute( reshape( res.U{mu}, [ res.rank(mu), res_sz(2), p, res.rank(mu+1)] ), [1 2 4 3]);
0284         end
0285         
0286         % ugly hack to calc the residual
0287         for j = 1:p
0288             tmp = res;
0289             tmp.U{mu} = tmp.U{mu}(:,:,:,j);
0290             resi_norm(j) = norm(tmp);
0291             residuums(j,2*(i-1)+2) = resi_norm(j);
0292         end
0293         micro_res = [micro_res, resi_norm];
0294 
0295         if precondition == 1
0296             res = apply(prec, res);
0297         end
0298         res = contract( X, res, [mu, mu+1]);
0299         res_left = permute( res{1}, [1 2 4 3] );
0300         res_left = reshape( res_left, [size(res{1},1)*size(res{1},2)*p, size(res{1},3)]);
0301         res_combined =  res_left * unfold(res{2},'right');
0302 
0303         if precondition == 2
0304             expBglobal = constr_precond_outer( A, X, mu, mu+1 );
0305             res_size = [size(res{1},1), size(res{1},2), p, size(res{2},2), size(res{2},3)];
0306             res_combined = reshape( res_combined, res_size );
0307             res_combined = permute( res_combined, [1 2 4 5 3]);  % move p to back
0308             res_combined = apply_global_precond( res_combined, res_size([1 2 4 5 3]), expBglobal );
0309             res_combined = ipermute( res_combined, [1 2 4 5 3]);  % move p to back
0310         end
0311         
0312         res_combined = reshape( res_combined, [size(res{1},1)*size(res{1},2)*p, size(res{2},2)*size(res{2},3)]); 
0313 
0314         [uu,ss,vv] = svd( res_combined, 'econ');
0315         s = trunc_singular( diag(ss), opts.tolOP, true, opts.maxrankRes );
0316 
0317         res{1} = reshape( uu(:,1:s)*ss(1:s,1:s), [size(res{1},1), size(res{1},2), p, s]);
0318         res{2} = reshape( vv(:,1:s)', [s, size(res{2},2), size(res{2},3)]);
0319 
0320         right = cat(1, X.U{mu+1}, res{2} );
0321         res{1} = permute( res{1}, [1 2 4 3]);
0322         U = cat(3, X.U{mu}, res{1});
0323 
0324         tmp = permute( U, [1, 2, 4, 3] );
0325         U = reshape( tmp, [size(tmp,1)*size(tmp,2), p*size(tmp,4)] );
0326 
0327         [U,S,V] = svd( U, 'econ' );
0328         s = trunc_singular( diag(S), opts.tol, true, opts.maxrank );
0329         U = U(:,1:s);
0330         X.U{mu} = reshape( U, [sz(1), sz(2), s] );
0331         W = S(1:s,1:s)*V(:,1:s)';
0332         W = reshape( W, [s, p, size(tmp,4)]);
0333         W = permute( W, [1, 3, 2]);
0334         for k = 1:p
0335             C{k} = tensorprod_ttemps( right, W(:,:,k), 1);
0336         end
0337         X.U{mu+1} = C{1};
0338 
0339         if opts.verbose
0340             disp( ['Augmented system of size (', num2str( [sz(1)*p, sz(2)*sz(3)]), '). Cut-off tol: ', num2str(opts.tol) ])
0341             disp( sprintf( 'Number of SVs: %i. Truncated to: %i => %g %%', length(diag(S)), s, s/length(diag(S))*100))
0342             disp(' ')
0343         end
0344 
0345     end
0346     %residuums = [residuums, norm(apply(A,X) - lambda(end)*X)];
0347 
0348     fprintf(1, '---------------    finshed sweep.    ---------------\n') 
0349     disp(['Current residuum: ', num2str(residuums(:,2*(i-1)+2).')])
0350     disp(' ')
0351     disp(' ')
0352 end
0353 
0354 evs = zeros(p,1);
0355 % Evaluate rayleigh quotient for the p TT/MPS tensors
0356 for i=1:p
0357     evec = X;
0358     evec.U{d} = C{i};
0359     evs(i) = innerprod( evec, apply(A, evec));%/ innerprod( evec, evec);
0360 end
0361 evalue = [evalue, evs];
0362 end
0363 
0364 function [left, right] = Afun_prepare( A, x, idx )
0365     y = A.apply(x); 
0366     if idx == 1
0367         right = innerprod( x, y, 'RL', idx+1 );
0368         left = [];
0369     elseif idx == x.order
0370         left = innerprod( x, y, 'LR', idx-1 );
0371         right = [];
0372     else
0373         left = innerprod( x, y, 'LR', idx-1 );
0374         right = innerprod( x, y, 'RL', idx+1 ); 
0375     end
0376 end
0377 
0378 function res = Afun_block_optim( A, U, sz, left, right, mu )
0379 
0380 
0381     p = size(U, 2);
0382     
0383     V = reshape( U, [sz, p] );
0384     V = A.apply( V, mu );
0385     
0386     if mu == 1
0387         tmp = reshape( permute( V, [3 1 2 4] ), [size(V, 3), sz(1)*sz(2)*p]);
0388         tmp = right * tmp;
0389         tmp = reshape( tmp, [size(right, 1), sz(1), sz(2), p]);
0390         tmp = ipermute( tmp, [3 1 2 4] );
0391     elseif mu == A.order
0392         tmp = reshape( V, [size(V,1), sz(2)*sz(3)*p]);
0393         tmp = left * tmp;
0394         tmp = reshape( tmp, [size(left, 1), sz(2), sz(3), p]);
0395     else
0396         tmp = reshape( permute( V, [3 1 2 4] ), [size(V, 3), size(V, 1)*sz(2)*p]);
0397         tmp = right * tmp;
0398         tmp = reshape( tmp, [size(right, 1), size(V, 1), sz(2), p]);
0399         tmp = ipermute( tmp, [3 1 2 4] );
0400 
0401         tmp = reshape( tmp, [size(V, 1), sz(2)*sz(3)*p]);
0402         tmp = left * tmp;
0403         tmp = reshape( tmp, [size(left, 1), sz(2), sz(3), p]);
0404 
0405     end
0406 
0407     res = reshape( tmp, [prod(sz), p] );
0408 end
0409      
0410 
0411 function res = apply_local_precond( U, sz, expB)
0412 
0413     p = size(U, 2);
0414 
0415     x = reshape( U, [sz, p] );
0416     res = zeros( [sz, p] );
0417 
0418     for i = 1:size( expB, 2)
0419         tmp = reshape( x, [sz(1), sz(2)*sz(3)*p] );
0420         tmp = reshape( expB{1,i}*tmp, [sz(1), sz(2), sz(3), p] );
0421 
0422         tmp = reshape( permute( tmp, [2 1 3 4] ), [sz(2), sz(1)*sz(3)*p] );
0423         tmp = ipermute( reshape( expB{2,i}*tmp, [sz(2), sz(1), sz(3), p] ), [2 1 3 4] );
0424 
0425         tmp = reshape( permute( tmp, [3 1 2 4] ), [sz(3), sz(1)*sz(2)*p] );
0426         tmp = ipermute( reshape( expB{3,i}*tmp, [sz(3), sz(1), sz(2), p] ), [3 1 2 4] );
0427 
0428         res = res + tmp;
0429     end
0430     res = reshape( res, [prod(sz), p] );
0431     
0432 end
0433 
0434 function res = apply_global_precond( x, res_sz, expB)
0435 
0436     res = zeros( res_sz );
0437 
0438     for i = 1:size( expB, 2)
0439         tmp = reshape( x, [res_sz(1), prod(res_sz(2:5))] );
0440         tmp = reshape( expB{1,i}*tmp, res_sz );
0441 
0442         tmp = reshape( permute( tmp, [2 1 3 4 5] ), [res_sz(2), res_sz(1)*prod(res_sz(3:5))] );
0443         tmp = ipermute( reshape( expB{2,i}*tmp, [res_sz(2), res_sz(1), res_sz(3:5)] ), [2 1 3 4 5] );
0444 
0445         tmp = reshape( permute( tmp, [3 1 2 4 5] ), [res_sz(3), prod(res_sz([1:2,4:5]))] );
0446         tmp = ipermute( reshape( expB{3,i}*tmp, [res_sz(3), res_sz([1:2,4:5])] ), [3 1 2 4 5] );
0447 
0448         tmp = reshape( permute( tmp, [4 1 2 3 5] ), [res_sz(4), prod(res_sz([1:3,5]))] );
0449         tmp = ipermute( reshape( expB{4,i}*tmp, [res_sz(4), res_sz([1:3,5])] ), [4 1 2 3 5] );
0450 
0451         res = res + tmp;
0452     end
0453     
0454 end