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.

 See also: checkdiff checkgradient checkretraction

CROSS-REFERENCE INFORMATION ^

This function calls: 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 % See also: checkdiff checkgradient checkretraction
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 %   April 3, 2015 (NB):
0028 %       Works with the new StoreDB class system.
0029 %
0030 %   Nov. 1, 2016 (NB):
0031 %       Issues a call to getGradient rather than getDirectionalDerivative.
0032 
0033         
0034     % Verify that the problem description is sufficient.
0035     if ~canGetCost(problem)
0036         error('It seems no cost was provided.');
0037     end
0038     if ~canGetGradient(problem)
0039         warning('manopt:checkhessian:nograd', ...
0040                 'It seems no gradient was provided.');
0041     end
0042     if ~canGetHessian(problem)
0043         warning('manopt:checkhessian:nohess', ...
0044                 'It seems no Hessian was provided.');
0045     end
0046     
0047     x_isprovided = exist('x', 'var') && ~isempty(x);
0048     d_isprovided = exist('d', 'var') && ~isempty(d);
0049     
0050     if ~x_isprovided && d_isprovided
0051         error('If d is provided, x must be too, since d is tangent at x.');
0052     end
0053     
0054     % If x and / or d are not specified, pick them at random.
0055     if ~x_isprovided
0056         x = problem.M.rand();
0057     end
0058     if ~d_isprovided
0059         d = problem.M.randvec(x);
0060     end
0061     
0062     %% Check that the directional derivative and the Hessian at x along d
0063     %% yield a second order model of the cost function.
0064     
0065     % Compute the value f0 at f, directional derivative df0 at x along d,
0066     % and Hessian along [d, d].
0067     storedb = StoreDB();
0068     xkey = storedb.getNewKey();
0069     f0 = getCost(problem, x, storedb, xkey);
0070     df0 = problem.M.inner(x, d, getGradient(problem, x, storedb, xkey));
0071     d2f0 = problem.M.inner(x, d, getHessian(problem, x, d, storedb, xkey));
0072     
0073     % Compute the value of f at points on the geodesic (or approximation
0074     % of it) originating from x, along direction d, for stepsizes in a
0075     % large range given by h.
0076     h = logspace(-8, 0, 51);
0077     value = zeros(size(h));
0078     for i = 1 : length(h)
0079         y = problem.M.exp(x, d, h(i));
0080         ykey = storedb.getNewKey();
0081         value(i) = getCost(problem, y, storedb, ykey);
0082     end
0083     
0084     % Compute the quadratic approximation of the cost function using f0,
0085     % df0 and d2f0 at the same points.
0086     model = polyval([.5*d2f0 df0 f0], h);
0087     
0088     % Compute the approximation error
0089     err = abs(model - value);
0090     
0091     % And plot it.
0092     loglog(h, err);
0093     title(sprintf(['Hessian check.\nThe slope of the continuous line ' ...
0094                    'should match that of the dashed\n(reference) line ' ...
0095                    'over at least a few orders of magnitude for h.']));
0096     xlabel('h');
0097     ylabel('Approximation error');
0098     
0099     line('xdata', [1e-8 1e0], 'ydata', [1e-16 1e8], ...
0100          'color', 'k', 'LineStyle', '--', ...
0101          'YLimInclude', 'off', 'XLimInclude', 'off');
0102     
0103     % In a numerically reasonable neighborhood, the error should decrease
0104     % as the cube of the stepsize, i.e., in loglog scale, the error
0105     % should have a slope of 3.
0106     window_len = 10;
0107     [range, poly] = identify_linear_piece(log10(h), log10(err), window_len);
0108     hold all;
0109     loglog(h(range), 10.^polyval(poly, log10(h(range))), 'LineWidth', 3);
0110     hold off;
0111     
0112     fprintf('The slope should be 3. It appears to be: %g.\n', poly(1));
0113     fprintf(['If it is far from 3, then directional derivatives or ' ...
0114              'the Hessian might be erroneous.\n']);
0115     fprintf(['Note: if the exponential map is only approximate, and it '...
0116              'is not a second-order approximation,\nthen it is normal ' ...
0117              'for the slope test to reach 2 instead of 3. Check the ' ...
0118              'factory for this.\n' ...
0119              'If tested at a critical point, then even for a first-order '...
0120              'retraction the slope test should yield 3.\n']);
0121 
0122     
0123     %% Check that the Hessian at x along direction d is a tangent vector.
0124     if isfield(problem.M, 'tangent')
0125         hess = getHessian(problem, x, d, storedb, xkey);
0126         phess = problem.M.tangent(x, hess);
0127         residual = problem.M.lincomb(x, 1, hess, -1, phess);
0128         err = problem.M.norm(x, residual);
0129         fprintf('The residual should be zero, or very close. ');
0130         fprintf('Residual: %g.\n', err);
0131         fprintf(['If it is far from 0, then the Hessian is not in the ' ...
0132                  'tangent plane.\n']);
0133     else
0134         fprintf(['Unfortunately, Manopt was unable to verify that the '...
0135                  'Hessian is indeed a tangent vector.\nPlease verify ' ...
0136                  'this manually.']);
0137     end    
0138     
0139     %% Check that the Hessian at x is symmetric.
0140     d1 = problem.M.randvec(x);
0141     d2 = problem.M.randvec(x);
0142     h1 = getHessian(problem, x, d1, storedb, xkey);
0143     h2 = getHessian(problem, x, d2, storedb, xkey);
0144     v1 = problem.M.inner(x, d1, h2);
0145     v2 = problem.M.inner(x, h1, d2);
0146     value = v1-v2;
0147     fprintf(['<d1, H[d2]> - <H[d1], d2> should be zero, or very close.' ...
0148              '\n\tValue: %g - %g = %g.\n'], v1, v2, value);
0149     fprintf('If it is far from 0, then the Hessian is not symmetric.\n');
0150     
0151 end

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005