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     % 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.
0143     if canGetGradient(problem)
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