Issues on the use of a modified Bunch and Kaufman decomposition for large scale Newton’s equation

In this work, we deal with Truncated Newton methods for solving large scale (possibly nonconvex) unconstrained optimization problems. In particular, we consider the use of a modified Bunch and Kaufman factorization for solving the Newton equation, at each (outer) iteration of the method. The Bunch and Kaufman factorization of a tridiagonal matrix is an effective and stable matrix decomposition, which is well exploited in the widely adopted SYMMBK (Bunch and Kaufman in Math Comput 31:163–179, 1977; Chandra in Conjugate gradient methods for partial differential equations, vol 129, 1978; Conn et al. in Trust-region methods. MPS-SIAM series on optimization, Society for Industrial Mathematics, Philadelphia, 2000; HSL, A collection of Fortran codes for large scale scientific computation, http://www.hsl.rl.ac.uk/; Marcia in Appl Numer Math 58:449–458, 2008) routine. It can be used to provide conjugate directions, both in the case of 1×1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1\times 1$$\end{document} and 2×2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\times 2$$\end{document} pivoting steps. The main drawback is that the resulting solution of Newton’s equation might not be gradient–related, in the case the objective function is nonconvex. Here we first focus on some theoretical properties, in order to ensure that at each iteration of the Truncated Newton method, the search direction obtained by using an adapted Bunch and Kaufman factorization is gradient–related. This allows to perform a standard Armijo-type linesearch procedure, using a bounded descent direction. Furthermore, the results of an extended numerical experience using large scale CUTEst problems is reported, showing the reliability and the efficiency of the proposed approach, both on convex and nonconvex problems.


Introduction
In this paper we consider the unconstrained optimization problem min f (x), (1.1) where f : R n −→ R is a possibly nonconvex real valued function, n is large and we solve (1.1) by means of an iterative method which generates the sequence {x h }. We consider the standard assumptions that f is twice continuously differentiable, and that {x h } ⊂ , where ⊂ R n is compact. Furthermore, no sparsity pattern is assumed for the Hessian matrix ∇ 2 f (x). Then, we approach the solution of (1.1) by means of the iterative process being p h ∈ IR n a search direction and α h > 0 a steplength. Dealing with large scale problems, Truncated Newton methods are often considered a good choice for their solution. Truncated Newton methods are also called Newton-Krylov methods, since a Krylov subspace method is usually used, for iteratively (approximately) solving the Newton equation at each iteration. In [21] a general description of a Truncated Newton method is reported. As well known, given an initial guess x h , a Truncated Newton method for problem (1.1) is based on two nested loops: • The outer iterations, where the current estimate of the solution x h+1 is updated using x h ; • The inner iterations, where an iterative algorithm is used for computing, at each outer iteration h, the approximate solution p h of the Newton equation The Conjugate Gradient (CG) method is among the most commonly used iterative algorithms (see e.g. [1,17]) for solving (1.3). The main drawback of CG method is that it is well-posed only in case ∇ 2 f (x h ) is positive definite. Conversely, if matrix ∇ 2 f (x h ) is indefinite, the CG algorithm could break down and some difficulties may arise to define an approximate solution p h of (1.3). In order to cope also with indefinite linear systems, some modified CG algorithms have been proposed in [8,9,11,12,14]. In particular, if the CG breaks down when solving (1.3), we can take as Newton-type direction p h in the recursion (1.2): • A scaled steepest descent direction (i.e., p h = −C∇ f (x h ), with C positive definite); • A suitable gradient-related adaptation of the current best approximate solution of (1.3); • A suitable (descent) negative curvature direction. Moreover, in [18] a non-standard approach is proposed in order to build a gradientrelated direction, by separately accumulating positive and negative curvature directions computed by the CG.
Considering a different standpoint, the CG algorithm may be also used in the indefinite case, as proposed in [13,18], where a suitable combination of conjugate directions proves to yield a gradient-related direction.
In this work we are interested to compare the quality of the search direction obtained by both CG and SYMMBK algorithms (the latter being proposed by Chandra in [5]), when solving (1.3), in order to provide eventually the gradient-related direction p h in (1.2). On this guideline, we first need to prove that both these solvers yield a gradient-related direction to the optimization framework. I.e., the direction p h in (1.2) is eventually bounded and of sufficient descent, namely for any h ≥h ≥ 1 the next definition holds. Definition 1.1 Given the direction p h in (1.2), which approximately solves (1.3), let h, σ 1 , σ 2 and M be finite positive values independent of h, withh ≥ 1. Then, we say that p h is gradient-related at x h if for any h ≥h the next couple of conditions hold: (1.5) For this purpose, as regards the CG, the reader may refer to [11,12,18]. On the other hand, we also need to analyze the Lanczos process, adopted inside SYMMBK algorithm. As well known, one of the advantages of using the Lanczos process in place of CG algorithm is that it does not break down in the indefinite case. We recall that in the case the matrix ∇ 2 f (x h ) is positive definite, from a theoretical point of view, using the Lanczos process or CG algorithm is to a large extent equivalent [23,24].
The paper is organized as follows. In Sect. 2 some preliminaries on Truncated Newton method and Lanczos process are reported; then, the Bunch and Kaufman decomposition, as well as some basics on SYMMBK (see also [6]), are detailed. In Sect. 3 we show how to compute a gradient-related direction by using the Bunch and Kaufman decomposition. Section 4 details numerical experiments. Finally, Sect. 4 reports some concluding remarks.
We indicate by · the Eclidean norm of real vectors and matrices. Moreover, λ (C) represents the th eigenvalue of the real symmetric matrix C. Finally, e i is the ith real unit vector.

Preliminaries
In this section we report some basics we will use in the sequel. In particular, first we recall a general scheme of a linesearch-based Truncated Newton algorithm. Afterwards, we report the well known Lanczos process. Finally, we detail the Bunch and Kaufman factorization in order to decompose the tridiagonal matrix provided by the Lanczos process.

Truncated Newton method
A practical linesearch-based Truncated Newton algorithm can be described as follows (see e.g. [21]).

INNER ITERATIONS
Compute p h which approximately solves equation (1.3) and satisfies the truncation rule Compute the steplength α h by using a linesearch procedure.

End Do
In large scale settings, iterative methods like the CG algorithm can be used in the inner iterations. In the case the Hessian matrix in (1.3) is indefinite, then the Lanczos algorithm may successfully be used as an alternative to the CG method, since it does not suffer for a possible pivot breakdown. We remark that within the Truncated Newton methods, the importance of an efficient truncation criterion for the inner iterations was pointed out in [8,9,22]. More recently, an adaptive truncation rule, in order to enhance the residual-based criterion, was proposed in [3,4].

Assumption 2.1 Given the function
For completeness we close this section by recalling a very general result of global convergence for Truncated Newton methods (see any textbook on Continuous Optimization).

Basics of the Lanczos process
As well known, the Lanczos process (see [6]) is an efficient tool to transform the symmetric indefinite linear system (1.3) into a symmetric tridiagonal one. On this guideline, consider the Newton equation (1.3) and rewrite it as follows (here we drop for simplicity the dependency on the subscript h): A practical implementation of the Lanczos process applied to (2.1), which generates the orthogonal Lanczos vectors q i , i ≥ 1, is reported in the scheme below.
Defining Q k ∈ IR n×k as the matrix whose columns are the Lanczos vectors, that is and the symmetric tridiagonal matrix T k ∈ IR k×k the Lanczos process generates T k from (2.1) in k ≤ n steps. On the other hand, after k ≤ n iterations the next additional relations also hold for the Lanczos process: In particular, if γ k+1 = 0 in (2.4), then from (2.4)-(2.8) we have (see also [5]) Relations (2.9) also reveal that, since T k is nonsingular, as long as γ k+1 = 0 in (2.4) then computing the vector y k from T k y k = b e 1 may easily yield d k = Q k y k , being d k an approximate solution of the linear system in (2.1).

Bunch and Kaufman decomposition
As reported in Sect. 2.1, unlike the CG method, the Lanczos process does not break down if matrix A in (2.1) is indefinite nonsingular. Furthermore, by (2.5)-(2.6) the nonsingularity of A also yields the nonsingularity of matrix T k . On the other hand, observe that a suitable decomposition of the tridiagonal matrix T k in (2.3) is mandatory, in order to easily compute y k by solving the tridiagonal system T k z = b e 1 . However, in the case the matrix A is indefinite, the Cholesky factorization of T k may not exist. In this regard, the Bunch and Kaufman decomposition of T k may be adopted, as in SYMMBK algorithm. We detail now for T k the Bunch and Kaufman decomposition (see [2]) where B k is a block diagonal matrix with 1×1 or 2×2 diagonal blocks, and S k is a unit lower triangular matrix, such that its non-zeroes are restricted to three main diagonals. The diagonal blocks of B k represent pivoting elements in the factorization (2.10). In particular, in order to control the growth factor when computing the decomposition (2.10), the scalar is selected (see [5]), and the choice of 1 × 1 pivot or 2 × 2 pivot for B k is based on the following rule (which only refers, without loss of generality, to the first step of the decomposition (2.10)): • If |δ 1 | > η(γ 2 ) 2 , then δ 1 in (2.3) is used as a 1×1 pivot to generate the factorization (2.10) of T k . Namely, we have T k = S k B k S T k with After a brief computation we have is used as a 2 × 2 pivot, to generate the factorization (2.10), so that Following a simple recursion, the tridiagonal matricesT k−1 in (2.12) andT k−2 in (2.15) can be similarly decomposed, so that the entries of S k in (2.10) can be fully computed after m pivoting steps of the Bunch and Kaufman decomposition, with m ≤ k (see [5] and relations (2.19)-(2.20) or (2.25)-(2.26)).

Basics of SYMMBK
Here we report some basics of the SYMMBK algorithm (first proposed by Chandra in [5]), that uses the Bunch and Kaufman factorization in (2.10) to compute d k in (2.9). Technicalities in this section are needed to provide formulae (2.22)-(2.24) and (2.28)-(2.30), which are essential to prove properties for our proposal in Sect. 3. To this aim, recalling (2.10) we introduce the n × k matrix By (2.10), (2.16) and (2.17), the equations in (2.9) can be rewritten as follows: (2.18) Then, to calculate the vector d k of (2.9), we equivalently need to compute both W k and ζ (k) . In this regard, we must provide in (2.18) the matrices S k and B k , distinguishing between the cases where either a 1 × 1 or a 2 × 2 pivoting step was performed for their computation. This can be easily done by considering that the first step is performed as in (2.11) (1 × 1 pivot) or (2.14) (2 × 2 pivot). Then, for any successive step we use information from the previous pivot, as reported hereafter (we assume that applying the rule ( j)-( j j) in Sect. 3 the Bunch and Kaufman decomposition performed m steps, with m ≤ k). To preserve the notation used by Chandra [5], here we simply focus on the last pivot; then, we infer results for the intermediate pivots, too.
If the last pivot performed is 1 × 1 (i.e. according with the rule ( j) we have δ k large enough), then from Sect. 2.3 the first relation 4 are computed from the second last step. In particular, if the second last pivot was 1 × 1, then After some computations, replacing the expressions of B k and S k in (2.18) we obtain the following formulae With an obvious generalization the formulae (2.22)-(2.24) hold replacing the step k with the step i ≤ k, whenever a 1 × 1 pivot is performed by the Bunch and Kaufman decomposition.
Similarly, if the last pivot performed is now 2 × 2 (i.e. according with the rule ( j j) we have δ k−1 relatively small and 4 are computed from the second last step. In particular, if the second last pivot was 1 × 1, then where now After some computations, replacing again the expressions of B k and S k in (2.18), we obtain the following formulae As for (2.22)-(2.24), the formulae (2.28)-(2.30) hold replacing the step k with the step i ≤ k, as long as at step i a 2 × 2 pivot is used. In the end, we conclude that SYMMBK algorithm generates the vectors z i , i = 1, . . . , k, which represent a sequence of approximations of the vector d k satisfying (2.18). Moreover, the direction d k = z k satisfies (2.18) (or equivalently (2.9)).

Our proposal
As detailed in Sect. 2.3, the Bunch and Kaufman decomposition used within the SYMMBK algorithm is able to cope with Newton's equation (2.1), also in the indefinite case. However, the search direction provided by (2.9) might not be gradient-related, at the hth outer iteration. In this section we aim at suitably adapting SYMMBK algorithm, in order to guarantee that, using the computed approximate solution of (1.3), we can yield a gradient-related direction. For this purpose, first of all we need to slightly modify the partial pivoting rule used for choosing 1 × 1 or 2× 2 pivot in the Bunch and Kaufman decomposition. In particular, given ω ∈ (0, 1) the choice for 1 × 1 pivot or 2 × 2 pivot in Sect. 2.3, in order to obtain (2.11) or (2.14), is replaced by the following rule (again for simplicity we refer only to the first step of the decomposition): ( j) if |δ 1 | > ωη(γ 2 ) 2 , then δ 1 is used as a 1 × 1 pivot; ( j j) if |δ 1 | ≤ ωη(γ 2 ) 2 , then the 2 × 2 block (2.13) is used as a 2 × 2 pivot. Now we show that the novel condition |δ 1 | > ωη(γ 2 ) 2 in ( j) (which yields (2.11)), along with the novel condition |δ 1 | ≤ ωη(γ 2 ) 2 in ( j j) (which yields (2.14)), allow to use an adapted SYMMBK algorithm for the construction of a gradient-related direction. To this aim, we recall that repeatedly performing the factorizations (2.11) and (2.14) we finally obtain for T k the decomposition (2.10), where (m ≤ k) Of course, since the matrix A in (2.1) is possibly indefinite, the direction d k in (2.9) might be of ascent and/or not gradient-related. However, we show how to possibly modify d k in order to obtain the gradient-related direction p h , to be used in (1.2).
Thus, in both the above subcases S 3,1 and S 3,2 in (3.6) are evidently bounded by a constant value, with (3.9) Finally, iterating the above procedure for any step, by (3.2), (3.6), (3.8), (3.9) and using the test ( j)-( j j) in the Bunch and Kaufman decomposition to decide pivoting 1 × 1 or 2 × 2, we obtain that (a) By (2.11)-(2.12) and (2.14)-(2.15), any pivoting element B j k in (3.1) has a determinant which is bounded away from zero; (b) All the entries of the matrix S k in (2.11) and (2.14) are bounded by constant values independent of k.
This implies (see Sect. 2.4) that any vector in the sequences ζ 1 w 1 , . . . , ζ k w k and z 1 , . . . , z k is bounded. Moreover, the vector d k = z k fulfilling (2.9) is unique and bounded.
Observe that by the previous result, the condition (i) in Sect. 3 is partly guaranteed to hold. However, using SYMMBK algorithm, the direction d k computed in Proposition 3.1 might not be a direction of sufficient descent. To guarantee also the fulfillment of (ii), we propose to suitably use the information collected by the sequence {ζ i w i }, in order to compute p h endowed with descent property.
From Proposition 3.1 and Sect. 2.4, the real vector is bounded. Thus, we propose to compute the search direction p h , at the outer iteration h of the Truncated Newton method, in accordance with the next scheme:

Reverse-Scheme
Data: Set the initial vector p h = z 0 = 0, along with the parameter φ > 0; If at step i of the Bunch and Kaufman decomposition in SYMMBK a 1 × 1 pivot is performed, when solving the linear system T k z = b e 1 , then If at step i of the Bunch and Kaufman decomposition in SYMMBK a 2 × 2 pivot is performed, when solving the linear system T k z = b e 1 , then set Set p h = p h + uu + vv.

End Do
We remark that in case the matrix A in (2.1) is positive definite, and the parameter ω in ( j)-( j j) is sufficiently small (i.e. no 2 × 2 pivot is performed), then the vector p h computed by the Reverse-Scheme coincides with Newton's direction. Now, in the next proposition we prove that the direction p h , obtained by Reverse-Scheme at any iterate of the sequence {x h } of the Truncated Newton method, is gradient-related. This result is not an immediate consequence of Proposition 3.1, since the latter proposition merely guarantees the simple boundedness of the direction d k , for a given k ≥ 1. Conversely, the next result proves that, under mild additional assumptions, all the vectors in the sequence { p h } are of sufficient descent and eventually are uniformly bounded by a positive finite constant. Proof The statement trivially holds if in Proposition 3.1 the Lanczos process performs just one iteration, being γ 2 < ε. In all the other cases (i.e. the Lanczos process has performed at least 2 iterations), we first prove that p h is of sufficient descent for f (x) at x h . In this regard, it suffices to prove the sufficient descent property for p h after the first step of the Bunch and Kaufman decomposition in SYMMBK, since the Reverse-Scheme accumulates only non-ascent directions in p h . To this aim, we distinguish between two cases, which may occur at iterate x h , when applying the Bunch and Kaufman decomposition to matrix T k in (2.9): • The first step corresponds to a 1 × 1 pivot. Then, by (2.16) and (2.22)-(2.24), we have for the vector uu in Reverse-Scheme (3.10) so that More explicitly, by (2.17) and relations (2.9)-(2.11) we have where the asterisks represent irrelevant entries. Then, ζ 1 = b /δ 1 . Moreover, by Assumption 2.1 the set is compact, so that (3.12) • The first step corresponds to a 2×2 pivot. Then, recalling (2.28) we have w 2 = q 2 , so that by (3.10) for the vectors uu and vv in Reverse-Scheme the next relations hold (3.13) Now, again by (2.17) and relations (2.9)-(2.10), (2.14) we get where again the asterisks represent irrelevant entries. Thus, it results Finally, similarly to (3.11), (3.13) and Reverse-Scheme the relation (3.14) Relations (3.12) and (3.14) prove that for any h ≥ 1 the condition (1.5) is eventually fulfilled, with As regards the boundedness of d k in (2.9) and p h , for any h ≥ 1, observe that by (2.9) we first have . Now, let m 1 be the number of steps where the test ( j) is fulfilled (i.e. m 1 is the number of 1 × 1 pivot steps), so that (k − m 1 )/2 is the number of 2 × 2 pivot steps. Also observe that the inverse S −1 k of S k in (2.11) is given by and the inverse S −1 k of S k in (2.14) is given by Thus, the bound (3.4) on S 2,1 , along with the bounds (3.8) and (3.9) on S 3,1 and S 3,2 yield Moreover, as regards the diagonal blocks of B −1 k (see (2.11) and (2.14)), for any i ≥ 1 we have by the proof of Proposition 3.1 and the compactness of the following bounds: where β is given in (3.15). Similarly, by (2.18) we have W k = Q k S −T k , so that any column w i of W k satisfies (see (3.15)) This implies, by (3.16) and the boundedness of vectors ζ 1 w 1 , . . . , ζ k w k (see Proposition 3.1), that also uu and vv in Reverse-Scheme are bounded by a positive constant, for any x h ∈ . Thus, p h is bounded by a positive constant independent of h.

Numerical experiments
In order to assess the effectiveness of the approach proposed in this paper, we performed an extensive numerical testing. We used the HSL_MI02 Fortran routine available from the HSL Mathematical Software Library [19] which implements the SYMMBK algorithm, where we modified the pivoting rule as described in Sect. 3, namely ( j)-( j j) on page 9. Then, we embedded this routine within a standard implementation of a linesearch based Truncated Newton method, namely the one proposed in [15] (we refer to this paper for any implementation details), without considering any preconditioning strategy. Moreover, we replaced the CG used as iterative solver by the Truncated Newton method in [15] with the modified SYMMBK algorithm. It is worthwhile to note that in [15], in order to handle the indefinite case, the inner CG iterations are not terminated whenever a negative curvature direction is encountered. Indeed, by drawing inspiration from [18], positive and negative curvature CG directions {d i } are separately accumulated, and then suitably combined to obtain the search direction. We report the following scheme of the CG method adopted in [15], in order to briefly describe the latter strategy.

End Do
In particular, at the kth CG inner iteration, the following vectors are defined where for a given ∈ (0, 1) Then, the vector p P k + p N k is used as search direction. This choice has several advantages for which we refer to [15]. In particular, it can be proved that p P k + p P k is a gradient-related direction. However, as shown in [7], the use of p P k and p N k obtained by separating the contribution of negative and positive curvature directions actually corresponds to considering a modified positive definite linear system, in place of the original indefinite one. For this reason, in our numerical experiments, in order to avoid unfair conclusions, we prefer not to use the modified search direction p P k + p N k and adopt a more standard approach in the literature. Therefore, in our implementation we use the approximate solution of the original linear system as search direction. An additional safeguard is needed whenever the Hessian matrix is nearly singular at the hth iterate of the Truncated Newton method. In this case, the steepest descent direction is considered as "alternative" direction (both when the CG or the modified SYMMBK algorithms are adopted). As regards the termination criterion of the inner iteration, the standard residual-based criterion is used (see e.g. [3,4,21]), which we make slightly tighter (with respect to that in [15]) by multiplying the tolerance of the criterion by 10 −2 . In this manner, a more accurate solution of the Newton system at each outer iteration of the methods is required, and possibly a more significant comparison between the use of the modified SYMMBK in place of the standard CG can be obtained. Finally, the parameter ξ in the proof of Proposition 3.1 is set as ξ = max{(1 − η|δ 2 |), 10 −1 }, so that we consequently chose ω = min{1, In addition, the parameter φ in Reverse-Scheme was set as φ = 10 −10 : in this regard we strongly remark that in our numerical experience we always observed As regards the test problems, we selected all the large scale unconstrained test problems in the CUTEst collection [16], including both convex and nonconvex functions. In particular, we considered 112 problems whose dimension ranges from 1000 to 10,000. We report the results obtained in terms of number inner iterations, number of function evaluations and CPU time. The comparison between the use of the CG algorithm and the modified SYMMBK, when performing the inner iterations, is reported by using performance profiles [10]. Here the abscissa axis details values of the positive parameter τ ≥ 1. In particular, for a given valueτ of the parameter τ , each profile represents the percentage of problems solved by the corresponding method, using at mostτ times the resource (inner iterations/function evaluations/CPU time) used by the best method. Figure 1 reports the comparison in terms of inner iterations. By observing these plots, the performance of the two approaches seems to be similar and the use of the modified SYMMBK is slightly preferable in terms of robustness. The profiles in terms of function evaluations are reported in Fig. 2. In this case, the superiority of using the modified SYMMBK algorithm is more evident. Finally, as regards the comparison in terms of CPU, Fig. 3 indicates that the use of the standard CG leads to a slight improvement in terms of efficiency, but it is definitely less performing in terms of robustness. We recall indeed that though the CG can easily provide a gradient-related search direction in convex problems, it might prematurely stop in nonconvex problems. Now, in order to better assess the quality of the search direction obtained by using the modified SYMMBK algorithm, we performed further numerical experiments. In particular we computed the normalized directional derivative of the search direction  .2), namely ∇ f (x h ) T p h / ∇ f (x h ) p h , at each outer iteration of the CG and the modified SYMMBK method respectively, for some test problems. As an example, we report in Fig. 4 the normalized directional derivative obtained as the number of outer iterations increases, for the test problems DIXMAANJ-1500 (upper-left), DIXMAANK-3000 (upper-right), DIXMAANL-1500 (mid-left), DIXMAANL-3000 (mid-right), DIXMAANJ-3000 (lower-left), SPARSINE-1000 (lower-right), where the number indicates the dimension of the test problem.
On these problems the method based on the modified SYMMBK algorithm outperforms the one based on the standard CG (a much smaller number of iterations is required to achieve convergence). Plots reveal that this is due to the fact that the quality of the Newton-type direction, obtained by means of the modified SYMMBK algorithm, is definitely better. Indeed, on one hand the use of the modified SYMMBK algorithm for computing the search direction enables to discard the steepest descent direction more frequently than the use of the standard CG. On the other hand, the introduction of the parameter ω in the modified SYMMBK algorithm (see ( j)-( j j)) allows to have values of the directional derivative which are more likely bounded away from zero. Similar behaviour is evidenced when solving many other test problems of the CUTEst collection. Note also that in the plots of Fig. 4, both the compared methods tend to initially select the steepest descent method as a search direction. This should not be surprising, inasmuch as in the early outer iterations the truncation criterion of the Truncated Newton method is less tight. Finally, we are persuaded that the modified test ( j)-( j j) in SYMMBK might be worth investigating, also to construct a suitable negative curvature direction (see [14]), in nonconvex problems. This issue plays a keynote role in case the adopted Truncated Newton method is asked to provide also convergence to stationary points, where second order optimality conditions hold. Furthermore, there is the chance that adaptively updating the parameter ω in ( j)-( j j) may further enhance the performance of our method, preserving the theory developed in Sect. 3.