Home > manopt > solvers > preconditioners > preconhessiansolve.m

preconhessiansolve

PURPOSE ^

Preconditioner based on the inverse Hessian, by solving linear systems.

SYNOPSIS ^

function preconfun = preconhessiansolve(problem, options)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005