Preconditioner based on the inverse Hessian, by solving linear systems. function preconfun = preconhessiansolve(problem) function preconfun = preconhessiansolve(problem, options) Input: A Manopt problem structure (already containing the manifold and enough information to compute the Hessian of the cost) and an options structure (optional, currently ignored). Notice that if the Hessian is not positive definite, then its inverse is not positive definite either and this preconditioner is not suitable. If the Hessian cannot be computed on 'problem', a warning is issued. An approximation of the Hessian will be used instead, and the present preconditioner will attempt to invert that (although it may not be a linear operator). If no approximate Hessian is provided either, a generic approximation is used. Behavior is unspecified. Output: Returns a function handle, encapsulating a generic preconditioner of the Hessian based on solving linear systems of the form: Hessian(x)[preconfun(x, xdot)] = xdot, where x is the point on the manifold, xdot is the input to the preconditioner (a tangent vector) and preconfun(x, xdot) is returned (also a tangent vector). The solve may be approximate. The returned preconfun has this calling pattern: function precxdot = preconfun(x, xdot) function precxdot = preconfun(x, xdot, storedb) function precxdot = preconfun(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, the gradient and the Hessian (typically, problem.cost, problem.grad and problem.hess, or problem.egrad and problem.ehess). Then, to use this generic purpose Hessian preconditioner: problem.precon = preconhessiansolve(problem, options); Passing that problem structure to the conjugategradients solver (which uses preconditioning) configured in steepest descent mode results in a type of Riemannian Newton method. See also: conjugategradients
0001 function preconfun = preconhessiansolve(problem, options) 0002 % Preconditioner based on the inverse Hessian, by solving linear systems. 0003 % 0004 % function preconfun = preconhessiansolve(problem) 0005 % function preconfun = preconhessiansolve(problem, options) 0006 % 0007 % Input: 0008 % 0009 % A Manopt problem structure (already containing the manifold and enough 0010 % information to compute the Hessian of the cost) and an options structure 0011 % (optional, currently ignored). Notice that if the Hessian is not positive 0012 % definite, then its inverse is not positive definite either and this 0013 % preconditioner is not suitable. 0014 % 0015 % If the Hessian cannot be computed on 'problem', a warning is issued. An 0016 % approximation of the Hessian will be used instead, and the present 0017 % preconditioner will attempt to invert that (although it may not be a 0018 % linear operator). If no approximate Hessian is provided either, a generic 0019 % approximation is used. Behavior is unspecified. 0020 % 0021 % Output: 0022 % 0023 % Returns a function handle, encapsulating a generic preconditioner of the 0024 % Hessian based on solving linear systems of the form: 0025 % Hessian(x)[preconfun(x, xdot)] = xdot, 0026 % where x is the point on the manifold, xdot is the input to the 0027 % preconditioner (a tangent vector) and preconfun(x, xdot) is returned 0028 % (also a tangent vector). The solve may be approximate. 0029 % 0030 % The returned preconfun has this calling pattern: 0031 % 0032 % function precxdot = preconfun(x, xdot) 0033 % function precxdot = preconfun(x, xdot, storedb) 0034 % function precxdot = preconfun(x, xdot, storedb, key) 0035 % 0036 % x is a point on the manifold problem.M, xdot is a tangent vector to that 0037 % manifold at x, storedb is a StoreDB object, and key is the StoreDB key to 0038 % point x. 0039 % 0040 % Usage: 0041 % 0042 % Typically, the user will set problem.M and other fields to define the 0043 % cost, the gradient and the Hessian (typically, problem.cost, problem.grad 0044 % and problem.hess, or problem.egrad and problem.ehess). Then, to use this 0045 % generic purpose Hessian preconditioner: 0046 % 0047 % problem.precon = preconhessiansolve(problem, options); 0048 % 0049 % Passing that problem structure to the conjugategradients solver 0050 % (which uses preconditioning) configured in steepest descent mode results 0051 % in a type of Riemannian Newton method. 0052 % 0053 % See also: conjugategradients 0054 0055 % This file is part of Manopt: www.manopt.org. 0056 % Original author: Nicolas Boumal, April 9, 2015. 0057 % Contributors: 0058 % Change log: 0059 0060 % Check availability of the Hessian, or at least of an approximation. 0061 if ~canGetHessian(problem) && ~canGetApproxHessian(problem) 0062 % Note: we do not give a warning if an approximate Hessian is 0063 % explicitly given in the problem description, as in that case the 0064 % user seems to be aware of the issue. 0065 warning('manopt:getHessian:approx', ... 0066 ['No Hessian provided. Using an FD approximation instead.\n' ... 0067 'To disable this warning: warning(''off'', ''manopt:getHessian:approx'')']); 0068 problem.approxhess = approxhessianFD(problem); 0069 end 0070 0071 % Set local defaults here, and merge with user options, if any. 0072 localdefaults = struct(); 0073 if ~exist('options', 'var') || isempty(options) 0074 options = struct(); 0075 end 0076 options = mergeOptions(localdefaults, options); 0077 0078 % Build and return the function handle here. This extra construct via 0079 % funhandle makes it possible to make storedb and key optional. 0080 preconfun = @funhandle; 0081 function precxdot = funhandle(x, xdot, storedb, key) 0082 % Allow omission of the key, and even of storedb. 0083 if ~exist('key', 'var') 0084 if ~exist('storedb', 'var') 0085 storedb = StoreDB(); 0086 end 0087 key = storedb.getNewKey(); 0088 end 0089 precxdot = hessiansolvehelper(options, problem, x, xdot, ... 0090 storedb, key); 0091 end 0092 0093 end 0094 0095 0096 function precxdot = hessiansolvehelper(options, problem, x, xdot, storedb, key) 0097 % This function does the actual work. 0098 0099 % Exclude the case where xdot is zero 0100 norm_xdot = problem.M.norm(x, xdot); 0101 if norm_xdot < eps 0102 precxdot = problem.M.zerovec(x); 0103 return; 0104 end 0105 0106 % Get a shorthand for the Hessian of the cost on M at x. 0107 hessian = @(u) getHessian(problem, x, u, storedb, key); 0108 0109 % Setup an optimization problem on the tangent space to problem.M at x. 0110 M = problem.M; 0111 tgtspace = tangentspacefactory(M, x); 0112 prblm.M = tgtspace; 0113 prblm.cost = @cost; 0114 prblm.grad = @grad; 0115 prblm.hess = @(u, udot) 2*hessian(hessian(udot))/norm_xdot; 0116 0117 function [f, store] = cost(u, store) 0118 if ~isfield(store, 'residue') 0119 Hu = hessian(u); 0120 store.residue = M.lincomb(x, 1, Hu, -1, xdot); 0121 end 0122 f = M.norm(x, store.residue).^2 / norm_xdot; 0123 end 0124 function [g, store] = grad(u, store) 0125 if ~isfield(store, 'residue') 0126 Hu = hessian(u); 0127 store.residue = M.lincomb(x, 1, Hu, -1, xdot); 0128 end 0129 g = 2 * hessian(store.residue) / norm_xdot; 0130 end 0131 0132 % checkgradient(prblm); pause; 0133 % checkhessian(prblm); pause; 0134 0135 localdefaults.solver = @trustregions; 0136 localdefaults.verbosity = 0; 0137 % Merge local defaults with user options, if any. 0138 if ~exist('options', 'var') || isempty(options) 0139 options = struct(); 0140 end 0141 options = mergeOptions(localdefaults, options); 0142 0143 % Solve the linear system by solving the optimization problem. 0144 precxdot = manoptsolve(prblm, M.zerovec(), options); 0145 0146 end