RiemannPrecondSteep

   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 
0006 function [xR, residuum, gradnorm, cost, times] = RiemannPrecondSteep( L, F, X0, Lh, Ph, opts )
0007 % L is the actual operator
0008 %    needs a apply(L, X) interface but can be anything
0009 % Lh is the operator that represents the (inexact) Euclidean Hessian
0010 %    needs a apply(L, X) interface but can be anything
0011 % Ph is the preconditioner for Lh; should be a TTeMPS_op_laplace operator
0012 %
0013 % When L is Laplace+perturbation, taking Lh and Ph both the Laplacian works
0014 % well.
0015 
0016 t_start = tic();
0017 % set default opts
0018 if ~exist( 'opts', 'var');       opts = struct();     end
0019 if ~isfield( opts, 'maxiter');   opts.maxiter = 500;  end
0020 if ~isfield( opts, 'tol');       opts.tol = 1e-16;     end
0021 if ~isfield( opts, 'safe_norm');       opts.safe_norm = false;     end
0022 
0023 d = X0.order;
0024 n = X0.size;
0025 
0026 
0027 
0028 [xL, xR, G] = gauge_matrices( X0 );
0029 
0030 %xL = orthogonalize(X, d);
0031 %xR = orthogonalize(X, 1);
0032 
0033 cost = zeros(opts.maxiter, 1);
0034 residuum = zeros(opts.maxiter, 1);
0035 gradnorm = zeros(opts.maxiter, 1);
0036 times = zeros(opts.maxiter, 1);
0037 normRHS = norm(F);
0038 
0039 for i = 1:opts.maxiter
0040     
0041     g = euclid_grad( L, xR, F );
0042     
0043     cost(i) = cost_function_res( xR, g );
0044     residuum(i) = norm(g, opts.safe_norm) / normRHS;
0045     times(i) = toc(t_start);
0046     
0047     % test for stopping criterion
0048     if abs(residuum(i)) < opts.tol
0049         sprintf( 'Current residual: %g', residuum(i))
0050         residuum = residuum(1:i);
0051         cost = cost(1:i);
0052         times = times(1:i);
0053         sprintf( 'RiemannLinsolve CONVERGED after %i iterations', i )
0054         break
0055     end
0056     
0057     grad = TTeMPS_tangent_orth( xL, xR, g );
0058     gradnorm(i) = norm( grad );
0059     
0060     sprintf('steepest descent step %i', i)
0061     %P_grad = solvePrecond( L, P, grad, xL, xR, opts );
0062     P_grad = solvePrecond_noSaddle( Lh, Ph, grad, xL, xR, opts, G );      
0063         
0064        
0065         
0066     %check_precond_laplace(P, grad, P_grad)
0067     %eta = -P_grad;
0068     
0069     %line search
0070     alpha = linesearch_linearized( L, P_grad, g )
0071     %alpha = linesearch_linearized2( L, P_grad, grad )
0072     %alpha = -1;
0073     xR = tangentAdd(  P_grad, alpha, true );
0074     
0075     [xL, G] = left_orth_with_gauge( xR );
0076     %xL = orthogonalize( X, d );
0077     %xR = orthogonalize( X, 1 );
0078 end
0079 
0080 end
0081 
0082 function res = cost_function( L, X, F )
0083 res = 0.5*innerprod( X, apply(L, X) ) - innerprod( X, F );
0084 end
0085 function res = cost_function_res( X, res )
0086 res = 0.5*innerprod( X, res );
0087 end
0088 
0089 function res = euclid_grad( L, X, F )
0090 res = apply(L, X) - F;
0091 end
0092 
0093 function alpha = linesearch_linearized( L, xi, g )
0094 eta = tangent_to_TTeMPS( xi );
0095 alpha = -innerprod( eta, g );
0096 alpha = alpha / innerprod( eta, apply(L, eta) );
0097 end
0098 
0099 function alpha = linesearch_linearized2( L, xi, grad )
0100 alpha = -innerprod( xi, grad );
0101 eta = tangent_to_TTeMPS( xi );
0102 alpha = alpha / innerprod( eta, apply(L, eta) );
0103 end