0001 function [X, C, evalue, residuums, micro_res, objective, elapsed_time] = amen_eigenvalue(A, prec, p, rr, opts)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
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
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
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 )))
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
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
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
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
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]);
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]);
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
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
0356 for i=1:p
0357 evec = X;
0358 evec.U{d} = C{i};
0359 evs(i) = innerprod( evec, apply(A, 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