function problem_critpt = criticalpointfinder(problem)
% Creates a Manopt problem whose optima are the critical points of another.
%
% problem_critpt = criticalpointfinder(problem)
%
% Given a Manopt problem structure 'problem', this tool returns a new
% problem structure, 'problem_critpt', such that the global optima of the
% new problem coincide with the critical points of the original problem.
% This can be useful notably in empirical studies of the properties of
% saddle points of a problem.
%
% Concretely, if f is the cost function of the given problem, grad f
% denotes its (Riemannian) gradient and Hess f denotes its (Riemannian)
% Hessian, then the new problem has a cost function g defined by:
%
%   g(x) = (1/2)*norm(grad f(x))^2,
%
% where x is a point on the manifold problem.M (the new problem lives on
% the same manifold), and norm(.) = problem.M.norm(x, .) is the Riemannian
% norm on the tangent space at x. The Riemannian gradient of g is elegantly
% obtained from knowledge of f:
%
%   grad g(x) = Hess f(x)[grad f(x)]
%
% If the Hessian of f is not available in the given problem, Manopt will
% approximate it automatically to compute an approximate gradient of g.
% If the Hessian of f is available, then an approximate Hessian of g is
% defined in the returned problem as
%
%  approxhess g(x)[u] = Hess f(x)[ Hess f(x)[u] ].
%
% This approximation is exact if x is a critical point of f, which is
% enough to ensure superlinear local convergence to critical points of f
% using the trustregions algorithm, for example.
%
% Once problem_critpt is obtained, it can be passed to any of the solvers
% of Manopt to compute critical points of the original problem. Supplying
% an initial point to the solver allows to aim for a critical point in a
% specific neighborhood of the search space.
%
%
% Usage example:
% 
% The code below creates a problem whose optima are dominant eigenvectors
% of a matrix A and whose critical points are any eigenvectors of A, then
% compute critical points using the present tool:
%
% n = 100; A = randn(n); A = .5*(A+A');
% problem.M = spherefactory(n);
% problem.cost  = @(x) -x'*(A*x);
% problem.egrad = @(x) -2*A*x;
% problem.ehess = @(x, xdot) -2*A*xdot;
% problem_critpt = criticalpointfinder(problem);
% opts.tolcost = .5*(1e-5)^2; % aim for a gradient smaller than 1e-5
% [x, fx] = trustregions(problem_critpt, [], opts); % random initial guess
% fprintf('Norm of the gradient at x: %g\n', sqrt(2*fx));
% fprintf('This is small if x is close to being an eigenvector: %g\n',...
%         norm((x'*A*x)*x - A*x));
% % The two displayed numbers are equal up to a factor 2.
%
%
% See also: trustregions

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Jan. 25, 2017.
% Contributors: 
% Change log: 

% TODO: Determine a safe way of using the caching functionalities of Manopt
%       with this tool. The issue in passing along storedb and key in the
%       costgrad and approxhess functions is that the storedb will be
%       associated to problem_critpt, not to problem. This may cause bugs
%       that would be very difficult to catch. To be on the safe side,
%       caching is not used at all here, but this may cause running times
%       to be longer than necessary. To create a local storedb associated
%       to problem and to only use the key seems to also not be a viable
%       solution, since there is no clear way of resetting it to zero
%       everytime a solver is called on problem_critpt.
%       -- Jan. 26, 2017 (NB)

    problem_critpt.M = problem.M;
    problem_critpt.costgrad = @costgrad;
    
    % If the Hessian is available for the problem, we build an approximate
    % Hessian based on it. Otherwise, there is no reason to believe that
    % this approximate Hessian would be better than the standard
    % approximate Hessian created by Manopt.
    if canGetHessian(problem)
        problem_critpt.approxhess = @approxhess;
    end
    
    function [g, gradg] = costgrad(x)
        
        gradf = getGradient(problem, x);
        Hessf_gradf = getHessian(problem, x, gradf);
        
        g = .5*problem.M.norm(x, gradf)^2;
        gradg = Hessf_gradf;
        
    end
    
    % This is not quite the Hessian because there should be a third-order
    % derivative term (which is inaccessible), but: at critical points
    % (where grad f(x) = 0 for the f of problem.cost) this Hessian is
    % exact, so it will allow for superlinear local convergence in
    % algorithms such as trustregions.
    function HHu = approxhess(x, u)
        
        Hu  = getHessian(problem, x, u);
        HHu = getHessian(problem, x, Hu);
        
    end

end
