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.

 See also: hessianextreme canGetPrecon canGetSqrtPrecon

CROSS-REFERENCE INFORMATION ^

This function calls: 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     if ~canGetHessian(problem)
0119         warning('manopt:hessianspectrum:nohessian', ...
0120                 ['The Hessian appears to be unavailable.\n' ...
0121                  'Will try to use an approximate Hessian instead.\n'...
0122                  'Since this approximation may not be linear or '...
0123                  'symmetric,\nthe computation might fail and the '...
0124                  'results (if any)\nmight make no sense.']);
0125     end
0126 
0127     vec = @(u_mat) problem.M.vec(x, u_mat);
0128     mat = @(u_vec) problem.M.mat(x, u_vec);
0129     tgt = @(u_mat) problem.M.tangent(x, u_mat);
0130     
0131     % n: size of a vectorized tangent vector
0132     % dim: dimension of the tangent space
0133     % necessarily, n >= dim.
0134     % The vectorized operators we build below will have at least n - dim
0135     % zero eigenvalues.
0136     n = length(vec(problem.M.zerovec(x)));
0137     dim = problem.M.dim();
0138     
0139     % It is usually a good idea to force a gradient computation to make
0140     % sure precomputable things are precomputed.
0141     if canGetGradient(problem)
0142         [unused1, unused2] = getCostGrad(problem, x, storedb, key); %#ok
0143     end
0144     
0145     hess = @(u_mat) tgt(getHessian(problem, x, tgt(u_mat), storedb, key));
0146     hess_vec = @(u_vec) vec(hess(mat(u_vec)));
0147     
0148     % Regardless of preconditioning, we can only have a symmetric
0149     % eigenvalue problem if the vec/mat pair of the manifold is an
0150     % isometry:
0151     vec_mat_are_isometries = problem.M.vecmatareisometries();
0152     
0153     
0154     if ~useprecon
0155 
0156         % No preconditioner to use: simply use the Hessian as is.
0157 
0158         eigs_opts.issym = vec_mat_are_isometries;
0159         eigs_opts.isreal = true;
0160         lambdas = eigs(hess_vec, n, dim, 'LM', eigs_opts);
0161             
0162     elseif canGetSqrtPrecon(problem)
0163 
0164         % There is a preconditioner, and we have its square root: deal with
0165         % the symmetric composition SqrtPrecon o Hessian o SqrtPrecon.
0166 
0167         sqrtprec = @(u_mat) tgt(getSqrtPrecon(problem, x, tgt(u_mat), storedb, key));
0168         sqrtprec_vec = @(u_vec) vec(sqrtprec(mat(u_vec)));
0169 
0170         eigs_opts.issym = vec_mat_are_isometries;
0171         eigs_opts.isreal = true;
0172         lambdas = eigs(@(u_vec) ...
0173                       sqrtprec_vec(hess_vec(sqrtprec_vec(u_vec))), ...
0174                       n, dim, 'LM', eigs_opts);
0175             
0176     elseif canGetPrecon(problem)
0177             
0178         % There is a preconditioner, but we don't have its square root:
0179         % deal with the non-symmetric composition Precon o Hessian.
0180 
0181         prec = @(u_mat) tgt(getPrecon(problem, x, tgt(u_mat), storedb, key));
0182         prec_vec = @(u_vec) vec(prec(mat(u_vec)));
0183         % prec_inv_vec = @(u_vec) pcg(prec_vec, u_vec);
0184 
0185         eigs_opts.issym = false;
0186         eigs_opts.isreal = true;
0187         lambdas = eigs(@(u_vec) prec_vec(hess_vec(u_vec)), ...
0188                        n, dim, 'LM', eigs_opts);
0189         
0190     else
0191         
0192         error('No preconditioner is available in the problem structure.');
0193         
0194     end
0195     
0196     lambdas = sort(lambdas);
0197 
0198 end

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