Solves trust-region subproblem with TRSgep in a subspace of tangent space.
function trsoutput = trs_gep(problem, trsinput, options, storedb, key)
minimize <eta, grad> + .5*<eta, Hess(eta)>
subject to <eta, eta> <= Delta^2
This is meant to be used by the trustregion solver.
To use this method, specify trs_gep as an option, and your chosen
subspace dimension in the problem structure, as follows:
options.subproblemsolver = @trs_gep;
options.gepsubspacedim = n; % Integer in the range 1:problem.M.dim().
% If omitted, default is problem.M.dim().
x = trustregions(problem, [], options);
Note: trs_gep does not use preconditioning.
In principle, trs_gep solves the trust-region subproblem exactly in a
subspace of the tangent space with dimension options.gepsubspacedim.
If that dimension is equal to the manifold dimension, then the solver
is meant to find a global optimum of the TRS, up to numerical issues.
This function achieves that goal as follows: it creates an orthonormal
basis for (a subspace of) the tangent space using tangentorthobasis,
it expresses the Hessian and gradient of the cost function at the
current point x (restricted to the subspace) in the chosen basis, and it
passes those objects to TRSgep.
The basis is constructed by tangentorthobasis with randomly tangent
vectors (linearly independent with probability 1) then orthonormalized.
Therefore, if gepsubspacedim is less than the manifold dimension, the
minimization is executed over a random subspace. In that scenario, if the
gradient is nonzero, the gradient is included in the basis to be
orthonormalized. This ensures that the point returned by this solver is
always as good as the Cauchy point.
Constructing the basis itself can be time consuming in high dimensions,
and aiming for an exact solve of the TRS as well. This subproblem solver
is meant mostly for research, not for efficiency.
Inputs:
problem: Manopt optimization problem structure
trsinput: structure with the following fields:
x: point on the manifold problem.M
fgradx: gradient of the cost function of the problem at x
Delta: trust-region radius
options: structure containing options for the subproblem solver
storedb, key: manopt's caching system for the point x
Options specific to this subproblem solver (default value):
gepsubspacedim (problem.M.dim())
A value between 1 and problem.M.dim() inclusive that specifies the
dimension of the subpsace over which we wish to solve the
trust-region subproblem.
Output: the structure trsoutput contains the following fields:
eta: approximate solution to the trust-region subproblem at x
Heta: Hess f(x)[eta] -- this is necessary in the outer loop, and it
is often naturally available to the subproblem solver at the
end of execution, so that it may be cheaper to return it here.
limitedbyTR: true if a boundary solution is returned
printstr: logged information to be printed by trustregions.
stats: structure with the following statistics:
hessvecevals: number of Hessian-vector calls issued
trs_gep can also be called in the following way (by trustregions) to
obtain part of the header to print and an initial stats structure:
function trsoutput = trs_gep([], [], options)
In this case, trsoutput contains the following fields:
printheader: subproblem header to be printed before the first loop of
trustregions
initstats: struct with initial values for stored stats in subsequent
calls to trs_gep. Used in the first call to savestats
in trustregions to initialize the info struct properly.
See also: TRSgep trs_tCG trustregions0001 function trsoutput = trs_gep(problem, trsinput, options, ~, ~) 0002 % Solves trust-region subproblem with TRSgep in a subspace of tangent space. 0003 % 0004 % function trsoutput = trs_gep(problem, trsinput, options, storedb, key) 0005 % 0006 % minimize <eta, grad> + .5*<eta, Hess(eta)> 0007 % subject to <eta, eta> <= Delta^2 0008 % 0009 % This is meant to be used by the trustregion solver. 0010 % To use this method, specify trs_gep as an option, and your chosen 0011 % subspace dimension in the problem structure, as follows: 0012 % 0013 % options.subproblemsolver = @trs_gep; 0014 % options.gepsubspacedim = n; % Integer in the range 1:problem.M.dim(). 0015 % % If omitted, default is problem.M.dim(). 0016 % x = trustregions(problem, [], options); 0017 % 0018 % Note: trs_gep does not use preconditioning. 0019 % 0020 % In principle, trs_gep solves the trust-region subproblem exactly in a 0021 % subspace of the tangent space with dimension options.gepsubspacedim. 0022 % If that dimension is equal to the manifold dimension, then the solver 0023 % is meant to find a global optimum of the TRS, up to numerical issues. 0024 % 0025 % This function achieves that goal as follows: it creates an orthonormal 0026 % basis for (a subspace of) the tangent space using tangentorthobasis, 0027 % it expresses the Hessian and gradient of the cost function at the 0028 % current point x (restricted to the subspace) in the chosen basis, and it 0029 % passes those objects to TRSgep. 0030 % 0031 % The basis is constructed by tangentorthobasis with randomly tangent 0032 % vectors (linearly independent with probability 1) then orthonormalized. 0033 % Therefore, if gepsubspacedim is less than the manifold dimension, the 0034 % minimization is executed over a random subspace. In that scenario, if the 0035 % gradient is nonzero, the gradient is included in the basis to be 0036 % orthonormalized. This ensures that the point returned by this solver is 0037 % always as good as the Cauchy point. 0038 % 0039 % Constructing the basis itself can be time consuming in high dimensions, 0040 % and aiming for an exact solve of the TRS as well. This subproblem solver 0041 % is meant mostly for research, not for efficiency. 0042 % 0043 % 0044 % Inputs: 0045 % problem: Manopt optimization problem structure 0046 % trsinput: structure with the following fields: 0047 % x: point on the manifold problem.M 0048 % fgradx: gradient of the cost function of the problem at x 0049 % Delta: trust-region radius 0050 % options: structure containing options for the subproblem solver 0051 % storedb, key: manopt's caching system for the point x 0052 % 0053 % Options specific to this subproblem solver (default value): 0054 % gepsubspacedim (problem.M.dim()) 0055 % A value between 1 and problem.M.dim() inclusive that specifies the 0056 % dimension of the subpsace over which we wish to solve the 0057 % trust-region subproblem. 0058 % 0059 % Output: the structure trsoutput contains the following fields: 0060 % eta: approximate solution to the trust-region subproblem at x 0061 % Heta: Hess f(x)[eta] -- this is necessary in the outer loop, and it 0062 % is often naturally available to the subproblem solver at the 0063 % end of execution, so that it may be cheaper to return it here. 0064 % limitedbyTR: true if a boundary solution is returned 0065 % printstr: logged information to be printed by trustregions. 0066 % stats: structure with the following statistics: 0067 % hessvecevals: number of Hessian-vector calls issued 0068 % 0069 % 0070 % trs_gep can also be called in the following way (by trustregions) to 0071 % obtain part of the header to print and an initial stats structure: 0072 % 0073 % function trsoutput = trs_gep([], [], options) 0074 % 0075 % In this case, trsoutput contains the following fields: 0076 % printheader: subproblem header to be printed before the first loop of 0077 % trustregions 0078 % initstats: struct with initial values for stored stats in subsequent 0079 % calls to trs_gep. Used in the first call to savestats 0080 % in trustregions to initialize the info struct properly. 0081 % 0082 % See also: TRSgep trs_tCG trustregions 0083 0084 % This file is part of Manopt: www.manopt.org. 0085 % Original author: Victor Liao, June 13, 2022. 0086 % Contributors: Nicolas Boumal 0087 % Change log: 0088 0089 % trustregions only wants header and default values for stats. 0090 if nargin == 3 && isempty(problem) && isempty(trsinput) 0091 trsoutput.printheader = sprintf('%9s %s', 'hessvec', 'stopreason'); 0092 trsoutput.initstats = struct('hessvecevals', 0); 0093 return; 0094 end 0095 0096 x = trsinput.x; 0097 grad = trsinput.fgradx; 0098 Delta = trsinput.Delta; 0099 0100 M = problem.M; 0101 0102 % Set local defaults here 0103 localdefaults.gepsubspacedim = M.dim(); 0104 0105 % Merge global and local defaults, then merge w/ user options, if any. 0106 localdefaults = mergeOptions(getGlobalDefaults(), localdefaults); 0107 if ~exist('options', 'var') || isempty(options) 0108 options = struct(); 0109 end 0110 options = mergeOptions(localdefaults, options); 0111 0112 % dimension of subspace over which to solve the TRS. 0113 n = options.gepsubspacedim; 0114 0115 assert(n >= 1 && n <= M.dim() && n == round(n), ... 0116 'options.gepsubspacedim must be an integer between 1 and M.dim().'); 0117 0118 % If gradient is nonzero, then even if n < M.dim() we 0119 % guarantee to do at least as well as the cauchy point 0120 % by including the gradient in the basis of the subspace. 0121 % The column vector grad_vec contains the coordinates of grad in the 0122 % basis B. 0123 grad_norm = M.norm(x, grad); 0124 if grad_norm ~= 0 0125 % Append n-1 random tangent vectors to the gradient, 0126 % then orthonormalize with Gram-Schmidt. 0127 B = tangentorthobasis(M, x, n, {grad}); 0128 grad_vec = zeros(n, 1); 0129 grad_vec(1) = grad_norm; 0130 else 0131 B = tangentorthobasis(M, x, n); 0132 grad_vec = zeros(n, 1); 0133 end 0134 0135 % Express the Hessian of the cost function f at x in the basis B. 0136 % If B is a basis for a subspace of T_x M rather than for the whole 0137 % tangent space, then H represents the Hessian restriced to that 0138 % subspace. 0139 H = hessianmatrix(problem, x, B); 0140 0141 % This is where the actual work happens. 0142 [eta_vec, limitedbyTR] = TRSgep(H, grad_vec, Delta); 0143 0144 % Construct the tangent vector eta using its coordinates eta_vec in the 0145 % basis B. 0146 eta = lincomb(M, x, B, eta_vec); 0147 0148 % We want to return Heta, defined by: 0149 % Heta = getHessian(problem, x, eta, storedb, key). 0150 % This however requires one Hessien-vector call. 0151 % We can avoid issuing that cal with the two lines below. 0152 % This is likely to be faster, but may be less accurate numerically. 0153 % Which is better may depend on the application. 0154 Heta_vec = H*eta_vec; 0155 Heta = lincomb(M, x, B, Heta_vec); 0156 0157 if limitedbyTR 0158 stopreason_str = 'Exact trs_gep boundary sol'; 0159 else 0160 stopreason_str = 'Exact trs_gep interior sol'; 0161 end 0162 0163 printstr = sprintf('%9d %s', n, stopreason_str); 0164 stats = struct('hessvecevals', n); 0165 0166 trsoutput.eta = eta; 0167 trsoutput.Heta = Heta; 0168 trsoutput.limitedbyTR = limitedbyTR; 0169 trsoutput.printstr = printstr; 0170 trsoutput.stats = stats; 0171 end