A comparison of eigenvalue-based algorithms and the generalized Lanczos trust-region algorithm for Solving the trust-region subproblem
AA COMPARISON OF EIGENVALUE-BASED ALGORITHMS ANDTHE GENERALIZED LANCZOS TRUST-REGION ALGORITHMFOR SOLVING THE TRUST-REGION SUBPROBLEM ∗ ZHONGXIAO JIA y AND
FA WANG z Abstract.
Solving the trust-region subproblem (TRS) plays a key role in numerical optimizationand many other applications. Based on a fundamental result that the solution of TRS of size n ismathematically equivalent to finding the rightmost eigenpair of a certain matrix pair of size 2 n ,eigenvalue-based methods are promising due to their simplicity. For n large, the implicitly restartedArnoldi (IRA) and refined Arnoldi (IRRA) algorithms are well suited for this eigenproblem. For areasonable comparison of overall efficiency of the algorithms for solving TRS directly and eigenvalue-based algorithms, a vital premise is that the two kinds of algorithms must compute the approximatesolutions of TRS with (almost) the same accuracy, but such premise has been ignored in the literature.To this end, we establish close relationships between the two kinds of residual norms, so that, givena stopping tolerance for IRA and IRRA, we are able to determine a reliable one that GLTR shoulduse so as to ensure that GLTR and IRA, IRRA deliver the converged approximate solutions withsimilar accuracy. We also make a convergence analysis on the residual norms by the GeneralizedLanczos Trust-Region (GLTR) algorithm for solving TRS directly, the Arnoldi method and therefined Arnoldi method for the equivalent eigenproblem. A number of numerical experiments arereported to illustrate that IRA and IRRA are competitive with GLTR and IRRA outperforms IRA. Key words. trust-region subproblem, eigenvalue problem, implicitly restarted Arnoldi algorith-m, implicitly restarted refined Arnoldi algorithm, GLTR algorithm, Krylov subspace, residual norm,stopping criterion
AMS subject classi(cid:12)cations.
1. Introduction.
The trust-region subproblem (TRS) is(1.1) min ∥ s ∥ B ≤ ∆ q ( s ) with q ( s ) = g T s + 12 s T As, where A ∈ R n × n is symmetric, B ∈ R n × n is symmetric positive definite, the nonzerovector g ∈ R n , ∆ > ∥ · ∥ B is the B -normdefined by ∥ s ∥ B = √ s T Bs . The matrix B is not necessarily equal to the identity I ,and it is often constructed to impose a smoothness condition on a solution to (1.1)for the ill-posed problem and to incorporate scaling of variables in optimization [22].For instance, it is argued in [2] that a good choice is B = J − T J − for some invertiblematrix J or the Hermitian polar factor of A [13].Problem (1.1) arises from many applications, e.g., from nonlinear numerical op-timization [2, 21], where q ( s ) is a quadratic model of min f ( s ) at the current ap-proximate solution, A is Hessian and g is the gradient of f at the current approxi-mate solution, and many others, e.g., Tikhonov regularization of ill-posed problems[22, 23], graph partitioning problems [12], the constrained eigenvalue problem [6], andthe Levenberg–Marquardt algorithm for solving nonlinear least squares problems [21].The following theorem [2, 20] provides a theoretical basis for many TRS solutionalgorithms and give necessary and sufficient conditions, called the optimal conditions, (cid:3) This work was supported in part by the National Natural Science Foundation of China (No.11771249) y Corresponding author. Department of Mathematical Sciences, Tsinghua University, 100084 Bei-jing, China. ( [email protected] ) z Department of Mathematical Sciences, Tsinghua University, 100084 Beijing, China.( [email protected] ) 1
Z. JIA AND F. WANG for the solution of TRS (1.1).
Theorem 1.1.
A vector s opt is a solution to (1.1) if and only if there exists theoptimal Lagrangian multiplier λ opt ≥ such that ∥ s opt ∥ B ∆ , (1.2) ( A + λ opt B ) s opt = − g, (1.3) λ opt (∆ − ∥ s opt ∥ B ) = 0 , (1.4) A + λ opt B ≽ , (1.5) where the notation ≽ indicates that a symmetric matrix is semi-positive de(cid:12)nite. TRS algorithms have been extensively studied for a few decades, and a number ofalgorithms have been available. On the basis of [1], the authors in [19] have classifiedthe known algorithms into four categories: accurate methods for small to mediumsized problems [20], accurate methods for large sparse problems [4, 10, 11, 24, 22, 23,27], approximate methods [28, 31] which are designed to approximately solve (1.1),and eigenvalue-based methods [1, 6], which are suitable for both small to mediumsized problems and large sparse problems.Over the past two decades, the GLTR method [8] has become a commonly usedalgorithm for large scale TRS’s. It consists of two phases: In the first phase, itstarts with the Truncated Conjugate Gradient (TCG) algorithm with the zero initialguess [28, 31]. When A is positive definite and ∥ A − g ∥ B ≤ ∆, which corresponds to λ opt = 0, the algorithm ultimately returns a converged approximate solution to s opt = − A − g and stops. If the TCG iterate exceeds the trust-region boundary or a negativecurvature is encountered, which means that the trust-region constraint is activatedand corresponds to a semi-positive definite or indefinite A but λ opt >
0, the GLTRmethod enters the second phase and switches to the Lanczos method that accuratelysolves a sequence of the projected problems restricted to some Krylov subspaces, andit proceeds in such a way until the approximate solution converges to s opt . The GLTRmethod may require many iterations, so that storing the basis vectors may becomeunacceptable. One may store them in secondary memory and calls them when formingGLTR iterates, but this way is time consuming. Besides, numerical difficulties such asloss of orthogonality among the Lanczos basis vectors make the algorithm intricate.To this end, Zhang et al. [33, 34] have recently proposed a complicated restartedGLTR algorithm, which needs further explorations and developments.A potentially great breakthrough has come to the solution of TRS’s when Adachiet al. [1] have proved that TRS (1.1) can be treated by solving the generalizedeigenvalue problem of a certain matrix pair of order 2 n if the solution s opt to the TRS(1.1) is in the easy case and reaches the trust-region boundary. Precisely, they haveshown that a solution to TRS (1.1) can be determined by the rightmost eigenvalueand the associated eigenvector of some 2 n × n matrix pair. Therefore, computingthe solution s opt of TRS (1.1) amounts to computing the rightmost eigenpair of ageneralized eigenvalue problem, for which a number of standard numerical methodsare applicable [7, 29], depending on the size of n . Concretely, for A large, Adachi etal. have used the implicitly restarted Arnoldi (IRA) algorithm proposed by Sorensen[26], i.e., the Matlab function eigs , to solve the generalized eigenvalue problem. IRAis simple to implement, and has the advantage of low storage and ingenious restartingstrategy. Numerical experiments in [1, 34] have indicated that IRA can be competitivewith the GLTR algorithm in terms of CPU time and accuracy and may thus be apromising algorithm for large scale TRS’s. However, the comparison standards used igenvalue-based algorithms for TRS s opt withsimilar accuracy. Therefore, we should find relationships between them, relate themand set them correctly in order to make two kinds of algorithms give almost the samesolution accuracy for TRS (1.1) when making an efficiency comparison.For large scale matrix eigenproblems, it has been known from, e.g., [14, 18, 29],that Arnoldi algorithms, such as IRA, may have convergence problems [14, 18] for com-puting eigenvectors. Precisely, the approximate eigenvectors, i.e., the Ritz vectors,may converge erratically and even may not converge even though the correspondingRitz values converge. Jia [15, 16] has corrected this problem by proposing the refinedRayleigh–Ritz method that uses new approximate eigenvectors, called refined Ritzvectors or refined eigenvector approximations, to approximate the desired eigenvec-tors; see [16, 29] for details. More generally, given a general projection subspace,the refined Ritz vector is more accurate than the Ritz vector [17]. Particularly, forthe refined Arnoldi method [15], Jia [16] has developed an implicitly restarted refinedArnoldi (IRRA) algorithm, which applies the same implicitly restarting scheme as inIRA to the refined Arnoldi method but uses the new refined shifts to construct a bet-ter projection subspace at each restart than IRA does, so that IRRA may outperformIRA substantially.In this paper, we extend the refined Arnoldi method [15] and IRRA [18] for thestandard eigenvalue problems to generalized ones. Then our main goal is to makeextensive experiments on TRS’s with B = I and B ̸ = I and to show whether or notIRA and IRRA are practically competitive with the GLTR algorithm without restart.For this purpose, we need to settle down two important issues. One of them is toselect a fair comparison measure. As is well known, in exact arithmetic, the Lanczosprocess in GLTR exploits elegant three term recurrences to compute an orthonormalor B -orthonormal basis of the underlying Krylov subspace with B = I or B ̸ = I , butthe Arnoldi process in IRA and IRRA uses long recurrences to do similar things. Forsubspaces of small dimension, due to algorithmic features, the number of matrix-vectorproducts (MVs) dominates the efficiency of the GLTR method, the Arnoldi methodand the refined Arnoldi method for subspaces of small dimension. As a result, thetotal number of MVs for GLTR, IRA and IRRA is a reasonable measure to judgetheir overall efficiency when subspace dimension for IRA and IRRA is small. Thereare certainly other measures, such as the CPU time. But the reliability of CPU time ismuch involved, and it strongly depends on programming optimization, programminglanguages and computing environment, etc. In any event, we must stress that MVsis only a rough measure of overall performance, and we ignore complicated effects ofseveral other factors, e.g., the sparsity of A and B , reorthogonalization in the Lanczosprocess and the Arnoldi process, and calls of the Lanczos basis vectors from secondarystorages when iterations used by GLTR are too many and secondary storages have tobe exploited.The other important issue that deserves enough attention is to establish correctand fair stopping criteria when comparing the overall efficiency of GLTR and IRA,IRRA. This issue has received no attention in the literature, and one commonly termi-nates them by simply setting user-prescribed tolerances for respective residual norms Z. JIA AND F. WANG in an irrelevant way [1, 34]. This will cause unreasonable or incorrect comparisonsbecause the residual norms themselves for GLTR and IRA, IRRA are completely dif-ferent and are generally not of similar or comparable size. That is, the same tolerancesfor the two kinds of algorithms do not mean that the TRS residual norm of approxi-mate solution recovered from the eigenvalue-based algorithms is the same as or verycomparable to the residual norm by algorithms for solving TRS directly, causing thattwo kinds of algorithms do not compute the converged approximate solutions of TRSwith very similar accuracy.In view of the above, as the premise of efficiency comparison of two kinds ofalgorithms, we should fix one stopping criterion for one kind of algorithm as thereference standard and then determine how small the other kind should be, so thatwe can finally compute approximate solutions with almost the same accuracy bythe two kinds of algorithms. In this paper, we shall establish a concise relationshipbetween the different stopping criteria for GLTR and IRA, IRRA. Based on it, weprescribe the stopping criteria for IRA and IRRA, then decide the one for GLTR,so that we can make a reasonable efficiency comparison. In the meantime, we willtheoretically investigate how fast GLTR, the Arnoldi method and the refined Arnoldimethod converge by analyzing the residual norms obtained by GLTR and by IRA andIRRA for the equivalent matrix eigenproblem.This paper is organized as follows. In Section 2, we give some preliminariesand review the GLTR method. In Section 3, we introduce the equivalence of thesolutions of (1.1) and the generalized eigenvalue problem of a certain 2 n × n matrixpair. We describe a general framework of an eigenvalue-based algorithm in [1], andintroduce IRA [26]. We then propose the refined Arnoldi method and develop thecorresponding IRRA algorithm [16] in the context of generalized eigenvalue problems.Section 4 is devoted to deriving relationships between the TRS residual norms andresidual norms of the equivalent generalized eigenvalue problem. In Section 5, weinvestigate the convergence of residual norms obtained by GLTR, the Arnoldi methodand the refined Arnoldi method for the equivalent eigenproblem of TRS. In Section6, we report numerical experiments to compare IRA and IRRA with GLTR, showingthat they are competitive with GLTR. Meanwhile, we compare IRRA with IRA anddemonstrate that IRRA is more efficient than IRA. Finally, we conclude the paper inSection 7.Throughout this paper, denote by the superscript T the transpose of a matrix orvector, by ∥ · ∥ the 2-norm of a matrix or vector, by I the identity matrix with orderclear from the context, and by e i the i th column of I . All vectors are column vectorsand are typeset in lowercase letters.
2. Preliminaraies.2.1. A solution to (1.1) . As we can see from Theorem 1.1, if (1.1) has no solu-tion with ∥ s opt ∥ B = ∆, then A is positive definite and s opt = − A − g with ∥ s opt ∥ B < ∆and λ opt = 0. Otherwise, ∥ s opt ∥ B = ∆ and λ opt >
0, the solution s opt to TRS (1.1)is unique and s opt = − ( A + λ opt B ) − g . All these correspond to the so-called “easycase” [2, 8, 20, 21] or “nondegenerate case” [2, 8, 20, 21] .If A is semi-positive definite or indefinite and g ⊥ N ( A − α n B ) , the null space of A − α n B , where α n is the leftmost eigenvalue of the matrix pair( A, B ), then we have the following definition [2, 8, 21].
Definition 2.1 (
Hard Case ). The solution of TRS (1.1) is a hard case if g is igenvalue-based algorithms for TRS orthogonal to the eigenspace corresponding to the eigenvalue α n of A − λB , and theoptimal Lagrangian multiplier is λ opt = − α n . In the hard case, TRS (1.1) may have multiple optimal solutions [21, p.87-88],which must reach the trust-region boundary and can be expressed by s opt = − ( A − α n B ) † g + ηu n , (2.1)where u n ∈ N ( A − α n B ), ∥ ( A − α n B ) † g ∥ B < ∆, and the superscript † denotes theMoore-Penrose generalized inverse.It is known that the scalar η satisfies η = ∆ − ∥ ( A − α n B ) † g ∥ B > ∥ u n ∥ B = 1 because of the B -orthogonality of ( A − α n B ) † g and u n .From (2.1), we see that there are infinitely many s opt whenever the eigenvalue α n of the matrix pair ( A, B ) is multiple, since N ( A − α n B ) is bigger than one. We notonly need to find the minimum B -norm solution to a symmetric semi-positive definitesingular system but also need to compute the eigenspace of the matrix pair ( A, B )associated with the eigenvalue α n . The hard case must be treated separately and hasbeen studied for years; see, e.g., [1, 5, 8, 20, 21, 24].As has been addressed in [2], the hard case rarely occurs in practice, as it requiresthat both A be semi-definite or indefinite and g be orthogonal to N ( A − α n B ). Inthe sequel, we are only concerned with the easy case. Withoutloss of generality, we always assume that ∥ s opt ∥ B = ∆ throughout the paper forotherwise A must be symmetric positive definite and the Conjugate Gradient (CG)method can solve (1.1). In this subsection, we briefly review the GLTR method [8].At iteration k , the GLTR method exploits the symmetric Lanczos process togenerate a B -orthonormal basis { p i } ki =0 of the underlying ( k + 1)-dimensional Krylovsubspace K k ( B − g, B − A ) = span { B − g, ( B − A ) B − g, . . . , ( B − A ) k B − g } , which can be written in matrix form AP k = BP k T k + β k +1 Bp k +1 e Tk +1 , (2.3) P Tk g = β e , β = ∥ B − g ∥ B = ∥ B − g ∥ , (2.4) g = β Bq , (2.5)where P k = ( p , p , . . . , p k ) is B -orthonormal, i.e., P Tk BP k = I , the matrix(2.6) T k = P Tk AP k = δ β β δ . . .. . . . . . . . .. . . δ k − β k β k δ k ∈ R ( k +1) × ( k +1) is symmetric tridiagonal, and e k +1 is the last column of the identity matrix of order k + 1. The GLTR method projects TRS (1.1) onto the ( k + 1)-dimensional Krylov Z. JIA AND F. WANG subspace K k ( B − g, B − A ) and yields the reduced TRS:(2.7) min ∥ h ∥≤ ∆ ϕ ( h ) = β e T h + 12 h T T k h. From Theorem 1.1, the vector h ( k ) is a solution to (2.7) if and only if there exists theoptimal Lagrangian multiplier λ ( k ) ≥ ∥ h ( k ) ∥ ∆ , ( T k + λ ( k ) I ) h ( k ) = − β e ,λ ( k ) (∆ − ∥ h ( k ) ∥ ) = 0 ,T k + λ ( k ) I ≽ . GLTR computes the iterate s ( k )gltr = P k h k to approximate s opt , which solves the pro-jected TRS(2.8) min s ∈K k ( B (cid:0) g;B (cid:0) A ) ; ∥ s ∥ B ≤ ∆ q ( s ) . The GLTR method exploits the Mor´e-Sorensen method [10] to solve the reduced TRS(2.7). It is shown in [8] that the B − -residual norm of the approximate solution( λ ( k ) , s ( k )gltr ) of (1.3) satisfies(2.9) ∥ ( A + λ ( k ) B ) s ( k )gltr ∥ B (cid:0) = β k +1 | e Tk +1 h ( k ) | , which is used as an efficient stopping criterion without forming s ( k )gltr explicitly untilthe residual norm drops below a user-prescribed tolerance.Importantly, if ∥ s opt ∥ B = ∆, we must have ∥ s ( k )gltr ∥ B = ∆ as k increases [8, 33]. Infact, it has been shown in [19] that ∥ s ( k )gltr ∥ B = ∆ after very few iterations. This hasbeen numerically confirmed to occur for k = 1 in [19, 33].In order to get an accurate approximation, one needs to expand the Krylov sub-space successively and store P k so as to form the iterate s ( k )gltr . In finite precisionarithmetic, numerical difficulties such as loss of orthogonality and excessive storagemay limit the applicability of the method. A complicated restarted variant of GLTRhas been proposed in [34] to reduce storage, and it has been numerically justified tobe comparable to GLTR in terms of MVs on some problems.
3. The equivalence of TRS and a certain matrix eigenvalue problemand the IRA, IRRA algorithms.3.1. The equivalence of TRS (1.1) and a matrix eigenvalue problem.
Define(3.1) M = ( − A gg T ∆ B − A ) ∈ R n × n , e B = ( B B ) ∈ R n × n . Adachi et al. [1] prove that if ∥ s opt ∥ B = ∆ then the solution of TRS (1.1) can beequivalently formulated as the rightmost eigenpair of the matrix pair ( M, e B ), as shownbelow. Theorem 3.1 ( [1] ). Let ( λ opt , s opt ) satisfy Theorem 1.1 with ∥ s opt ∥ B = ∆ . Thenthe rightmost eigenvalue µ of ( M, e B ) is real and simple, and the optimal Lagrangian igenvalue-based algorithms for TRS multiplier λ opt = µ . Let y = ( y T , y T ) T be an eigenvector of ( M, e B ) associated with µ , and suppose that g T y ̸ = 0 . Then the unique TRS solution is (3.2) s opt = − ∆ g T y y . Remark 3.1.
In (cid:12)nite precision arithmetic, it is shown in [1] that s opt can becomputed stably by the formula (3.3) s opt = − sign ( g T y )∆ y ∥ y ∥ B . Adachi et al. [1] have proved that g T y = 0 corresponds to the hard case, i.e., λ opt = − α n and g ⊥ N ( A − α n B ) . Therefore, in the easy case, g T y ̸ = 0 is guaran-teed, and (3.2) holds; furthermore, they have shown that g T y = 0 implies y = 0.Numerically, in order to distinguish the hard and easy cases, they design a threshold τ < ∥ y ∥ < τ is recognized as the hard case, where τ is chosen asthe level of square root of the machine precision.Algorithm 1 is a modified form of Algorithm 2 in [1] that forms a framework ofany eigenvalue-based algorithm for computing the solution s opt to (1.1) in the easycase. Algorithm 1
The eigenvalue-based algorithm in [1] (GEP)1. (Consider the case λ opt = 0) Solve As = − g for s by the CG algorithm, and s opt = s if ∥ s ∥ B < ∆. Otherwise, go to step 2.2. Use an eigensolver to compute the rightmost eigenvalue µ of the matrix pair( M, e B ) and the corresponding eigenvector y = ( y T , y T ) T .3. If ∥ y ∥ ≤ τ for a user-prescribed tolerance τ , then treat TRS as hard case andsolve it separately. Otherwise, s opt = − sign ( g T y )∆ y ∥ y ∥ B .Adachi et al. have exploited IRA, i.e., the Matlab function eigs , in step 2 of thisalgorithm and made numerical experiments on varying sized TRS’s to compare theoverall efficiency of IRA with that of the non-restarted GLTR. In the experiments,they set stopping criteria for the relative residual norms in both eigs and GLTR to bethe level of machine precision. The experiments have indicated that GLTR is a littlemore efficient. Comparisons have also been carried out in [34] for eigs , GLTR and theirrestarted GLTR as well as some others, where the stopping criterion for eigs was set tobe the default value, i.e., the level of machine precision, and the ones for B − -relativeresidual norms of GLTR and others applied to (1.1) were 10 − . However, as we havecommented in the introduction, such a comparison is problematic and unreasonablesince the residual norms obtained by GLTR and an eigenvalue-based algorithm aredefined differently, and for given user-prescribed stopping tolerances such as thosein [1, 34], two kinds of algorithms may compute converged approximate solutions ofTRS (1.1) with greatly different accuracy, so that the efficiency comparison is unfairand may be misleading. Therefore, it is vital to seek intimate relationships betweenthese two kinds of residual norms and fix a user-prescribed stopping criterion forone kind and then decide the correct one for the other kind. Only in this way canwe compute the approximate solution with very comparable accuracy and make areasonable efficiency comparison. We shall solve this problem in Section 4. Z. JIA AND F. WANG
In step 2 of Algorithm 1, one commonlyuses the IRA algorithm to carry out the task [1, 34]. In this subsection, we shallbriefly introduce some properties of the underlying Arnoldi method and the basedIRA algorithm that are suited for computing Ritz approximations to the rightmosteigenpairs of a large matrix (pair).Since B is symmetric positive definite, the matrix e B in (3.1) is too. We considerthe eigenproblem of ( M, e B ):(3.4) M x i = µ i e Bx i , where µ i , i = 1 , , . . . , n , are the eigenvalues labeled in descending order of their realparts: Re ( µ ) ≥ Re ( µ ) ≥ · · · ≥ Re ( µ n ) , and the x i are corresponding eigenvectors of e B -norm one. By Theorem 3.1, we areconcerned with ( µ , y ) := ( µ , x ), the rightmost eigenpair of ( M, e B ), when solvingTRS (1.1).Given an initial v of e B -norm one, the k -step Arnoldi process can be written inmatrix form: M V k = e BV k H k + h k +1 ;k e Bv k +1 e Tk = e BV k +1 e H k , (3.5)where V k = ( v , v , . . . , v k ) ∈ R n × k and V k +1 = ( V k , v k +1 ) are e B -orthonormal, H k ∈ R k × k is upper Hessenberg, and e H k is the ( k + 1) × k upper Hessenberg matrix whichis the same as H k except for an additional row whose only nonzero entry is h k +1 ;k inposition ( k + 1 , k ).Denote by ( µ ( k ) i , z ( k ) i ), i = 1 , , . . . , k , the eigenpairs of H k with ∥ z ( k ) i ∥ = 1, whichare labeled as Re ( µ ( k )1 ) ≥ Re ( µ ( k )2 ) ≥ · · · ≥ Re ( µ ( k ) k ) . Then the Arnoldi method uses some of µ ( k ) i , i = 1 , , . . . , k , called the Ritz values ofthe matrix pair ( M, e B ) with respect to the k -dimensional Krylov subspace(3.6) K k ( v , e B − M ) . = span { v , e B − M v , ( e B − M ) v , . . . , ( e B − M ) k − v } , to approximate some exterior µ i , e.g., the rightmost one(s), and the correspondingeigenvectors x i of ( M, e B ) are approximated by the Ritz vectors of e B -norm one:(3.7) x ( k ) i = V k z ( k ) i . The e B − -residual norm of any Ritz pair ( µ ( k ) i , x ( k ) i ) satisfies(3.8) ∥ ( M − µ ( k ) i e B ) x ( k ) i ∥ e B (cid:0) = h k +1 ;k | e Tk z ( k ) i | , which is used as a stopping criterion without forming x ( k ) i until the convergence occurs.However, an unfavorable aspect of the Arnoldi method is that one cannot knowin advance how many steps will be required before the eigenpairs of interest are wellapproximated by the Ritz pairs. In order to recover Ritz vectors, one is obliged either igenvalue-based algorithms for TRS V k or to recompute them. Therefore, the possible excessive storage andexpensive computational cost limit the efficiency and applicability of the method. Tothis end, restarting is generally necessary.Sorensen [26] has proposed an elegant implicit restarting technique: Supposethat some k specific eigenpairs are of interest and the ( k + p )-step Arnoldi process hasbeen performed and one has computed H k + p . Sorensen’s means does not restart theArnoldi process from the scratch. Rather, it implicitly restarts the Arnoldi processand realizes a new ( k + p )-step Arnoldi process with an updated initial vector. It saves k MVs with the matrix e B − M at much less cost of p implicitly shifted QR iterationson the small matrix H k + p and naturally generates a new k -step Arnoldi process bytruncating the implicitly shifted QR form. Then starting with iteration k + 1, thestandard Arnoldi process is performed to generate a new ( k + p )-step Arnoldi processon the updated ( k + p )-dimensional Krylov subspace. The resulting algorithm is theIRA (implicitly restarted Arnoldi) algorithm. One of the crucial points for the overallsuccess of IRA is proper selection of the p shifts involved. As a thumb, the shiftsshould be the best possible approximations to some unwanted eigenvalues of ( M, e B ).Within the framework of the Arnoldi method, Sorensen has suggested to use the p unwanted eigenvalues of H k + p as shifts, called exact shifts [26]. Since then exactshifts have been used in all public IRA-based softwares, e.g., the Matlab function eigs . Algorithm 2 sketches IRA for computing the k rightmost eigenpairs of ( M, e B ),where the superscript H denotes the conjugate transpose of a matrix. We commentthat Step 6 is performed using implicitly shifted QR iterations in real arithmetic inpractical computations even if some shifts are complex conjugates. All the Arnoldi based algorithms, such asIRA, may have convergence problems [14, 18]: The Ritz vectors obtained by themmay converge erratically and even may not converge even if the subspace contains suf-ficiently accurate approximations to the desired eigenvectors, while the correspondingRitz values converge. To correct this deficiency, Jia [15, 16] has proposed a refinedArnoldi method and developed an IRRA (implicitly restarted refined Arnoldi) al-gorithm. The refined Arnoldi method is mathematically different from the Arnoldimethod for extracting approximate eigenvectors from the underlying Krylov subspace,and the new eigenvector approximations, called the refined Ritz vectors, can be en-sured to converge unconditionally provided that the subspace is sufficiently accurate.Moreover, the refined Ritz vectors are more accurate than the Ritz vectors [17].In this subsection, we extend the refined Arnoldi method in [15] and the IRRAalgorithm in [16] for the standard eigenvalue problem to the generalized eigenvalueproblem of the matrix pair ( M, e B ). When adapting the refined Arnoldi method [15] tothe generalized one, rather than using Ritz vectors x ( k ) i as approximate eigenvectors,for each µ ( k ) i of interest we seek e B -normalized e x ( k ) i ∈ K k ( v , e B − M ) satisfying theoptimal condition(3.9) (cid:13)(cid:13)(cid:13)( M − µ ( k ) i e B ) e x ( k ) i (cid:13)(cid:13)(cid:13) e B (cid:0) = min x ∈K k ( v ; e B (cid:0) M ) ; ∥ x ∥ e B =1 (cid:13)(cid:13)(cid:13)( M − µ ( k ) i e B ) x (cid:13)(cid:13)(cid:13) e B (cid:0) , and use it to approximate the desired eigenvector x i . The approximation ˜ x ( k ) i is calledthe refined Ritz vector, or more generally, the refined approximate eigenvector, fromthe underlying subspace associated with µ ( k ) i .The following results are from [15] for e B = I , and can be straightforwardlyextended to the refined Arnoldi method applied to ( M, e B ) by making use of the0 Z. JIA AND F. WANG
Algorithm 2
IRA1. Choose an initial vector v with ∥ v ∥ e B (cid:0) = 1, the number k of desired rightmosteigenpairs of ( M, e B ), the number p of the shifts, and a stopping tolerance ϵ > k -step Arnoldi process for ( M, e B ) with the starting vector v : M V k = e BV k H k + r k e Tk .
3. Apply p additional steps of the Arnoldi process to obtain the ( k + p )-step Arnoldidecomposition: M V k + p = e BV k + p H k + p + e Br k + p e Tk + p ;4. Compute the eigenpairs of H k + p , and select its k eigenvalues with the largest realparts to approximate the desired eigenvalues µ i , i = 1 , , . . . , k of M and theremaining p ones as shifts σ , σ , . . . , σ p ;5. Q = I ;6. for j = 1 , , . . . , pQR factorization: H k + p − σ j I = Q j R j ; H k + p = Q Hj H k + p Q j ; Q = QQ j ;end7. H k = H k + p (1 : k, k );8. V k = V k + p Q (: , k );9. r k = H k + p ( k + 1 , k ) V k + p Q (: , k + 1) + Q ( k + p, k ) r k ;10. If the right-hand sides of (3.8) drop below ϵ for i = 1 , , . . . , k , stop, and take( µ ( k ) i , x ( k ) i ) in (3.7) as approximations to ( µ i , x i ); otherwise, go to Step 3.Arnoldi process (3.5). Theorem 3.2 ( [15] ). Let e z ( k ) i be the right singular vector of e H k − µ ( k ) i ˜ I associatedwith σ min ( e H k − µ ( k ) i ˜ I ) , the smallest singular value of e H k − µ ( k ) i ˜ I , where ˜ I is the sameas I except for an additional zero row. Then e x ( k ) i = V k e z ( k ) i , (3.10) (cid:13)(cid:13)(cid:13)( M − µ ( k ) i I ) e x ( k ) i (cid:13)(cid:13)(cid:13) e B (cid:0) = σ min ( e H k − µ ( k ) i ˜ I ) . (3.11)This theorem indicates that one can reliably and efficiently compute the refinedRitz vectors e x ( k ) i by computing the singular value decomposition of a small sized e H k − µ ( k ) i ˜ I . Particularly, (3.11) is used as the stopping criterion without explicitlyforming e x ( k ) i until the convergence occurs.The IRRA algorithm [16] applies the implicit restarting technique by Sorensen[26] to the Arnoldi process straightforwardly. A key contribution in [16] consists inproviding a selection of better shifts, called the refined shifts, based on the informationavailable during the refined Arnoldi method. The refined shifts have been theoreticallyshown and numerically justified to be better than the exact shifts used in IRA in thesense that they generate a better updated ( k + p )-dimensional Krylov subspace duringeach cycle. Algorithm 3 sketches IRRA that computes the k rightmost eigenpairs of igenvalue-based algorithms for TRS M, e B ). We point out that Steps 5–6 are managed to be run only in real arithmeticand Step 8 is performed by using implicitly shifted QR iterations in real arithmeticeven if there are complex conjugate shifts [18]. Algorithm 3
IRRA1. Choose an initial vector v with ∥ v ∥ e B (cid:0) = 1, the number k of desired rightmosteigenpairs of ( M, e B ), the number p of the shifts, and a stopping tolerance ϵ > k -step Arnoldi process for ( M, e B ) with the starting vector v : M V k = e BV k H k + r k e Tk .
3. Apply p additional steps of the Arnoldi process to obtain the ( k + p )-step Arnoldidecomposition: M V k + p = e BV k + p H k + p + e Br k + p e Tk + p ;4. Compute the eigenvalues µ ( k ) i of H k + p and the vectors e z ( k + p ) i in Theorem 3.2corresponding to the k eigenvalues µ ( k ) i with the largest real parts, and define thematrix Z k = ( e z ( k + p )1 , e z ( k + p )2 , . . . , e z ( k + p ) k );5. Compute the full QR factorization: Z k = ( U k , b U k ) ( R k × O b R k ) ;6. Compute the eigenvalues of the matrix b U Hk H k + p b U k as shifts σ , σ , . . . , σ p ;7. Q = I ;8. for j = 1 , , . . . , pQR factorization: H k + p − σ j I = Q j R j ; H k + p = Q Hj H k + p Q j ; Q = QQ j ;end9. H k = H k + p (1 : k, k );10. V k = V k + p Q (: , k );11. r k = H k + p ( k + 1 , k ) V k + p Q (: , k + 1) + Q ( k + p, k ) r k ;12. If the right-hand sides of (3.11) drop below ϵ for i = 1 , , . . . , k , stop, and take( µ ( k ) i , e x ( k ) i ) in (3.10) as approximations to ( µ i , x i ); otherwise, go to Step 3.For large scale TRS’s, in Step 2 of Algorithm 1, we can use IRA and IRRAto compute only one rightmost eigenpair ( µ , x ) of ( M, e B ). Define y = x . Wethen recover s opt by the formula (3.3). The resulting algorithms are abbreviated asTRS IRA and TRS IRRA, respectively.
4. Relationships between the TRS residual norms and the eigenvalue-based ones.
For B ̸ = I , by a change of variables: b A ← B − AB − , b g ← B − g, and c ← B s, TRS (1.1) can be equivalently transformed into the standard TRS with B = I :(4.1) min ∥ c ∥≤ ∆ b q ( c ) with b q ( c ) = b g T c + 12 c T b Ac,
Let c opt be the solution to (4.1) and s ( k )gltr and c ( k )gltr be the GLTR iterates appliedto the subspaces K k ( g, B − A ) and K k ( b g, b A ), respectively. Then it is shown in [19, 33]2 Z. JIA AND F. WANG that GLTR generates the identical Lagrangian multipliers λ ( k ) for (1.1) and (4.1) and s ( k )gltr = B − c ( k )gltr , (4.2) ∥ s ( k )gltr − s opt ∥ B = ∥ c ( k )gltr − c opt ∥ , (4.3) q ( s ( k )gltr ) − q ( s opt ) = b q ( c ( k )gltr ) − b q ( c opt ) , (4.4) ∥ ( A + λ ( k ) B ) s ( k )gltr + g ∥ B (cid:0) = ∥ ( b A + λ ( k ) I ) c ( k )gltr + b g ∥ , (4.5)More generally, independently of a specific solver for TRS (1.1) directly, let ˆ c bean approximate solution of (4.1), and define ˆ s = B − ˆ c to be the approximate solutionof (1.1). Then for a given ˆ λ ≥ ∥ ( A + ˆ λB )ˆ s + g ∥ B (cid:0) = ∥ ( b A + ˆ λI )ˆ c + b g ∥ . Define c M = e B − M e B − . Then the eigenvalues of c M are identical to the eigenvalues µ i , i = 1 , , . . . , n of( M, e B ), and its eigenvectors w i are related to the eigenvectors x i of the matrix pairby x i = e B − w i , from which it follows that ∥ x i ∥ e B = 1 if ∥ w i ∥ = 1 and vice versa.Mathematically, the Arnoldi method applied to K k ( v , e B − M ) and K k ( e B − v , c M )compute the same Ritz values µ ( k ) and the Ritz vectors x ( k ) and w ( k ) , respectively,where x ( k ) = e B − w ( k ) , and their residual norms satisfy(4.7) ∥ ( M − µ ( k ) e B ) x ( k ) ∥ e B (cid:0) = ∥ ( c M − µ ( k ) I ) w ( k ) ∥ . More generally, independently of any specific eigensolver for the eigenproblem of( M, e B ), let (ˆ λ, ˆ w ) be an approximation to the rightmost eigenpair ( λ opt , w ) of c M , anddefine ˆ y = e B − ˆ w to be the corresponding approximate eigenvector of ( M, e B ). Then(4.8) ∥ ( M − ˆ λ e B )ˆ y ∥ e B (cid:0) = ∥ ( c M − ˆ λI ) ˆ w ∥ . Based on (3.3), we recover the approximation ˆ s of the solution s opt from ˆ y as follows:Let ˆ y = (ˆ y T , ˆ y T ) T . Then the approximation ˆ s of s opt is recovered by(4.9) ˆ s = − sign ( g T ˆ y )∆ ˆ y ∥ ˆ y ∥ B . In terms of (4.5), (4.6) and (4.7), (4.8), clearly it suffices to consider the standardTRS (1.1) with B = I and its mathematically equivalent eigenproblem when seekingpossible intimate relationships between the residual norms of approximate solutionsof TRS and those of approximate eigenpairs of M . Based on them, we are able todesign reliable and proper stopping tolerances for GLTR and TRS IRA, TRS IRRAin order that they deliver converged approximate solutions of TRS (1.1) with verycomparable accuracy, so that we can make a reasonable and fair comparison of theoverall efficiency of the three algorithms.Notice that a common stopping criterion for any iterative algorithm that solves(1.1) requires that the TRS residual norm satisfy(4.10) ∥ ( A + ˆ λI )ˆ s + g ∥ ≤ tol igenvalue-based algorithms for TRS tol , where ˆ λ and ˆ s are obtained approximations to λ opt and s opt ,respectively. Correspondingly, for a large M , a stopping criterion for any eigenvalue-based iterative projection algorithm is(4.11) (cid:13)(cid:13)(cid:13) ( M − ˆ λI )ˆ y (cid:13)(cid:13)(cid:13) ≤ tol , with ∥ ˆ y ∥ = 1 and tol a given tolerance. Here we exclude those numerically backwardstable algorithms such as the QR algorithm for small to moderate sized matrix eigen-value problems, which automatically compute the eigenpair approximations at thelevel of machine precision in the sense of backward errors, i.e., tol = ∥ M ∥O ( ϵ mach )with ϵ mach being the machine precision. Whenever (4.11) is met, we terminate thealgorithm used, and recover the converged solution of TRS (1.1) as follows: FromTheorem 3.1 and (3.3), ˆ λ is an approximation of λ opt , and the converged approximatesolution(4.12) ˆ s = − ∆ g T ˆ y ˆ y = − sign ( g T ˆ y )∆ ˆ y ∥ ˆ y ∥ . As we have stressed previously, in order to fairly compare the efficiency of an iter-ative algorithm for solving (1.1) directly and an eigenvalue-based one, a key premiseis that two kinds of algorithms must compute the converged solutions with almostthe same accuracy. But this crucial issue has been ignored in the literature, e.g.,[1, 34]. Since one commonly measures solution accuracy by the residual norm of anapproximate solution, the comparison premise becomes the following: Given an ap-proximation (ˆ λ, ˆ y ), suppose that it satisfies criterion (4.11) and ˆ s is recovered fromˆ y by the formula (4.12). Then how small is the corresponding TRS residual norm ∥ ( A + ˆ λI )ˆ s + g ∥ , that is, how should tol in (4.10) be chosen? The other side of theproblem is: Suppose that (ˆ λ, ˆ s ) satisfies (4.10). Then what should tol be chosen?Note that an eigenvalue-based algorithm can recover ˆ s from ˆ y but an algorithm, e.g.,GLTR, that solves (1.1) directly does not provide an approximation ˆ y of the eigen-vector y of ( M, e B ). Therefore, it is only possible to decide tol from tol , but notthe other way around. To this end, we naturally attempt to select an appropriatevalue tol for a given tol . This goal can be achieved, provided that we can find closerelationships between ∥ ( A + ˆ λI )ˆ s + g ∥ and ∥ ( M − ˆ λI )ˆ y ∥ .We can present the following basic result. Theorem 4.1.
With the previous notation, we have ∥ ( A + ˆ λI )ˆ s + g ∥ ≤ ∆ ∥ ˆ y ∥ ∥ ( M − ˆ λI )ˆ y ∥ . (4.13) Proof.
Making use of (4.12), we obtain (cid:13)(cid:13)(cid:13)(cid:13) ( M − ˆ λI ) ( ˆ y ˆ y )(cid:13)(cid:13)(cid:13)(cid:13) = ( ∥ ˆ y ∥ ∆ ) ∥ ( A + ˆ λI )ˆ s + g ∥ + (cid:13)(cid:13)(cid:13) ˆ y − ( A + ˆ λI )ˆ y (cid:13)(cid:13)(cid:13) . (4.14)Therefore, we have ( ∥ ˆ y ∥ ∆ ) ∥ ( A + ˆ λI )ˆ s + g ∥ ≤ ∥ ( M − ˆ λI )ˆ y ∥ , Z. JIA AND F. WANG meaning that ∥ ( A + ˆ λI )ˆ s + g ∥ ≤ ∆ ∥ ˆ y ∥ ∥ ( M − ˆ λI )ˆ y ∥ . Remark 4.1.
For B ̸ = I , based on (4.5), (4.6) and (4.7), (4.8), relation (4.13) becomes ∥ ( A + ˆ λB )ˆ s + g ∥ B (cid:0) ≤ ∆ ∥ ˆ y ∥ B ∥ ( M − ˆ λB )ˆ y ∥ B (cid:0) . (4.15) Remark 4.2.
This theorem indicates that, given the stopping tolerance tol := tol in (4.11) for an eigenvalue-based algorithm, the TRS stopping tolerance tol in (4.10)for an algorithm that solves (1.1) directly should be chosen as (4.16) tol = ∆ ∥ ˆ y ∥ tol. Notice that the two terms of the right-hand side in (4.14) are the squared norms ofupper and lower parts of the left-hand side residual. They should be comparable insize. Therefore, the bound (4.13) is very sharp and almost an equality within roughlya multiple of two, showing that tol de(cid:12)ned by (4.16) is very reliable. This suggestsus to use such a choice in GLTR for a given stopping tolerance in IRA and IRRA. Remark 4.3.
For ˆ y very small, TRS is numerically in the hard case. (4.13)and the above remark show that the residual norm ∥ ( A + ˆ λI )ˆ s + g ∥ may not be small,meaning that ˆ s is not a meaningful approximation of s opt . For example, if ∥ ˆ y ∥ is assmall as ∥ ( M − ˆ λI )ˆ y ∥ , ∥ ( A + ˆ λI )ˆ s + g ∥ = O (∆) . This implies that the eigenvalue-basedalgorithms for computing the rightmost eigenpair of ( M, e B ) is not capable of solvingthe TRS in the hard case, as expected. Remark 4.4.
In (cid:12)nite precision arithmetic, stopping criteria for relative residualnorms must be used for reliability and general purpose. We terminate an eigenvalue-based algorithm whenever ∥ ( M − ˆ λI )ˆ y ∥∥ M ∥ ≤ tol is met for a user-prescribed tolerance tol , where ∥ · ∥ is the 1-norm of a matrix.Such criterion is independent of the scaling of M and allows tol to be at the level of ϵ mach and the premature or too late termination cannot occur. Correspondingly, thestopping criterion for algorithms that solve TRS directly becomes (4.17) ∥ ( A + ˆ λI )ˆ s + g ∥∥ M ∥ ≤ ∆ ∥ ˆ y ∥ tol. Whenever ∥ ˆ y ∥ is fairly small for a modest ∆ , the relative residual norm of the recov-ered ˆ s cannot attain the level of ϵ mach . This may be a limitation of eigenvalue-basedalgorithms for TRS.
5. The convergence of the residual norms obtained by GLTR, theArnoldi method and the re(cid:12)ned Arnoldi method for solving (1.1) with B = I . In this section, we investigate the residual convergence of non-restartedGLTR, the Arnoldi method and the refined Arnoldi method for solving TRS (1.1)with B = I . The results are straightforwardly applicable to the TRS with B ̸ = I byexploiting (4.6), (4.8) and the fact that the eigenvalues of c M and the pair ( M, e B ) arethe same. igenvalue-based algorithms for TRS The authors have studied theresidual convergence of the GLTR method in [19], where the following result is estab-lished.
Theorem 5.1.
Suppose ∥ s opt ∥ = ∥ s ( k )gltr ∥ = ∆ , and de(cid:12)ne η k = ∥ g ∥ ( α + λ opt ) ∥ g ∥ + ( α + λ opt ) ∆ , η k = 2( α + λ opt ) ∥ g ∥ + ( α + λ opt ) ∆ , where α and α n are the rightmost and leftmost eigenvalues of A , respectively. Then ∥ ( A + λ ( k ) I ) s ( k )gltr + g ∥ ≤ ( η k ∆ ∥ g ∥ + 8( α + λ opt ) η k ∆ ) ( √ κ − √ κ + 1 ) k +1) (5.1) + 4 √ κ ∆( α + λ opt ) ( √ κ − √ κ + 1 ) k +1 , where κ = (cid:11) + (cid:21) opt (cid:11) n + (cid:21) opt . Remark 5.1.
The second term of the right-hand side in (5.1) dominates thebound soon as k increases, and ∥ ( A + λ ( k ) I ) s ( k )gltr + g ∥ decreases at least like ( √ (cid:20) − √ (cid:20) +1 ) k +1 .The numerical experiments have con(cid:12)rmed that the bound is very sharp and can bealmost attainable. Recallthat ( λ opt , y ) := ( µ , x ) is the rightmost eigenpair of M and y = ( y T , y T ) T . Inthe Arnoldi method, λ ( k ) = µ ( k )1 is the rightmost Ritz value of M with respect to K k ( v , M ) that is used to approximate λ opt = µ , and y ( k ) := x ( k )1 is the correspondingRitz vector of M that approximates y . Let y ( k ) = ( y ( k )1 y ( k )2 ) . Then in terms of (4.12), an approximation s ( k )A of s opt is obtained by(5.2) s ( k )A = − ∆ g T y ( k )2 y ( k )1 = − sign ( g T y ( k )2 )∆ y ( k )1 ∥ y ( k )1 ∥ . Let π k = V k V Tk be the orthogonal projector onto K k ( v , M ). Then π k M π k is therestriction of M to K k ( v , M ) and its matrix representation is H k defined in (3.5). Adirect application of Theorem 3.8 in [14] to our context gives the following result onthe convergence of µ ( k )1 to µ = λ opt . Lemma 5.2.
Suppose that ∥ s opt ∥ = ∥ s ( k )A ∥ = ∆ , and let ε k = sin \ ( y, K k ( v , M )) with y the eigenvector of M associated with µ . Then for ε k small it holds that (5.3) | µ − µ ( k )1 | ≤ κ ( µ ( k )1 ) γ k ε k + O ( ε k ) , where κ ( µ ( k )1 ) is the spectral condition number of µ ( k )1 and γ k = ∥ π k M ( I − π k ) ∥ . In Theorem 3.8 of [14], sin \ ( y; K k ( v ; M )) in the right-hand side of (5.3) is tan \ ( y; K k ( v ; M )),but it is obvious that the sine and tangent can be replaced each other in the right-hand side whensin \ ( y; K k ( v ; M )) becomes small. Z. JIA AND F. WANG
Recall that ( µ ( k )1 , z ( k )1 ) is the rightmost eigenpair of H k and y ( k ) = V k z ( k )1 . Let thecolumns of Z ( k ) ⊥ be an orthonormal basis of the orthogonal complement of span { z ( k )1 } so that ( z ( k )1 , Z ( k ) ⊥ ) is orthogonal. Then we have(5.4) ( ( z ( k )1 ) T ( Z ( k ) ⊥ ) T ) M k ( z ( k )1 , Z ( k ) ⊥ ) = ( µ ( k )1 f Tk C k ) , where f Tk = ( z ( k )1 ) T M k Z ( k ) ⊥ and C k = ( Z ( k ) ⊥ ) T M k Z ( k ) ⊥ . Note that the eigenvalues of C k are the Ritz values other than µ ( k )1 of M with respect to K k ( v , M ).In our notation, it is direct from Theorem 3.2 in [18] to obtain the following result. Lemma 5.3 ( [18] ). With the previous notation, assume that sep( µ , C k ) > .Then sin \ ( y ( k ) , y ) ≤ ( ∥ M ∥ √ − ε k sep( µ , C k ) ) ε k , (5.5) where sep( µ , C k ) = σ min ( C k − µ I ) is the separation of µ and C k . Now we can establish the main result.
Theorem 5.4.
With the previous notation, assume that ∥ s opt ∥ = ∥ s ( k )A ∥ = ∆ ,and sep( µ , C k ) > . Then ∥ ( A + µ ( k )1 I ) s ( k )A + g ∥ ≤ c k ε k + O ( ε k ) , (5.6) where c k = ∆ ∥ y ( k )1 ∥ ( ∥ M − µ ( k )1 I ∥ ( ∥ M ∥ √ − ε k sep( µ , C k ) ) + | α | κ ( µ ( k )1 ) γ k ) with α = y T y ( k ) .Proof. Note that y and y ( k ) are real with ∥ y ∥ = ∥ y ( k ) ∥ = 1. Sincesin \ ( y ( k ) , y ) = ∥ ( I − yy T ) y ( k ) ∥ = ∥ y ( k ) − ( y T y ( k ) ) y ∥ = ∥ y ( k ) − αy ∥ and M y = µ y , we obtain ∥ ( M − µ ( k )1 I ) y ( k ) ∥ = ∥ ( M − µ ( k )1 I )( y ( k ) − αy ) + α ( µ − µ ( k )1 ) y ∥≤ ∥ ( M − µ ( k )1 I )( y ( k ) − αy ) ∥ + ∥ α ( µ − µ ( k )1 ) y ∥≤ ∥ M − µ ( k )1 I ∥∥ y ( k ) − αy ∥ + | α || µ − µ ( k )1 | = ∥ M − µ ( k )1 I ∥ sin \ ( y ( k ) , y ) + | α || µ − µ ( k )1 | . Substituting bound (5.3) for | µ − µ ( k )1 | and bound (5.5) for sin \ ( y ( k ) , y ) into theabove relation, we obtain ∥ ( M − (cid:22) ( k )1 I ) y ( k ) ∥ ≤ ( ∥ M − (cid:22) ( k )1 I ∥ ( ∥ M ∥ √ − " k sep( (cid:22) ; C k ) ) + | (cid:11) | (cid:20) ( (cid:22) ( k )1 ) (cid:13) k ) " k + O ( " k ) : From Theorem 4.1, ∥ ( A + µ ( k )1 I ) s ( k )A + g ∥ ≤ ( ∆ ∥ y ( k )1 ∥ ) ∥ ( M − µ ( k )1 I ) y ( k ) ∥ . igenvalue-based algorithms for TRS c k and ε k . As we can see, c k involves a few factors,where sep( µ , C k ) ≈ sep( µ ( k )1 , C k ) and | α | ≈ y ( k ) → y , and ∥ y ( k )1 ∥ = ∥ z ( k )1 ∥ . Weshall discuss ε k later on. Recall that e y ( k ) := e x ( k )1 is the refined Ritz vector of M from K k ( v , M ) to approximate y . Regarding the residual norms of the Ritz pair ( µ ( k )1 , y ( k ) ) and the refined Ritz pair( µ ( k )1 , e y ( k ) ), Theorem 4.1 of [17] proves that the strict inequality holds: ∥ ( M − µ ( k )1 I ) e y ( k ) ∥ < ∥ ( M − µ ( k )1 I ) y ( k ) ∥ , provided that ( µ ( k )1 , y ( k ) ) ̸ = ( µ , y ), the exact rightmost eigenpair of M . This meansthat, in terms of residual norms, e y ( k ) is generally more accurate than y ( k ) as anapproximation to y .Let e y ( k ) = ( e y ( k )1 e y ( k )2 ) . Then the refined Arnoldi method computes the approximation s ( k )RA of s opt by theformula(5.7) s ( k )RA = − ∆ g T e y ( k )2 e y ( k )1 = − sign ( g T e y ( k )2 )∆ e y ( k )1 ∥ e y ( k )1 ∥ . The following result from [18, Theorem 7.1] estimates the residual norm ∥ ( M − µ ( k )1 I ) e y ( k ) ∥ . Lemma 5.5.
With the previous notation, we have ∥ ( M − µ ( k )1 I ) e y ( k ) ∥ ≤ ∥ M − µ ( k )1 I ∥ ε k + | µ − µ ( k )1 | √ − ε k , (5.8) where ε k = sin \ ( y, K k ( v , M )) . Combining (4.13) with (5.8) yields the following result on the TRS residual normimmediately.
Theorem 5.6.
With the previous notation, we have ∥ ( A + µ ( k )1 I ) s ( k )RA + g ∥ ≤ ( ∥ M − µ ( k )1 I ∥ + κ ( µ ( k )1 ) γ k ) ∆ ∥ e y ( k )1 ∥ √ − ε k ε k + O ( ε k ) . (5.9)Theorem 5.1 indicates that ∥ ( A + λ ( k ) I ) s ( k )gltr + g ∥ decreases at least as fast as ( √ (cid:20) − √ (cid:20) +1 ) k +1 , that is to say, κ = (cid:11) + (cid:21) opt (cid:11) n + (cid:21) opt is the main factor affecting the convergenceof the residual norm in the GLTR method. Theorem 5.4 and Theorem 5.6 show thatconvergence of the TRS residual norms by the Arnoldi method and the refined Arnoldimethod strongly relies on ε k = sin \ ( y, K k ( v , M )).A number of estimates on ε k have been available in [14, 15] for a general ma-trix unsymmetric M , which is directly applicable to our specific context and thus8 Z. JIA AND F. WANG omitted. Unfortunately, it is impossible to draw any definitive conclusion on whichof ( √ (cid:20) − √ (cid:20) +1 ) k +1 and ε k decays faster. As is seen from the results in [14, 15], amongothers, how fast ε k decays critically depends on the eigenvalue distribution of M , andthe better separated is µ from µ , the faster ε k decays. Except the fact that therightmost eigenvalue µ = λ opt , nothing is known on the other eigenvalues of M andthe relationships between them and the eigenvalues of A . Theoretically, it is wellpossible that the GLTR method converges faster in some cases but the converse maybe true in other cases for the Arnoldi method and the refined Arnoldi method. Ourlater numerical experiments will confirm this. Table 1
Test matrices. name n nonzeros sparsity application bcsstk29 13,992 619,488 0.3164% Structural Problembcsstk33 8,738 591,904 0.7752% Structural Problempcrystk02 13,965 968,583 0.4967% Duplicate Materials Problemgaron2 13,535 373,235 0.2037% Computational Fluid Dynamics Problemlhr11 10,964 231,806 0.1928% Chemical Process Simulation Problemlhr11c 10,964 233,741 0.1944% Chemical Process Simulation Problembarth5 15,606 61,484 0.0253% Duplicate Structural Problemnemeth01 9,506 725,054 0.8024% Theoretical/Quantum Chemistry Problem Sequencenemeth17 9,506 629,620 0.6968% Subsequent Theoretical/Quantum Chemistry Problempkustk02 10,800 810,000 0.6944% Structural Problemnd3k 9,000 3,279,690 4.0490% 2D/3D Problemcoupled 11,341 97,193 0.0756% Circuit Simulation ProblemPres Poisson 14,822 715,804 0.3258% Computational Fluid Dynamics Problemstokes64s 12,546 140,034 0.0890% Computational Fluid Dynamics Problemtuma2 12,992 49,365 0.0292% 2D/3D Problemopt1 15,449 1,930,655 0.8089% Structural Problemramage02 16,830 2,866,352 1.012% Computational Fluid Dynamics ProblemSi5H12 19,896 738,598 0.1866% Theoretical/Quantum Chemistry Problemnet25 9,520 401,200 0.4427% Optimization Problempf2177 9,728 725,144 0.7663% Optimization Problemc-41 9,769 101,635 0.1065% Optimization Problemc-44 10,728 85,000 0.0739% Optimization ProblemPGPgiantcompo 10,680 48,632 0.0426% Undirected Multigraphwing nodal 10,937 150,976 0.1262% Undirected Graph
6. Numerical experiments.
We now report numerical experiments to comparethe overall performance of the three algorithms: TRS IRA, TRS IRRA and GLTR.All the experiments were performed on an Intel Core (TM) i7 with CPU 3.6GHz,and 8 GB memory under the Microsoft Windows 10 64-bit system with the machineprecision ϵ mach = 2 . × − . We used the Matlab 2020a function eigs of IRA[26] and Matlab code of IRRA [18], called reigs , and the routine galahad gltr of theGALAHAD library [9], a Fortran 90 code of GLTR. igenvalue-based algorithms for TRS is available on the internet, wherethe Matlab function eigs is used to solve the eigenvalue problem. In our numericalexperiments, we have taken different subspace dimensions m = k + p = 20 , , , eigs and reigs . We take the stopping tolerance tol = 10 − for the relative residualnorms of approximate eigenpairs of ( M, e B ) and the maximum number of restarts tobe 600 in eigs and reigs . We use (4.17) to set the stopping tolerance of GLTR. Due tothe paper length, we only report the results for the case that the subspace dimensionis 30 since we have observed similar phenomena for the dimension m = 20 , ,
70. Itshould be noted that m − ( k + 3) shifts are used in the codes eigs and reigs . In ourcontext, k = 1, and each cycle of eigs and reigs costs m − galahad gltr , we have used the double precision arithmetic, set the maximumnumber of iterations to be the problem dimension n , and performed the Lanczosprocess without reorthogonalization.The matrices A in our test problems are A = G + G T , where the matrices G are taken from the Suite Sparse Matrix Collection [3] with nine examples for testinga Lanczos method for the extreme Lorentz eigenvalue problem and fifteen matricesused in other applications [34], such as the fast factorization pivoting methods forsparse symmetric indefinite systems [25]. Some characteristics of these matrices arelisted in Table 1, where nonzeros is the number of nonzero entries, and the sparsityis nonzeros/n . We take the matrix B = I and the symmetric positive definitetridiagonal matrix B = tridiag(1 , , g to be a unit length vector generated by the Matlabfunction randn(n,1) , and the trust region radii are ∆ = 1 and ∆ = 100, respectively.We compute the relative residual normRes . = ∥ ( A + ˆ λB )ˆ s + g ∥ B (cid:0) ∥ g ∥ B (cid:0) , (6.1)where (ˆ λ, ˆ s ) is a converged approximation obtained by each of GLTR, TRS IRA andTRS IRRA using our stopping criteria.The overall efficiency of each algorithm is measured by MVs. We should pointout that, making use of the structures of M and B defined by (3.1), a matrix-vectorproduct in eigs and reigs is two MVs in GLTR. That is, taking MVs in GLTR asreference, each iteration of eigs and reigs from k + 3 = 4 upwards costs two MVs. Inthe later tables, two MVs used in eigs and reigs precisely means one iteration. Weshall report ratio . = MVs IRA − MVs
IRRA
MVs
IRA , where MVs IRA and MVs
IRRA denote the MVs used by IRA and IRRA, respective-ly. This quantity measures the percentage of the overall efficiency improvement ofTRS IRRA over TRS IRA.Tables 2–3 report the results on the test problems with B = I and ∆ = 1 , (cid:24) nakatsukasa/codes/TRSgep.m. Z. JIA AND F. WANG
Table 2
Numerical results on sparse A = G + G T with G in Table 1 and B = I , ∆ = 1 . Matrix GLTR TRS IRA TRS IRRAMVs Res MVs Res MVs Res ratiobcsstk29 371 9 . e −
08 1226 2 . e −
08 1112 1 . e −
08 9 . . e −
08 528 1 . e −
08 456 9 . e −
09 13.6%pcrystk02 245 8 . e −
08 560 1 . e −
08 500 3 . e −
08 10.7%garon2 213 1 . e −
10 316 2 . e −
11 292 2 . e −
11 7.6%lhr11 137 2 . e −
09 320 1 . e −
10 296 1 . e −
10 7.5%lhr11c 137 2 . e −
09 320 3 . e −
10 296 1 . e −
10 7.5%barth5 189 5 . e −
09 320 1 . e −
10 296 3 . e −
10 7.5%nemeth01 247 1 . e −
10 488 2 . e −
11 396 5 . e −
11 23.0%nemeth17 47 3 . e −
12 86 1 . e −
12 80 3 . e −
12 7.0%pkustk02 181 2 . e −
06 422 1 . e −
06 368 1 . e −
06 12.8%nd3k 193 5 . e −
10 442 1 . e −
10 424 8 . e −
11 4.1%coupled 345 2 . e −
08 794 4 . e −
09 664 4 . e −
09 16.4%Pres Poisson 127 8 . e −
11 264 1 . e −
11 232 7 . e −
11 12.1%stokes64s 169 1 . e −
10 404 7 . e −
12 360 1 . e −
11 10.9%tuma2 127 2 . e −
11 200 2 . e −
12 180 7 . e −
12 9.0%opt1 167 7 . e −
08 296 5 . e −
08 296 2 . e −
08 0.0%ramage02 175 8 . e −
08 734 1 . e −
08 664 1 . e −
08 9.5%Si5H12 215 3 . e −
10 372 4 . e −
11 360 6 . e −
11 3.2%net25 31 3 . e −
11 128 3 . e −
11 108 2 . e −
11 15.6%pf2177 285 1 . e −
08 660 4 . e −
09 620 7 . e −
09 6.1%c-41 363 2 . e −
07 1088 1 . e −
07 1000 4 . e −
07 8.1%c-44 267 2 . e −
07 788 4 . e −
08 758 2 . e −
08 3.8%PGPgiantcompo 101 6 . e −
09 216 6 . e −
10 194 9 . e −
10 10.2%wing nodal 299 8 . e −
09 584 4 . e −
10 560 1 . e −
10 4.1%nemeth01 and tuma2 with ∆ = 1, TRS IRA and TRS IRRA use considerably feweriterations than GLTR does. These results have confirmed that there is no consid-erable winner between GLTR and the eigenvalue-based algorithms TRS IRA andTRS IRRA. Regarding TRS IRA and TRS IRRA, it is clear that TRS IRRA is moreand can be considerably more efficient than TRS IRA. A careful calculation shows thatthe average improvements are 9.15% for ∆ = 1 and 10.88% for ∆ = 100, respectively.Furthermore, the recorded Res’s have important implications. When it is roughlyequal to our stopping tolerance 10 − for TRS IRA and TRS IRRA, ∆ / ∥ ˆ y ∥ = O (1),meaning that ∥ ˆ y ∥ is not small and TRS is indeed in the easy case. If Res is muchless than 10 − , then ∥ ˆ y ∥ is quite small, meaning that the solution accuracy of TRSis impaired though TRS is still in the easy case. igenvalue-based algorithms for TRS Table 3
Numerical results on sparse A = G + G T with G in Table 1 and B = I , ∆ = 100 . Matrix GLTR TRS IRA TRS IRRAMVs Res MVs Res MVs Res ratiobcsstk29 329 9 . e −
04 1224 8 . e −
05 1204 3 . e −
05 1.6%bcsstk33 217 9 . e −
05 468 1 . e −
05 376 1 . e −
05 19.7%pcrystk02 185 1 . e −
03 560 1 . e −
04 500 2 . e −
04 10.7%garon2 1687 1 . e −
06 7194 3 . e −
07 5340 1 . e −
07 25.8%lhr11 125 2 . e −
05 320 3 . e −
06 296 1 . e −
06 7.5%lhr11c 125 2 . e −
05 320 5 . e −
06 296 3 . e −
06 7.5%barth5 151 5 . e −
05 320 2 . e −
06 296 1 . e −
06 7.5%nemeth01 781 5 . e −
06 7132 1 . e −
06 6156 2 . e −
06 13.7%nemeth17 261 6 . e −
09 1080 2 . e −
09 652 1 . e −
09 39.6%pkustk02 149 2 . e −
02 378 1 . e −
02 368 3 . e −
02 2.7%nd3k 1297 3 . e −
06 10668 1 . e −
07 8328 3 . e −
07 21.9%coupled 281 4 . e −
04 722 2 . e −
05 652 7 . e −
05 9.7%Pres Poisson 911 5 . e −
07 3218 1 . e −
07 2918 8 . e −
07 9.3%stokes64s 1991 1 . e −
04 32144 3 . e −
05 31260 3 . e −
05 2.8%tuma2 487 4 . e −
07 1204 3 . e −
08 1196 2 . e −
08 0.7%opt1 133 1 . e −
03 296 4 . e −
04 268 9 . e −
04 9.5%ramage02 137 9 . e −
04 628 3 . e −
04 582 1 . e −
04 7.3%Si5H12 197 3 . e −
06 420 5 . e −
07 416 2 . e −
07 1.0%net25 29 3 . e −
06 118 1 . e −
06 108 1 . e −
06 8.5%pf2177 245 1 . e −
04 724 3 . e −
05 620 2 . e −
05 14.4%c-41 297 5 . e −
03 1088 1 . e −
03 1000 4 . e −
03 8.1%c-44 217 1 . e −
03 788 3 . e −
04 690 2 . e −
04 12.4%PGPgiantcompo 89 4 . e −
05 230 5 . e −
06 200 5 . e −
06 13.1%wing nodal 247 9 . e −
05 584 8 . e −
06 548 5 . e −
06 6.2%Since GLTR is non-restarted, its memory requirement becomes harsh and evenunacceptable if too many iterations are needed unless the Lanczos basis vectors arestored in second memory, so that they can be called to produce the final solution. Asa result, considerable communication time is needed whenever GLTR produces thefinal solution by calling the previously saved basis vectors. In contrast, TRS IRA andTRS IRRA have low fast storage requirement for fairly small projection subspaces,and it is easy and simple to implement them.For B = tridiag(1 , , B = I , and the previous comments apply here. Besidesgaron2, nemeth01 and tuma2, the algorithms TRS IRA and TRS IRRA with ∆ = 1also use much fewer iterations than GLTR for wing nodal. In the meantime, a careful2 Z. JIA AND F. WANG calculation shows that the average improvements are 8.73% for ∆ = 1 and 8.90% for∆ = 100, respectively. A remarkable difference from Tables 2–3 is that all the Res’sin Table 5 are much larger than 10 − , which implies that all the ∆ / ∥ ˆ y ∥ B are quitelarge, or equivalently the ∥ ˆ y ∥ B are small. Table 4
Numerical results on sparse A = G + G T with G in Table 1 and B = tridiag(1 ; ; , ∆ = 1 . Matrix GLTR TRS IRA TRS IRRAMVs Res MVs Res MVs Res ratiobcsstk29 209 2 . e −
09 448 8 . e −
10 422 1 . e −
09 5.8%bcsstk33 185 7 . e −
09 468 5 . e −
09 392 1 . e −
09 16.2%pcrystk02 171 3 . e −
09 580 4 . e −
09 472 2 . e −
09 18.6%garon2 255 4 . e −
11 424 1 . e −
11 372 2 . e −
11 12.3%lhr11 147 1 . e −
09 392 5 . e −
10 344 1 . e −
09 7.5%lhr11c 147 1 . e −
08 422 9 . e −
09 368 4 . e −
09 12.8%barth5 455 2 . e −
08 1080 3 . e −
08 1020 1 . e −
08 5.6%nemeth01 199 1 . e −
08 360 5 . e −
09 320 7 . e −
09 11.1%nemeth17 129 1 . e −
11 250 4 . e −
12 224 3 . e −
11 10.4%pkustk02 187 2 . e −
07 378 1 . e −
07 368 1 . e −
07 2.7%nd3k 383 1 . e −
10 740 1 . e −
10 728 5 . e −
11 1.6%coupled 307 1 . e −
06 908 4 . e −
06 746 3 . e −
06 17.8%Pres Poisson 143 1 . e −
11 298 4 . e −
12 268 1 . e −
11 10.1%stokes64s 243 8 . e −
11 512 9 . e −
11 476 7 . e −
11 7.0%tuma2 299 4 . e −
08 630 3 . e −
08 580 3 . e −
08 7.9%opt1 141 1 . e −
08 308 1 . e −
09 268 2 . e −
09 13.0%ramage02 135 3 . e −
08 260 2 . e −
08 260 2 . e −
08 0.0%Si5H12 455 4 . e −
11 854 3 . e −
11 782 4 . e −
11 8.4%net25 63 3 . e −
08 116 7 . e −
08 108 6 . e −
08 6.9%pf2177 115 6 . e −
05 260 2 . e −
05 236 3 . e −
05 9.2%c-41 409 7 . e −
07 1932 5 . e −
07 1744 2 . e −
07 9.7%c-44 421 1 . e −
07 1360 1 . e −
07 1310 1 . e −
07 3.7%PGPgiantcompo 105 6 . e −
09 238 2 . e −
09 228 5 . e −
09 4.2%wing nodal 363 4 . e −
08 624 4 . e −
09 580 2 . e −
09 7.1%
7. Conclusion.
We have proposed and developed an eigenvalue-based formula-tion the refined Arnoldi method and IRRA algorithm to solve a large scale TRS’s.The refined Arnoldi method can overcome the potential non-convergence deficien-cy of the Arnoldi method when computing eigenvectors. In order to compare anyeigenvalue-based algorithm with the GLTR algorithm correctly and fairly, we haveestablished the relationship between ∥ ( A + ˆ λB )ˆ s + g ∥ B (cid:0) and ∥ ( M − ˆ λI )ˆ y ∥ B (cid:0) , where igenvalue-based algorithms for TRS Table 5
Numerical results on sparse A = G + G T with G in Table 1 and B = tridiag(1 ; ; , ∆ = 100 . Matrix GLTR TRS IRA TRS IRRAMvs Res Mvs Res Mvs Res ratiobcsstk29 175 2 . e −
05 448 2 . e −
05 422 1 . e −
05 5.8%bcsstk33 165 5 . e −
05 422 1 . e −
05 372 4 . e −
05 11.8%pcrystk02 139 8 . e −
05 632 1 . e −
05 472 1 . e −
05 25.3%garon2 2121 4 . e −
07 12384 1 . e −
07 11188 9 . e −
08 9.7%lhr11 133 1 . e −
05 370 5 . e −
06 360 2 . e −
06 2.7%lhr11c 133 1 . e −
05 422 7 . e −
05 380 3 . e −
05 10.0%barth5 383 2 . e −
04 1008 4 . e −
04 928 1 . e −
04 7.9%nemeth01 157 2 . e −
04 390 2 . e −
05 320 8 . e −
05 18.0%nemeth17 157 7 . e −
06 384 4 . e −
06 368 3 . e −
06 4.2%pkustk02 157 7 . e −
04 378 8 . e −
04 368 1 . e −
04 2.7%nd3k 2471 1 . e −
06 16544 9 . e −
07 16128 4 . e −
08 2.5%coupled 263 1 . e −
02 1230 1 . e −
02 828 1 . e −
02 32.7%Pres Poisson 1025 2 . e −
07 4318 4 . e −
08 3648 2 . e −
08 15.5%stokes64s 111 3 . e −
01 9316 4 . e −
01 8640 4 . e −
01 7.3%tuma2 223 4 . e −
04 630 4 . e −
04 584 4 . e −
05 7.3%opt1 117 1 . e −
04 261 1 . e −
04 260 1 . e −
04 0.4%ramage02 85 3 . e −
04 260 2 . e −
04 260 3 . e −
04 0.0%Si5H12 483 4 . e −
07 854 3 . e −
07 782 1 . e −
07 8.4%net25 51 8 . e −
04 116 9 . e −
04 108 7 . e −
04 6.9%pf2177 55 2 . e −
02 260 7 . e −
02 236 3 . e −
02 9.2%c-41 305 1 . e −
02 1744 7 . e −
03 1632 9 . e −
03 6.4%c-44 331 2 . e −
03 1412 2 . e −
03 1308 4 . e −
03 7.4%PGPgiantcompo 89 5 . e −
05 256 1 . e −
05 232 6 . e −
05 9.4%wing nodal 273 8 . e −
04 642 1 . e −
04 628 9 . e −
05 2.2%ˆ s is the approximation to s opt recovered from the converged eigenvector of the matrixpair ( M, B ). This issue is important but has been ignored in the literature whencomparing the efficiency and accuracy of an eigenvalue-based algorithm, such as IRA,and an iterative algorithm that solves TRS directly. Based on the obtained result, wecan first preset a stopping tolerance for eigenvalue-based algorithms, e.g., IRA andIRRA, and then determine a reliable one for the algorithms that solve TRS directly,e.g., GLTR, so that the two kinds of algorithms can compute converged approximatesolutions with very comparable accuracy. Therefore, we have provided a correct andfair comparison standard for the two kinds of algorithms. In the meantime, we haveanalyzed the residual norms obtained by GLTR and the Arnoldi method, the refinedArnoldi method for the equivalent eigenproblem of TRS, showing how they converge4
Z. JIA AND F. WANG as the subspace dimension increases.We have made extensive numerical experiments and confirmed that TRS IRAand TRS IRRA is comparable to GLTR. The results have also demonstrated thatTRS IRRA is more efficient than TRS IRA and improves the latter by roughly 10%in terms of matrix-vector products.
REFERENCES[1]
S. Adachi, S. Iwata, Y. Nakatsukasa, and A. Takeda , Solving the trust-region subproblemby a generalized eigenvalue problem , SIAM J. Optim., 27 (2017), pp. 269–291.[2]
A. R. Conn, N. I. M. Gould, and P. L. Toint , Trust-region Methods , SIAM, Philadelphia,2000.[3]
T. Davis, and Y. Hu , The University of Florida sparse matrix collection , ACM Trans. Math.Software, 38 (2011), pp. 1–25.[4]
J. B. Erway, P. E. Gill, and J. D. Griffin , Iterative methods for finding a trust-region step ,SIAM J. Optim., 20 (2009), pp. 1110–1131.[5]
C. Fortin, and H. Wolkowicz , The trust region subproblem and semidefinite programming ,Optim. Methods Softw., 19 (2004), pp. 41–67.[6]
W. Gander, C. H. Golub, and U. von Matt , A constrained eigenvalue problem , LinearAlgebra Appl., 114 (1989), pp. 815–839.[7]
G. H. Golub and C. F. van Loan , Matrix Computations , Fourth Edition, The John HopkinsUnversity, 2013.[8]
N. I. M. Gould, S. Lucidi, M. Roma, and P. L. Toint , Solving the trust-region subproblemusing the Lanczos method , SIAM J. Optim., 9 (1999), pp. 504–525.[9]
N. I. M. Gould, D. Orban, and P. L. Toint , GALAHAD, a library of thread-safe Fortran90 packages for large-scale nonlinear optimization , ACM Trans. Math. Softw., 29 (2004),pp. 353–372.[10]
N. I. M. Gould, D. P. Robinson, and H. S. Thorne , On solving trust-region and otherregularised subproblems in optimization , Math. Prog. Comput., 2 (2010), pp. 21–57.[11]
W. W. Hager , Minimizing a quadratic over a sphere , SIAM J. Optim., 12 (2001), pp. 188–208.[12]
W. W. Hager, and Y. Krylyuk , Graph partitioning and continuous quadratic programming ,SIAM J. Alg. Discrete Methods, 12( 1999), pp. 500–523.[13]
N. J. Higham , Functions of Matrices: Theory and Computation , SIAM, Philadelphia, PA,2008.[14]
Z. Jia , The convergence of generalized Lanczos methods for large unsymmetric eigenproblems ,SIAM J. Matrix Anal. Appl., 16 (1995), pp. 843–862.[15]
Z. Jia , Refined iterative algorithms based on Arnoldi’s process for large unsymmetric eigen-problems , Linear Algebra Appl., 259 (1997), pp. 1–23.[16]
Z. Jia , Polynomial characterizations of the approximate eigenvectors by the refined Arnoldimethod and an implicitly restarted refined Arnoldi algorithm , Linear Algebra Appl., 287(1999), pp. 191–214.[17]
Z. Jia , Some theoretical comparisons of refined Ritz vectors and Ritz vectors , Science in ChinaSer. A Math., 47 (2004), Supp. 222–233.[18]
Z. Jia, and G. W. Stewart , An analysis of the Rayleigh–Ritz method for approximatingeigenspaces , Math. Comput., 70 (2001), pp. 637–648.[19]
Z. Jia, and F. Wang , The convergence of the generalized Lanczos trust-region method for thetrust-region subproblem , SIAM J. Optim., accepted, 2021.[20]
J. J. Mor(cid:19)e, and D. C. Sorensen , Computing a trust region step , SIAM J. Sci. Statist. Com-put., 4 (1983), pp. 553–572.[21]
J. Nocedal and S. J. Wright , Numerical Optimization , Second Edition, Springer, 2006.[22]
M. Rojas, S. A. Santos, and D. C. Sorensen , A new matrix-free algorithm for the large-scaletrust-region subproblem , SIAM J. Optim., 11 (2001), pp. 611–646.[23]
M. Rojas, S. A. Santos, and D. C. Sorensen , Algorithm 873: LSTRS: MATLAB softwarefor large-scale trust-region subproblems and regularization , ACM Trans. Math. Softw., 34(2008), pp. 1–28.[24]
F. Rendl, and H. Wolkowicz , A semidefinite framework for trust region subproblems withapplications to large scale minimization , Math. Prog., 77 (1997), pp. 273–299.[25]
O. Schenk and K. Gartner , On fast factorization pivoting methods for symmetric indefinitesystems , Electron. Trans. Numer. Anal., 23 (2006), pp. 158–179.[26]
D. C. Sorensen , Implicit application of polynomial filters in a k-step Arnoldi Method , SIAMigenvalue-based algorithms for TRS J. Matrix Anal. Appl., 13 (1992), pp. 357–385.[27]
D. C. Sorensen , Minimization of a large-scale quadratic function subject to a spherical con-straint , SIAM J. Optim., 7 (1997), pp. 141–161.[28]
T. Steihaug , The conjugate gradient method and trust regions in large scale optimization ,SIAM J. Numer. Anal., 20 (1983), pp. 626–637.[29]
G. W. Stewart , Matrix Algorithms II: Eigensystems, SIAM, Philadelphia, PA, 2001.[30]
G. W. Stewart and J.-G. Sun , Matrix Perturbation Theory , Academic Press, INC., Boston,1990.[31]
P. L. Toint , Towards an efficient sparsity exploiting Newton method for minimization , SparseMatrices and Their Uses, Academic Press, 1981, pp. 57–88.[32]
Y. Yuan , On the truncated conjugate gradient method , Math. Prog., 87 (2000), pp. 561–573.[33]
L. H. Zhang, C. G. Shen, and R. C. Li , On the generalized Lanczos trust-region method ,SIAM J. Optim., 27 (2017), pp. 2110–2142.[34]