0001
0002
0003
0004
0005 function [X, residuum, cost] = alsLinsolve( L, F, X, opts )
0006
0007
0008 if ~exist( 'opts', 'var'); opts = struct(); end
0009 if ~isfield( opts, 'nSweeps'); opts.nSweeps = 4; end
0010
0011 d = X.order;
0012 n = X.size;
0013
0014
0015 normF = norm(F);
0016 cost = cost_function( L, X, F );
0017 residuum = norm( apply(L, X) - F ) / normF;
0018
0019 for sweep = 1:opts.nSweeps
0020 X = orthogonalize(X, 1);
0021 for i = 1:d-1
0022 disp( ['Current core: ', num2str(i)] )
0023
0024 Li = contract( L, X, i );
0025 Fi = contract( X, F, i );
0026
0027 Ui = Li \ Fi(:);
0028 X.U{i} = reshape( Ui, size(X.U{i}) );
0029 X = orth_at( X, i, 'left', true );
0030
0031 residuum = [residuum; norm( apply(L, X) - F ) / normF];
0032 cost = [cost; cost_function( L, X, F )];
0033 end
0034 for i = d:-1:2
0035 disp( ['Current core: ', num2str(i)] )
0036
0037 Li = contract( L, X, i );
0038 Fi = contract( X, F, i );
0039
0040 Ui = Li \ Fi(:);
0041 X.U{i} = reshape( Ui, size(X.U{i}) );
0042 X = orth_at( X, i, 'right', true );
0043 residuum = [residuum; norm( apply(L, X) - F ) / normF];
0044 cost = [cost; cost_function( L, X, F )];
0045 end
0046
0047 end
0048
0049
0050 end
0051
0052 function res = cost_function( L, X, F )
0053 res = 0.5*innerprod( X, apply(L, X) ) - innerprod( X, F );
0054 end
0055
0056 function res = euclid_grad( L, X, F )
0057 res = apply(L, X) - F;
0058 end
0059