linearsystem_compare

 Example code for the algorithms described in 

   D. Kressner, M. Steinlechner, and B. Vandereycken.
   Preconditioned low-rank Riemannian optimization for linear systems with tensor product structure.
   Technical report, July 2015. Revised February 2016. To appear in SIAM J. Sci. Comput.

0001 % Example code for the algorithms described in
0002 %
0003 %   D. Kressner, M. Steinlechner, and B. Vandereycken.
0004 %   Preconditioned low-rank Riemannian optimization for linear systems with tensor product structure.
0005 %   Technical report, July 2015. Revised February 2016. To appear in SIAM J. Sci. Comput.
0006 
0007 %   TTeMPS Toolbox.
0008 %   Michael Steinlechner, 2013-2016
0009 %   Questions and contact: michael.steinlechner@epfl.ch
0010 %   BSD 2-clause license, see LICENSE.txt
0011 
0012 clear all;
0013 close all;
0014 clc;
0015 
0016 rng(11);
0017 n = 20;
0018 d = 5;
0019 r = 5;
0020 
0021 nn = n*ones(1,d);
0022 % rr = [1, r*ones(1,d-1), 1];
0023 rr = [1 15 20 20 15 1];
0024 
0025 % Choose one.
0026 % operatorname = 'diffusion';
0027 operatorname = 'newton';
0028 
0029 if strcmpi(operatorname,'diffusion')
0030     A = anisotropicdiffusion( n, d );
0031 elseif strcmpi(operatorname,'newton')
0032     A = newton_potential( n, d );
0033 else
0034     error('Unknown operator requested.')
0035 end
0036 
0037 filename = [operatorname, '_d', num2str(d), '_n', num2str(n), '_r', num2str(r)];
0038 
0039 
0040 if ~exist(['amen/', filename, '.mat'],'file')
0041 
0042     % rng(11);
0043 
0044     % Lapl_op   : pure dD Laplacian as TTeMPS_op_laplace for preconditioning
0045     Riem_CG_TOL = 1e-5; % can be low
0046     ALS_CG_TOL = 1e-10; % cannot be too low, stagnates for d = 10; n = 100; r = 40; otherwise
0047     
0048     Lapl_op = TTeMPS_op_laplace( A.L0, d ); 
0049     Lapl_op = initialize_precond( Lapl_op );
0050   
0051     % rhs is rank one of all ones, scaled to Frobenius norm one
0052     % setup rhs as rank 1 of all ones
0053     rr1 = [1, 1*ones(1,d-1), 1];
0054     F = TTeMPS_rand( rr1, nn );
0055     for i = 1:d
0056         F.U{i} = ones(size(F.U{i}));
0057     end
0058     F = 1/norm(F) * F;
0059 
0060     % Choose operators:
0061     %    A is actual system
0062     %    prec_L_op is the approximate Hessian that will be solved
0063     %    prec_P_op is the preconditioner for prec_L_op. HAS TO BE of Laplace structure
0064     prec_L_op = Lapl_op;
0065     prec_P_op = Lapl_op;
0066 
0067     % One micro step of ALS decreases the error by a huge
0068     % factor. Make the result of one micro step the initial guess.
0069     % X0 = construct_initial_guess(A, F, rr, nn);
0070     X0 = TTeMPS_rand(rr, nn);
0071     X0_amen = X0;...construct_initial_guess(A, F, [1,1*ones(1,d-1),1], nn);
0072 
0073     %%%%%%%%%%%%%%%%%%%%%%% MANOPT NEW TEST CASE %%%%%%%%%%%%%%%%%%%%%%%%%%%
0074     disp('Now testing: Trust Regions!!')
0075 
0076     % Initialize TT fixed rank factory
0077     TT = fixedTTrankfactory(nn, rr);
0078 
0079     % Set up problem in Manopt
0080     problem.M = TT;
0081 
0082     % Setting up the original tensor competion problem for Manopt
0083     problem.cost = @(x) cost_lin(x, A, F);
0084     problem.egrad = @(x) euclidgrad(x, A, F);
0085     problem.ehess = @(x, u) euclidhess(u, A);
0086     % problem.precon = @(x, u) precond_laplace_noSaddle(prec_P_op, u);
0087     
0088     alpha_0 = 1;
0089 
0090     options.Delta0 = alpha_0;
0091     options.Delta_bar = 1e50;
0092     options.maxiter = 400;
0093     options.maxinner = 1000;
0094     options.maxtime = inf;
0095     options.tolgradnorm = 1e-10;
0096     % Collects cost on test points throughout iterations
0097     % options.statsfun = @(problem, x, stats) stats_test(problem, x, stats, A_Gamma, Gamma);
0098 
0099     checkgradient(problem);
0100     checkhessian(problem);
0101     pause;
0102 
0103     % [xL, costManopt, stats] = trustregions(problem, X0, options);
0104 
0105     % norm(xL)
0106     pause;
0107 
0108     % alpha_0 = 1;
0109 
0110     % options.Delta0 = alpha_0;
0111     % options.Delta_bar = 1e50;
0112     % options.maxiter = 400;
0113     % options.maxinner = 1000;
0114     % options.maxtime = inf;
0115     % options.tolgradnorm = 1e-10;
0116 
0117     % problem.ehess = @(x, u) euclidhess(u, A);
0118 
0119     % [xL2, costManopt, stats] = trustregions(problem, xL, options);
0120 
0121     % ==============================
0122     % TEST CASE 1: SIMPLE AMEn
0123     % ==============================
0124     disp('TEST CASE 1: Prec. AMEn, max rank res = 4');
0125     % !!! ALS is always preconditioned with Laplacian part
0126     opts = struct( 'nSweeps', 3, ...
0127                    'solver', 'pcg', ...
0128                    'maxrankRes', 4, ...
0129                    'prec', true)
0130                    
0131     tic
0132     [X_amen_prec1, residuum_amen_prec1, cost_amen_prec1, times_amen_prec1] = amen_fast( A, F, X0_amen, opts )
0133     t_amen_prec1 = toc;
0134     norm(X_amen_prec1)
0135     pause;
0136 
0137     [xL, costManopt, stats] = trustregions(problem, X_amen_prec1, options);
0138 
0139     norm(xL)
0140     pause;
0141 
0142     % ==============================
0143     % TEST CASE 1: SIMPLE AMEn
0144     % ==============================
0145     disp('TEST CASE 1: Prec. AMEn, max rank res = 8');
0146     % !!! ALS is always preconditioned with Laplacian part
0147     opts = struct( 'nSweeps', 3, ...
0148                    'solver', 'pcg', ...
0149                    'maxrankRes', 8, ...
0150                    'prec', true)
0151                    
0152     tic
0153     [X_amen_prec2, residuum_amen_prec2, cost_amen_prec2, times_amen_prec2] = amen_fast( A, F, X0_amen, opts )
0154     t_amen_prec2 = toc;
0155     norm(X_amen_prec2)
0156     pause;
0157     
0158     
0159     % ==============================
0160     % TEST CASE 1: SIMPLE ALS
0161     % ==============================
0162     disp('TEST CASE 1: Simple ALS');
0163     % !!! ALS is always preconditioned with Laplacian part
0164     opts = struct( 'nSweeps', 3, ...
0165                    'solver', 'pcg', ...
0166                    'pcg_accuracy', ALS_CG_TOL);
0167                    
0168     tic
0169     [X_als, residuum_als, cost_als, times_als] = alsLinsolve_fast( A, F, X0, opts )
0170     t_als = toc;
0171 
0172     % ==============================
0173     % TEST CASE 2: Riemannian SD with 1 application of Prec.
0174     % ==============================
0175     disp('TEST CASE 2: Riemannian SD with 1 steps of the approx. prec.');
0176     opts = struct('maxiter', 30, ...
0177                   'precond_tol', Riem_CG_TOL, ...
0178                   'precond_maxit', 1)
0179     
0180     opts.precond_tol = Riem_CG_TOL;
0181     opts.precond_maxit = 1; % 1 application of prec_P_op
0182         
0183     tic;
0184     [X_SD1, residuum_SD1, gradnorm_SD1, cost_SD1, times_SD1] = RiemannLinsolve( A, F, X0, prec_L_op, prec_P_op, opts )
0185     t_pcg = toc;
0186     
0187     % ==============================
0188     %% TEST CASE 4: Truncated (Gauss-)Newton with approx. Prec
0189     %% ==============================
0190     disp('TEST CASE 4: Truncated (Gauss-)Newton.');
0191     opts = struct('maxiter', 20, ...
0192                   'truncatedNewton', true);
0193     tic
0194     [X_tn, residuum_tn, gradnorm_tn, cost_tn, times_tn] = RiemannLinsolve( A, F, X0, prec_L_op, prec_P_op, opts )
0195     t_tn = toc;
0196 
0197     save([filename, '.mat'], 'X_als','residuum_als', 'cost_als', 'times_als', ...
0198                    'X_amen_prec1','residuum_amen_prec1', 'cost_amen_prec1', 'times_amen_prec1', ...
0199                    'X_amen_prec2','residuum_amen_prec2', 'cost_amen_prec2', 'times_amen_prec2', ...
0200                    'X_SD1','residuum_SD1', 'gradnorm_SD1', 'cost_SD1', 'times_SD1', ...
0201                    'X_tn', 'residuum_tn', 'gradnorm_tn', 'cost_tn', 'times_tn')
0202 
0203 
0204     
0205 else
0206     load(['amen/', filename, '.mat'])
0207 end
0208 
0209 
0210 % setup plotting
0211 col = lines(5);
0212 leg = {};
0213 
0214 % dummy plot for legend entry
0215 figure(1)
0216 semilogy(1,residuum_als(1), '-o', 'Color', col(1,:),'linewidth',2)
0217 hold on
0218 figure(2)
0219 semilogy(times_als(1), residuum_als(1), '-o', 'Color', col(1,:),'linewidth',2)
0220 hold on
0221 drawnow
0222 
0223 figure(1)
0224 semilogy(1,residuum_amen_prec1(1), '-o', 'Color', col(3,:),'linewidth',2)
0225 hold on
0226 figure(2)
0227 semilogy(times_amen_prec1(1), residuum_amen_prec1(1), '-o', 'Color', col(3,:),'linewidth',2)
0228 hold on
0229 drawnow
0230 
0231 figure(1)
0232 semilogy(1,residuum_amen_prec2(1), '-o', 'Color', col(2,:),'linewidth',2)
0233 hold on
0234 figure(2)
0235 semilogy(times_amen_prec2(1), residuum_amen_prec2(1), '-o', 'Color', col(2,:),'linewidth',2)
0236 hold on
0237 drawnow
0238 
0239 % Start drawing the graphs
0240 leg = [leg, 'ALS', 'AMEn, max. rank res. = 4', 'AMEn, max. rank res. = 8'];
0241 
0242 figure(1)
0243 semilogy(residuum_SD1, '--', 'Color', col(5,:),'linewidth',2)
0244 figure(2)
0245 semilogy(times_SD1, residuum_SD1, '--', 'Color', col(5,:),'linewidth',2)
0246 drawnow
0247 
0248 leg = [leg, 'Prec. steepest descent'];
0249 
0250 figure(1)
0251 semilogy(residuum_tn, '-k','linewidth',2)
0252 figure(2)
0253 semilogy(times_tn, residuum_tn, '-k','linewidth',2)
0254 leg = [leg, 'Approx. Newton'];
0255 drawnow
0256 
0257 figure(1)
0258 semilogy(residuum_als, '-', 'Color', col(1,:),'linewidth',2)
0259 iterations = 1:length(residuum_als);
0260 semilogy(iterations(1:d:end),residuum_als(1:d:end), 'o', 'Color', col(1,:),'linewidth',2)
0261 hold on
0262 figure(2)
0263 semilogy(times_als, residuum_als, '-', 'Color', col(1,:),'linewidth',2)
0264 semilogy(times_als(1:d:end), residuum_als(1:d:end), 'o', 'Color', col(1,:),'linewidth',2)
0265 hold on
0266 drawnow
0267 
0268 figure(1)
0269 semilogy(residuum_amen_prec1, '-', 'Color', col(3,:),'linewidth',2)
0270 iterations = 1:length(residuum_amen_prec1);
0271 semilogy(iterations(1:d:end),residuum_amen_prec1(1:d:end), 'o', 'Color', col(3,:),'linewidth',2)
0272 hold on
0273 figure(2)
0274 semilogy(times_amen_prec1, residuum_amen_prec1, '-', 'Color', col(3,:),'linewidth',2)
0275 semilogy(times_amen_prec1(1:d:end), residuum_amen_prec1(1:d:end), 'o', 'Color', col(3,:),'linewidth',2)
0276 hold on
0277 drawnow
0278 
0279 figure(1)
0280 semilogy(residuum_amen_prec2, '-', 'Color', col(2,:),'linewidth',2)
0281 iterations = 1:length(residuum_amen_prec2);
0282 semilogy(iterations(1:d:end),residuum_amen_prec2(1:d:end), 'o', 'Color', col(2,:),'linewidth',2)
0283 hold on
0284 figure(2)
0285 semilogy(times_amen_prec2, residuum_amen_prec2, '-', 'Color', col(2,:),'linewidth',2)
0286 semilogy(times_amen_prec2(1:d:end), residuum_amen_prec2(1:d:end), 'o', 'Color', col(2,:),'linewidth',2)
0287 hold on
0288 drawnow
0289 
0290 figure(1)
0291 legend(leg);
0292 set(gca,'fontsize',16)
0293 xlabel('Iterations')
0294 ylabel('Relative residual')
0295 set(gca,'fontsize',16)
0296 ylim([1e-10,1e1])
0297 
0298 figure(2)
0299 legend(leg);
0300 set(gca,'fontsize',16)
0301 xlabel('Time [s]')
0302 ylabel('Relative residual')
0303 set(gca,'fontsize',16)
0304 ylim([1e-10,1e1])
0305 
0306 
0307 function c = cost_lin(X, A, F)
0308     grad = euclidgrad(X, A, F);
0309     c = 0.5 * innerprod( X, grad);
0310 end
0311 
0312 
0313 function g = euclidgrad(X, A, F)
0314     g = (1/ (20^5))*(apply(A, X) - F);
0315 end
0316 
0317 function h = euclidhess(u, A)
0318     uTT = tangent_to_TTeMPS(u);
0319     h = (1/ (20^5))*apply(A, uTT);
0320 end