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:
• doubly_stochastic_denoising Find a doubly stochastic matrix closest to a given matrix, in Frobenius norm.
• basic_example_AD A basic example that shows how to apply automatic differentiation to
• complex_example_AD A basic example that shows how to define the cost funtion for
• complextest_AD1 Test AD for a complex optimization problem on a product manifold (struct)
• complextest_AD2 Test AD for a complex optimization problem on a power manifold (cell)
• complextest_AD3 Test AD for a complex optimization problem on a manifold which is stored
• realtest_AD1 Test AD for a real optimization problem on a product manifold (struct)
• realtest_AD2 Test AD for a real optimization problem on a power manifold (cell)
• realtest_AD3 Test AD for a real optimization problem on a manifold which is stored in
• linearsystem_compare Example code for the algorithms described in

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