Home > manopt > solvers > hessianapproximations > approxhessianFD.m

# approxhessianFD

## PURPOSE Hessian approx. fnctn handle based on finite differences of the gradient.

## SYNOPSIS function hessfun = approxhessianFD(problem, options)

## DESCRIPTION ``` Hessian approx. fnctn handle based on finite differences of the gradient.

function hessfun = approxhessianFD(problem)
function hessfun = approxhessianFD(problem, options)

Input:

A Manopt problem structure (already containing the manifold and enough
information to compute the cost gradient) and an options structure
(optional), containing one option:
options.stepsize (positive double; default: 2^-14).

If the gradient cannot be computed or approximated on 'problem',
a warning is issued.

Output:

Returns a function handle, encapsulating a generic finite difference
approximation of the Hessian of the problem cost. The finite difference
is based on computations of the gradient.

The returned hessfun has this calling pattern:

function hessfd = hessfun(x, xdot)
function hessfd = hessfun(x, xdot, storedb)
function hessfd = hessfun(x, xdot, storedb, key)

x is a point on the manifold problem.M, xdot is a tangent vector to that
manifold at x, storedb is a StoreDB object, and key is the StoreDB key to
point x.

Usage:

Typically, the user will set problem.M and other fields to define the
problem.egrad). Then, to use this generic purpose Hessian approximation:

problem.approxhess = approxhessianFD(problem, options);

## CROSS-REFERENCE INFORMATION This function calls:
• StoreDB
• canGetApproxGradient Checks whether an approximate gradient can be computed for this problem.
• canGetGradient Checks whether the gradient can be computed for a problem structure.
• mergeOptions Merges two options structures with one having precedence over the other.
This function is called by:
• positive_definite_karcher_mean Computes a Karcher mean of a collection of positive definite matrices.
• arc Adaptive regularization by cubics (ARC) minimization algorithm for Manopt
• preconhessiansolve Preconditioner based on the inverse Hessian, by solving linear systems.
• trustregions Riemannian trust-regions solver for optimization on manifolds.

## SOURCE CODE ```0001 function hessfun = approxhessianFD(problem, options)
0002 % Hessian approx. fnctn handle based on finite differences of the gradient.
0003 %
0004 % function hessfun = approxhessianFD(problem)
0005 % function hessfun = approxhessianFD(problem, options)
0006 %
0007 % Input:
0008 %
0009 % A Manopt problem structure (already containing the manifold and enough
0010 % information to compute the cost gradient) and an options structure
0011 % (optional), containing one option:
0012 %    options.stepsize (positive double; default: 2^-14).
0013 %
0014 % If the gradient cannot be computed or approximated on 'problem',
0015 % a warning is issued.
0016 %
0017 % Output:
0018 %
0019 % Returns a function handle, encapsulating a generic finite difference
0020 % approximation of the Hessian of the problem cost. The finite difference
0021 % is based on computations of the gradient.
0022 %
0023 % The returned hessfun has this calling pattern:
0024 %
0025 %   function hessfd = hessfun(x, xdot)
0026 %   function hessfd = hessfun(x, xdot, storedb)
0027 %   function hessfd = hessfun(x, xdot, storedb, key)
0028 %
0029 % x is a point on the manifold problem.M, xdot is a tangent vector to that
0030 % manifold at x, storedb is a StoreDB object, and key is the StoreDB key to
0031 % point x.
0032 %
0033 % Usage:
0034 %
0035 % Typically, the user will set problem.M and other fields to define the
0037 % problem.egrad). Then, to use this generic purpose Hessian approximation:
0038 %
0039 %   problem.approxhess = approxhessianFD(problem, options);
0040 %
0042
0043 % The Riemannian Trust-Region method, used in combination with the present
0044 % Hessian approximation, is called RTR-FD. Some convergence theory for it
0045 % is available in this paper:
0046 %
0047 % @incollection{boumal2015rtrfd
0048 %   author={Boumal, N.},
0049 %   title={Riemannian trust regions with finite-difference Hessian approximations are globally convergent},
0050 %   year={2015},
0051 %   booktitle={Geometric Science of Information}
0052 % }
0053
0054 % This file is part of Manopt: www.manopt.org.
0055 % Original author: Nicolas Boumal, April 8, 2015.
0056 % Contributors:
0057 % Change log:
0058 %
0059 %   Feb. 19, 2015 (NB):
0060 %       It is sufficient to ensure positive radial linearity to guarantee
0061 %       (together with other assumptions) that this approximation of the
0062 %       Hessian will confer global convergence to the trust-regions method.
0063 %       Formerly, in-code comments referred to the necessity of having
0064 %       complete radial linearity, and that this was harder to achieve.
0065 %       This appears not to be necessary after all, which simplifies the
0066 %       code.
0067 %
0068 %   April 3, 2015 (NB):
0069 %       Works with the new StoreDB class system.
0070 %
0071 %   April 8, 2015 (NB):
0072 %       Changed to approxhessianFD, which now returns a function handle
0073 %       that encapsulates the getHessianFD functionality. Will be better
0074 %       aligned with the other Hessian approximations to come (which may
0075 %       want to use storedb.internal), and now allows specifying the step
0076 %       size.
0077 %
0078 %   March 17, 2020 (NB):
0079 %       Following a bug report by Marco Sutti, added the instruction
0080 %           storedb.remove(key1);
0081 %       to avoid memory usage ramping up when many inner iterations are
0082 %       needed inside of tCG for trustregions. The deciding factor is that
0083 %       there is no need to cache the gradient at the artificially produced
0084 %       point used here for finite differencing, as this point is not
0085 %       visible outside of this function: there is no reason we would visit
0086 %       it again.
0087
0088     % This Hessian approximation is based on the gradient:
0089     % check availability.
0092                 'approxhessianFD requires the gradient to be computable.');
0093     end
0094
0095     % Set local defaults here, and merge with user options, if any.
0096     localdefaults.stepsize = 2^-14;
0097     if ~exist('options', 'var') || isempty(options)
0098         options = struct();
0099     end
0100     options = mergeOptions(localdefaults, options);
0101
0102     % Finite-difference parameter: how far do we look?
0103     stepsize = options.stepsize;
0104
0105     % Build and return the function handle here. This extra construct via
0106     % funhandle makes it possible to make storedb and key optional.
0107     hessfun = @funhandle;
0108     function hessfd = funhandle(x, xdot, storedb, key)
0109         % Allow omission of the key, and even of storedb.
0110         if ~exist('key', 'var')
0111             if ~exist('storedb', 'var')
0112                 storedb = StoreDB();
0113             end
0114             key = storedb.getNewKey();
0115         end
0116         hessfd = hessianFD(stepsize, problem, x, xdot, storedb, key);
0117     end
0118
0119 end
0120
0121
0122 function hessfd = hessianFD(stepsize, problem, x, xdot, storedb, key)
0123 % This function does the actual work.
0124 %
0125 % Original code: Dec. 30, 2012 (NB).
0126
0127     % Extract the input vector norm.
0128     norm_xdot = problem.M.norm(x, xdot);
0129
0130     % First, check whether the step xdot is not too small.
0131     if norm_xdot < eps
0132         hessfd = problem.M.zerovec(x);
0133         return;
0134     end
0135
0136     % Determine how far to retract xdot, so that the point reached does not
0137     % depend on the norm of xdot. This is what ensures radial linearity of
0138     % this present Hessian approximation.
0139     c = stepsize / norm_xdot;
0140
0141     % Compute the gradient at the current point.
0143
0144     % Compute a point a little further along xdot, and the gradient there.
0145     % Since this is a new point, we need a new key for it, for storedb.
0146     x1 = problem.M.retr(x, xdot, c);
0147     key1 = storedb.getNewKey();
0149
0150     % Remove any caching associated to that new point, since there is no
0151     % reason we would be visiting it again.
0152     storedb.remove(key1);
0153
0154     % Transport grad1 back from x1 to x.