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
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.

## CROSS-REFERENCE INFORMATION This function calls:
• StoreDB
• canGetApproxHessian Checks whether an approximate Hessian can be computed for this problem.
• canGetHessian Checks whether the Hessian can be computed for a problem structure.
• getHessian Computes the Hessian of the cost function at x along d.
• mergeOptions Merges two options structures with one having precedence over the other.
• approxhessianFD Hessian approx. fnctn handle based on finite differences of the gradient.
• trustregions Riemannian trust-regions solver for optimization on manifolds.
• manoptsolve Gateway helper function to call a Manopt solver, chosen in the options.
• tangentspacefactory Returns a manifold structure representing the tangent space to M at x.
This function is called by:

## 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
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 %
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;
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
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 Tue 19-May-2020 18:46:12 by m2html © 2005