Home > manopt > solvers > barzilaiborwein > barzilaiborwein.m

# barzilaiborwein

## PURPOSE

Riemannian Barzilai-Borwein solver with non-monotone line-search.

## SYNOPSIS

function [x, cost, info, options] = barzilaiborwein(problem, x, options)

## DESCRIPTION

``` Riemannian Barzilai-Borwein solver with non-monotone line-search.

function [x, cost, info, options] = barzilaiborwein(problem)
function [x, cost, info, options] = barzilaiborwein(problem, x0)
function [x, cost, info, options] = barzilaiborwein(problem, x0, options)
function [x, cost, info, options] = barzilaiborwein(problem, [], options)

Apply the Barzilai-Borwein minimization algorithm to the problem defined
in the problem structure, starting at x0 if it is provided (otherwise, at
a random point on the manifold). To specify options whilst not specifying
an initial guess, give x0 as [] (the empty matrix).

The algorithm uses its own special non-monotone line-search strategy.
Therefore, no lin-search algorithm should be specified in the problem
structure or in the options structure.

In most of the examples bundled with the toolbox (see link below), the
solver can be replaced by the present one if need be.

The outputs x and cost are the last reached point on the manifold and its
cost. This is not necessarily the best point generated since the method
is not monotone. The struct-array info contains information about the
iterations:
iter : the iteration number (0 for the initial guess)
cost : cost value
time : elapsed time in seconds
stepsize : norm of the last tangent vector retracted
linesearch : information logged by the line-search algorithm
And possibly additional information logged by options.statsfun.
For example, type [info.gradnorm] to obtain a vector of the successive

The options structure is used to overwrite the default values. All
options have a default value and are hence optional. To force an option
value, pass an options structure with a field options.optionname, where
optionname is one of the following and the default value is indicated
between parentheses:

The algorithm terminates if the norm of the gradient drops below this.
maxiter (1000)
The algorithm terminates if maxiter iterations have been executed.
maxtime (Inf)
The algorithm terminates if maxtime seconds elapsed.
minstepsize (1e-10)
The algorithm terminates if the linesearch returns a displacement
vector (to be retracted) smaller in norm than this value.
linesearch (@linesearch_hint)
This option should not be changed, as the present solver has its
own dedicated line-search strategy.
strategy ('direct')
The strategy used for the Barzilai-Borwein stepsize
'direct', compute the direct step <s_k,s_k>/<s_k,y_k>
'inverse', compute the inverse step <s_k,y_k>/<y_k,y_k>
'alternate', alternates between direct and inverse step
lambdamax (1e3)
The maximum stepsize allowed by the Barzilai-Borwein method
lambdamin (1e-3)
The minimum stepsize allowed by the Barzilai-Borwein method
lambda0 (1/10)
The initial stepsize of the Barzilai-Borwein method
ls_nmsteps (10)
The non-monotone line-search checks a sufficient decrease with respect
to the previous ls_nmsteps objective function values.
statsfun (none)
Function handle to a function that will be called after each
iteration to provide the opportunity to log additional statistics.
They will be returned in the info struct. See the generic Manopt
documentation about solvers for further information.
stopfun (none)
Function handle to a function that will be called at each iteration
to provide the opportunity to specify additional stopping criteria.
See the generic Manopt documentation about solvers for further
information.
verbosity (3)
Integer number used to tune the amount of output the algorithm
generates during execution (mostly as text in the command window).
The higher, the more output. 0 means silent.
storedepth (2)
Maximum number of different points x of the manifold for which a
store structure will be kept in memory in the storedb. If the
caching features of Manopt are not used, this is irrelevant. For
this algorithm, a store depth of 2 should always be sufficient.

The implementation of the Barzilai-Borwein method is based on the paper:

B. Iannazzo, M. Porcelli, "The Riemannian Barzilai-Borwein method with
nonmonotone line-search and the matrix geometric mean computation",
IMA Journal of Numerical Analysis, to appear, https://doi.org/10.1093/imanum/drx015.

## CROSS-REFERENCE INFORMATION

This function calls:
• StoreDB
• applyStatsfun Apply the statsfun function to a stats structure (for solvers).
• canGetApproxGradient Checks whether an approximate gradient can be computed for this problem.
• canGetCost Checks whether the cost function can be computed for a problem structure.
• canGetGradient Checks whether the gradient can be computed for a problem structure.
• canGetLinesearch Checks whether the problem structure can give a line-search a hint.
• getCostGrad Computes the cost function and the gradient at x in one call if possible.
• getGlobalDefaults Returns a structure with default option values for Manopt.
• mergeOptions Merges two options structures with one having precedence over the other.
• stoppingcriterion Checks for standard stopping criteria, as a helper to solvers.
• approxgradientFD Gradient approx. fnctn handle based on finite differences of the cost.
• linesearch_hint Armijo line-search based on the line-search hint in the problem structure.
This function is called by:

## SOURCE CODE

```0001 function [x, cost, info, options] = barzilaiborwein(problem, x, options)
0002 % Riemannian Barzilai-Borwein solver with non-monotone line-search.
0003 %
0004 % function [x, cost, info, options] = barzilaiborwein(problem)
0005 % function [x, cost, info, options] = barzilaiborwein(problem, x0)
0006 % function [x, cost, info, options] = barzilaiborwein(problem, x0, options)
0007 % function [x, cost, info, options] = barzilaiborwein(problem, [], options)
0008 %
0009 % Apply the Barzilai-Borwein minimization algorithm to the problem defined
0010 % in the problem structure, starting at x0 if it is provided (otherwise, at
0011 % a random point on the manifold). To specify options whilst not specifying
0012 % an initial guess, give x0 as [] (the empty matrix).
0013 %
0014 % The algorithm uses its own special non-monotone line-search strategy.
0015 % Therefore, no lin-search algorithm should be specified in the problem
0016 % structure or in the options structure.
0017 %
0018 % In most of the examples bundled with the toolbox (see link below), the
0019 % solver can be replaced by the present one if need be.
0020 %
0021 % The outputs x and cost are the last reached point on the manifold and its
0022 % cost. This is not necessarily the best point generated since the method
0023 % is not monotone. The struct-array info contains information about the
0024 % iterations:
0025 %   iter : the iteration number (0 for the initial guess)
0026 %   cost : cost value
0027 %   time : elapsed time in seconds
0029 %   stepsize : norm of the last tangent vector retracted
0030 %   linesearch : information logged by the line-search algorithm
0031 %   And possibly additional information logged by options.statsfun.
0032 % For example, type [info.gradnorm] to obtain a vector of the successive
0034 %
0035 % The options structure is used to overwrite the default values. All
0036 % options have a default value and are hence optional. To force an option
0037 % value, pass an options structure with a field options.optionname, where
0038 % optionname is one of the following and the default value is indicated
0039 % between parentheses:
0040 %
0042 %       The algorithm terminates if the norm of the gradient drops below this.
0043 %   maxiter (1000)
0044 %       The algorithm terminates if maxiter iterations have been executed.
0045 %   maxtime (Inf)
0046 %       The algorithm terminates if maxtime seconds elapsed.
0047 %   minstepsize (1e-10)
0048 %       The algorithm terminates if the linesearch returns a displacement
0049 %       vector (to be retracted) smaller in norm than this value.
0050 %   linesearch (@linesearch_hint)
0051 %       This option should not be changed, as the present solver has its
0052 %       own dedicated line-search strategy.
0053 %   strategy ('direct')
0054 %       The strategy used for the Barzilai-Borwein stepsize
0055 %       'direct', compute the direct step <s_k,s_k>/<s_k,y_k>
0056 %       'inverse', compute the inverse step <s_k,y_k>/<y_k,y_k>
0057 %       'alternate', alternates between direct and inverse step
0058 %   lambdamax (1e3)
0059 %       The maximum stepsize allowed by the Barzilai-Borwein method
0060 %   lambdamin (1e-3)
0061 %       The minimum stepsize allowed by the Barzilai-Borwein method
0062 %   lambda0 (1/10)
0063 %       The initial stepsize of the Barzilai-Borwein method
0064 %   ls_nmsteps (10)
0065 %       The non-monotone line-search checks a sufficient decrease with respect
0066 %       to the previous ls_nmsteps objective function values.
0067 %   statsfun (none)
0068 %       Function handle to a function that will be called after each
0069 %       iteration to provide the opportunity to log additional statistics.
0070 %       They will be returned in the info struct. See the generic Manopt
0071 %       documentation about solvers for further information.
0072 %   stopfun (none)
0073 %       Function handle to a function that will be called at each iteration
0074 %       to provide the opportunity to specify additional stopping criteria.
0075 %       See the generic Manopt documentation about solvers for further
0076 %       information.
0077 %   verbosity (3)
0078 %       Integer number used to tune the amount of output the algorithm
0079 %       generates during execution (mostly as text in the command window).
0080 %       The higher, the more output. 0 means silent.
0081 %   storedepth (2)
0082 %       Maximum number of different points x of the manifold for which a
0083 %       store structure will be kept in memory in the storedb. If the
0084 %       caching features of Manopt are not used, this is irrelevant. For
0085 %       this algorithm, a store depth of 2 should always be sufficient.
0086 %
0087 %
0088 % The implementation of the Barzilai-Borwein method is based on the paper:
0089 %
0090 % B. Iannazzo, M. Porcelli, "The Riemannian Barzilai-Borwein method with
0091 % nonmonotone line-search and the matrix geometric mean computation",
0092 % IMA Journal of Numerical Analysis, to appear, https://doi.org/10.1093/imanum/drx015.
0093 %
0095
0096 % This file is part of Manopt: www.manopt.org.
0097 % Original author: Margherita Porcelli, May 31, 2017
0098 % Contributors: Nicolas Boumal, Bruno Iannazzo
0099 % Change log:
0100 %
0101 %   Aug. 2, 2018 (NB):
0102 %       Now using storedb.remove() to keep the cache lean.
0103
0104
0105     % Verify that the problem description is sufficient for the solver.
0106     if ~canGetCost(problem)
0107         warning('manopt:getCost', ...
0108                 'No cost provided. The algorithm will likely abort.');
0109     end
0111         % Note: we do not give a warning if an approximate gradient is
0112         % explicitly given in the problem description, as in that case the
0113         % user seems to be aware of the issue.
0116                 'It may be necessary to increase options.tolgradnorm.\n' ...
0117                 'To disable this warning: warning(''off'', ''manopt:getGradient:approx'')']);
0119     end
0120
0121     % Ensure options exists as a structure
0122     if ~exist('options', 'var') || isempty(options)
0123         options = struct();
0124     end
0125
0126     % Set local defaults here
0127     localdefaults.minstepsize = 1e-10;
0128     localdefaults.maxiter = 1000;
0130
0131     % Upper and lower bound for the Barzilai-Borwein stepsize
0132     localdefaults.lambdamax = 1e3;
0133     localdefaults.lambdamin = 1e-3;
0134     % Initial Barzilai-Borwein stepsize
0135     localdefaults.lambda0 = 1/10;
0136
0137     % Barzilai-Borwein strategy (direct, inverse or alternate)
0138     localdefaults.strategy = 'direct';
0139
0140     % Line-search parameters
0141     % 1) Make sure the user didn't try to define a line search
0142     if canGetLinesearch(problem) || isfield(options, 'linesearch')
0143         error('manopt:BB:ls', ...
0144               ['The problem structure may not specify a line-search ' ...
0145                'hint for the BB solver,\nand the options structure ' ...
0146                'may not specify a line-search algorithm for BB.']);
0147     end
0148     % 2) Define the line-search parameters
0149     problem.linesearch = @(x, d, storedb, key) 1;
0150     options.linesearch = @linesearch_hint;
0151     % The Armijo sufficient decrease parameter
0152     localdefaults.ls_suff_decr = 1e-4;
0153     % The previous steps checked in the non-monotone line-search strategy
0154     localdefaults.ls_nmsteps = 10;
0155
0156
0157     % Merge global and local defaults, then merge w/ user options, if any.
0158     localdefaults = mergeOptions(getGlobalDefaults(), localdefaults);
0159     options = mergeOptions(localdefaults, options);
0160
0161
0162     % Shorthands for some parameters
0163     strategy = options.strategy;
0164     lambdamax = options.lambdamax;
0165     lambdamin = options.lambdamin;
0166     lambda0 = options.lambda0;
0167
0168     timetic = tic();
0169
0170
0171     % If no initial point x is given by the user, generate one at random.
0172     if ~exist('x', 'var') || isempty(x)
0173         x = problem.M.rand();
0174     end
0175
0176     % Create a store database and get a key for the current x
0177     storedb = StoreDB(options.storedepth);
0178     key = storedb.getNewKey();
0179
0180     % Compute objective-related quantities for x
0183
0184     % Some variables below need to store information about iterations. We
0185     % preallocate for a reasonable amount of intended iterations to avoid
0186     % memory re-allocations.
0187     mem_init_size = min(10000, options.maxiter+1);
0188
0189     % Store the cost value
0190     f = zeros(mem_init_size, 1);
0191     f(1) = cost;
0192
0193     % Iteration counter (at any point, iter is the number of fully executed
0194     % iterations so far)
0195     iter = 0;
0196
0197     % Save stats in a struct array info, and preallocate.
0198     stats = savestats();
0199     info(1) = stats;
0200     info(mem_init_size).iter = [];
0201
0202     if options.verbosity >= 2
0203         fprintf(' iter\t                cost val\t     grad. norm\n');
0204     end
0205
0206     % Set the initial Barzilai-Borwein stepsize
0207     lambda = lambda0;
0208
0209     % Start iterating until stopping criterion triggers
0210     while true
0211
0212         % Display iteration information
0213         if options.verbosity >= 2
0215         end
0216
0217         % Start timing this iteration
0218         timetic = tic();
0219
0220         % Run standard stopping criterion checks
0221         [stop, reason] = stoppingcriterion(problem, x, options, ...
0222                                                              info, iter+1);
0223
0224         % If none triggered, run specific stopping criterion check
0225         if ~stop && stats.stepsize < options.minstepsize
0226             stop = true;
0227             reason = sprintf(['Last stepsize smaller than minimum '  ...
0228                               'allowed; options.minstepsize = %g.'], ...
0229                               options.minstepsize);
0230         end
0231
0232         if stop
0233             if options.verbosity >= 1
0234                 fprintf([reason '\n']);
0235             end
0236             break;
0237         end
0238
0239         % Pick the descent direction as minus the gradient (scaled)
0240         desc_dir = problem.M.lincomb(x, -lambda, grad);
0241
0242         % Execute the nonmonotone line search
0243         k = iter + 1;
0244         start = max(1, k - options.ls_nmsteps + 1);
0245
0246         [stepsize, newx, newkey, lsstats] = ...
0247             options.linesearch(problem, x, desc_dir, max(f(start:k)), ...
0248                             -lambda * gradnorm^2, options, storedb, key);
0249
0250         % Updates the value of lambda
0251         lambda = lambda * lsstats.alpha;
0252
0253         % Compute the new cost-related quantities for newx
0256
0257         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0258         % BARZILAI-BORWEIN STRATEGY
0259
0260         % Store the cost value
0261         f(iter+2) = newcost;
0262
0263         % Transport the old gradient to newx
0265
0266         % Compute the difference between grandients
0268
0269         % Compute the transported step
0270         Stransp =  problem.M.lincomb(x, -lambda, grad_transp);
0271
0272         % Compute the new Barzilai-Borwein step following the strategy
0273         % direct strategy
0274         if strcmp(strategy, 'direct')
0275           num = problem.M.norm(newx, Stransp)^2;
0276           den = problem.M.inner(newx, Stransp, Y);
0277           if den > 0
0278             lambda = min( lambdamax, max(lambdamin, num/den) );
0279           else
0280             lambda = lambdamax;
0281           end
0282         end
0283
0284         % inverse strategy
0285         if strcmp(strategy, 'inverse')
0286           num = problem.M.inner(newx, Stransp, Y);
0287           den = problem.M.norm(newx, Y)^2;
0288
0289           if num > 0
0290             lambda = min( lambdamax, max(lambdamin, num/den) );
0291           else
0292             lambda = lambdamax;
0293           end
0294         end
0295
0296         % alternate strategy
0297         if strcmp(strategy, 'alternate')
0298           num = problem.M.norm(newx, Stransp)^2;
0299           den = problem.M.inner(newx, Stransp, Y);
0300           den2 = problem.M.norm(newx, Y)^2;
0301           if (den > 0)
0302             if mod(iter,2)==0
0303                 lambda = min( lambdamax, max(lambdamin, num/den) );
0304         else
0305                 lambda = min( lambdamax, max(lambdamin, den/den2) );
0306             end
0307           else
0308             lambda = lambdamax;
0309           end
0310         end
0311
0312         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0313
0314         % Update iterate info
0315         storedb.removefirstifdifferent(key, newkey);
0316         x = newx;
0317         key = newkey;
0318         cost = newcost;
0321
0322         % iter is the number of iterations we have accomplished.
0323         iter = iter + 1;
0324
0325         % Make sure we don't use too much memory for the store database.
0326         storedb.purge();
0327
0328         % Log statistics for freshly executed iteration
0329         stats = savestats();
0330         info(iter+1) = stats;
0331
0332     end
0333
0334     info = info(1:iter+1);
0335
0336     if options.verbosity >= 1
0337         fprintf('Total time is %f [s] (excludes statsfun)\n', ...
0338                 info(end).time);
0339     end
0340
0341
0342
0343     % Routine in charge of collecting the current iteration stats
0344     function stats = savestats()
0345         stats.iter = iter;
0346         stats.cost = cost;
0348         if iter == 0
0349             stats.stepsize = NaN;
0350             stats.time = toc(timetic);
0351             stats.linesearch = [];
0352         else
0353             stats.stepsize = stepsize;
0354             stats.time = info(iter).time + toc(timetic);
0355             stats.linesearch = lsstats;
0356         end
0357         stats = applyStatsfun(problem, x, storedb, key, options, stats);
0358     end
0359
0360 end```

Generated on Mon 10-Sep-2018 11:48:06 by m2html © 2005