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