checkretraction

 Check the order of agreement of a retraction with an exponential.
 
 function checkretraction(M)
 function checkretraction(M, x)
 function checkretraction(M, x, v)

 checkretraction performs a numerical test to check the order of agreement
 between the retraction and the exponential map in a given Manopt
 manifold structure M. The test is performed at the point x if it is
 provided (otherwise, the point is picked at random) and along the tangent
 vector v at x if one is provided (otherwise, a tangent vector at x is
 picked at random.)

 See also: checkmanifold checkdiff checkgradient checkhessian

0001 function checkretraction(M, x, v)
0002 % Check the order of agreement of a retraction with an exponential.
0003 %
0004 % function checkretraction(M)
0005 % function checkretraction(M, x)
0006 % function checkretraction(M, x, v)
0007 %
0008 % checkretraction performs a numerical test to check the order of agreement
0009 % between the retraction and the exponential map in a given Manopt
0010 % manifold structure M. The test is performed at the point x if it is
0011 % provided (otherwise, the point is picked at random) and along the tangent
0012 % vector v at x if one is provided (otherwise, a tangent vector at x is
0013 % picked at random.)
0014 %
0015 % See also: checkmanifold checkdiff checkgradient checkhessian
0016 
0017 % This file is part of Manopt: www.manopt.org.
0018 % Original author: Nicolas Boumal, Oct. 21, 2016.
0019 % Contributors:
0020 % Change log:
0021 
0022     if ~isfield(M, 'exp')
0023         error(['This manifold has no exponential (M.exp): ' ...
0024                'no reference to compare the retraction.']);
0025     end
0026     if ~isfield(M, 'dist')
0027         error(['This manifold has no distance (M.dist): ' ...
0028                'this is required to run this check.']);
0029     end
0030 
0031     if ~exist('x', 'var') || isempty(x)
0032         x = M.rand();
0033         v = M.randvec(x);
0034     end
0035     
0036     if ~exist('v', 'var') || isempty(v)
0037         v = M.randvec(x);
0038     end
0039     
0040     % Compare the retraction and the exponential over steps of varying
0041     % length, on a wide log-scale.
0042     tt = logspace(-12, 0, 251);
0043     ee = zeros(size(tt));
0044     for k = 1 : numel(tt)
0045         t = tt(k);
0046         ee(k) = M.dist(M.exp(x, v, t), M.retr(x, v, t));
0047     end
0048     
0049     % Plot the difference between the exponential and the retration over
0050     % that span of steps, in log-log scale.
0051     loglog(tt, ee);
0052     
0053     % We hope to see a slope of 3, to confirm a second-order retraction. If
0054     % the slope is only 2, we have a first-order retration. If the slope is
0055     % less than 2, this is not a retraction.
0056     % Slope 3
0057     line('xdata', [1e-12 1e0], 'ydata', [1e-30 1e6], ...
0058          'color', 'k', 'LineStyle', '--', ...
0059          'YLimInclude', 'off', 'XLimInclude', 'off');
0060     % Slope 2
0061     line('xdata', [1e-14 1e0], 'ydata', [1e-20 1e8], ...
0062          'color', 'k', 'LineStyle', ':', ...
0063          'YLimInclude', 'off', 'XLimInclude', 'off');
0064      
0065 
0066     % Figure out the slope of the error in log-log, by identifying a piece
0067     % of the error curve which is mostly linear.
0068     window_len = 10;
0069     [range, poly] = identify_linear_piece(log10(tt), log10(ee), window_len);
0070     hold all;
0071     loglog(tt(range), 10.^polyval(poly, log10(tt(range))), 'LineWidth', 3);
0072     hold off;
0073     
0074     xlabel('Step size multiplier t');
0075     ylabel('Distance between Exp(x, v, t) and Retr(x, v, t)');
0076     title(sprintf('Retraction check.\nA slope of 2 is required, 3 is desired.'));
0077     
0078     fprintf('Check agreement between M.exp and M.retr. Please check the\n');
0079     fprintf('factory file of M to ensure M.exp is a proper exponential.\n');
0080     fprintf('The slope must be at least 2 to have a proper retraction.\n');
0081     fprintf('For the retraction to be second order, the slope should be 3.\n');
0082     fprintf('It appears the slope is: %g.\n', poly(1));
0083     fprintf('Note: if exp and retr are identical, this is about zero: %g.\n', norm(ee));
0084     fprintf('In the latter case, the slope test is irrelevant.\n');
0085 
0086 end