Home > manopt > tools > checkhessian.m

# checkhessian

## PURPOSE

Checks the consistency of the cost function and the Hessian.

## SYNOPSIS

function checkhessian(problem, x, d)

## DESCRIPTION

Checks the consistency of the cost function and the Hessian.

function checkhessian(problem)
function checkhessian(problem, x)
function checkhessian(problem, x, d)

checkhessian performs a numerical test to check that the directional
derivatives and Hessian defined in the problem structure agree up to
second 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).

It is also tested that the result of applying the Hessian along that
direction is indeed a tangent vector, and that the Hessian operator is
symmetric w.r.t. the Riemannian metric.

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

## CROSS-REFERENCE INFORMATION

This function calls:
• StoreDB
• canGetCost Checks whether the cost function can be computed for a problem structure.
• canGetGradient Checks whether the gradient can be computed for a problem structure.
• canGetHessian Checks whether the Hessian can be computed for a problem structure.
• getCost Computes the cost function at x.
• getHessian Computes the Hessian of the cost function at x along d.
• identify_linear_piece Identify a segment of the curve (x, y) that appears to be linear.
This function is called by:

## SOURCE CODE

0001 function checkhessian(problem, x, d)
0002 % Checks the consistency of the cost function and the Hessian.
0003 %
0004 % function checkhessian(problem)
0005 % function checkhessian(problem, x)
0006 % function checkhessian(problem, x, d)
0007 %
0008 % checkhessian performs a numerical test to check that the directional
0009 % derivatives and Hessian defined in the problem structure agree up to
0010 % second order with the cost function at some point x, along some direction
0011 % d. The test is based on a truncated Taylor series (see online Manopt
0012 % documentation).
0013 %
0014 % It is also tested that the result of applying the Hessian along that
0015 % direction is indeed a tangent vector, and that the Hessian operator is
0016 % symmetric w.r.t. the Riemannian metric.
0017 %
0018 % Both x and d are optional and will be sampled at random if omitted.
0019 %
0021
0022 % This file is part of Manopt: www.manopt.org.
0023 % Original author: Nicolas Boumal, Dec. 30, 2012.
0024 % Contributors:
0025 % Change log:
0026 %
0027 %   March 26, 2017 (JB):
0028 %       Detects if the approximated quadratic model is exact
0029 %       and provides the user with the corresponding feedback.
0030 %
0031 %   April 3, 2015 (NB):
0032 %       Works with the new StoreDB class system.
0033 %
0034 %   Nov. 1, 2016 (NB):
0035 %       Issues a call to getGradient rather than getDirectionalDerivative.
0036
0037
0038     % Verify that the problem description is sufficient.
0039     if ~canGetCost(problem)
0040         error('It seems no cost was provided.');
0041     end
0044                 'It seems no gradient was provided.');
0045     end
0046     if ~canGetHessian(problem)
0047         warning('manopt:checkhessian:nohess', ...
0048                 'It seems no Hessian was provided.');
0049     end
0050
0051     x_isprovided = exist('x', 'var') && ~isempty(x);
0052     d_isprovided = exist('d', 'var') && ~isempty(d);
0053
0054     if ~x_isprovided && d_isprovided
0055         error('If d is provided, x must be too, since d is tangent at x.');
0056     end
0057
0058     % If x and / or d are not specified, pick them at random.
0059     if ~x_isprovided
0060         x = problem.M.rand();
0061     end
0062     if ~d_isprovided
0063         d = problem.M.randvec(x);
0064     end
0065
0066     %% Check that the directional derivative and the Hessian at x along d
0067     %% yield a second order model of the cost function.
0068
0069     % Compute the value f0 at f, directional derivative df0 at x along d,
0070     % and Hessian along [d, d].
0071     storedb = StoreDB();
0072     xkey = storedb.getNewKey();
0073     f0 = getCost(problem, x, storedb, xkey);
0074     df0 = problem.M.inner(x, d, getGradient(problem, x, storedb, xkey));
0075     d2f0 = problem.M.inner(x, d, getHessian(problem, x, d, storedb, xkey));
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 quadratic approximation of the cost function using f0,
0089     % df0 and d2f0 at the same points.
0090     model = polyval([.5*d2f0 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(['Hessian check.\nThe slope of the continuous line ' ...
0098                    'should match that of the dashed\n(reference) line ' ...
0099                    'over at least a few orders of magnitude for h.']));
0100     xlabel('h');
0101     ylabel('Approximation error');
0102
0103     line('xdata', [1e-8 1e0], 'ydata', [1e-16 1e8], ...
0104          'color', 'k', 'LineStyle', '--', ...
0105          'YLimInclude', 'off', 'XLimInclude', 'off');
0106
0107
0108     if ~all( err < 1e-12 )
0109         % In a numerically reasonable neighborhood, the error should
0110         % decrease as the cube of the stepsize, i.e., in loglog scale, the
0111         % error should have a slope of 3.
0112         isModelExact = false;
0113         window_len = 10;
0114         [range, poly] = identify_linear_piece(log10(h), log10(err), window_len);
0115     else
0116         % The 2nd order model is exact: all errors are (numerically) zero
0117         % Fit line from all points, use log scale only in h.
0118         isModelExact = true;
0119         range = 1:numel(h);
0120         poly = polyfit(log10(h), err, 1);
0121         % Set mean error in log scale for plot
0122         poly(end) = log10(poly(end));
0123         % Change title to something more descriptive for this special case.
0124         title(sprintf(...
0125               ['Hessian check.\n'...
0126                'It seems the quadratic model is exact:\n'...
0127                'Model error is numerically zero for all h.']));
0128     end
0129     hold all;
0130     loglog(h(range), 10.^polyval(poly, log10(h(range))), 'LineWidth', 3);
0131     hold off;
0132
0133     if ~isModelExact
0134         fprintf('The slope should be 3. It appears to be: %g.\n', poly(1));
0135         fprintf(['If it is far from 3, then directional derivatives or ' ...
0136                  'the Hessian might be erroneous.\n']);
0137         fprintf(['Note: if the exponential map is only approximate, and it '...
0138                  'is not a second-order approximation,\nthen it is normal ' ...
0139                  'for the slope test to reach 2 instead of 3. Check the ' ...
0140                  'factory for this.\n' ...
0141                  'If tested at a critical point, then even for a first-order '...
0142                  'retraction the slope test should yield 3.\n']);
0143     else
0144         fprintf(['The quadratic model appears to be exact ' ...
0145                  '(within numerical precision),\n'...
0146                  'hence the slope computation is irrelevant.\n']);
0147     end
0148
0149
0150     %% Check that the Hessian at x along direction d is a tangent vector.
0151     if isfield(problem.M, 'tangent')
0152         hess = getHessian(problem, x, d, storedb, xkey);
0153         phess = problem.M.tangent(x, hess);
0154         residual = problem.M.lincomb(x, 1, hess, -1, phess);
0155         err = problem.M.norm(x, residual);
0156         fprintf('The residual should be zero, or very close. ');
0157         fprintf('Residual: %g.\n', err);
0158         fprintf(['If it is far from 0, then the Hessian is not in the ' ...
0159                  'tangent plane.\n']);
0160     else
0161         fprintf(['Unfortunately, Manopt was unable to verify that the '...
0162                  'Hessian is indeed a tangent vector.\nPlease verify ' ...
0163                  'this manually.']);
0164     end
0165
0166     %% Check that the Hessian at x is symmetric.
0167     d1 = problem.M.randvec(x);
0168     d2 = problem.M.randvec(x);
0169     h1 = getHessian(problem, x, d1, storedb, xkey);
0170     h2 = getHessian(problem, x, d2, storedb, xkey);
0171     v1 = problem.M.inner(x, d1, h2);
0172     v2 = problem.M.inner(x, h1, d2);
0173     value = v1-v2;
0174     fprintf(['<d1, H[d2]> - <H[d1], d2> should be zero, or very close.' ...
0175              '\n\tValue: %g - %g = %g.\n'], v1, v2, value);
0176     fprintf('If it is far from 0, then the Hessian is not symmetric.\n');
0177
0178 end

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