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
0001 function [H, basis] = hessianmatrix(problem, x, basis) 0002 % Computes a matrix which represents the Hessian in some tangent basis. 0003 % 0004 % [H, basis] = hessianmatrix(problem, x) 0005 % [H, basis] = hessianmatrix(problem, x, basis) 0006 % 0007 % problem is a Manopt problem structure with a manifold and cost function. 0008 % x is a point on the manifold problem.M. 0009 % basis (optional) is an orthonormal basis for the tangent space to the 0010 % manifold at x. If no basis is supplied, one will be generated at random. 0011 % If the basis spans only a subspace of the tangent space at x, 0012 % then the returned matrix represents the Hessian restricted to that subspace. 0013 % 0014 % H is an n-by-n symmetric matrix (with n the number of vectors in the basis) 0015 % such that H(i, j) is the inner product between basis{i} 0016 % and Hess(basis{j}), with respect to the metric on the tangent space to 0017 % problem.M at x, where Hess(basis{j}) is the vector obtained after 0018 % applying the Hessian at x to basis{j}. 0019 % 0020 % For optimization, it is usually not useful to compute the Hessian matrix, 0021 % as this quickly becomes expensive. This tool is provided mostly for 0022 % exploration and debugging rather than to be used algorithmically in 0023 % solvers. To access the spectrum of the Hessian, it may be more practical 0024 % to call hessianextreme or hessianspectrum. This should coincide with eig(H). 0025 % 0026 % 0027 % Example of equivalence: 0028 % 0029 % Hu = getHessian(problem, x, u) 0030 % 0031 % is equivalent to (but much faster than): 0032 % 0033 % B = tangentorthobasis(M, x); 0034 % H = hessianmatrix(problem, x, B); 0035 % u_vec = tangent2vec(M, x, B, u); 0036 % Hu_vec = H*u_vec; 0037 % Hu = lincomb(M, x, B, Hu_vec); 0038 % 0039 % Note that there will be some error due to numerical round-off. 0040 % 0041 % 0042 % See also: hessianspectrum hessianextreme tangentorthobasis orthogonalize tangent2vec 0043 0044 % This file is part of Manopt: www.manopt.org. 0045 % Original author: Nicolas Boumal, July 14, 2016. 0046 % Contributors: 0047 % Change log: 0048 0049 % TODO: refactor using operator2matrix 0050 0051 0052 % No warning if an approximate Hessian is available, as then the user 0053 % is presumably aware of what they are doing. 0054 if ~canGetHessian(problem) && ~canGetApproxHessian(problem) 0055 warning('manopt:hessianmatrix:nohessian', ... 0056 ['The Hessian appears to be unavailable.\n' ... 0057 'Will try to use an approximate Hessian instead.\n'... 0058 'Since this approximation may not be linear or '... 0059 'symmetric,\nthe computation might fail and the '... 0060 'results (if any)\nmight make no sense.']); 0061 end 0062 0063 0064 % Unless an orthonormal basis for the tangent space at x is provided, 0065 % pick a random one. 0066 if ~exist('basis', 'var') || isempty(basis) 0067 n = problem.M.dim(); 0068 basis = tangentorthobasis(problem.M, x, n); 0069 else 0070 n = numel(basis); 0071 end 0072 0073 % Create a store database and get a key for x 0074 storedb = StoreDB(1); 0075 key = storedb.getNewKey(); 0076 0077 % Apply the Hessian at x to each basis vector 0078 Hbasis = cell(n, 1); 0079 for k = 1 : numel(Hbasis) 0080 Hbasis{k} = getHessian(problem, x, basis{k}, storedb, key); 0081 end 0082 0083 % H is the matrix which contains the inner products of 0084 % the ((basis vectors)) with the ((Hessian applied to basis vectors)). 0085 H = zeros(n); 0086 for i = 1 : n 0087 H(i, i) = problem.M.inner(x, basis{i}, Hbasis{i}); 0088 for j = (i+1) : n 0089 H(i, j) = problem.M.inner(x, basis{i}, Hbasis{j}); 0090 H(j, i) = H(i, j); 0091 end 0092 end 0093 0094 end