approxhessianFD

 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
 cost and the gradient (typically, problem.cost and problem.grad or
 problem.egrad). Then, to use this generic purpose Hessian approximation:

   problem.approxhess = approxhessianFD(problem, options);

 See also: trustregions

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
0036 % cost and the gradient (typically, problem.cost and problem.grad or
0037 % problem.egrad). Then, to use this generic purpose Hessian approximation:
0038 %
0039 %   problem.approxhess = approxhessianFD(problem, options);
0040 %
0041 % See also: trustregions
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.
0090     if ~canGetGradient(problem) && ~canGetApproxGradient(problem)
0091         warning('manopt:approxhessianFD:nogradient', ...
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.
0142     grad = getGradient(problem, x, storedb, key);
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();
0148     grad1 = getGradient(problem, x1, storedb, key1);
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.
0155     grad1 = problem.M.transp(x1, x, grad1);
0156     
0157     % Return the finite difference of them: (grad1 - grad)/c.
0158     hessfd = problem.M.lincomb(x, 1/c, grad1, -1/c, grad);
0159     
0160 end