alsLinsolve

   TTeMPS Toolbox. 
   Michael Steinlechner, 2013-2016
   Questions and contact: michael.steinlechner@epfl.ch
   BSD 2-clause license, see LICENSE.txt

0001 %   TTeMPS Toolbox.
0002 %   Michael Steinlechner, 2013-2016
0003 %   Questions and contact: michael.steinlechner@epfl.ch
0004 %   BSD 2-clause license, see LICENSE.txt
0005 function [X, residuum, cost] = alsLinsolve( L, F, X, opts )
0006 
0007 % set default opts
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