Home > manopt > tools > checkdiff.m

checkdiff

PURPOSE ^

Checks the consistency of the cost function and directional derivatives.

SYNOPSIS ^

function checkdiff(problem, x, d, force_gradient)

DESCRIPTION ^

 Checks the consistency of the cost function and directional derivatives.

 function checkdiff(problem)
 function checkdiff(problem, x)
 function checkdiff(problem, x, d)

 checkdiff performs a numerical test to check that the directional
 derivatives defined in the problem structure agree up to first order with
 the cost function at some point x, along some direction d. The test is
 based on a truncated Taylor series (see online Manopt documentation).

 Both x and d are optional and will be sampled at random if omitted.

 See also: checkgradient checkhessian

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function checkdiff(problem, x, d, force_gradient)
0002 % Checks the consistency of the cost function and directional derivatives.
0003 %
0004 % function checkdiff(problem)
0005 % function checkdiff(problem, x)
0006 % function checkdiff(problem, x, d)
0007 %
0008 % checkdiff performs a numerical test to check that the directional
0009 % derivatives defined in the problem structure agree up to first order with
0010 % the cost function at some point x, along some direction d. The test is
0011 % based on a truncated Taylor series (see online Manopt documentation).
0012 %
0013 % Both x and d are optional and will be sampled at random if omitted.
0014 %
0015 % See also: checkgradient checkhessian
0016 
0017 % If force_gradient = true (hidden parameter), then the function will call
0018 % getGradient and infer the directional derivative, rather than call
0019 % getDirectionalDerivative directly. This is used by checkgradient.
0020 
0021 % This file is part of Manopt: www.manopt.org.
0022 % Original author: Nicolas Boumal, Dec. 30, 2012.
0023 % Contributors:
0024 % Change log:
0025 %
0026 %   March 26, 2017 (JB):
0027 %       Detects if the approximated linear model is exact
0028 %       and provides the user with the corresponding feedback.
0029 %
0030 %   April 3, 2015 (NB):
0031 %       Works with the new StoreDB class system.
0032 
0033     if ~exist('force_gradient', 'var')
0034         force_gradient = false;
0035     end
0036         
0037     % Verify that the problem description is sufficient.
0038     if ~canGetCost(problem)
0039         error('It seems no cost was provided.');
0040     end
0041     if ~force_gradient && ~canGetDirectionalDerivative(problem)
0042         error('It seems no directional derivatives were provided.');
0043     end
0044     if force_gradient && ~canGetGradient(problem)
0045         % Would normally issue a warning, but this function should only be
0046         % called with force_gradient on by checkgradient, which will
0047         % already have issued a warning.
0048     end
0049         
0050     x_isprovided = exist('x', 'var') && ~isempty(x);
0051     d_isprovided = exist('d', 'var') && ~isempty(d);
0052     
0053     if ~x_isprovided && d_isprovided
0054         error('If d is provided, x must be too, since d is tangent at x.');
0055     end
0056     
0057     % If x and / or d are not specified, pick them at random.
0058     if ~x_isprovided
0059         x = problem.M.rand();
0060     end
0061     if ~d_isprovided
0062         d = problem.M.randvec(x);
0063     end
0064 
0065     % Compute the value f0 at f and directional derivative at x along d.
0066     storedb = StoreDB();
0067     xkey = storedb.getNewKey();
0068     f0 = getCost(problem, x, storedb, xkey);
0069     
0070     if ~force_gradient
0071         df0 = getDirectionalDerivative(problem, x, d, storedb, xkey);
0072     else
0073         grad = getGradient(problem, x, storedb, xkey);
0074         df0 = problem.M.inner(x, grad, d);
0075     end
0076     
0077     % Compute the value of f at points on the geodesic (or approximation
0078     % of it) originating from x, along direction d, for stepsizes in a
0079     % large range given by h.
0080     h = logspace(-8, 0, 51);
0081     value = zeros(size(h));
0082     for i = 1 : length(h)
0083         y = problem.M.exp(x, d, h(i));
0084         ykey = storedb.getNewKey();
0085         value(i) = getCost(problem, y, storedb, ykey);
0086     end
0087     
0088     % Compute the linear approximation of the cost function using f0 and
0089     % df0 at the same points.
0090     model = polyval([df0 f0], h);
0091     
0092     % Compute the approximation error
0093     err = abs(model - value);
0094     
0095     % And plot it.
0096     loglog(h, err);
0097     title(sprintf(['Directional derivative check.\nThe slope of the '...
0098                    'continuous line should match that of the dashed\n'...
0099                    '(reference) line over at least a few orders of '...
0100                    'magnitude for h.']));
0101     xlabel('h');
0102     ylabel('Approximation error');
0103     
0104     line('xdata', [1e-8 1e0], 'ydata', [1e-8 1e8], ...
0105          'color', 'k', 'LineStyle', '--', ...
0106          'YLimInclude', 'off', 'XLimInclude', 'off');
0107     
0108      
0109     if ~all( err < 1e-12 )
0110         % In a numerically reasonable neighborhood, the error should
0111         % decrease as the square of the stepsize, i.e., in loglog scale,
0112         % the error should have a slope of 2.
0113         isModelExact = false;
0114         window_len = 10;
0115         [range, poly] = identify_linear_piece(log10(h), log10(err), window_len);
0116     else
0117         % The 1st order model is exact: all errors are (numerically) zero
0118         % Fit line from all points, use log scale only in h.
0119         isModelExact = true;
0120         range = 1:numel(h);
0121         poly = polyfit(log10(h), err, 1);
0122         % Set mean error in log scale for plot.
0123         poly(end) = log10(poly(end));
0124         % Change title to something more descriptive for this special case.
0125         title(sprintf(...
0126               ['Directional derivative check.\n'...
0127                'It seems the linear model is exact:\n'...
0128                'Model error is numerically zero for all h.']));
0129     end
0130     hold all;
0131     loglog(h(range), 10.^polyval(poly, log10(h(range))), 'LineWidth', 3);
0132     hold off;
0133     
0134     if ~isModelExact
0135         fprintf('The slope should be 2. It appears to be: %g.\n', poly(1));
0136         fprintf(['If it is far from 2, then directional derivatives ' ...
0137                  'might be erroneous.\n']);
0138     else
0139         fprintf(['The linear model appears to be exact ' ...
0140                  '(within numerical precision),\n'...
0141                  'hence the slope computation is irrelevant.\n']);
0142     end
0143     
0144 end

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005