Home > manopt > solvers > trustregions > trs_gep.m

# trs_gep

## PURPOSE

Solves trust-region subproblem with TRSgep in a subspace of tangent space.

## SYNOPSIS

function trsoutput = trs_gep(problem, trsinput, options, ~, ~)

## DESCRIPTION

``` 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
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.

## CROSS-REFERENCE INFORMATION

This function calls:
• getGlobalDefaults Returns a structure with default option values for Manopt.
• mergeOptions Merges two options structures with one having precedence over the other.
• round ROUND Approximate TTeMPS tensor within a prescribed tolerance.
• round ROUND Approximate TTeMPS tensor within a prescribed tolerance.
• round ROUND Approximate TTeMPS operator within a prescribed tolerance.
• TRSgep Solves trust-region subproblem by a generalized eigenvalue problem.
• hessianmatrix Computes a matrix which represents the Hessian in some tangent basis.
• lincomb Computes a linear combination of tangent vectors in the Manopt framework.
• tangentorthobasis Returns an orthonormal basis of tangent vectors in the Manopt framework.
This function is called by:

## SOURCE CODE

```0001 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;
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.
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);
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```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005