Home > manopt > manifolds > ttfixedrank > TTeMPS_1.1 > examples > linearsystem_compare.m

# linearsystem_compare

## PURPOSE

Example code for the algorithms described in

## SYNOPSIS

This is a script file.

## DESCRIPTION

``` 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.```

## CROSS-REFERENCE INFORMATION

This function calls:
• disp DISP Display TT/MPS tensor.
• innerprod INNERPROD Inner product between two TT/MPS tensors.
• norm NORM Norm of a TT/MPS tensor.
• disp DISP Display TT/MPS block-mu tensor.
• innerprod INNERPROD Inner product between two TT/MPS tensors.
• norm NORM Norm of a TT/MPS block-mu tensor.
• apply APPLY Application of TT/MPS operator to a TT/MPS tensor
• disp DISP Display TT/MPS operator.
• TTeMPS_op_laplace
• apply APPLY Application of TT/MPS Laplace-like operator to a TT/MPS tensor
• disp DISP Display TT/MPS operator.
• TTeMPS_rand TTEMPS_RAND Create random TTeMPS tensor
• RiemannLinsolve Riemannian approx. Newton for linear systems. For more information, we refer to the report
• alsLinsolve_fast TTeMPS Toolbox.
• amen_fast TTeMPS Toolbox.
• anisotropicdiffusion
• newton_potential
• fixedTTrankfactory Manifold of tensors of fixed Tensor Train (TT) rank, embedded geometry
• trustregions Riemannian trust-regions solver for optimization on manifolds.
• checkgradient Checks the consistency of the cost function and the gradient.
• checkhessian Checks the consistency of the cost function and the Hessian.
This function is called by:

## SOURCE CODE

```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
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);
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;
0096     % Collects cost on test points throughout iterations
0097     % options.statsfun = @(problem, x, stats) stats_test(problem, x, stats, A_Gamma, Gamma);
0098
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;
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
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)
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```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005