Documents Live, a web authoring and publishing system
If you see this, something is wrong
Table of contents
First published on Tuesday, Jul 7, 2026 and last modified on Tuesday, Jul 7, 2026 by François Chaplais.
CNRS, ENS de Lyon, Inria, Université Claude Bernard Lyon 1, LIP, UMR 5668, Lyon, France Email
ENS de Lyon, CNRS, Inria, Université Claude Bernard Lyon 1, LIP, UMR 5668, Lyon, France Email
Sorbonne Université, CNRS, LIP6, F-75005 Paris, France Email
ENS de Lyon, CNRS, Inria, Université Claude Bernard Lyon 1, LIP, UMR 5668, Lyon, France Email
Newton's method, mixed precision, error analysis, floating-point arithmetic, inexact Newton, quasi-Newton.
Abstract
Second-order optimization methods, such as Newton’s algorithm, achieve fast local convergence and high accuracy, but their practical use is often limited by high computational costs. To mitigate this issue, variants such as inexact and quasi-Newton methods are widely used. A complementary and promising approach to improve the efficiency of the method is to employ mixed precision arithmetic, using different floating-point precisions for different operations, based on their impact on the convergence and accuracy of the method. In this work, we perform an error analysis of Newton’s method accounting for different sources of inexactness, including approximations and rounding errors. We present a convergence analysis for the generated sequence, establishing bounds on the convergence rate and attainable accuracy. This theoretical framework covers quasi-Newton and inexact Newton methods, and is leveraged to propose mixed precision algorithms. We present a wide set of numerical experiments to illustrate our theoretical results and the behavior of Newton’s method and its approximate variants in mixed precision floating-point arithmetic.
1 Introduction
Modern computational science, including fields such as machine learning, inverse problems, image restoration, and physical simulations, relies on the efficient solution of complex large-scale optimization problems. These problems are often ill-conditioned or ill-posed, making them particularly challenging to solve efficiently. Second-order methods, most notably Newton’s method, address this challenge by exploiting curvature information through the Hessian matrix of second derivatives. By accounting for the local geometry of the objective function, they generate search directions that are both appropriately scaled and well oriented. As a result, second-order methods are often the method of choice in these settings, offering robust performance together with quadratic convergence in a neighborhood of the optima [1, 2].
Despite these advantages, the adoption of Newton’s method is hindered by its high computational cost per iteration, which is due to the need of forming the Hessian and solving the resulting (potentially large) linear system needed to compute the search direction. To mitigate these costs, various approximations of Newton’s method have been developed, including quasi-Newton methods that replace the Hessian by some approximations [3], and inexact Newton methods that solve the linear systems approximately [4]. Nevertheless, even approximate Newton methods remain quite computationally intensive and so, in this work, we are interested in another promising and complementary direction to reduce the computational burden of second-order methods: the use of low precision floating-point arithmetic. The rise of specialized hardware accelerators, such as NVIDIA’s tensor cores [5], has driven the successful use of low precision arithmetic to reduce computational and memory costs. However, the naive application of low precision can lead to a significant loss of accuracy, and, in the case of iterative methods, to a slower (or lack of) convergence. In order to minimize computational cost while preserving convergence and solution accuracy, the use of mixed precision algorithms has spread in different fields. Notably, many such algorithms have been successfully developed in numerical linear algebra; see [6] for a survey. In contrast, the use of mixed precision arithmetic for general nonlinear optimization algorithms remains largely unexplored and lacks a unifying theoretical framework.
In this work, we address this theoretical gap by proposing a rigorous error analysis for Newton’s method for optimization, accounting for various sources of inexactness. Our central contribution is to perform a convergence study of Newton’s method under the assumption that the three main steps of the algorithm (gradient computation, Hessian-related steps, and iterate update) are computed inexactly, with possibly different levels of precision. This study allows for assessing the impact of the perturbed operations on the final solution accuracy and on the convergence rate. In turn, this allows for deriving precise guidelines for assigning a different precision to the various operations.
Importantly, our error analysis is general enough to encompass a wide range of approximations, including not only floating-point arithmetic, but also inexact Newton and quasi-Newton methods (for the latter, we specifically focus on the Gauss–Newton method). This allows us not only to recover known convergence results for such variants of Newton’s method when all the operations are performed in exact arithmetic, but also to incorporate rounding errors in such variants. This leads us to analyze and propose mixed precision implementations of inexact and quasi-Newton methods, whose convergence and accuracy is covered by our theoretical results. Our analysis highlights the interplay between approximation and rounding errors, and shows how to balance both sources of error.
We illustrate our theoretical findings through extensive numerical experiments, demonstrating the soundness and generality of our error analysis in practice and providing useful insights into the behavior of Newton’s method and its variants in mixed precision arithmetic.
1.1 Related work
Mixed precision optimization (meaning optimization methods in which different computations can have different accuracy, or different quantities can be subject to errors of different magnitudes) has been widely studied in the literature. The predominant approach consists in analyzing the convergence of optimization algorithms in the presence of inexact function or derivative evaluations [7].
Many works consider frameworks in which the accuracy of the estimates increases with time, imposing a decreasing absolute or relative error on the function and gradient approximations to ensure the convergence properties of the methods [8, 9, 10]. Most often these methods are studied in the finite sum context, suited to machine learning applications, where the approximations are built by subsampling techniques [11, 12, 13, 14]. Another line of work inherits techniques from derivative-free optimization and adopts fully probabilistic frameworks based on the fully linear assumption on the models employed [9, 10, 15]. Most of the works on inexact optimization consider globally convergent methods and propose a worst-case analysis [16, 17, 18], which counts the number of iterations necessary to drive the norm of the gradient below a given threshold; this leads to bounds that are usually pessimistic and rarely observed in practice.
We adopt a different perspective here, coming from numerical and error analysis. We consider Newton’s method without a globalization strategy and we focus on a deterministic local convergence analysis, to target sharper and more informative bounds on the convergence of the iterates. From this perspective many works have considered root-finding Newton’s method, see for example [19, 20, 21, 22, 23, 24, 25, 26]. These tighter analyses allow in particular for assessing much more precisely the effect of rounding errors on Newton’s methods. Not many works in optimization are concerned with this aspect; we mention [27, 28], which focus on first-order methods, and [29], which, although studying a second-order optimizer, does not address Hessian approximations.
Our work is closest to Tisseur’s [24], which is focused on root-finding Newton’s method. Our analysis is an adaptation and extension of that in [24] to the distinct setting of unconstrained optimization. While the two methods are intimately related, the analysis for nonlinear systems does not directly apply to optimization problems, and we believe that the latter deserves a proper dedicated analysis. In particular, specializing the analysis to optimization problems allows for taking into account sources of errors that are specific to this context, such as gradient or Hessian approximations. Moreover, while adapting Tisseur’s analysis to our context, we have also made some small changes and improvements to make the bounds more readable and slightly sharper.
1.2 Organization of the paper
The article is organized as follows. Section 2 presents the error analysis of mixed precision Newton’s method in a general framework. Section 3 discusses the application of our framework to floating-point arithmetic, and inexact and Gauss–Newton’s methods. Section 4 presents numerical experiments with all these variants. Section 5 provides concluding remarks.
1.3 Notations
All computed quantities are denoted by a hat. We denote by \( \lVert \cdot \rVert \) any vector norm and the corresponding operator norm, unless otherwise specified, and by \( \kappa(A) = \lVert A \rVert \lVert A^{-1} \rVert \) the condition number of a matrix \( A\) . Finally, given a point \( x\in \mathbb{R}^n\) and a radius \( \rho >0\) , we denote by \( B_{\rho}(x) = \{y \in \mathbb{R}^n : \lVert y - x \rVert < \rho\}\) the open ball of radius \( \rho\) centered in \( x\) .
2 Mixed precision Newton’s method for optimization
Given a twice continuously differentiable function \( f: \mathbb{R}^n \to \mathbb{R}\) , a general minimization problem can be stated as follows:
(1)
We denote by \( g:\mathbb{R}^n\rightarrow \mathbb{R}^n\) the gradient of \( f\) and by \( H:\mathbb{R}^n\rightarrow \mathbb{R}^{n\times n}\) its Hessian. Newton’s method for optimization is an iterative algorithm that generates a sequence of approximations \( \{x_{i}\}\) to a minimizer \( x^*\) of \( f\) using the following update rule, assuming that \( H(x_i)\) is positive definite:
(2)
The behavior of Newton’s method in exact arithmetic is well understood [1, Thm. 3.5]. If initialized close enough to a solution, Newton’s method generates a well-defined sequence that converges to that solution at a quadratic convergence rate. In this section, our goal is to study the behavior of the method when the operations in Equation 2 are subject to errors. These errors may simply be rounding errors due to finite precision or approximation errors that arise when inexact or quasi-Newton variants are used. Moreover, we allow different operations to be affected by errors of different size, and so we refer to this method as mixed precision Newton’s method.
Specifically, we consider the following error model, where computed quantities affected by errors are marked by a hat.
Error Model 1
At each iteration \( i\) , mixed precision Newton’s step satisfies
(3.a)
(3.b)
(3.c)
(3.d)
(3.e)
This model depends on three error terms \( e_i^g\) , \( E_i^H\) , and \( e_i^+\) , which are each bounded as follows.
- We bound the norm of the gradient error \( \lVert e_i^g \rVert \) with an absolute error of size \( \epsilon^g_i\) , which is a general term that accounts for inexactness in the gradient evaluation.
- We bound the norm of the Hessian error \( \lVert E_i^H \rVert \) with a relative normwise error of size \( \epsilon^H_i\) . This accounts for both the error incurred in forming the Hessian and the backward error for solving the associated linear system.
- We bound the norm of the update error \( \lVert e_i^+ \rVert \) with a relative normwise error of size \( \epsilon_i\) . This accounts for the error in updating the iterate and rounding it to the working precision.
By using the subscript \( i\) , we account for the nonconstant behavior of rounding errors, and for potential adaptive precision strategies that vary the precisions across iterations—even though exploring this setting is outside the scope of this article, it remains a field of interest for future work.
We will assume Lipschitz continuity on the Hessian, as defined in Definition 1, which is quite standard to ensure quadratic convergence [1].
Definition 1
A function \( \varphi\) is said to be Lipschitz continuous in an open set \( \Omega \subseteq \mathbb{R}^n\) if there exists a constant \( L \geq 0\) , called the Lipschitz constant, such that
(4)
Our main result is stated in the following theorem.
Theorem 1
Let \( f:\mathbb{R}^n\rightarrow\mathbb{R}\) be twice continuously differentiable. Let \( x^* \in \mathbb{R}^n\) be a minimizer for \( f\) , and assume that \( H(x^*)\) is nonsingular. Let \( H\) be Lipschitz continuous with Lipschitz constant \( L_H\) , as defined in Definition 1, in an open neighborhood \( \Omega\) of \( x^*\) . Let \( \{\widehat{x}_{i}\}\) be the sequence generated by mixed precision Newton’s method under Error Model 1. If, at some iteration \( i\) ,
(5)
then
(6)
where
(7)
(8)
(9)
Moreover, if Equation 5 holds for all \( i\) and
(10)
for \( \theta_{\max} \in [0, 1)\) , then there exists \( \rho > 0\) such that if \( \widehat{x}_0 \in B_{\rho}(x^*)\) , the sequence \( \{\widehat{x}_{i}\}\) is well defined and satisfies Equation 6 for all \( i\) until
(11)
Proof
See Section A.1.
Theorem 1 shows that the error decreases until the first iteration for which it becomes smaller than \( \gamma_i/(1-\theta_{\max})\) . Assuming \( \theta_{\max}\) is sufficiently less than \( 1\) , which is satisfied when \( L_H\lVert H(\widehat{x}_i)^{-1} \rVert \lVert \widehat{x}_i - x^* \rVert \) and \( \nu_i\) are safely below one, this means that the quality of the possible approximation of \( x^*\) is mainly determined by \( \gamma_i\) , which we call the limiting accuracy. To first order,
(12)
(13)
Hence, the limiting accuracy depends, to first order, on the gradient error \( \epsilon^g_i\) and on the working precision \( \epsilon_i\) , but not on the Hessian error \( \epsilon^H_i\) . This shows that we may tolerate some errors in forming the Hessian and solving the associated linear system without impacting the final solution accuracy. Moreover, since \( \epsilon^g_i\) depends on the error in evaluating the gradient, which may be large, and since \( \epsilon^g_i\) is multiplied by \( \lVert H(\widehat{x}_i)^{-1} \rVert \) , Equation 13 suggests that the error incurred in the gradient evaluation is the one that impacts the limiting accuracy the most. Therefore, if we wish to obtain the highest possible solution quality, the gradient should be evaluated as accurately as possible. Finally, note that the behavior of \( \gamma_i\) over the iterations is not necessarily monotone, especially due to \( \epsilon^g_i\) , which can change over time depending on the conditioning of the gradient at the current iterate \( \widehat{x}_i\) .
Let us now turn to the convergence rate of the method. Inequality Equation 6 shows that it contains both a quadratic term \( \alpha_i\) (see Equation 7) and a linear term \( \beta_i\) (see Equation 8). To first order, these terms behave as
(14)
First, in exact arithmetic, and in absence of approximations on the gradient and the Hessian, we have \( \epsilon_i = \epsilon^g_i = \epsilon^H_i = \nu_i = 0\) for all \( i\) , and so \( \beta_i = \gamma_i=0\) . Using the same argument as in [1, Thm. 3.5], we have that, if \( \widehat{x}_i\) is close enough to the solution, \( \lVert H(\widehat{x}_i)^{-1} \rVert \leq 2\lVert H(x^*)^{-1} \rVert \) , recovering the standard Newton’s quadratic convergence, under the same assumptions. Moreover, even in presence of errors, Newton’s method may still converge quadratically if the linear term \( \beta_i\) is small compared with the quadratic one \( \alpha_i\) ; conversely, if \( \beta_i\) becomes dominant, then this will deteriorate the convergence to a linear rate, or may even prevent convergence if \( \beta_i \ge 1\) . Importantly, \( \beta_i\) depends on the Hessian error, which shows that forming and solving the Hessian system approximately will affect the convergence rate of the method. Finally, in order to preserve convergence, the level of error \( \epsilon^H_i\) introduced in the Hessian should be chosen to be inversely proportional to the condition number of the Hessian \( \kappa(H(\widehat{x}_i))\) : the more ill-conditioned the Hessian, the smaller the tolerated error.
2.1 Analogous result for the gradient norm convergence
We conclude this section by proving an analogous result for the gradient norm.
Theorem 2
Let \( f:\mathbb{R}^n\rightarrow\mathbb{R}\) be twice continuously differentiable. Let \( x^* \in \mathbb{R}^n\) be a minimizer for \( f\) , and assume that \( H(x^*)\) is nonsingular. Let \( H\) be Lipschitz continuous with Lipschitz constant \( L_H\) , as defined in Definition 1, in an open neighborhood \( \Omega\) of \( x^*\) . Let \( \{\widehat{x}_{i}\}\) be the sequence generated by mixed precision Newton’s method under Error Model 1. If Equation 5 holds for some iteration \( i\) , then
(15)
where
(16)
(17)
where \( \mu_i := L_H\lVert H(\widehat{x}_i)^{-1} \rVert \lVert \widehat{x}_i - x^* \rVert \) , \( \tau_i := L_H \lVert H(\widehat{x}_i)^{-1} \rVert \gamma_i\) , and \( \gamma_i\) is defined in Equation 9.
Moreover, if Equation 5 holds for all \( i\) and \( \phi_i < \phi_{\max} \in [0, 1)\) , then there exists \( \rho > 0\) such that if \( \widehat{x}_0 \in B_{\rho}(x^*)\) , the sequence \( \{g(\widehat{x}_{i})\}\) is well defined and satisfies Equation 15 for all \( i\) until
(18)
Proof
See Section A.2.
As for Theorem 1, the gradient norm decreases until the first iteration for which it becomes smaller than \( \psi_i /(1-\phi_{\max})\) . To first order, \( \phi_i\) and \( \psi_i\) behave as
(19)
(20)
The assumptions of this second theorem are stronger than the ones of Theorem 1, since they also require a bound on the term \( \gamma_i\) . This assumption is necessary to ensure that the absolute error on the solution decreases enough to show a decrease in the gradient norm too. Moreover, in classical Newton’s method, the gradient norm decreases quadratically. In this theorem the term responsible for the quadratic convergence is hidden in \( \phi_i\) , which contains a factor \( \lVert \widehat{x}_i - x^* \rVert \) , which multiplies \( \|g(\widehat{x}_i)\|\) .
Additionally, we see that the error \( \epsilon^g_i\) on the gradient evaluation impacts not only \( \psi_i\) , but also the convergence rate \( \phi_i\) of the gradient norm, through the term \( \tau_i\) . This contrasts with Theorem 1, in which \( \epsilon^g_i\) does not affect the convergence rate of the error on the solution \( \lVert \widehat{x}_i - x^* \rVert \) , but only its limiting accuracy.
3 Newton’s approximations
Error Model 1 is quite general and encompasses different sources of errors. First of all, it covers the rounding errors arising from the use of finite precision floating-point arithmetic; we consider these errors in Section 3.1. However, it can account for more general sources of errors, typically arising from approximations introduced to make classical Newton’s method more suitable for large-scale problems. Specifically, we consider inexact Newton (Section 3.2) and Gauss–Newton (Section 3.3) methods; for both of these variants, we not only discuss the errors introduced by their approximations, but also their interplay with rounding errors, that is, we consider mixed precision inexact Newton and Gauss–Newton in floating-point arithmetic.
3.1 Floating-point Newton
In any floating-point arithmetic compliant with the IEEE 754 standard [30], the elementary operations satisfy the following model [31, sect. 2.2]:
(21)
where \( u\) is the unit roundoff of the precision used and \( \rm{\operatorname{f.2ptl}}(\cdot)\) represents the results computed in floating-point arithmetic. Hence, floating-point arithmetic introduces relative errors proportional to the unit roundoff of the arithmetic.
In our context, we consider a mixed precision approach for Newton’s method, outlined in Algorithm 1, which uses three floating-point arithmetics with different unit roundoffs:
- \( u\) is the unit roundoff of the working precision, used for storing and updating the iterates;
- \( u_g\) is the unit roundoff of the arithmetic used for evaluating the gradient;
- \( u_H\) is the unit roundoff of the arithmetic used for forming and solving the Hessian system.
Algorithm
Input: initial guess \( x_0\) , Hessian \( H\) , gradient \( g\)
Output: an approximation \( x_{i+1}\) to the minimizer \( x^*\)
Algorithm 1 Mixed precision Newton
Let us now discuss how the unit roundoffs \( u\) , \( u_H\) , \( u_g\) relate to the corresponding error terms \( \epsilon_i\) , \( \epsilon^H_i\) , \( \epsilon^g_i\) in Error Model 1. By Equation 21, we readily have \( \epsilon_i = u\) . For the Hessian system, a backward stable solver will deliver an error \( \epsilon^H_i\) of order \( u_H\) ; for example, for a direct solver based on Cholesky factorization, \( \epsilon^H_i = O(n^2)u_H\) [31, eq. (10.7)], where the dimensional constant in \( O(n^2)\) is known to be pessimistic [32]. Finally, \( \epsilon^g_i\) will be a (potentially large) multiple of \( u_g\) , but its precise value is very much dependent on the expression of the gradient, the point at which it is evaluated, and the method of evaluation. We will discuss in Section 4.1 how to measure these errors in practice.
Given the discussion in the previous section, the setting of interest is \( u_g \le u \le u_H\) : we consider the use of a potentially higher precision to evaluate the gradient (to improve the limiting accuracy) and of a potentially lower precision to form and solve the Hessian system (to reduce the computational cost, while preserving high limiting accuracy, at the price of potentially deteriorating the convergence rate).
3.2 Inexact Newton
Newton’s method requires solving a linear system of the form \( H(\widehat{x}_i)\widehat{d}_i = -g(\widehat{x}_i)\) at each iteration. For large-scale problems, solving this system exactly with a direct method can be quite expensive. Instead, inexact Newton’s methods [1, chap. 7.1] solve this system approximately by an iterative solver such as the conjugate gradient (CG) method [33]. The standard criterion to stop such an iterative method is to stop whenever the computed \( \widehat{d}_i\) satisfies, for given tolerances \( 0\leq \eta_i<1\) ,
(22)
Our framework can be applied to inexact Newton’s method using the Rigal–Gaches theorem [31, Thm. 7.1], which shows that the following two statements are equivalent [31, eq. (17.33b)]:
- \( \exists E_i^H : (H(\widehat{x}_i)+E_i^H)\widehat{d}_i = -g(\widehat{x}_i), \lVert E_i^H \rVert \le \epsilon^H_i\|H(\widehat{x}_i)\|\) ;
- \( \lVert H(\widehat{x}_i)\widehat{d}_i + g(\widehat{x}_i) \rVert \le \epsilon^H_i\lVert H(\widehat{x}_i) \rVert \lVert \widehat{d}_i \rVert \) .
Hence inexact Newton’s method in exact arithmetic and with no approximations on the gradient satisfies Error Model 1 with
(23)
where \( \zeta_i = \lVert H(\widehat{x}_i) \rVert \lVert \widehat{d}_i \rVert /\lVert g(\widehat{x}_i) \rVert \) satisfies
(24)
In the regime where, due to the inexactness in the linear system solution, the term \( \beta_i\) dominates over the term \( \alpha_i\) , Theorem 1 proves a convergence rate of the form
(25)
This can be compared with the standard convergence theory of inexact Newton’s method, as in [4], [1, chap. 7.1], where a convergence rate is proved in the energy norm induced by \( H(x^*)^2\) with \( \beta_i = \eta_\mathrm{max}\) for any \( \eta_\mathrm{max}\) such that \( \forall i, \eta_i < \eta_\mathrm{max}\) , which implies \( \lVert \widehat{x}_{i+1} - x^* \rVert \le \eta_{\max} \kappa(H(\widehat{x}_i))\lVert \widehat{x}_i -x^* \rVert \) .
Exploiting the generality of our framework, we can combine inexact Newton’s method with mixed precision floating-point arithmetic. This amounts to modifying Algorithm 1 so that the Hessian system on line 3 is solved by an iterative solver with stopping tolerance \( \eta_i\) . Then Error Model 1 is satisfied with \( \epsilon^H_i \leq c_\mathrm{solver} u_H + \eta_i/\zeta_i\) , where the constant \( c_\mathrm{solver}\) depends on the specifics of the iterative solver used; for CG, \( c_\mathrm{solver}=O(nk^2)\) , where \( k\) is the number of iterations [34]. This gives some indication on how to choose \( \eta_i\) based on the unit roundoff \( u_H\) of the arithmetic used to solve the Hessian system (or vice versa): in order to equilibrate both sources of inexactness, we should set \( \eta_i \approx \zeta_i c_\mathrm{solver} u_H\) . We will illustrate this rule of thumb experimentally in Section 4.3.
3.3 Gauss–Newton
Our error model Error Model 1 can potentially encompass quasi-Newton methods. Indeed, the term \( E_i^H\) can be used to represent the errors arising from the approximation of the Hessian matrix. In this section, we present as an example the Gauss–Newton method, a quasi-Newton method specifically designed for nonlinear least-squares problems, that is, problems of the form
where \( r: \mathbb{R}^n \to \mathbb{R}^m\) is referred to as the residual. Exploiting the problem’s structure, the method builds a Hessian approximation only using first-order information of \( r\) [19]. In fact, the gradient of \( f\) can be expressed as \( g(x) = J(x)^Tr(x)\) , with \( J: \mathbb{R}^n \to \mathbb{R}^{m \times n}\) the Jacobian matrix of \( r\) , and the second-order derivatives as
(26)
where \( r_i(x)\) is the \( i\) th component of the residual \( r(x)\) . The Gauss–Newton method approximates the Hessian \( H(x)\) with \( J(x)^TJ(x)\) , thus neglecting the term \( S(x)\) , which contains the second-order derivatives of the residual. Its iterations thus read:
(27)
The convergence of the Gauss–Newton method depends on the relative importance of the discarded term \( S(\widehat{x}_i)\) with respect to \( J(\widehat{x}_i)^TJ(\widehat{x}_i)\) . If this term is negligible, we can recover the fast quadratic convergence of Newton’s method, but, if it is not, the convergence can degrade to a linear one or the method may not converge at all. This result, for exact arithmetic, can for instance be found in [19, Thm. 10.2.1].
It is possible to apply Error Model 1 to the Gauss–Newton method by interpreting the discarded term \( S(\widehat{x}_i)\) as a perturbation matrix \( E_i^H\) . Defining \( E_i^H = -S(\widehat{x}_i)\) , we have \( \epsilon^H_i = \lVert S(\widehat{x}_i) \rVert /\lVert H(\widehat{x}_i) \rVert \) and Theorem 1 applies with \( \alpha_i\approx L_H \|H(\widehat{x}_i)^{-1}\|\) and \( \beta_i \approx \epsilon^H_i\kappa(H(\widehat{x}_i)) = \lVert S(\widehat{x}_i) \rVert \lVert H(\widehat{x}_i)^{-1} \rVert \) . We can compare this convergence result to the one in [19, Thm. 10.2.1]. The latter is based on the key assumption that there exists a \( \sigma\) such that
(28)
where \( \lambda_\mathrm{min}(\cdot)\) denotes the smallest eigenvalue of a matrix. Since the first inequality is certainly satisfied with \( \sigma=\lVert S(x^*) \rVert \) , a sufficient condition for this assumption to hold is \( \lVert S(x^*) \rVert < \lambda_\mathrm{min}(J(x^*)^TJ(x^*))\) .
We can relate this assumption to ours as follows, using the 2-norm, denoted as \( \|\cdot\|_2\) . Our assumptions in Error Model 1 and Theorem 1 require
and thus
By Weyl’s inequality [35, Thm. 4.3.1], for all \( x\) such that \( S(x)\) is positive semidefinite, it holds
Thus if \( S(x^*)\) is positive semidefinite and \( \|S(x^*)\|_2<\lambda_{\min}(J(x^*)^TJ(x^*))\) , then also \( \|S(x^*)\|_2< \lambda_{\min}(H(x^*))\) .
Concerning the convergence rate, in [19, Thm. 10.2.1] we have for all \( i\)
assuming \( J(x)\) to be \( L_J\) -Lipschitz, as defined in Definition 1. Assuming again \( S(x^*)\) to be positive semidefinite, and that it is \( L_S\) -Lipschitz, our rates satisfy
We thus obtain a convergence rate bound similar to [19] when \( \widehat{x}_i\) approaches \( x^*\) .
Once again, note that our framework accounts for both the Gauss–Newton approximation of the Hessian and any other source of inexactness, in particular the use of mixed precision floating-point arithmetic. We can obtain such a mixed precision Gauss–Newton method by modifying Algorithm 1 by replacing the Hessian matrix \( H(\widehat{x}_i)\) on line 3 by \( J(\widehat{x}_i)^TJ(\widehat{x}_i)\) . Then Error Model 1 is satisfied with \( \epsilon^H_i = c_\mathrm{solver} u_H + \lVert S(\widehat{x}_i) \rVert /\lVert H(\widehat{x}_i) \rVert \) , where \( c_\mathrm{solver}\) is a constant depending on the method used for solving the Gauss–Newton linear system. This shows that Gauss–Newton can be quite resilient to the use of low precision for the Hessian, since the \( u_H\) term will only impact the convergence rate if it is dominant compared to \( \lVert S(\widehat{x}_i) \rVert \) . Conversely, Gauss–Newton may perform just as well as Newton when using low precision if the rounding errors dominate. We will illustrate this observation experimentally in Section 4.3.
4 Numerical Experiments
In this section, we present numerical experiments to validate our bounds and illustrate the conclusions that we can draw from them. The code used to perform these experiments is available online . After describing our experimental setting in Section 4.1, we focus first on standard Newton’s method in Section 4.2, and then consider the inexact Newton and Gauss–Newton variants in Section 4.3 and Section 4.4, respectively.
4.1 Experimental setting
We now describe the setup used throughout this section. We first outline the implementation details, including the floating-point arithmetics and precision combinations considered, and then introduce the test problems used to validate our theoretical bounds.
4.1.1 Implementation
The algorithms have been implemented using python, leveraging the NumPy library to use different floating-point arithmetics. We consider the standard double (fp64) and single (fp32) precisions, as well as bfloat16 (abbreviated bf16 in the charts) arithmetic, simulated via the library ml_dtypes . We also use extended precision in order to compute the reference solution \( x^*\) ; we use np.float128, though that does not provide quadruple precision as the name suggests, but rather an 80-bit “long double” precision. We use the notation \( u\equiv\) fp32 to indicate that the precision with unit roundoff \( u\) has been set to fp32 (for example).
We will consider different precision combinations, denoted as tuples \( (\cdot, \cdot, \cdot)\) , indicating, respectively, the precisions with unit roundoff \( u_g\) , \( u\) , and \( u_H\) in Algorithm 1.
For almost all the experiments, we plot the convergence history of both the relative error \( \lVert \widehat{x}_i - x^* \rVert /\lVert x^* \rVert \) and the gradient norm \( \lVert g(\widehat{x}_i) \rVert \) (always on the left and right part of the figures, respectively).
We also plot the bound on the relative error and the gradient norm, as derived in Theorem 1 and Theorem 2 respectively. We use slightly more transparent, dashed curves for these bounds and only plot them when the required assumptions are satisfied. In order to compute these quantities, we must compute both \( \epsilon^g_i\) and \( \epsilon^H_i\) . For \( \epsilon^g_i\) , we compute \( g(\widehat{x}_i)\) in precision np.float128, and set \( \epsilon^g_i\) to the norm of the difference between \( g(\widehat{x}_i)\) and the gradient computed in the chosen precision with unit roundoff \( u_g\) . The error on the Hessian \( \epsilon^H_i\) is computed using the Rigal–Gaches formula [31, Thm. 7.1] for the backward error, evaluated in precision np.float128.
4.1.2 Test problems
For most experiments, we will use the following two test problems. Note that, in this section only, \( x_k\) refers to the \( k\) th component of the vector \( x\) , not the \( k\) th iterate of the algorithm.
-
ENGVAL1 : a standard minimization problem from the CUTEst dataset [36]; the function to be minimized is
\[ \begin{equation} f(x) = 3+\sum_{k=0}^{n-2} (x_k^2+x_{k+1}^2)^2 -4x_k, \end{equation} \](29)
where \( x \in \mathbb{R}^n\) and \( n=100\) in our setting.
-
SINREG: a least-squares regression problem on some syntethic vectors of datapoints \( z, y \in \mathbb{R}^m\) , where the function to be minimized is:
\[ \begin{equation} f(x) = \frac{1}{2}\lVert F(x, z) - y \rVert ^2 = \frac{1}{2}\lVert r(x) \rVert ^2, \end{equation} \](30)
and \( F : \mathbb{R}^{n} \times \mathbb{R}^m \to \mathbb{R}^m\) is a model parametrized by \( x \in \mathbb{R}^n\) . In our experiments, we use \( m=50\) and \( F\) defined as
\[ \begin{equation} F(x, z) = x_0 z + \sum_{k=1}^{\left\lceil\frac{n-1}{2}\right\rceil} x_{2k-1} z^{k+1} + \sum_{k=1}^{\left\lfloor\frac{n-1}{2}\right\rfloor}x_{2k} \sin(x_{2k} z), \end{equation} \](31)
where the power in \( z^{k+1}\) is applied componentwise. This test problem has a diagonal Hessian whose condition number can be easily controlled. In each set of experiments, we define a reference solution \( \bar{x}^*\) and \( y\) is computed as \( F(\bar{x}^*,z) + \xi\) , where \( \xi \in \mathbb{R}^m\) is noise randomly sampled from a uniform \( [-\delta, \delta]\) distribution. The value of \( \delta\) is \( 10^{-1}\) , unless otherwise specified.
For both problems, since no solution \( x^*\) is available beforehand, \( x^*\) is set to the solution found by standard Newton’s method in uniform extended precision (that is, with all operations performed in np.float128) after at most \( 500\) iterations.
We will also perform experiments on a wider range of CUTEst problems in Section 4.2.3.
4.2 Standard Newton in mixed precision
In this section we consider standard Newton’s method in mixed precision floating-point arithmetic, as described in Algorithm 1. In this first set of experiments, we do not want the Hessian system solution to be affected by any inexactness other than floating-point errors; thus the linear systems are solved directly using LU factorization. Hence, in this setting we have \( \epsilon^H_i \approx u_H\) .
4.2.1 Floating-point errors
We consider the SINREG problem as defined in Equation 30–Equation 31, because it allows us to study the impact of the condition number of the Hessian. Indeed, choosing \( n=4\) and \( \bar{x}^* = (1.0123, 2.01234, 1.01231, 2.01234)\) , we have \( \kappa(H(x^*)) \approx 2 \times 10^1\) .
Figure 1 shows the convergence of mixed precision Newton for this problem. The assumptions of Theorem 1 and Theorem 2 are easily satisfied in this case, even for low precisions. Therefore, the method converges for all precision combinations and respects the theoretical bounds. Moreover, these theoretical bounds (transparent dashed curves) are quite descriptive as they capture well the actual convergence behavior (solid curves), for all precision combinations.
Comparing the purple and blue curves (for which \( u_H\equiv\) fp32) with the orange and green ones (for which \( u_H\equiv\) bfloat16), we can see that \( u_H\) does not impact the limiting accuracy and only slightly impacts the convergence rate, since the problem is well conditioned.
Moreover, the purple and blue curves mostly overlap and, similarly, the orange and green curves also mostly overlap. This shows that both the convergence rate and limiting accuracy are essentially unchanged whether we set \( u_g\equiv\) fp32 or \( u_g\equiv\) fp64. In particular, the attainable relative error is mainly determined by the working precision \( u\) in these charts. Further experiments, that we omit for brevity, however confirm that it is actually regulated by \( \max\left(\epsilon^g_i \lVert H(x^*)^{-1} \rVert , u \lVert x^* \rVert \right)\) , as predicted by the theory through \( \gamma_i\) in Equation 13. In this case, the gradient is computed analytically, its evaluation is not significantly affected by propagation of rounding errors, and the Hessian matrix is well conditioned, and so the contribution of \( u\) dominates. In this context, then, we can just use two different precisions, computing the gradient in the same precision as the target working precision (\( u_g=u\) ) and using a lower precision for the Hessian (\( u_H \gg u\) ).
This behavior changes when considering different problems where the Hessian is ill conditioned or the gradient is strongly approximated, as we will see in subsequent examples. For instance, setting \( n = 8\) in Equation 31 and taking the new solution point to be the concatenation of two copies of \( \bar{x}^*\) , we now have \( \kappa(H(x^*)) \approx 4 \times 10^5\) . Figure 4 shows the convergence of Newton’s method on this new problem with the same four precisions sets as previously.
In this setting, the Hessian is so ill-conditioned that the assumptions of Theorem 1, Theorem 2 may not be satisfied if the precisions are too low or if the starting point is too far from the solution. We only plot the theoretical bounds (transparent curves) when these assumptions are satisfied. For example, when \( u_H\equiv\) bfloat16 (orange and green curves), the assumptions are never satisfied and in fact the method diverges. When \( u_H\equiv\) fp32 (purple and blue curves), the method does converge, although the assumptions are only satisfied for the purple curve (when \( u_g\equiv\) fp64) and only starting at iteration \( i=8\) , when the current iterate is sufficiently near the exact solution. This shows that in such extreme cases the theory can be too conservative and unable to guarantee convergence even though it is empirically observed.
Moreover, comparing the purple and blue curves reveals the impact of the conditioning on the relative error, which makes \( \epsilon^g_i\) the dominant error source when \( u_g = u\) . This problem instance thus illustrates that evaluating the gradient in a higher precision than the working precision (\( u_g \ll u\) ) can be beneficial to improve the limiting accuracy.
4.2.2 Finite differences
In this section we illustrate the effect of the errors coming from the gradient approximation and their interplay with rounding errors. Specifically, we consider forward finite differences [37]:
where \( e_k\) is the \( k\) th canonical basis vector, \( \mathbb{R}^+ \ni h \ll 1\) , and the computation is done in precision with unit roundoff \( u_g\) . The choice of \( h\) is of fundamental importance. The error on the gradient coming from forward finite differences is indeed of order \( \epsilon^g_i \approx h + u_g/h\) [1, eq. (8.5)]. The common approach is thus to choose \( h \approx \sqrt{u_g}\) , which balances approximation and finite precision errors and leads to a total error of order \( \epsilon^g_i=\sqrt{u_g}\) .
Figure 7 compares the convergence of mixed precision Newton on the ENGVAL1 problem (defined in Equation 29) with the gradient approximated with finite differences, for different steps \( h\) and for \( u_g\equiv\) fp64 or fp32. We also report the baseline convergence using the exact gradient to show the impact of these approximations.
The figure confirms that the choice \( h \approx \sqrt{u_g}\) leads to the best limiting accuracy, itself of order \( \sqrt{u_g}\) . Hence, for \( u_g = u\) (Figure 8), the limiting accuracy is significantly worse when the gradient is approximated than when using the analytical gradient, due to \( \epsilon^g_i\) dominating the \( \gamma_i\) term. The figure also shows that, as expected, this approximation does not impact the convergence rate of the method. Our framework can thus correctly predict the behavior of the method even in this case, showing that it can cover gradient errors that are not solely due to finite precision. Finally, this case provides another example where using higher precision for evaluating the gradient (\( u_g \ll u\) ) significantly improves the final accuracy.
4.2.3 Performance profile on CUTEst problems
We now evaluate the behavior of mixed precision Newton on a wider range of test problems coming from the optimization problems dataset CUTEst [36]. We consider 30 problems for which second-order derivatives are available in the python library, and whose Hessian is not too ill-conditioned (\( \kappa(H(x^*)) < 10^{15}\) ); their dimension ranges from 2 to 120. We compare the uniform fp32 method, in which all computations are performed in fp32, with three mixed precision variants using different precision sets for \( (u_g, u, u_H)\) : (fp64, fp32, fp32) uses higher precision for the gradient, (fp32, fp32, bfloat16) uses lower precision for the Hessian, and (fp64, fp32, bfloat16) does both at the same time, thereby using three different arithmetics. All the methods use the same working precision, \( u\equiv\) fp32. According to our theory (see Equation 17), the limiting accuracy on the gradient norm at a given iteration \( i\) will be, at best, of order \( \lVert H(\widehat{x}_i) \rVert \lVert \widehat{x}_i \rVert u\) . We therefore use this quantity as stopping criterion: all the methods run for at most \( 1000\) iterations and stop earlier if the gradient norm becomes smaller than \( \lVert H(\widehat{x}_i) \rVert \lVert \widehat{x}_i \rVert u\) .
In order to compare the behavior of the different precision sets on such a large number of problems, we use a performance profile [38], displayed in Figure 10. The performance profile is built considering the number of iterations as a performance metric. For each problem, the best performing method is the one that converges in the least number of iterations, and the others are compared to it by computing the ratio between their number of iterations and the best one. The performance profile then reports, for each method, the percentage of problems for which this ratio is below a certain threshold \( \tau\) .
The performance profile shows that using a higher precision for the gradient (fp64 instead of fp32) barely has any impact on either the number of successfully solved problems or the number of iterations. There is only one problem for which using \( u_g\equiv\) fp64 instead of \( u_g\equiv\) fp32 allows convergence, namely HAIRY, whose Hessian is not ill-conditioned (\( \kappa(H(x^*)) \approx 10^3\) ), but which shows a high error on the gradient when using lower precision (\( \epsilon^g_i \approx 0.1\) for many iterations when using fp32 everywhere). This suggests that, at least for these CUTEst problems, the evaluation of gradient in floating-point arithmetic is typically well behaved and does not require the use of a precision higher than the working precision.
On the other hand, the use of lower precision for the Hessian (bfloat16 instead of fp32) does have a more visible impact. The number of solved problems is reduced from 23 to 21: the two problems for which setting \( u_H\equiv\) bfloat16 prevents convergence are MEXHAT and PENALTY1. For these problems, the assumption on the Hessian conditioning of Theorem 1 is in fact not satisfied, which shows that our theory correctly detects the possible lack of convergence . Moreover, the average number of iterations is increased from about \( 10\) to about \( 18\) ; this increase is limited to at most a factor \( \tau=2\) for 60% of the problems. This suggests that many of these CUTEst problems present a sufficiently well-conditioned Hessian to be solved at a reduced computational cost by using lower precision than the working precision for the Hessian.
4.3 Mixed precision inexact Newton
We now consider the mixed precision inexact Newton method proposed in Section 3.2. We use the CG method to approximately solve the Hessian linear systems. We use Equation 22 as stopping criterion, with a fixed tolerance \( \eta_i=\eta\) across all nonlinear iterations \( i\) . The CG solver is initialized with the zero vector, and the maximum number of linear iterations is set to \( 100\) . We consider the ENGVAL1 problem defined in Equation 29, which leads to Hessian matrices with moderate condition numbers between \( 10^2\) and \( 10^3\) depending on the nonlinear iteration.
Recall that, according to our theory and as discussed in Section 3.2, we expect the convergence rate of the method to be mainly driven by \( \epsilon_i^H\approx u_H+\eta/\zeta_i\) , with \( \zeta_i\) defined in Equation 23. Moreover, for this problem, we have \( \zeta_i \approx 1\) for all \( i\) . Hence, we expect the convergence rate to be determined by the maximum of \( u_H\) and \( \eta\) .
We confirm this experimentally in Figure 12. We compare the number of nonlinear iterations required by mixed precision inexact Newton’s method, with precisions \( u_g=u\equiv\) fp64 and varying \( u_H\) and \( \eta\) . The figure shows that, as long as the stopping tolerance is large (\( \eta \ge10^{-3}\) ), there is no difference between using \( u_H\equiv\) fp64, fp32, or bfloat16, because the error is dominated by the solver’s inexactness. Thus, for such highly approximate solvers, we may safely use low precision arithmetic for the Hessian without impacting the convergence rate. For smaller tolerances \( \eta \le 10^{-4}\) , a difference between \( u_H\equiv\) bfloat16 and \( u_H\equiv\) fp32 or fp64 appears, with the former stagnating at 14–15 nonlinear iterations regardless of \( \eta\) , since the error is then dominated by \( u_H\) .
One could expect to see a similar difference appear between \( u_H\equiv\) fp64 and \( u_H\equiv\) fp32 when \( \eta\) becomes smaller than the fp32 unit roundoff, but this is not the case. The reason for this behavior is explained in Figure 15, Figure 14, which compare the relative error convergence with \( u_H\equiv\) fp64 and \( u_H\equiv\) fp32, and plots the three terms composing its theoretical bound separately. It shows that the dominant term is first \( \alpha_i\lVert \widehat{x}_i-x^* \rVert ^2\) , which decreases with \( i\) , until it becomes smaller than \( \beta_i\lVert \widehat{x}_i-x^* \rVert \) ; however, at that point, the relative error is already below the limiting accuracy \( \gamma_i\) . Hence, the unit roundoff \( u_H\) , which only appears in the term \( \beta_i\) , does not have any impact on the convergence rate in this situation. On the other hand, looking at the convergence with \( u_H\equiv\) bfloat16 in Figure 13, we see that the term \( \beta_i\lVert \widehat{x}_i-x^* \rVert \) becomes dominant earlier, and so delays convergence by a few iterations.
Finally, we want to validate the theoretical bounds when \( \zeta_i \gg 1\) for some \( i\) , in which case \( \eta_i\) should not be fixed across iterations, but rather depend on the current iterate through \( \zeta_i\) . In this experiment, we use the same test case as in Figure 4, for which \( \zeta_i\) ranges from \( 10^1\) to \( 10^3\) , with no particular pattern across iterations. In this setting, our bounds predict that the stopping condition derived in Section 3.2 (\( \eta_i = \zeta_i u_H\) ) should be as effective as \( \eta_i = u_H\) , that is, should guarantee the same convergence rate as if using a direct solver, while being less computationally expensive.
In Figure 16 we confirm that, for this problem, setting \( \eta_i = u_H\) for all \( i\) guarantees the same converge rate as using a direct solver. This is visible by comparing the purple curve with the one of Figure 4. Moreover, we see that setting \( \eta_i = \zeta_i u_H\) slightly increases the number of nonlinear iterations needed for the method to converge, but significantly reduces the number of CG steps overall, from 150 (across 8 nonlinear iterations) to 100 (across 9 nonlinear iterations).
4.4 Mixed precision Gauss–Newton
We conclude our experiments with the mixed precision Gauss–Newton method proposed in Section 3.3. We consider the least-squares problem SINREG defined in Equation 30–Equation 31, again using \( n=4\) and \( \bar{x}^* = (1.0123, 2.01234, 1.01231, 2.01234)\) , and with starting point \( x_0 = \bar{x}^* + 10^{-1}\) . We solve the Gauss–Newton system by LU factorization in precision with unit roundoff \( u_H\) . Thus, in this case, we expect the error term \( \epsilon^H_i\) to be driven by the maximum of \( u_H\) and the relative norm of the discarded term \( S(\widehat{x}_i)\) . We recall that the noise \( \xi\) added to the model outputs \( y_i\) is sampled from a uniform \( [-\delta, \delta]\) distribution. To make the norm of the discarded term \( S(\widehat{x}_i)\) vary, we test different noise sizes \( \delta\) .
In Figure 19, we compare the convergence in gradient norm of mixed precision Newton and Gauss–Newton; the relative error behaves similarly. We set \( u=u_g\equiv\) fp64 and we compare \( u_H\equiv\) fp64 (left plot) with \( u_H\equiv\) bfloat16 (right plot). The noise size is here set to \( \delta = 5\) , which leads to a discarded term of quite large relative norm \( \lVert S(\widehat{x}_i) \rVert /\lVert H(\widehat{x}_i) \rVert \) (starting around \( 3\times 10^{-2}\) for the first iterations and decreasing to \( 3\times 10^{-4}\) at convergence). As a result, in the left plot (with \( u_H\equiv\) fp64), the Hessian error \( \epsilon^H_i\) is dominated by this term and is much larger with Gauss–Newton than with standard Newton, so the former converges at a significantly slower rate than the latter. On the other hand, in the right plot (with \( u_H\equiv\) bfloat16), the Hessian error is now dominated by the rounding errors in bfloat16 arithmetic, and the two methods achieve a similar convergence rate. This illustrates that lower precision arithmetic can be safely used for the Hessian when it is already approximated by Gauss–Newton or, conversely, that Gauss–Newton can safely replace standard Newton if the Hessian precision is low.
In Figure 22, we plot the number of nonlinear iterations required by the two methods for different noise sizes \( \delta\) . When \( \delta\) is small, there is no significant difference between the two methods, even with \( u_H\equiv\) fp64, because the residual term discarded when using Gauss–Newton is small and does not dominate the Hessian error. As \( \delta\) increases, when using \( u_H\equiv\) fp64, the convergence of Gauss–Newton becomes slower than that of Newton. On the other hand, when using \( u_H\equiv\) bfloat16, the two methods converge in a similar number of iterations even for large \( \delta\) ; interestingly, Gauss–Newton may sometimes converge slightly faster than Newton.
5 Conclusion
We have presented a general framework for mixed precision Newton’s method for optimization, where the three main operations of the algorithm (gradient evaluation, Hessian system, and solution update) are affected by inexactness. Our error model is generic and can be applied to various sources of inexactness, including floating-point arithmetic, approximate linear solvers (inexact Newton), Hessian matrix approximations (Gauss–Newton), and any combination thereof. The main results of our analysis are Theorem 1, Theorem 2, which show how the convergence rate and attainable accuracy of Newton’s method are affected by these different sources of inexactness and provide guidelines for choosing the precisions of the different operations. For the inexact Newton and Gauss–Newton variants, we link our theory with known convergence results for exact arithmetic in the literature and discuss what changes in floating-point arithmetic, highlighting the interplay between these approximate variants and rounding errors. We have performed an extensive set of numerical experiments to validate the theoretical analysis and illustrate the interesting behavior of mixed precision Newton. Our results show that the empirically observed convergence rate and attainable accuracy match their theoretically predicted behavior when the assumptions underlying our theory are satisfied, and that the bounds are quite descriptive. Tests on a broad range of problems from the CUTEst dataset highlight the robustness and wide applicability of mixed precision Newton’s method, and suggest a significant potential for using lower precisions while only marginally sacrifing the accuracy and convergence rate.
This work opens the way to many promising perspectives:
- the error model and convergence theory could be extended to other quasi-Newton methods, such as BFGS [1, chap. 6.1], where Hessian inexactness depends on the gradient accuracy;
- the analysis could be extended to stochastic optimization [39], where inexactness in gradient and Hessian evaluations may also arise from data sampling;
- mixed precision Newton could be compared with perturbed first-order methods, such as the gradient descent method analysed in [28], to assess whether the use of Hessian information continues to provide advantages when operating under low or mixed precision arithmetic;
- a practical high-performance implementation of mixed precision Newton and its variants could be developed, to assess quantitatively the performance gains in memory, time, and energy that can be achieved by the use of lower precisions on modern hardware.
Acknowledgments
Funding for the PhD thesis of G. C. was provided by the Graduate+ MATHINFI Programme. This work was also partially supported by the Fondation Simone et Cino Del Duca and by projects managed by the French National Research Agency (ANR): France 2030 NumPEx Exa-MA (ANR-22-EXNU-0002), PEPR IA SHARP (ANR-23-PEIA-0008), MixHPC (ANR-23-CE46-0005-01), FPT-4 (ANR-24-CE46-7572), and MEPHISTO (ANR-24-CE23-7039).
A Proofs
A.1 Proof of Theorem 1
For the proof of this theorem, we will need the following two lemmas.
Lemma 1 (Lem. 4.1.12, [19])
Let \( g: \mathbb{R}^n \to \mathbb{R}^m\) be continuously differentiable in the open convex set \( \Omega \subseteq \mathbb{R}^n\) . Assume, for \( v \in \Omega\) , the matrix \( H\) of the first-order derivatives of \( g\) to be \( L_H\) -Lipschitz continuous at \( v\) in the neighborhood \( \Omega\) . Then, for any \( w \in \Omega\) , it holds
Lemma 2 (Thm. 3.1.4, [19])
Let \( M \in \mathbb{R}^{n \times n}\) . If \( \lVert M \rVert < 1\) , then \( (I+M)^{-1}\) exists and
Also, let \( A, B \in \mathbb{R}^{n \times n}\) and \( \lVert A^{-1}(B-A) \rVert < 1\) , then \( B\) is nonsingular and
Proof (Proof of Theorem 1)
Since \( H\) is continuous and \( H(x^*)\) is positive definite, there exists \( \rho > 0\) such that \( \forall x \in B_{\rho}(x^*)\) , \( H(x)\) is positive definite, and thus invertible.
First we prove that if \( \widehat{x}_i \in B_{\rho}(x^*)\) , then \( \widehat{x}_{i+1}\) is well defined, that is, that \( H(\widehat{x}_i) + E_i^H\) is invertible. To do so, we rewrite the perturbed Hessian:
Using Lemma 2 with \( M\equiv H(\widehat{x}_i)^{-1}E_i^H\) shows that \( I + H(\widehat{x}_i)^{-1}E_i^H\) is invertible since, by the definition of \( E_i^H\) in Equation 3.d and by Equation 5, we have
(32)
As the product of two invertible matrices, \( H(\widehat{x}_i) + E_i^H\) is then invertible, with inverse given by
(33)
Now, we will bound the error \( \lVert \widehat{x}_{i+1}-x^* \rVert \) . By Equation 3.b, we have
Taking norms yields
(34)
We now reformulate the first term of the right-hand side:
and using Lemma 2 and Equation 32 we can bound its norm by
(35)
Then, using Lemma 2 again with \( M\equiv H(\widehat{x}_i)^{-1}E_i^H\) and Equation 32, we obtain
(36)
Since \( g(x^*)=0\) , Lemma 1 with \( w\equiv \widehat{x}_i\) and \( v\equiv x^*\) yields
(37)
Moreover, by Equation 3.a, Equation 3.e, and Equation 36, we have
(38)
We then use Equation 37 to bound the norm of \( g(\widehat{x}_i)\) :
(39)
Finally, collecting inequalities from Equation 35 to Equation 39, and Equation 3.c to bound each of the terms in Equation 34 yields
(40)
where \( \alpha_i\) , \( \beta_i\) , and \( \gamma_i\) are defined in Equation 7, Equation 8, and Equation 9. We have thus proved that if \( \widehat{x}_i \in B_{\rho}(x^*)\) and \( \nu_i < 1\) , then Equation 40 holds.
Assume now that \( \widehat{x}_0\in B_{\rho}(x^*)\) . Let \( i_1\) be the first \( i\) for which the relative error does not decrease, that is, \( \lVert \widehat{x}_{i+1} - x^* \rVert \ge \lVert \widehat{x}_i - x^* \rVert \) . Then, for all \( i < i_1\) , since \( \lVert \widehat{x}_{i+1} - x^* \rVert < \lVert \widehat{x}_i - x^* \rVert \) , \( \widehat{x}_{i+1} \in B_{\rho}(x^*)\) remains in the ball. If, moreover, \( \nu_i<1\) for all \( i\le i_1\) , then we can inductively apply Equation 40 until \( i=i_1\) . Let us assume that there exists a \( \theta_{\max} \in [0, 1)\) such that \( \theta_i := \alpha_i\lVert \widehat{x}_i - x^* \rVert + \beta_i < \theta_{\max}\) for all \( i < i_1\) . Then, if \( \gamma_i \le (1-\theta_{\max})\lVert \widehat{x}_i - x^* \rVert \) , we have
which shows that at \( i=i_1\) , for the error to stop decreasing, we must necessarily have \( \gamma_i > (1-\theta_{\max})\lVert \widehat{x}_i - x^* \rVert \) . Therefore Equation 40 holds for all \( i\) until the first \( i_0 \le i_1\) for which \( \lVert \widehat{x}_i-x^* \rVert < \gamma_i /(1-\theta_{\max})\) .
A.2 Proof of Theorem 2
Proof (Proof of Theorem 2)
If Equation 5 holds at iteration \( i\) , by Theorem 1 \( \widehat{x}_{i+1}\) is well defined and the error on the iterates is bounded as in Equation 6. We want to relate the norm of the gradient at iteration \( i+1\) with that at the previous iteration.
To do that, let us define
Note that by Equation 3.b and Equation 3.a, we have
which yields, by Equation 3.c, Equation 3.d and Equation 3.e, that
(41)
By Error Model 1, using Lemma 2 with \( M\equiv H(\widehat{x}_i)^{-1}E_i^H\) and Equation 32, Equation 33, we have
(42)
which gives
(43)
By Lemma 1 we have \( \lVert \omega_i \rVert \leq \frac{L_H}{2}\lVert \widehat{x}_{i+1} - \widehat{x}_i \rVert ^2.\) By Equation 3.b and using Equation 3.e and Equation 42, it follows that
(44)
By the triangle inequality and Theorem 1 we also have
(45)
with \( \theta_i\) defined in Equation 10. By multiplying Equation 44 and Equation 45 term by term, we obtain:
(46)
where the penultimate and the last terms on the right-hand side of the inequality are obtained using \( \lVert H(\widehat{x}_i) \rVert \lVert H(\widehat{x}_i)^{-1} \rVert = \kappa(H(\widehat{x}_i)) \geq 1\) . Substituting Equation 43 and Equation 46 into Equation 41 we have
with
(47)
and
(48)
where \( \mu_i = L_H\lVert H(\widehat{x}_i)^{-1} \rVert \lVert \widehat{x}_i - x^* \rVert \) and \( \tau_i = L_H\gamma_i\lVert H(\widehat{x}_i)^{-1} \rVert \) , with \( \gamma_i\) defined in Theorem 1. Rearranging the terms in Equation 47 and Equation 48, we prove the first result of the theorem.
For the second part of the theorem, we proceed as for the proof of Theorem 1 in Section 7.1, to prove that the gradient norm decreases linearly at a rate at least \( \phi_i\) until it reaches \( \psi_i /(1-\phi_{\max})\) .
References
[1] Jorge Nocedal and Stephen J. Wright. Numerical optimization. Springer Series in Operations Research and Financial Engineering. Springer Nature, 2006. doi:10.1007/978-0-387-40065-5 .
[2] Stephen Boyd and Lieven Vandenberghe. Convex Optimization. Cambridge University Press, 2004.
[3] José Mario Martínez. Practical quasi-newton methods for solving nonlinear systems. J. Comput. Appl. Math., 124(1-2):97–121, December 2000. doi:10.1016/S0377-0427(00)00434-9 .
[4] Ron S. Dembo, Stanley C. Eisenstat, and Trond Steihaug. Inexact Newton methods. SIAM J. Numer. Anal., 19(2):400–408, 1982. doi:10.1137/0719025 .
[5] Cuda ptx isa. nvidia, may 2024. release 8.5.
[6] Nicholas J. Higham and Theo Mary. Mixed precision algorithms in numerical linear algebra. Acta Numer., 31:347–414, 2022. doi:10.1017/S0962492922000022 .
[7] Coralia Cartis, Nicholas IM Gould, and Philippe L Toint. Evaluation Complexity of Algorithms for Nonconvex Optimization: Theory, Computation and Perspectives. SIAM, 2022.
[8] Stefania Bellavia, Serge Gratton, and Elisa Riccietti. A Levenberg–Marquardt method for large nonlinear least-squares problems with dynamic accuracy in functions and gradients. Numer. Math., 140(3):791–825, November 2018. doi:10.1007/s00211-018-0977-z .
[9] A. S. Bandeira, K. Scheinberg, and L. N. Vicente. Convergence of trust-region methods based on probabilistic models. SIAM J. Optim., 24(3):1238–1264, 2014. doi:10.1137/130915984 .
[10] R. Chen, M. Menickelly, and K. Scheinberg. Stochastic optimization using a trust-region method and random models. Math. Program., 169(2):447–487, June 2018. doi:10.1007/s10107-017-1141-8 .
[11] Peng Xu, Fred Roosta, and Michael W. Mahoney. Newton-type methods for non-convex optimization under inexact Hessian information. Mathematical Programming, 184(1):35–70, November 2020. doi:10.1007/s10107-019-01405-z .
[12] Stefania Bellavia, Greta Malaspina, and Benedetta Morini. Inexact Newton methods with matrix approximation by sampling for nonlinear least-squares and systems, August 2023. URL: https://optimization-online.org/?p=24079 .
[13] Stefania Bellavia, Nataša Krejić, and Nataša Krklec Jerinkić. Subsampled inexact Newton methods for minimizing large sums of convex functions. IMA J. Numer. Anal., 40(4):2309–2341, 10 2020. doi:10.1093/imanum/drz027 .
[14] Farbod Roosta-Khorasani and Michael W. Mahoney. Sub-sampled Newton methods. Math. Program., 174(1-2):293–326, March 2019. doi:10.1007/s10107-018-1346-5 .
[15] El Houcine Bergou, Youssef Diouane, Vyacheslav Kungurtsev, and Clément W. Royer. A stochastic Levenberg–Marquardt method using random models with complexity results. SIAM-ASA J. Uncertain. Quantif., 10(1):507–536, 2022. doi:10.1137/20M1366253 .
[16] Z. Yao, Peng Xu, Fred Roosta, Stephen J. Wright, and Michael W. Mahoney. Inexact Newton-CG algorithms with complexity guarantees. IMA J. Numer. Anal., 2021. doi:10.1093/imanum/drac043 .
[17] Stefania Bellavia, Gianmarco Gurioli, and Benedetta Morini. Adaptive cubic regularization methods with dynamic inexact hessian information and applications to finite-sum minimization. IMA J. Numer. Anal., 41(1):764–799, 01 2021. doi:10.1093/imanum/drz076 .
[18] Zhewei Yao, Peng Xu, Fred Roosta, and Michael W Mahoney. Inexact nonconvex Newton-type methods. INFORMS J. Optim., 3(2):154–182, 2021. doi:10.1287/ijoo.2019.0043 .
[19] J. E. Dennis, Jr. and Robert B. Schnabel. Numerical methods for unconstrained optimization and nonlinear equations, volume 16 of Classics in Applied Mathematics. Society for Industrial and Applied Mathematics, Philadelphia, PA, 1996. Corrected reprint of the 1983 original. doi:10.1137/1.9781611971200 .
[20] P. Lancaster. Error analysis for the Newton-Raphson method. Numer. Math., 9(1):55–68, November 1966. doi:10.1007/BF02165230 .
[21] H. Woźniakowski. Numerical stability for solving nonlinear equations. Numer. Math., 27(4):373–390, December 1977. doi:10.1007/BF01399601 .
[22] T. J. Ypma. The effect of rounding errors on Newton-like methods. IMA J. Numer. Anal., 1983. doi:10.1093/imanum/3.1.109 .
[23] T. J. Ypma. Local convergence of inexact newton methods. SIAM J. Numer. Anal., 21(3):583–590, 1984. doi:10.1137/0721040 .
[24] Françoise Tisseur. Newton's method in floating point arithmetic and iterative refinement of generalized eigenvalue problems. SIAM J. Matrix Anal. Appl., 22(4):1038–1057, 2001. doi:10.1137/S0895479899359837 .
[25] C. T. Kelley. Newton's method in mixed precision. SIAM Rev., 64(1):191–211, 2022. doi:10.1137/20M1342902 .
[26] C. T. Kelley. Newton's method in three precisions. Pac. J. Optim., 20(3):461–474, 2024. URL: https://arxiv.org/abs/2307.16051 .
[27] Dominique Monnet and Dominique Orban. A multi-precision quadratic regularization method for unconstrained optimization with rounding error analysis. Comput. Optim. Appl., 91(3):997–1031, 2025. doi:10.1007/s10589-025-00676-x .
[28] Artem Vasin and Alexander Gasnikov. Lower and upper bounds of the convergence rate of gradient methods with composite noise in gradient, 2026. arXiv:2603.12376 .
[29] S. Gratton and Ph. L. Toint. A note on solving nonlinear optimization problems in variable precision. Comput. Optim. Appl., 76(3):917–933, 2020. doi:10.1007/s10589-020-00190-2 .
[30] IEEE standard for floating-point arithmetic. IEEE Std 754-2019 (Revision of IEEE 754-2008), pages 1–84, July 2019. doi:10.1109/IEEESTD.2019.8766229 .
[31] Nicholas J. Higham. Accuracy and Stability of Numerical Algorithms. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, second edition, 2002. doi:10.1137/1.9780898718027 .
[32] Nicholas J. Higham and Theo Mary. A new approach to probabilistic rounding error analysis. SIAM J. Sci. Comput., 41(5):A2815–A2835, 2019. doi:10.1137/18M1226312 .
[33] Magnus R. Hestenes and Eduard Stiefel. Methods of conjugate gradients for solving linear systems. J. Research Nat. Bur. Standards, 49:409–436, 1952. URL: https://nvlpubs.nist.gov/nistpubs/jres/049/jresv49n6p409_A1b.pdf .
[34] Thomas Bake, Erin Carson, and Yuxin Ma. Forward and backward error bounds for a mixed precision preconditioned conjugate gradient algorithm, 2025. arXiv:2510.11379 .
[35] Roger A. Horn and Charles R. Johnson. Matrix Analysis. Cambridge University Press, 1985.
[36] S. Gratton and Ph. L. Toint. S2MPJ and CUTEst optimization problems for Matlab, Python and Julia. Optim. Methods Softw., 40(4):871–903, 2025. doi:10.1080/10556788.2025.2490640 .
[37] Marco Rando, Cesare Molinari, Lorenzo Rosasco, and Silvia Villa. A structured tour of optimization with finite differences, 2025. arXiv:2505.19720 .
[38] Elizabeth D. Dolan and Jorge J. Moré. Benchmarking optimization software with performance profiles. Math. Program., 91(2):201–213, January 2002. doi:10.1007/s101070100263 .
[39] R. H. Byrd, S. L. Hansen, Jorge Nocedal, and Y. Singer. A stochastic quasi-Newton method for large-scale optimization. SIAM J. Optim., 26(2):1008–1031, 2016. doi:10.1137/140954362 .