function [H, basis] = hessianmatrix(problem, x, basis)
% Computes a matrix which represents the Hessian in some tangent basis.
%
% [H, basis] = hessianmatrix(problem, x)
% [H, basis] = hessianmatrix(problem, x, basis)
%
% problem is a Manopt problem structure with a manifold and cost function.
% x is a point on the manifold problem.M.
% basis (optional) is an orthonormal basis for the tangent space to the
% manifold at x. If no basis is supplied, one will be generated at random.
% If the basis spans only a subspace of the tangent space at x,
% then the returned matrix represents the Hessian restricted to that subspace.
%
% H is an n-by-n symmetric matrix (with n the number of vectors in the basis)
% such that H(i, j) is the inner product between basis{i}
% and Hess(basis{j}), with respect to the metric on the tangent space to
% problem.M at x, where Hess(basis{j}) is the vector obtained after
% applying the Hessian at x to basis{j}.
%
% For optimization, it is usually not useful to compute the Hessian matrix,
% as this quickly becomes expensive. This tool is provided mostly for
% exploration and debugging rather than to be used algorithmically in
% solvers. To access the spectrum of the Hessian, it may be more practical
% to call hessianextreme or hessianspectrum. This should coincide with eig(H).
%
%
% Example of equivalence:
%
%     Hu = getHessian(problem, x, u)
%
% is equivalent to (but much faster than):
%
%     B = tangentorthobasis(M, x);
%     H = hessianmatrix(problem, x, B);
%     u_vec = tangent2vec(M, x, B, u);
%     Hu_vec = H*u_vec;
%     Hu = lincomb(M, x, B, Hu_vec);
%
% Note that there will be some error due to numerical round-off.
% 
%
% See also: hessianspectrum hessianextreme tangentorthobasis orthogonalize tangent2vec

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, July 14, 2016.
% Contributors: 
% Change log: 

% TODO: refactor using operator2matrix


    % No warning if an approximate Hessian is available, as then the user
    % is presumably aware of what they are doing.
    if ~canGetHessian(problem) && ~canGetApproxHessian(problem)
        warning('manopt:hessianmatrix:nohessian', ...
                ['The Hessian appears to be unavailable.\n' ...
                 'Will try to use an approximate Hessian instead.\n'...
                 'Since this approximation may not be linear or '...
                 'symmetric,\nthe computation might fail and the '...
                 'results (if any)\nmight make no sense.']);
    end
    

    % Unless an orthonormal basis for the tangent space at x is provided,
    % pick a random one.
    if ~exist('basis', 'var') || isempty(basis)
        n = problem.M.dim();
        basis = tangentorthobasis(problem.M, x, n);
    else
        n = numel(basis);
    end
    
    % Create a store database and get a key for x
    storedb = StoreDB(1);
    key = storedb.getNewKey();
    
    % Apply the Hessian at x to each basis vector
    Hbasis = cell(n, 1);
    for k = 1 : numel(Hbasis)
        Hbasis{k} = getHessian(problem, x, basis{k}, storedb, key);
    end
    
    % H is the matrix which contains the inner products of
    % the ((basis vectors)) with the ((Hessian applied to basis vectors)).
    H = zeros(n);
    for i = 1 : n
        H(i, i) = problem.M.inner(x, basis{i}, Hbasis{i});
        for j = (i+1) : n
            H(i, j) = problem.M.inner(x, basis{i}, Hbasis{j});
            H(j, i) = H(i, j);
        end
    end

end
