Iterative Grossone-Based Computation of Negative Curvature Directions in Large-Scale Optimization

We consider an iterative computation of negative curvature directions, in large-scale unconstrained optimization frameworks, needed for ensuring the convergence toward stationary points which satisfy second-order necessary optimality conditions. We show that to the latter purpose, we can fruitfully couple the conjugate gradient (CG) method with a recently introduced approach involving the use of the numeral called Grossone . In particular, recalling that in principle the CG method is well posed only when solving positivedeﬁnitelinearsystems,ourproposalexploitstheuseofgrossonetoenhancethe performanceoftheCG,allowingthecomputationofnegativecurvaturedirectionsin theindeﬁnitecase,too.Ouroverallmethodcouldbeusedtosigniﬁcantlygeneralizethe theoryinstate-of-the-artliterature.Moreover,itstraightforwardlyallowsthesolution ofNewton’sequationinoptimizationframeworks,eveninnonconvexproblems.We remarkthatouriterativeproceduretocomputeanegativecurvaturedirectiondoesnot requirethestorageofanymatrix,simplyneedingtostoreacoupleofvectors.This deﬁnitelyrepresentsanadvancewithrespecttocurrentresultsintheliterature.


Introduction
We consider the solution of the nonconvex unconstrained optimization problem min x∈R n f (x), where f : R n → R is a nonlinear smooth function and n is large.
Despite the use of the term 'minimization' in the last problem, most of the methods proposed in the literature (for its solution) generate a sequence of points {x k }, which is only guaranteed to converge to stationary points.Thus, specific methods need to be applied in case stationary points for the above problem, satisfying also second-order necessary optimality conditions, are sought (see, for instance, the seminal papers [1][2][3][4][5][6][7] in the framework of truncated Newton methods).Observe that additional care when using the latter methods is definitely mandatory, since imposing standard first-order stationarity conditions may not in general ensure convexity of the quadratic model of the objective function, in a neighborhood of the solution points.In this regard, the computation of so-called negative curvature directions for the objective function is an essential tool (see also the recent papers [4,8]), to guarantee convergence to stationary points which satisfy second-order necessary conditions.
Here, we want to provide a framework for the computation of negative curvature directions of quadratic functions, to be used within globally convergent iterative methods for large-scale nonlinear programming.Observe that the asymptotic convergence of iterative methods, toward second-order stationary points, implies that the Hessian matrix at limit points must be positive semidefinite.This fact requires that in principle the iterative methods adopted must be able to fully explore the eigenspaces of the Hessian matrix, at the current iterate, at least in a neighborhood of the stationary points.Equivalently, the optimization method adopted will have to efficiently cope also with nonconvexities of the objective function.The latter fact issues specific concerns in case n is large, since the computational effort to solve a nonlinear programming problem can be strongly affected by the scale.
In particular, as shown in [3,5,9], exploiting the nonconvexities of f (x) can be accomplished by suitable Newton-Krylov methods (in the context of Hessian-free truncated Newton methods), such that at each outer iteration j, a pair of search directions (s j , d j ) is computed, satisfying specific properties.Namely, the vector s j must be a direction which approximately solves Newton's equation ∇ 2 f (x j ) s = −∇ f (x j ) at x j .Its purpose is essentially that of ensuring the efficient convergence of the sequence {x j } to stationary points.On the other hand, the nonascent direction d j is a negative curvature direction (if any) for the objective function at x j .That is, d j is a nonascent direction such that d T j ∇ 2 f (x j )d j ≤ 0, satisfying suited conditions in order to force convergence to those stationary points where second-order necessary optimality conditions hold.In particular, d j should resemble an eigenvector corresponding to the least negative eigenvalue of the Hessian matrix ∇ 2 f (x j ).In [3,5] the direction d j is obtained as by-product of the Krylov-subspace method applied for solving Newton's equation, though an expensive storage is required in [5] and a heavy computational burden is necessary in the approach proposed in [3].
To overcome these drawbacks, in [9] an important novelty is introduced, namely the iterative computation of the sequence {d j }, not requiring neither expensive computations nor excessive storage.This novel approach is based on the use of the so-called Planar-CG method which is a modification of the conjugate gradient method in [11].In particular, in [9], at any iterate x j ∈ R n , the vector d j is computed through the linear combination of a few n-real vectors, generated by the used Krylov subspace method.It was proved that on nonconvex problems, the overall exact computation of d j simply requires at iterate x j the storage of at most four n-real vectors.Even if this approach revealed effective, it presents some drawbacks related to a cumbersome analysis.
In this paper, partially starting from the drawbacks of the proposal in [9], we aim at describing a strong simplification in the computation of the directions {d j }, by using a novel approach which extends some ideas in [12].Namely, we adopt the Krylovbased method CG ① defined in [12], being ① the symbol for the grossone (see [13]), in order to generate a suitable matrix factorization which allows the computation of {d j }.Similarly to [12], we first show that the CG is an ideal candidate to generate the latter matrix factorization.However, it may reveal serious disadvantages on nonconvex problems.In this regard, CG ① represents a natural generalization of the CG, and with some care allows to extend CG properties on indefinite problems.Then, the CG ① is used to generate directions in eigenspaces of the Hessian matrix associated with negative eigenvalues and to provide a suitable matrix factorization, which allows to exploit the results in [14].
We also propose a numerical experience, where we assess the effectiveness of the negative curvature directions computed in the current paper.We prefer to skip a numerical comparison between our proposal and those in [3,5,10], the latter requiring an expensive matrix storage or the need for recomputing some quantities/vectors.This risks to make the comparison unfair, inasmuch as in [9] and here we prove that an inexpensive iterative computation of d j is obtained, by storing at most two (four in [9]) working vectors.
To sum up, considering [9] as a reference paper with respect to our analysis, the main enhancements of the approach in the current paper can be summarized as follows: -in [9] the computation of the negative curvature direction d j , at iterate x j , requires the storage of up to four vectors, while here we propose a method requiring the storage of only two vectors; -the theory in [9] heavily relies on complicate matrix factorizations, due to the structure of the Planar-CG method therein adopted.Conversely, here the analysis through grossone only indirectly uses matrix factorizations provided by a Planar-CG method.Moreover, the Planar-CG method indirectly adopted here is definitely simpler that the one adopted in [9].Hence, here the theoretical analysis to prove convergence of the sequence {x j } to stationary points, satisfying second-order conditions, is drastically simplified; -the strategy adopted in [9] to compute the search directions is definitely more computationally expensive than the one proposed here; -as regards numerical results, we do not claim that our proposal is always more efficient than the one in [9], depending on the problem in hand.
The paper is organized as follows: Sect. 2 reports some preliminaries on the use of negative curvature directions within truncated Newton methods.In Sect.3, we give some basics on grossone, in order to motivate its use within the algorithm CG ① .Section 4 emphasizes the importance of certain matrix factorizations, in order to iteratively compute the final negative curvature direction.Section 5 stresses the importance of pairing CG ① with the approach in [14].Moreover, the last section explicitly yields our formula for determining the negative curvature direction.Finally, Sect.6 contains a numerical experience on our proposal, and Sect.7 reports some conclusions and future research perspectives.
In this paper, we use standard notations for vectors and matrices.With • , we indicate the Euclidean norm.λ[A] is a general eigenvalue of matrix A ∈ R n×n , and A 0 [A 0] indicates that A is positive definite [semidefinite].e k ∈ R n represents the kth unit vector, while the symbol ① represents the numeral grossone (see also [13]).

Negative Curvature Directions in Truncated Newton Methods
Hereafter, we will use the following scheme with f ∈ C 2 (R n ), as a general reference for an unconstrained optimization problem.Moreover, the equation represents Newton's equation associated with problem (1).The use of negative curvature directions in the framework of truncated Newton methods was introduced in the early papers [6,7], in order to define algorithms converging toward second-order critical points, namely stationary points where the Hessian matrix is positive semidefinite.Following the approach in [7], the sequence of negative curvature directions {d j } is expected to satisfy the conditions in the next assumption.Assumption 2.1 Given problem (1), with f ∈ C 2 (R n ), the nonascent directions in the sequence {d j } are bounded and satisfy the conditions where λ min ∇ 2 f (x j ) is the smallest eigenvalue of the Hessian matrix ∇ 2 f (x j ).
The approach adopted in [7] may be generalized to some extent (see, for instance, the proposal in [5]), by suitably weakening conditions the directions {d j } are subject to.Roughly speaking, the condition (a) in Assumption 2.1 implies that at any iterate x j , the nonascent vector d j must be a nonpositive curvature direction.Moreover, as in condition (b), when the quantity d T j ∇ 2 f (x j )d j approaches zero, then the sequence {x j } is approaching a region of convexity for the function f (x).Indeed, in such a case there will be no more chances to compute a negative curvature direction satisfying d T j ∇ 2 f (x j )d j < 0, so that eventually the condition d T j ∇ 2 f (x j )d j → 0 must hold.Of course, on convex problems, negative curvature directions are not present, so that points provided by Newton-Krylov methods eventually satisfy also second-order stationarity conditions.
We recall that the main purpose of Newton-Krylov methods is to compute the (possibly) infinite sequence {x j }, such that at least one of its subsequences is convergent to a stationary point of f (x).In this regard, Assumption 2.1 does not imply a unique choice of d j at iteration j.In fact, in order to fulfill (b) of Assumption 2.1, it suffices to compute d j so that d T j ∇ 2 f (x j )d j ≤ v T j ∇ 2 f (x j )v j , being v j an eigenvector associated with the smallest eigenvalue of ∇ 2 f (x j ).In addition, d j becomes essential only eventually; i.e., far from a stationary point, it might be unnecessary to force convergence toward regions of convexity for f (x).Nevertheless, using information associated with the negative curvature direction d j , also when far from a solution point, may considerably enhance efficiency.The latter fact was evidenced, for instance, in [3,5,10] and intuitively follows from the next reasoning.Given the local quadratic expansion at x j , the directional derivative of q j (d) along d, may strongly decrease when d is not only a descent vector, but also a negative curvature direction.This fact is explicitly used in [2,3,5,10] when a negative curvature direction d j is adopted, at any iteration j, even if it is possibly not associated with an eigenvector z corresponding to the smallest negative eigenvalue of ∇ 2 f (x j ).Further developments on the computation of negative directions, to be used within truncated Newton methods, have been introduced in the already mentioned papers [3,5,9].In particular, in [3,5] the direction d j is obtained as by-product, when applying a Krylov-subspace method to solve Newton's equation (2).However, since the Krylov method may perform, at iteration j, a number k of steps considerably smaller than n, not all the eigenspaces associated with the Hessian matrix ∇ 2 f (x j ) will be explored, vanishing the search of z.Hence, only an approximation of z may be available after k steps of the Krylov-based method.
In [9], the iterative computation of the directions {d j } is proposed.The innovative contribution in [9] consisted of explicitly providing the iterative computation of the sequence {d j }, without requiring burdensome re-computing (as in [3]) or any expensive storage (as in [5]).In particular, the Krylov-based procedure adopted in [9] to compute d j involves the use of the Planar-CG method, which represents an extension of the CG to nonconvex quadratic functions, where the Hessian matrix is possibly indefinite.The approach using this Planar-CG method surely proved to be effective, but it has a major disadvantage, requiring a fairly complex analysis which involves considering different and articulated subcases.element allowing one to express finite quantities.)From the foundational point of view, grossone has been introduced as an infinite unit of measure equal to the number of elements of the set N of natural numbers.(Notice that the noncontradictoriness of the ①-based computational methodology has been studied in depth in [15][16][17].)From the practical point of view, this methodology has given rise both to a new supercomputer patented in several countries (see [18]) and called Infinity Computer and to a variety of applications starting from optimization (see [12,[19][20][21][22][23][24]) and going through infinite series (see [13,[25][26][27][28]), fractals and cellular automata (see [25,[29][30][31][32]), hyperbolic geometry and percolation (see [33,34]), the first Hilbert problem and Turing machines (see [13,35,36]), infinite decision making processes and probability (see [13,[37][38][39]), numerical differentiation and ordinary differential equations (see [40][41][42][43]), etc.This methodology does not contradict traditional views on infinity and infinitesimals (Cantor, Leibnitz, Robinson, etc.) and proposes just another, more computationally oriented, way to deal with these objects.In particular, in order to avoid misunderstanding it should be stressed that there exist several differences (see [44] for a detailed discussion) that distinguish the numerical ①-based methodology from symbolically oriented nonstandard analysis of Robinson.Another important preliminary remark is that symbols traditionally used to work with infinities and infinitesimals (∞ introduced by Wallis, Cantor's ω, ℵ 0 , ℵ 1 , ..., etc.) are not used together with ①.Similarly, when the positional numeral system and the numeral 0 expressing zero had been introduced, symbols V, X and other symbols from the Roman numeral system had not been used in the positional numeral system.
The numeral ① allows one to construct different numerals involving infinite, finite, and infinitesimal parts and to execute numerical computations with all of them in a unique computational framework.As a result, it becomes possible to execute arithmetical operations with a variety of different infinities and infinitesimals.As a remarkable result, indeterminate forms such as ∞ − ∞ or ∞ • 0 are not present when one works with numbers expressed in the ①-based numeral system.Traditionally existing kinds of divergences do not appear, as well.They are substituted by expressions that can contain also finite, infinite and infinitesimal parts.
In order to give some examples of arithmetical operations that can be executed in the ①-based numeral system, let us consider the following numbers: ① and ① 16.3  (that are examples of infinities), along with ① −1 and ① − 16.3 (that are examples of infinitesimals).Then, we can compute, for instance, the following expressions: Table 1 A practical implementation of the CG algorithm for the symmetric positive definite linear system Ax = b, with A ∈ R n×n The Conjugate Gradient (CG) method Data: In general, in the ①-based numeral system the simplest infinitesimal numbers are represented by numerals having only negative finite powers of ① (e.g., 40.17① −13.26 +87.32① −25.7 , see also examples above).The simplest infinite numbers are represented by numerals having at least one positive power of ①.Then, it can be seen in ( 3) that ① 0 = 1; therefore, a finite number a can be represented in the new numeral system simply as a① 0 = a, where the numeral a itself can be written down by any convenient numeral system used to express finite numbers.These numbers are called purely finite because they do not contain infinitesimal parts.For instance, number 5 is purely finite and 5 + 3① −16.3 is finite but not purely finite, because it contains the infinitesimal part 3① −16.3 .Notice that all infinitesimals are not equal to zero.In particular, 1  ① > 0 because it is a result of division of two positive numbers.

The Matrix Factorizations We Need
In the current and the next section, we describe how to use some Krylov-subspace methods, in order to gain advantage of a suitable factorization for the (possibly) indefinite Hessian matrix ∇ 2 f (x j ).We strongly remark that we never explicitly compute here any Hessian decomposition, since our final achievements definitely rely on implicit decompositions, induced by Krylov-based methods.
As a general result, we highlight that computing negative curvature directions for f (x) at x j , which match the requirements in Assumption 2.1, may reduce to a simple task when suitable factorizations of ∇ 2 f (x j ) are available.To give an intuition of the latter fact, suppose both the relations are available at iterate x j , being M j ∈ R n×k , where C j , Q j , B j ∈ R k×k are nonsingular.In this regard, the CG (Table 1) is an example of a Krylov-based method, satisfying the following properties: -it provides the decompositions in (4) (see also (6)) when ∇ 2 f (x j ) 0, with j ≥ 1; -the matrices M j , C j , Q j , B j have a special structure, inasmuch as (see also [45]): The columns of M j are unit orthogonal vectors, C j is tridiagonal, Q j is unit lower bidiagonal, and B j is diagonal.
Note that the Lanczos process, which represents another renowned Krylov-based method in the literature, iteratively provides the left decomposition in (4), with C j tridiagonal, but not the right one.It is indeed necessary to couple the Lanczos process with a suitable factorization of C j , in order to obtain usable negative curvature directions or solvers for Newton's equation (see, e.g., SYMMLQ/MINRES [46], SYMMBK [47]).Now, given (4) suppose the vector w ∈ R k is an eigenvector of B j , associated with its negative eigenvalue λ, whose computation is 'relatively simple.'Moreover, suppose the vector y ∈ R k is easily available, such that Q T j y = w.Then, by ( 4) the equalities/inequalities In particular, thanks to the chain of equalities above, if λ is the smallest negative eigenvalue of B j , then M j y is also an eigenvector of ∇ 2 f (x j ), associated with the smallest eigenvalue of ∇ 2 f (x j ).The most renowned Krylov-subspace methods for symmetric linear systems (i.e., SYMMLQ, SYMMBK, CG, Planar-CG methods [48][49][50]) can all provide the factorizations (4) when applied to solve Newton's equation at the iterate x j .Hence, generating a negative curvature which satisfies (a) in Assumption 2.1, may not be a difficult goal.However, fulfilling also (b) and the boundedness of the latter negative curvature direction, is a less trivial task.Indeed, the counterexample in Sect. 4 of [7] issues such a drawback, when a modified Cholesky factorization of the Hessian matrix is possibly adopted.
We strongly remark this point, since our main effort here is that of coupling a Krylov-subspace method with the novel tool in the literature given by grossone.In particular, we want to show that by the use of a subset of properties which hold for grossone, we can yield an implicit matrix factorization as in (4), fulfilling also (b) and the boundedness of the final negative curvature direction d j in Assumption 2.1.
On this purpose, let us first state a general formal result for Krylov-subspace methods, which possibly summarizes the above considerations.The proof of the next lemma easily follows from Lemma 4.3 in [7] and Theorem 3.2 in [9].Lemma 4.1 Let problem (1) be given with f ∈ C 2 (R n ), and consider an iterative method for solving (1), which generates the sequence {x j }.Let the level set L 0 = {x ∈ R n : f (x) ≤ f (x 0 )} be compact, being any limit point x of {x j } a stationary point for (1), with λ[∇ 2 f ( x)] > λ > 0. Suppose n iterations of a Newton-Krylov method are performed to solve Newton's equation (2) at iterate x j , for a given j ≥ 0, so that the decompositions are available.Moreover, suppose R j ∈ R n×n is orthogonal, T j ∈ R n×n has the same eigenvalues of ∇ 2 f (x j ), with at least one negative eigenvalue, and L j , B j ∈ R n×n are nonsingular.Let z be the unit eigenvector corresponding to the smallest eigenvalue of B j , and let ȳ ∈ R n be the (bounded) solution of the linear system L T j y = z.Then, the vector d j = R j ȳ is bounded and satisfies Assumption 2.1.
The vector d j computed in Lemma 4.1 may be used to guarantee the satisfaction of Assumption 2.1, i.e., the sequence {d j } can guarantee convergence to second-order critical points.However, three main drawbacks of the approach in Lemma 4.1 are that (α) the eigenvector z of B j and the solution of the linear system L T j y = z should be of easy computation; (β) the corresponding vector ȳ should be provably bounded; (γ ) at iterate j the Newton-Krylov method adopted to solve (2) possibly does not perform n iterations.
Observe that according to the requirements in Assumption 2.1, after a careful consideration, the issue at item (γ ) is not really so relevant.Indeed, in any case (see also [9]) when j → ∞ the convergence of the Newton-Krylov method imposes that it eventually performs n iterations [51].On the other hand, in case at iterate x j , for a finite j, ∇ 2 f (x j ) 0 or a vector v ∈ R n such that v T ∇ 2 f (x j )v < 0 is unavailable, then the factorization (5) yet exists and we can simply set d j = 0, which satisfies (a) in Assumption 2.1 along with the boundedness requirement.Though the CG is not well posed when ∇ 2 f (x j ) 0, in [9] the authors reported that, in case n CG steps are performed without stopping when solving Newton's equation, the above items (α) and (β) can be relatively easily fulfilled, even in case ∇ 2 f (x j ) is indefinite.In particular, these results are obtained exploiting the factorizations in Lemma 4.1, for which the CG specifically yields (Table 1) Thanks to the above expressions of R j , B j and L j , in [9] the authors proved that after n steps the CG straightforwardly yields also the bounded negative curvature direction Moreover, d j in (7) satisfies Lemma 4.1, thanks to the fact that B j is diagonal (i.e., its eigenvectors coincide with the canonical basis), L j is unit lower bidiagonal (so that the solution of L T j y = e m is straightforwardly available by backtracking) and d j is provably bounded.
Our goal is that of possibly replicating an analogous reasoning, with other Krylovbased methods for indefinite linear systems, following similar guidelines.In this regard, observe that both the tasks (α) and (β) might be hardly guaranteed only using, for instance, the instruments in [14], essentially because comparing with the CG, the structure of the matrices L j and B j generated by the Planar-CG method in [14] is more cumbersome.Nevertheless, in the next sections, starting from the structure of matrices L j and B j , as computed by the algorithm in [14], we will show how to use CG ① in [12] in order to fulfill the hypotheses of Lemma 4.1.

Our Proposal: Preliminaries
To fill the gap outlined in the previous section, and recalling that in Lemma 4.1 we focus on the case where j → +∞, let us set for the sake of simplicity . This allows us to drop the dependency on the subscript j.Consider the method CG ① in [12] (which is also reported in Table 2, for the sake of completeness.Observe that the practical implementation of Step k in CG ① currently allows the test p T k Ap k = 0 to be replaced by the inequality The CG ① substantially coincides with the CG, as long as p T k Ap k = 0.Moreover, in case at Step k we have p T k Ap k = 0, from Section 5.1 of [12] the CG ① generates both the vectors r k+1 and p k+1 , such that they depend on ①.Furthermore, we have (after a simple computation, and using the standard Landau-Lifsits notation O(•)) Table 2 The CG ① algorithm for solving the symmetric indefinite linear system Ax = b, A ∈ R n×n .In a practical implementation of Step k of CG ① , the test p T k Ap k = 0 may be replaced by the inequality CG ① : Conjugate Gradient method coupled with grossone ① , and Recalling that by definition s = O(① −2 ), i.e., s① = O(① −1 ), then neglecting in α k+1 the terms containing negative powers of ① (corresponding indeed to infinitesimals with respect to the value ① 0 = 1-see Sect.3) we have This immediately implies that Note that by using CG ① , similarly to the CG to solve (2), we can recover the structure of B j and L j in (6), so that (7) formally applies.However, since s① is infinitesimal in (9), after some computation we have (see also Sect.5.1 of [12]) Then, in case m = k + 1 in (7), to compute the direction d j we would have which implies that d j is not bounded (being s① = O(① −1 )) and Lemma 4.1 cannot be fulfilled.This consideration should not be surprising.Indeed, it basically summarizes the fact that similarly to the CG, CG ① is unable to provide the diagonal matrix ( 6) with finite entries, in case the tridiagonal matrix T j in ( 5) is indefinite.Nevertheless, to overcome the latter limitation, now we show how to properly couple CG ① with the Planar-CG method in [14], in order to obtain a suitable bounded negative curvature direction which fulfills Lemma 4.1.

Coupling CG ① with the Algorithm in [14]
Following the taxonomy in Sect.4, assume without loss of generality that the Krylovsubspace method detailed in [14] is applied to solve Newton's equation and n steps are performed. 1Again this allows us to drop the dependency on the subscript j.After some computation, the following matrices are generated (see also [52]) where and where (see also [52]) the matrix R ∈ R n×n has n unit orthogonal columns and is given by with r k+1 = Ap k , while T ∈ R n×n is tridiagonal.Moreover, {α i }, {β i }, e k+1 are suitable scalars (being e k+1 = (Ap k ) T A(Ap k )/ Ap k 2 ).We also recall that β i > 0, for any i ≥ 1.Finally, in the above matrices L and B we assume (for the sake of simplicity) that the Krylov-based method in [14] has performed all CG steps, with the exception of only one planar iteration (namely the kth iteration-see [14] and [48]), corresponding to have p T k Ap k ≈ 0.Then, our novel approach proposes to introduce the numeral grossone, as in [13,[22][23][24], and follows some guidelines from [12], in order to exploit a suitable matrix factorization from (16), such that Lemma 4.1 is fulfilled.In this regard, consider matrix B in ( 16) and the next technical result.(16), and let
, along with λ k > 0 and λ k+1 < 0. Finally, if r i ≥ ε, for any i ≤ k, then Proof The first part of the proof follows after a short computation and observing that det As regards (20), since we have Then, replacing the factorization (18) into the expression of B in (16), we obtain the equivalent factorization T = L B L T = L B LT , where where L 11 , L 21 are defined in (13), L 33 in ( 14), B 11 , B 33 in (15) and We remark that unlike the matrix B, now B is a diagonal matrix, though L has now a slightly more complex structure than the matrix L. Note also that after an easy computation, we have in L Table 3 Correspondence between quantities/vectors computed by the algorithm in [14] and the algorithm CG ① in [12] Algorithm in [14] C G ① in [12] where (see [14]) Now, let us consider again the algorithm CG ① in [12], and assume that at Steps k and k + 1 it generated the coefficients α k and α k+1 in (9), when solving the linear system (12), being p T k Ap k ≈ 0 at Step k.In [12], we have already detailed the one-to-one relationship between the quantities generated by the algorithms in [14] and in Table 2, showing how CG ① can be considered, to large extent, an extension of the CG to the indefinite case.Table 3 specifically reports this relationship, showing how it can be possible to compute all the quantities in (23) using CG ① , in place of the algorithm in [14].Thus, similarly to the result obtained in (6), applying the CG, after n steps of CG ① we want to define an implicit matrix factorization for A as in (16), where now the 2 × 2 matrix on the left-hand side of ( 18) is suitably replaced by the matrix diag{1/α k , 1/α k+1 }.Now we establish a full correspondence between the matrix k in ( 18) and ( 23), obtained by the algorithm in [14], and the matrix diag{1/α k , 1/α k+1 } from [12].Since both α k s① = 0 and α k+1 /(s①) = 0, and by Lemma 5.1 λ k > 0 with λ k+1 < 0, we can always find the 2 × 2 (diagonal) positive definite matrix C k such that where k is defined in Lemma 5.1 and In practice, using CG ① we would like to rearrange the matrices L and B in ( 23), obtained by applying the algorithm in [14], so that the equalities T = L B L T = L B LT hold and the block k in B is suitably replaced by the left side of (25).Note that the diagonal matrix on the left side of ( 25) is scaled with respect to the matrix diag{1/α k , 1/α k+1 }, by using terms containing ①.Moreover, it is worth mentioning that by ( 8) and ( 9), both the diagonal entries of the matrix on the left side of ( 25) are finite and not infinitesimal.The rationale behind this choice is suggested by (11) and Lemma 4.1, where the easy computation of the vectors z and ȳ is sought.Indeed, we shortly show that the scaling in (25) both allows to easily find the final negative curvature direction d j in Lemma 4.1, and ensures that for any j the norm d j is suitably bounded.This finally implies that applying CG ① and exploiting Table 3, we can fulfill Assumption 2.1 without recurring first to the algorithm in [14].Now, from ( 26) and ( 9) we obtain showing that, apart from infinitesimals we ignored when writing (9), the diagonal entries of C k are independent of ①.Finally, by (25) and considering the matrix k in Lemma 5.1, we obtain This also implies that we can now equivalently modify the nonsingular matrices in (23) as where L 11 and L 21 are defined in ( 13), L 33 in ( 14), B 11 , B 33 in ( 15) and so that in Lemma 4.1 we have for matrix T j the expression We strongly remark that using CG ① and relation (25), we have simplified the expression of D, replacing it with D. This is obtained at the cost of a slight modification of matrix L into L: We shortly prove that this arrangement can easily allow the computation of a bounded negative curvature direction d j at x j .Once more we urge to remark that the computation of L and D can be completely carried on replacing the algorithm in [14] with CG ① , as the equivalence/correspondence in Table 3 reveals.(We highlight indeed that the iterate x k+2 in [14] and the iterate y k+2 in [12] coincide, when neglecting the infinitesimal terms containing s①.)The next lemma proves that L in ( 28) is nonsingular under the assumptions in Lemma 4.1.

Lemma 5.2 Let the assumptions in Lemma 4.1 hold, with T j = L D LT and L, D
defined in (28).Then, we have Proof The first two relations follow immediately from (27), Table 3 and recalling that C k is nonsingular.Moreover, since Now we are ready to compute at iterate x j the negative curvature direction d j which complies with Assumption 2.1, exploiting the decomposition T j = L D LT from Lemma 4.1.
Proposition 5.1 Suppose n iterations of CG ① algorithm are performed to solve Newton's equation (2), at iterate x j , so that the decompositions exist, where R is defined in (17), and L along with D is defined in (28).In the hypotheses of Lemma 4.1, let z be the unit eigenvector corresponding to the (negative) smallest eigenvalue of D and let ŷ be the solution of the linear system LT y = z.Then, the vector d j = R ŷ is bounded and satisfies Assumption 2.1.In addition, the computation of d j requires the storage of at most two n-real vectors.
Proof First observe that by [53], also in case at the iterate x j the Hessian matrix ∇ 2 f (x j ) is indefinite, there exists at most one step k, with 0 ≤ k ≤ n, such that in CG ① we might have p T k ∇ 2 f (x j ) p k = 0. Thus, similarly to the rest of the paper, without loss of generality in this proof we assume that possibly the equality p T k ∇ 2 f (x j ) p k = 0 only holds at step k.Moreover, the matrix D is diagonal, which implies that the unit vector associated with its ith eigenvalue μ i ( D) is by e i .
To fulfill Assumption 2.1, we first need to compute the vector ŷ in Lemma 4.1, i.e., we have to solve the linear system being z ∈ R n the unit eigenvector associated with the (negative) smallest eigenvalue of D. To this purpose, by Lemma 5.2 the vector ŷ exists and is bounded.Now, we distinguish among the next four subcases, where we use the notation k ∈ arg min i {μ i ( D)}, i.e., k is an index corresponding to the smallest eigenvalue μ k ( D) of D.
(I) In this subcase, we assume k / ∈ {k, k + 1} along with k < k.In particular, since D is diagonal, then (29) reduces to LT y = e k , i.e., by Lemma 5.2 and Table 3 whose solution ŷ ∈ R n can be explicitly computed recalling that, as in Table 3, 2 / r i 2 and backtracking from the value of ŷn to ŷ1 , we have Finally, as in Lemma 4.1 and recalling that for the algorithm CG ① we have p i = r i +β i−1 p i−1 , for any i ≥ 1, the corresponding negative curvature direction d j is given by which exactly coincides with the proposal in [9], when k / ∈ {k, k + 1} along with k < k.Finally, it is easily seen that the conditions r i ≥ ε, from algorithm CG ① , the quantity d j is bounded and the computation of d j simply requires the storage of the unique vector p k / r k .(II) In this subcase, we assume k / ∈ {k, k + 1} along with k > k + 1.Since again D is diagonal, then (29) reduces to Thus, again backtracking from the value of ŷn to ŷk +1 we first obtain Then, we have also while for ŷi , i ∈ {k + 1, k}, we have from above the relations Observing that by Lemma 5.1 λ k = λ k+1 , and recalling that in Table 3 √ which allow to backtrack and compute also the remaining entries ŷk−1 , . . ., ŷ1 of vector ŷ, being On the overall, the final computation of the negative curvature direction d j yields for this subcase Finally, following the guidelines in Table 2 of [9], the conditions r i ≥ ε from algorithm CG ① yield that d j is bounded.Moreover, with a similar analysis in [9] the computation of d j requires the storage of just two vectors.(III) In this subcase, we assume k = k.However, note that this subcase can never occur, since by ( 25) and therefore no negative curvature direction can be provided from the current step k. (IV) As a final subcase, we assume k = k + 1, i.e., by (25) Again, since D is diagonal, then the linear system (29) reduces to LT y = e k (or equivalently LT y = e k+1 ), with Now, we have for the last n − k entries of vector ŷ the expression On the other hand, the condition ŷk +1 = ŷk+2 = 0 and the above relation Recalling that now k = k + 1 and in Table 3 Ap k = r k+1 = r k , then As a consequence, and for ŷk −2 , . . ., ŷ1 , we have Finally, the overall negative curvature direction d j becomes now whose computation is well posed, since λ k+1 < 0. Again, by ( 20)-( 21), the fact that λ > 0 in Lemma 4.1 and the other hypotheses, the quantity d j is bounded.
In addition, the computation of d j evidently needs the storage of just two n-real vectors.
Observation 5.1 We remark that the computation of the negative curvature direction d j requires at most the additional storage of a couple of vectors, which confirms the competitiveness of the storage proposed in [9].Thus, the approach in this paper does not only prove to be applicable to large-scale problems, but it also simplifies the theory in [9], which is currently in our knowledge the only proposal of iterative computation of negative curvatures for large-scale problems, which does not need any recomputing (as in [3]), and which requires neither a full matrix factorization nor any matrix storage.
in [9], where we replaced the Krylov-based iterative procedure therein by the CG ① procedure.The codes were written in Fortran compiled with Gfortran 6 under Linux Ubuntu 18.04, and the runs were performed on a PC with Intel Core i7-4790K quadcore 4.00 GHz Processor and 32 GB RAM.Now we strongly remark the guidelines and the limits of the numerical experience reported in this section: -we show how to detect and assess negative curvatures for the Hessian matrix ∇ 2 f (x j ), at the current iterate x j ; -we compute negative curvatures which could be able to guarantee the overall convergence of the optimization method toward second-order critical points; -we do not claim that our proposal shows better numerical results with respect to [9], being the main focus of this paper on theoretical issues.Thus, our numerical experience only tests the reliability and the effectiveness of our method, rather than proposing a numerical comparison with the current literature; -we also intend to check for the quality of the stationary points detected by our approach.
In particular, we considered all the 112 large-scale unconstrained test problems in CUTEst [54] suite.The algorithm performs a classic nested loop of outer-inner iterations.Thus, at the current jth outer iteration the algorithm iteratively solves Newton's equation ∇ 2 f (x j )s = −∇ f (x j ), performing a certain number of inner iterations.To build an approximate solution s j of Newton's equation, and possibly a negative curvature direction, inner iterations are stopped whenever the following truncation rule is satisfied being {η j } a forcing sequence, with η j → 0. The condition η j → 0 guarantees superlinear convergence of the overall method, when close enough to the final stationary point.As regards settings and parameters of the linesearch procedure we adopted, as well as the overall stopping criterion, the reader can refer to [9].(We also recall that, unlike in [9], here we preferred not to include any nonmonotonicity in the used algorithm, in order to clearly distinguish the contribution from our idea.)At each inner iteration k ≥ 1, the algorithm in Table 2 detects a curvature of the objective function, by computing the term p T k ∇ 2 f (x j ) p k , a negative value of the last quantity indicating a negative curvature direction.
We compared two truncated Newton methods: the first (i) not including the use of negative curvature directions (namely NoNegCurv), so that convergence to simple stationary points could be guaranteed; the second (ii) including negative curvatures (namely NegCurv) which satisfy Assumption 2.1, implying convergence to stationary points where second-order necessary optimality conditions are fulfilled.Thus, by a comparison between them, we might have expected: (a) (ii) to be more efficient than (i) in terms of the computational effort (in our largescale setting we measured the computational effort through the number of inner iterations, which are representative of the overall computational burden, including CPU time); (b) the quality of the solutions detected by (ii), i.e., the value of the objective function at the solution, is expected to be on average not worse than in the case of (i), since for (ii) the solution points satisfy additional theoretical properties; (c) the stationarity (measured by ∇ f (x * ) ) of the final solution detected using (ii) is possibly expected to be competitive with respect to (i).This is because in a neighborhood of the solution point, our proposal is expected to collect more information on the objective function.
The above considerations are to large extent confirmed by our numerical experience as detailed in the following.First, note that using (ii), we detected negative curvatures on 40 test problems out of 112; of course, this does not imply that the remaining 72 test problems only include convex functions.It rather implies that on 72 problems, no regions of concavity for the objective function were encountered.For these 40 test problems, the obtained results in terms of number of (outer) iterations (it), number of function evaluations (nf), number of inner iterations (inner-it), optimal function value ( f (x * )), gradient norm at the optimal point ( g(x * ) ), solution time in seconds (time), are reported in Tables 4 and 5.In particular, for each test problem, we report results using both the NoNegCurv method (top row) and the NegCurv method (bottom row).
By observing these results, we first note that on two test problems both the algorithms fail to converge within the maximum CPU time of 900 s.The comparison on the remaining test problems shows that in most cases algorithm NegCurv performs the best in terms of solution time and inner iterations, confirming expectation (a) which is our main goal.The results highlight only one test problem (GENHUPMS 1000) where the use of NegCurv yields a significant worsening of the performance.We easily realize that including our procedure to compute negative curvature directions allows to both speed up the overall convergence and decrease the number of inner iterations.
The detailed results only partially validate also (b) and (c).As regards (b), since in a few test problems the algorithms converge to different points, a sound statistical analysis cannot be given, though a better optimal value is sometimes observed by using NegCurv algorithm.Similarly, as concerns (c), the values of ∇ f (x * ) provided by NoNegCurv and NegCurv seem to a large extent comparable on this test set.
To have an overview of the effectiveness and the robustness of the approach we propose in this paper, we now consider summary results by using performance profiles [55].The performance profiles represent a popular and widely used tool for providing objective information when benchmarking optimization algorithms.Their meaning can be summarized as follows: Suppose you have a set of solvers S to be compared on a set of test problems P.
For each problem p ∈ P and solver s ∈ S, define t ps the statistic obtained by running solver s on problem p. Namely, t ps ≥ 0 is a performance measure of interest, e.g., solution time, number of function evaluations, etc.The performance on problem p by solver s is compared with the best performance by any solver on this problem by means of the performance ratio r ps = t ps {t ps | s ∈ S} .Moreover, an upper bound r is chosen such that r ps ≤ r for all s ∈ S and p ∈ P and if a solver s fails to solve problem p then r ps is set to r .The performance profile of the solver s is the function namely the cumulative distribution function for the performance ratio.
In particular, we report in Figs. 1 (full profile) and 2 (detail profile) the performance profiles comparing NoNegCurv and NegCurv algorithms in terms of inner iterations.The test set considered includes the 40 test problems reported in Tables 4 and 5.
The detailed plot reported in Fig. 2 clearly shows the effectiveness of NegCurv algorithm with respect to NoNegCurv.Indeed, as an example, let us consider the value of abscissa 1.2 in Fig. 2. The plots show that the NegCurv algorithm is able to solve about 62% of the test problems within 1.2 times the number of inner iterations of the best algorithm.Conversely, the NoNegCurv algorithm is able to solve up to 78% of the test problems within the same number of inner iterations.On the other hand, in terms of robustness the algorithms can be considered comparable with a slight preference for NegCurv algorithm as evidenced by Fig. 1.The last consideration follows from the observation that for values of the abscissa parameter larger than 3.5, the two plots basically tend to overlap.

Conclusions
We proposed a novel approach for the efficient solution of large-scale unconstrained optimization problems, where the detected solutions likely are endowed with strong theoretical properties.Our proposal exploits the simplicity of an algebra associated with the numeral grossone, which was recently introduced in the literature to handle infinite and infinitesimal quantities.
We were able to extend the results in [9] in view of a theoretical simplification, avoiding to make reference to Planar-CG methods which require a more complex analysis.The theory in this paper allows us to guarantee that the iterative computation of negative curvatures does not need any matrix storage, while preserving convergence toward points satisfying second-order necessary optimality conditions.
Then, we also provided numerical results, which show the efficiency of our proposal.We remark that the focus of this paper is not on a numerical comparison among different algorithms which exploit negative curvature directions.Rather, we paired the approach in [9] with a novel paradigm provided by grossone, in light of preserving numerical efficiency within a sound theoretical framework, in dealing with nonconvex problems.This is the first stage toward a more complete numerical experience where the iterative algorithm CG ① can be fully tested, including even more challenging problems from real-world applications.
Observe that the proposed approach is independent under multiplication of the function by a positive scaling constant or adding a shifting constant.This is an important property that is specially exploited in the global optimization framework (see, e.g., [24]), since strongly homogeneous algorithms are definitely appealing.Furthermore,

(
top row) and NegCurv (bottom row) are reported.

(
top row) and NegCurv (bottom row) are reported.

Fig. 1 Fig. 2
Fig. 1Performance profile relative to a comparison between (ii) and (i), with respect to the number of inner iterations, where negative curvatures are (red line)/are not (green

Table 4
Complete results for CUTEst test problems where negative curvature directions are encountered.For each test problem, results for NoNegCurv

Table 5
Complete results for CUTEst test problems where negative curvature directions are encountered.For each test problem, results for NoNegCurv