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.

The slope test requires the exponential map of the manifold, or at least
a second-order retraction. If M.exp() is not available, the retraction is
used and a message is issued to instruct the user to check whether M.retr
is second-order or not. If it is not, the slope test is only valid at
critical points of the cost function (which can be computed by
optimization.)

## 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 %
0020 % The slope test requires the exponential map of the manifold, or at least
0021 % a second-order retraction. If M.exp() is not available, the retraction is
0022 % used and a message is issued to instruct the user to check whether M.retr
0023 % is second-order or not. If it is not, the slope test is only valid at
0024 % critical points of the cost function (which can be computed by
0025 % optimization.)
0026 %
0028
0029 % This file is part of Manopt: www.manopt.org.
0030 % Original author: Nicolas Boumal, Dec. 30, 2012.
0031 % Contributors:
0032 % Change log:
0033 %
0034 %   March 26, 2017 (JB):
0035 %       Detects if the approximated quadratic model is exact
0036 %       and provides the user with the corresponding feedback.
0037 %
0038 %   April 3, 2015 (NB):
0039 %       Works with the new StoreDB class system.
0040 %
0041 %   Nov. 1, 2016 (NB):
0042 %       Issues a call to getGradient rather than getDirectionalDerivative.
0043 %
0044 %   Dec. 6, 2017 (NB):
0045 %       Added message in case tangent2ambient might be necessary in
0046 %       defining ehess (this was a common difficulty for users.)
0047 %
0048 %   Aug. 2, 2018 (NB):
0049 %       Using storedb.remove() to avoid unnecessary cache build-up.
0050 %
0051 %   Sep. 6, 2018 (NB):
0052 %       Now checks whether M.exp() is available; uses retraction otherwise
0053 %       and issues a message that the user should check whether the
0054 %       retraction is second-order or not.
0055
0056
0057     % Verify that the problem description is sufficient.
0058     if ~canGetCost(problem)
0059         error('It seems no cost was provided.');
0060     end
0063                 'It seems no gradient was provided.');
0064     end
0065     if ~canGetHessian(problem)
0066         warning('manopt:checkhessian:nohess', ...
0067                 'It seems no Hessian was provided.');
0068     end
0069
0070     x_isprovided = exist('x', 'var') && ~isempty(x);
0071     d_isprovided = exist('d', 'var') && ~isempty(d);
0072
0073     if ~x_isprovided && d_isprovided
0074         error('If d is provided, x must be too, since d is tangent at x.');
0075     end
0076
0077     % If x and / or d are not specified, pick them at random.
0078     if ~x_isprovided
0079         x = problem.M.rand();
0080     end
0081     if ~d_isprovided
0082         d = problem.M.randvec(x);
0083     end
0084
0085     %% Check that the directional derivative and the Hessian at x along d
0086     %% yield a second order model of the cost function.
0087
0088     % Compute the value f0 at f, directional derivative df0 at x along d,
0089     % and Hessian along [d, d].
0090     storedb = StoreDB();
0091     xkey = storedb.getNewKey();
0092     f0 = getCost(problem, x, storedb, xkey);
0094     df0 = problem.M.inner(x, d, gradx);
0095     hessxd = getHessian(problem, x, d, storedb, xkey);
0096     d2f0 = problem.M.inner(x, d, hessxd);
0097
0098
0099     % Pick a stepping function: exponential or retraction?
0100     if isfield(problem.M, 'exp')
0101         stepper = problem.M.exp;
0102         extra_message = '';
0103     else
0104         stepper = problem.M.retr;
0105         fprintf(['* M.exp() is not available: using M.retr() instead.\n' ...
0106                  '* Please check the manifold documentation to see if\n' ...
0107                  '* the retraction is second order. If not, the slope\n' ...
0108                  '* test is allowed to fail at non-critical x.\n']);
0109         extra_message = ['(But do mind the message above: the slope is\n' ...
0110                          'allowed to be 2 at non-critical points x.)\n'];
0111     end
0112
0113
0114     % Compute the value of f at points on the geodesic (or approximation
0115     % of it) originating from x, along direction d, for stepsizes in a
0116     % large range given by h.
0117     h = logspace(-8, 0, 51);
0118     value = zeros(size(h));
0119     for i = 1 : length(h)
0120         y = stepper(x, d, h(i));
0121         ykey = storedb.getNewKey();
0122         value(i) = getCost(problem, y, storedb, ykey);
0123         storedb.remove(ykey); % no need to keep it in memory
0124     end
0125
0126     % Compute the quadratic approximation of the cost function using f0,
0127     % df0 and d2f0 at the same points.
0128     model = polyval([.5*d2f0 df0 f0], h);
0129
0130     % Compute the approximation error
0131     err = abs(model - value);
0132
0133     % And plot it.
0134     loglog(h, err);
0135     title(sprintf(['Hessian check.\nThe slope of the continuous line ' ...
0136                    'should match that of the dashed\n(reference) line ' ...
0137                    'over at least a few orders of magnitude for h.']));
0138     xlabel('h');
0139     ylabel('Approximation error');
0140
0141     line('xdata', [1e-8 1e0], 'ydata', [1e-16 1e8], ...
0142          'color', 'k', 'LineStyle', '--', ...
0143          'YLimInclude', 'off', 'XLimInclude', 'off');
0144
0145
0146     if ~all( err < 1e-12 )
0147         % In a numerically reasonable neighborhood, the error should
0148         % decrease as the cube of the stepsize, i.e., in loglog scale, the
0149         % error should have a slope of 3.
0150         isModelExact = false;
0151         window_len = 10;
0152         [range, poly] = identify_linear_piece(log10(h), log10(err), window_len);
0153     else
0154         % The 2nd order model is exact: all errors are (numerically) zero
0155         % Fit line from all points, use log scale only in h.
0156         isModelExact = true;
0157         range = 1:numel(h);
0158         poly = polyfit(log10(h), err, 1);
0159         % Set mean error in log scale for plot
0160         poly(end) = log10(poly(end));
0161         % Change title to something more descriptive for this special case.
0162         title(sprintf(...
0163               ['Hessian check.\n'...
0164                'It seems the quadratic model is exact:\n'...
0165                'Model error is numerically zero for all h.']));
0166     end
0167     hold all;
0168     loglog(h(range), 10.^polyval(poly, log10(h(range))), 'LineWidth', 3);
0169     hold off;
0170
0171     if ~isModelExact
0172         fprintf('The slope should be 3. It appears to be: %g.\n', poly(1));
0173         fprintf(['If it is far from 3, then directional derivatives,\n' ...
0174                  'the gradient or the Hessian might be erroneous.\n', ...
0175                  extra_message]);
0176     else
0177         fprintf(['The quadratic model appears to be exact ' ...
0178                  '(within numerical precision),\n'...
0179                  'hence the slope computation is irrelevant.\n']);
0180     end
0181
0182
0183     %% Check that the Hessian at x along direction d is a tangent vector.
0184     if isfield(problem.M, 'tangent')
0185         hess = getHessian(problem, x, d, storedb, xkey);
0186         phess = problem.M.tangent(x, hess);
0187         residual = problem.M.lincomb(x, 1, hess, -1, phess);
0188         err = problem.M.norm(x, residual);
0189         fprintf('Tangency residual should be zero, or very close; ');
0190         fprintf('residual: %g.\n', err);
0191         fprintf(['If it is far from 0, then the Hessian is not in the ' ...
0192                  'tangent space.\n']);
0193     else
0194         fprintf(['Unfortunately, Manopt was unable to verify that the '...
0195                  'output of the Hessian call is indeed a tangent ' ...
0197     end
0198
0199     %% Check that the Hessian at x is symmetric.
0200     d1 = problem.M.randvec(x);
0201     d2 = problem.M.randvec(x);
0202     h1 = getHessian(problem, x, d1, storedb, xkey);
0203     h2 = getHessian(problem, x, d2, storedb, xkey);
0204     v1 = problem.M.inner(x, d1, h2);
0205     v2 = problem.M.inner(x, h1, d2);
0206     value = v1-v2;
0207     fprintf(['<d1, H[d2]> - <H[d1], d2> should be zero, or very close.' ...
0208              '\n\tValue: %g - %g = %g.\n'], v1, v2, value);
0209     fprintf('If it is far from 0, then the Hessian is not symmetric.\n');
0210
0211     %% Check if the manifold at hand is one of those for which there should
0212     %  be a call to M.tangent2ambient, as this is a common mistake. If so,
0213     %  issue an additional message. Ideally, one would just need to check
0214     %  for the presence of tangent2ambient, but productmanifold (for
0215     %  example) generates one of those in all cases, even if it is just an
0216     %  identity map.
0217     if isfield(problem.M, 'tangent2ambient_is_identity') && ...
0218                                      ~problem.M.tangent2ambient_is_identity
0219
0220         fprintf('\n\n=== Special message ===\n');
0221         fprintf(['For this manifold, tangent vectors are represented\n' ...
0222                  'differently from their ambient space representation.\n' ...
0223                  'In practice, this means that when defining\n' ...
0224                  'v = problem.ehess(x, u), one may need to call\n' ...
0225                  'u = problem.M.tangent2ambient(x, u) first, so as to\n'...
0226                  'transform u into an ambient vector, if this is more\n' ...
0227                  'convenient. The output of ehess should be an ambient\n' ...
0228                  'vector (it will be transformed to a tangent vector\n' ...
0229                  'automatically).\n']);
0230
0231     end
0232
0233     if ~canGetHessian(problem)