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
Del: trust-region radius (positive 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
same license as Manopt.
See also: trs_gep trs_tCG_cached trs_tCG trustregions0001 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 % 0043 % See also: trs_gep trs_tCG_cached trs_tCG trustregions 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