Home > manopt > tools > hessianspectrum.m

# hessianspectrum

## PURPOSE

Returns the eigenvalues of the (preconditioned) Hessian at x.

## SYNOPSIS

function lambdas = hessianspectrum(problem, x, usepreconstr, storedb, key)

## DESCRIPTION

``` Returns the eigenvalues of the (preconditioned) Hessian at x.

function lambdas = hessianspectrum(problem, x)
function lambdas = hessianspectrum(problem, x, useprecon)
function lambdas = hessianspectrum(problem, x, useprecon, storedb)
function lambdas = hessianspectrum(problem, x, useprecon, storedb, key)

If useprecon is not set, or if it is set to 'noprecon' (default), this
computes and returns the eigenvalues of the Hessian operator (which needs
to be symmetric but not necessarily definite) on the tangent space at x.
There are problem.M.dim() eigenvalues. Matlab's eigs is used internally.

If useprecon is set to 'precon', the eigenvalues of the composition of
the Hessian with the preconditioner at x are computed: Precon o Hessian.
The preconditioner must have been defined in the problem structure and
has to be symmetric, positive definite. It is supposed to approximate the
inverse of the (Riemannian) Hessian. Ideally, the preconditioned Hessian
is better conditioned (smaller ratio of largest to smallest eigenvalue in
magnitude) than the non-preconditioned spectrum. The present tool can
help assess that.

The typical ways to define a preconditioner are via problem.precon or
problem.sqrtprecon (see comment below). These should be function handles
with the same input/output system as problem.hess for the Hessian.

If the Hessian is not available from the problem structure, an
approximate Hessian will be used. There are no guarantees of
interpretability, but this may nevertheless be useful at times.

Even though the Hessian and the preconditioner are both symmetric, their
composition is not symmetric. This can slow down the call to 'eigs'
substantially. If possible, you may specify the square root of the
preconditioner in the problem structure, as sqrtprecon. This operator on
the tangent space at x must also be symmetric, positive definite, and
such that SqrtPrecon o SqrtPrecon = Precon. Then, the spectrum of the
symmetric operator SqrtPrecon o Hessian o SqrtPrecon is computed: it is
the same as the spectrum of Precon o Hessian, but is usually faster to
compute. If both Precon and SqrtPrecon are provided, only SqrtPrecon will
be used.

The input and the output of the Hessian and of the preconditioner are
projected on the tangent space to avoid undesired contributions of the
ambient space.

storedb is a StoreDB object, key is the StoreDB key to point x.

Requires the manifold description in problem.M to have these functions:

u_vec = vec(x, u_mat) :
Returns a column vector representation of the normal (usually
matrix) representation of the tangent vector u_mat. vec must be an
isometry between the tangent space (with its Riemannian metric) and
a subspace of R^n where n = length(u_vec), with the 2-norm on R^n.
In other words: it is an orthogonal projector.

u_mat = mat(x, u_vec) :
The inverse of vec (its adjoint).

u_mat_clean = tangent(x, u_mat) :
Subtracts from the tangent vector u_mat any component that would
make it "not really tangent", by projection.

answer = vecmatareisometries() :
Returns true if the linear maps encoded by vec and mat are
isometries, false otherwise. It is better if the answer is yes.

## CROSS-REFERENCE INFORMATION

This function calls:
• StoreDB
• canGetApproxHessian Checks whether an approximate Hessian can be computed for this problem.
• canGetGradient Checks whether the gradient can be computed for a problem structure.
• canGetHessian Checks whether the Hessian can be computed for a problem structure.
• canGetPrecon Checks whether a preconditioner was specified in the problem description.
• canGetSqrtPrecon Checks whether a square root of preconditioner was specified in problem.
• getCostGrad Computes the cost function and the gradient at x in one call if possible.
• getHessian Computes the Hessian of the cost function at x along d.
• getPrecon Applies the preconditioner for the Hessian of the cost at x along d.
• getSqrtPrecon Applies the square root of the Hessian preconditioner at x along d.
This function is called by:

## SOURCE CODE

```0001 function lambdas = hessianspectrum(problem, x, usepreconstr, storedb, key)
0002 % Returns the eigenvalues of the (preconditioned) Hessian at x.
0003 %
0004 % function lambdas = hessianspectrum(problem, x)
0005 % function lambdas = hessianspectrum(problem, x, useprecon)
0006 % function lambdas = hessianspectrum(problem, x, useprecon, storedb)
0007 % function lambdas = hessianspectrum(problem, x, useprecon, storedb, key)
0008 %
0009 % If useprecon is not set, or if it is set to 'noprecon' (default), this
0010 % computes and returns the eigenvalues of the Hessian operator (which needs
0011 % to be symmetric but not necessarily definite) on the tangent space at x.
0012 % There are problem.M.dim() eigenvalues. Matlab's eigs is used internally.
0013 %
0014 % If useprecon is set to 'precon', the eigenvalues of the composition of
0015 % the Hessian with the preconditioner at x are computed: Precon o Hessian.
0016 % The preconditioner must have been defined in the problem structure and
0017 % has to be symmetric, positive definite. It is supposed to approximate the
0018 % inverse of the (Riemannian) Hessian. Ideally, the preconditioned Hessian
0019 % is better conditioned (smaller ratio of largest to smallest eigenvalue in
0020 % magnitude) than the non-preconditioned spectrum. The present tool can
0021 % help assess that.
0022 %
0023 % The typical ways to define a preconditioner are via problem.precon or
0024 % problem.sqrtprecon (see comment below). These should be function handles
0025 % with the same input/output system as problem.hess for the Hessian.
0026 %
0027 % If the Hessian is not available from the problem structure, an
0028 % approximate Hessian will be used. There are no guarantees of
0029 % interpretability, but this may nevertheless be useful at times.
0030 %
0031 % Even though the Hessian and the preconditioner are both symmetric, their
0032 % composition is not symmetric. This can slow down the call to 'eigs'
0033 % substantially. If possible, you may specify the square root of the
0034 % preconditioner in the problem structure, as sqrtprecon. This operator on
0035 % the tangent space at x must also be symmetric, positive definite, and
0036 % such that SqrtPrecon o SqrtPrecon = Precon. Then, the spectrum of the
0037 % symmetric operator SqrtPrecon o Hessian o SqrtPrecon is computed: it is
0038 % the same as the spectrum of Precon o Hessian, but is usually faster to
0039 % compute. If both Precon and SqrtPrecon are provided, only SqrtPrecon will
0040 % be used.
0041 %
0042 % The input and the output of the Hessian and of the preconditioner are
0043 % projected on the tangent space to avoid undesired contributions of the
0044 % ambient space.
0045 %
0046 % storedb is a StoreDB object, key is the StoreDB key to point x.
0047 %
0048 % Requires the manifold description in problem.M to have these functions:
0049 %
0050 %   u_vec = vec(x, u_mat) :
0051 %       Returns a column vector representation of the normal (usually
0052 %       matrix) representation of the tangent vector u_mat. vec must be an
0053 %       isometry between the tangent space (with its Riemannian metric) and
0054 %       a subspace of R^n where n = length(u_vec), with the 2-norm on R^n.
0055 %       In other words: it is an orthogonal projector.
0056 %
0057 %   u_mat = mat(x, u_vec) :
0058 %       The inverse of vec (its adjoint).
0059 %
0060 %   u_mat_clean = tangent(x, u_mat) :
0061 %       Subtracts from the tangent vector u_mat any component that would
0062 %       make it "not really tangent", by projection.
0063 %
0064 %   answer = vecmatareisometries() :
0065 %       Returns true if the linear maps encoded by vec and mat are
0066 %       isometries, false otherwise. It is better if the answer is yes.
0067 %
0068 % See also: hessianextreme canGetPrecon canGetSqrtPrecon
0069
0070 % This file is part of Manopt: www.manopt.org.
0071 % Original author: Nicolas Boumal, July 3, 2013.
0072 % Contributors:
0073 % Change log:
0074 %
0075 %   Dec. 18, 2014 (NB):
0076 %       The lambdas are now sorted when they are returned.
0077 %
0078 %   April 3, 2015 (NB):
0079 %       Works with the new StoreDB class system.
0080 %       Does no longer accept sqrtprecon as an input: the square root of
0081 %       the preconditioner may now be specified directly in the problem
0082 %       structure, following the same syntax as the preconditioner precon.
0083 %
0084 %   April 4, 2015 (NB):
0085 %       By default, the spectrum is computed without the preconditioner's
0086 %       effect, even if it is available. A new input option allows to
0087 %       switch this behavior without the need to change the problem
0088 %       structure.
0089
0090     % Allow omission of the key, and even of storedb.
0091     if ~exist('key', 'var')
0092         if ~exist('storedb', 'var')
0093             storedb = StoreDB();
0094         end
0095         key = storedb.getNewKey();
0096     end
0097
0098     % Manage the option to use or not use a preconditioner.
0099     % The input is a string. It is here transformed into a Boolean.
0100     if ~exist('usepreconstr', 'var') || isempty(usepreconstr)
0101         usepreconstr = 'noprecon';
0102     end
0103     switch lower(usepreconstr)
0104         case 'noprecon'
0105             useprecon = false;
0106         case 'precon'
0107             useprecon = true;
0108         otherwise
0109             % A bit of legacy code heads up.
0110             if isa(usepreconstr, 'function_handle')
0111                 warning('manopt:hessianspectrum:oldsyntax', ...
0112                         ['This function no longer expects sqrtprecon ' ...
0113                          'as input. Place it in the problem structure.']);
0114             end
0115             error('Input useprecon must be either ''precon'' or ''noprecon''.');
0116     end
0117
0118     % No warning if an approximate Hessian is available, as then the user
0119     % is presumably aware of what they are doing.
0120     if ~canGetHessian(problem) && ~canGetApproxHessian(problem)
0121         warning('manopt:hessianspectrum:nohessian', ...
0122                 ['The Hessian appears to be unavailable.\n' ...
0123                  'Will try to use an approximate Hessian instead.\n'...
0124                  'Since this approximation may not be linear or '...
0125                  'symmetric,\nthe computation might fail and the '...
0126                  'results (if any)\nmight make no sense.']);
0127     end
0128
0129     vec = @(u_mat) problem.M.vec(x, u_mat);
0130     mat = @(u_vec) problem.M.mat(x, u_vec);
0131     tgt = @(u_mat) problem.M.tangent(x, u_mat);
0132
0133     % n: size of a vectorized tangent vector
0134     % dim: dimension of the tangent space
0135     % necessarily, n >= dim.
0136     % The vectorized operators we build below will have at least n - dim
0137     % zero eigenvalues.
0138     n = length(vec(problem.M.zerovec(x)));
0139     dim = problem.M.dim();
0140
0141     % It is usually a good idea to force a gradient computation to make
0142     % sure precomputable things are precomputed.
0144         [unused1, unused2] = getCostGrad(problem, x, storedb, key); %#ok
0145     end
0146
0147     hess = @(u_mat) tgt(getHessian(problem, x, tgt(u_mat), storedb, key));
0148     hess_vec = @(u_vec) vec(hess(mat(u_vec)));
0149
0150     % Regardless of preconditioning, we can only have a symmetric
0151     % eigenvalue problem if the vec/mat pair of the manifold is an
0152     % isometry:
0153     vec_mat_are_isometries = problem.M.vecmatareisometries();
0154
0155
0156     if ~useprecon
0157
0158         % No preconditioner to use: simply use the Hessian as is.
0159
0160         eigs_opts.issym = vec_mat_are_isometries;
0161         eigs_opts.isreal = true;
0162         lambdas = eigs(hess_vec, n, dim, 'LM', eigs_opts);
0163
0164     elseif canGetSqrtPrecon(problem)
0165
0166         % There is a preconditioner, and we have its square root: deal with
0167         % the symmetric composition SqrtPrecon o Hessian o SqrtPrecon.
0168
0169         sqrtprec = @(u_mat) tgt(getSqrtPrecon(problem, x, tgt(u_mat), storedb, key));
0170         sqrtprec_vec = @(u_vec) vec(sqrtprec(mat(u_vec)));
0171
0172         eigs_opts.issym = vec_mat_are_isometries;
0173         eigs_opts.isreal = true;
0174         lambdas = eigs(@(u_vec) ...
0175                       sqrtprec_vec(hess_vec(sqrtprec_vec(u_vec))), ...
0176                       n, dim, 'LM', eigs_opts);
0177
0178     elseif canGetPrecon(problem)
0179
0180         % There is a preconditioner, but we don't have its square root:
0181         % deal with the non-symmetric composition Precon o Hessian.
0182
0183         prec = @(u_mat) tgt(getPrecon(problem, x, tgt(u_mat), storedb, key));
0184         prec_vec = @(u_vec) vec(prec(mat(u_vec)));
0185         % prec_inv_vec = @(u_vec) pcg(prec_vec, u_vec);
0186
0187         eigs_opts.issym = false;
0188         eigs_opts.isreal = true;
0189         lambdas = eigs(@(u_vec) prec_vec(hess_vec(u_vec)), ...
0190                        n, dim, 'LM', eigs_opts);
0191
0192     else
0193
0194         error('No preconditioner is available in the problem structure.');
0195
0196     end
0197
0198     lambdas = sort(lambdas);
0199
0200 end```

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005