Home > manopt > solvers > trustregions > TRSgep.m

# TRSgep

## PURPOSE

Solves trust-region subproblem by a generalized eigenvalue problem.

## SYNOPSIS

function [x, limitedbyTR, eigstrouble, accurate] = TRSgep(A, a, Del)

## DESCRIPTION

``` Solves trust-region subproblem by a generalized eigenvalue problem.

function [x, limitedbyTR, eigstrouble] = TRSgep(A, a, Del)
function [x, limitedbyTR, eigstrouble, accurate] = TRSgep(A, a, Del)

This function returns a solution x to the following optimization problem:

minimize .5*(x.'*A*x) + a.'*x
subject to x.'*x <= Del^2

The boolean 'limitedbyTR' is true if the solution would have been
different absent the norm constraint. In that case, the norm of x is Del.

The boolean 'eigstrouble' is false if the call to eigs went as planned
(no NaN values returned). If is true otherwise. Even when true, the
output x may still be good, but it is not certain.

Inputs:
A: nxn symmetric
a: nx1 vector, both real

If called with three outputs, then computationally expensive checks are
run to verify the accuracy of the output. If the output appears to be
globally optimal (as expected) within some demanding numerical
tolerances, then 'accurate' is true; otherwise it is false.

The iterative solver pcg sometimes issues a warning.
It can be disabled with:
warning('off', 'MATLAB:pcg:tooSmallTolerance');
It is adviseable to re-enable it (with 'on' instead of 'off') when done.

Code adapted from Yuji Nakatsukasa's code for the
paper by Satoru Adachi, Satoru Iwata, Yuji Nakatsukasa, and Akiko Takeda

Original code: https://people.maths.ox.ac.uk/nakatsukasa/codes/TRSgep.m
Reference paper: https://epubs.siam.org/doi/abs/10.1137/16M1058200

The authors kindly allowed us to include their code in Manopt under the

## CROSS-REFERENCE INFORMATION

This function calls:
• norm NORM Norm of a TT/MPS tensor.
• norm NORM Norm of a TT/MPS block-mu tensor.
This function is called by:
• trs_gep Solves trust-region subproblem with TRSgep in a subspace of tangent space.

## SOURCE CODE

```0001 function [x, limitedbyTR, eigstrouble, accurate] = TRSgep(A, a, Del)
0002 % Solves trust-region subproblem by a generalized eigenvalue problem.
0003 %
0004 % function [x, limitedbyTR, eigstrouble] = TRSgep(A, a, Del)
0005 % function [x, limitedbyTR, eigstrouble, accurate] = TRSgep(A, a, Del)
0006 %
0007 % This function returns a solution x to the following optimization problem:
0008 %
0009 %     minimize .5*(x.'*A*x) + a.'*x
0010 %     subject to x.'*x <= Del^2
0011 %
0012 % The boolean 'limitedbyTR' is true if the solution would have been
0013 % different absent the norm constraint. In that case, the norm of x is Del.
0014 %
0015 % The boolean 'eigstrouble' is false if the call to eigs went as planned
0016 % (no NaN values returned). If is true otherwise. Even when true, the
0017 % output x may still be good, but it is not certain.
0018 %
0019 % Inputs:
0020 %   A: nxn symmetric
0021 %   a: nx1 vector, both real
0022 %   Del: trust-region radius (positive real)
0023 %
0024 % If called with three outputs, then computationally expensive checks are
0025 % run to verify the accuracy of the output. If the output appears to be
0026 % globally optimal (as expected) within some demanding numerical
0027 % tolerances, then 'accurate' is true; otherwise it is false.
0028 %
0029 % The iterative solver pcg sometimes issues a warning.
0030 % It can be disabled with:
0031 %   warning('off', 'MATLAB:pcg:tooSmallTolerance');
0032 % It is adviseable to re-enable it (with 'on' instead of 'off') when done.
0033 %
0034 % Code adapted from Yuji Nakatsukasa's code for the
0035 % paper by Satoru Adachi, Satoru Iwata, Yuji Nakatsukasa, and Akiko Takeda
0036 %
0037 % Original code: https://people.maths.ox.ac.uk/nakatsukasa/codes/TRSgep.m
0038 % Reference paper: https://epubs.siam.org/doi/abs/10.1137/16M1058200
0039 %
0040 % The authors kindly allowed us to include their code in Manopt under the
0041 % same license as Manopt.
0042 %
0044
0045 % This file is part of Manopt: www.manopt.org.
0046 % Original author: Yuji Nakatsukasa, 2015.
0047 % Contributors: Revised by Nikitas Rontsis, December 2018
0048 % Change log:
0049 %   VL June 29, 2022:
0050 %       Modified original code to return limitedbyTR boolean and change
0051 %       ellipsoid norm constraint to unweighted norm.
0052 %   NB Aug. 19, 2022:
0053 %       Comments + cosmetic changes.
0054 %       Corrected determination of limitedbyTR.
0055 %       Added support for input a = 0.
0056 %       Clarified the logic around picking the Newton step or not.
0057 %   NB Aug. 26, 2022:
0058 %       Added optional accuracy checks.
0059
0060     n = size(A, 1);
0061
0062     eigstrouble = false;
0063
0064     % We set this flag to true iff the solution x we eventually return is
0065     % limited by the trust-region boundary.
0066     limitedbyTR = true;
0067
0068     % Tolerance for hard-case.
0069     % If this triggers, then the solver works harder to check itself.
0070     tolhardcase = 1e-4;
0071
0072     % If a is exactly zero, pcg (called below) abandons on the first
0073     % iteration. Instead, we give it a small input and re-check at the end.
0074     a_is_zero = (norm(a) == 0);
0075     if a_is_zero
0076         a = eps*randn(n, 1);
0077     end
0078
0079     % Compute the Newton step p1 up to some accuracy.
0080     % pcg sometimes issues a warning: see help above.
0081     [p1, ~, relres, ~] = pcg(A, -a, 1e-12, 500);
0082
0083     % If the Newton step is computed accurately and it is in the trust
0084     % region, then it may very well be the solution to the TRS.
0085     % We make a note of it, and will re-check at the end.
0086     newton_step_may_be_solution = (relres < 1e-5 && (p1'*p1 <= Del^2));
0087
0088     % This is the core of the code.
0089     MM1 = [sparse(n, n) speye(n) ; speye(n) sparse(n, n)];
0090     [V, lam1] = eigs(@(x) MM0timesx(A, a, Del, x), 2*n, -MM1, 1, 'lr');
0091
0092     % Sometimes the output is complex.
0093     if norm(real(V)) < 1e-3
0094         V = imag(V);
0095     else
0096         V = real(V);
0097     end
0098     lam1 = real(lam1);
0099
0100     % This is parallel to the solution:
0101     x = V(1:n);
0102     normx = norm(x);
0103
0104     % In the easy case, this naive normalization improves accuracy.
0105     x = x/normx*Del;
0106     % Take the correct sign.
0107     if x'*a > 0
0108         x = -x;
0109     end
0110
0111     % If we suspect a (numerically) hard case, work harder.
0112     if normx < tolhardcase
0113
0114         % Sometimes, eigs fails: it issues a warning and sets lam1 to NaN.
0115         % The vector V may still be defined and partially useable though.
0116         % Here, we estimate what lam1 ought to be (see also the code below
0117         % that checks accuracy) and use to proceed.
0118         if isnan(lam1)
0119             lam1 = -(x'*(A*x + a))/(x'*x);
0120             eigstrouble = true;
0121         end
0122
0123         x1 = V(n+1:end);
0124         alpha1 = lam1;
0125         Pvect = x1;
0126         % First try only k = 1, that is almost always enough.
0127         % pcg sometimes issues a warning: see help above.
0128         Afun = @(x) pcgforAtilde(A, lam1, Pvect, alpha1, x);
0129         [x2, ~] = pcg(Afun, -a, 1e-12, 500);
0130         % If large residual, repeat
0131         if norm((A+lam1)*x2 + a) > tolhardcase*norm(a)
0132             for ii = [3, 6, 9]
0133                 [Pvect, ~] = eigs(A, speye(n), ii, 'sa');
0134                 % pcg sometimes issues a warning: see help above.
0135                 Afun = @(x) pcgforAtilde(A, lam1, Pvect, alpha1, x);
0136                 [x2, ~] = pcg(Afun, -a, 1e-8, 500);
0137                 if norm((A+lam1)*x2 + a) < tolhardcase*norm(a)
0138                     break;
0139                 end
0140             end
0141         end
0142
0143         aa = x1'*x1;
0144         bb = 2*(x2'*x1);
0145         cc = x2'*x2 - Del^2;
0146         % Move to the boundary: set alp such that norm(x2+alp*x1) = Delta.
0147         % alp = (-bb + sqrt(bb^2 - 4*aa*cc))/(2*aa);
0148         alp = max(real(roots([aa, bb, cc])));
0149         x = x2 + alp*x1;
0150     end
0151
0152     % If we suspected that the Newton step might be the solution to the
0153     % TRS, we compare it to the boundary solution we just computed and pick
0154     % the best one.
0155     if newton_step_may_be_solution
0156         if (p1'*A*p1)/2 + a'*p1 < (x'*A*x)/2 + a'*x
0157             x = p1;
0158             limitedbyTR = false;
0159         end
0160     end
0161
0162     % If the input a was zero, then earlier in the code we replaced it with
0163     % a tiny random vector. Two things may have happened afterwards:
0164     % If A is positive definite, then the solution x is also a tiny vector.
0165     % In all likelihood, it did not hit the TR: then we know to replace x
0166     % with zero.
0167     % Otherwise, at least one eigenvalue of A is <= 0, and there exists a
0168     % solution on the boundary of the trust-region: that is what should
0169     % have been computed already, hence we do nothing.
0170     if a_is_zero && ~limitedbyTR
0171         x = zeros(n, 1);
0172     end
0173
0174
0175     % This is for debugging purposes only: it is expensive to run.
0176     % The code checks via a dual certificate that x is a global optimum,
0177     % up to some numerical tolerances. It also checks limitedbyTR.
0178     if nargout >= 3
0179         tol = 1e-13;
0180         mineig = min(eig(A));
0181         if norm(x) ~= 0
0182             % Estimate the dual variable for the norm constraint.
0183             mu = -(x'*(A*x + a))/(x'*x);
0184             % The vector x is optimal iff:
0185             %   norm(x) <= Del,
0186             %   M = A + mu*I is psd and mu >= 0,
0187             %   M*x + b = 0, and
0188             %   mu = 0 whenever we are not limited by TR.
0189             % We also need that limitedbyTR => norm(x) == Del.
0190             reltol = @(c) c + tol*max(1, c); % to check a <~ c with c >= 0.
0191             accurate = (norm(x) <= reltol(Del) && ...
0192                         max(0, -mineig) <= reltol(mu) && ...
0193                         all(abs(A*x+a + mu*x) <= reltol(abs(mu*x))) && ...
0194                         ( limitedbyTR || mu <= reltol(0)) && ...
0195                         (~limitedbyTR || Del <= reltol(norm(x))));
0196         else
0197             % The zero vector x is optimal iff a = 0 and A is psd.
0198             % Moreover, a solution x = 0 is clearly not limited by the TR.
0199             accurate = (norm(a) <= tol && mineig >= -tol && ~limitedbyTR);
0200         end
0201     end
0202
0203 end
0204
0205
0206
0207 function y = MM0timesx(A, g, Delta, x)
0208     % MM0 = [-Id A;
0209     %         A -g*g'/Delta^2];
0210     n = size(A, 1);
0211     x1 = x(1:n);
0212     x2 = x(n+1:end);
0213     y1 = -x1 + A*x2;
0214     y2 = A*x1 - g*(g'*x2)/Delta^2;
0215     y = [y1 ; y2];
0216 end
0217
0218 function y = pcgforAtilde(A, lamA, Pvect, alpha1, x)
0219     m = size(Pvect, 2);
0220     y = A*x + lamA*x;
0221     for ii = 1:m
0222         y = y + (alpha1*(x'*(Pvect(:, ii))))*(Pvect(:, ii));
0223     end
0224 end```

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