0001 function checkretraction(M, x, v)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
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
0041
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
0050
0051 loglog(tt, ee);
0052
0053
0054
0055
0056
0057 line('xdata', [1e-12 1e0], 'ydata', [1e-30 1e6], ...
0058 'color', 'k', 'LineStyle', '--', ...
0059 'YLimInclude', 'off', 'XLimInclude', 'off');
0060
0061 line('xdata', [1e-14 1e0], 'ydata', [1e-20 1e8], ...
0062 'color', 'k', 'LineStyle', ':', ...
0063 'YLimInclude', 'off', 'XLimInclude', 'off');
0064
0065
0066
0067
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