Stefania Bellavia

Serge Gratton

Elisa Riccietti
email: elisa.riccietti@unifi.it

A Levenberg-Marquardt method for large nonlinear least-squares problems with dynamic accuracy in functions and gradients

Keywords: Levenberg-Marquardt method, Dynamic accuracy, Large-scale, Nonlinear least-squares

In this paper we consider large scale nonlinear least-squares problems for which function and gradient are evaluated with dynamic accuracy and propose a Levenberg-Marquardt method for solving such problems. More precisely, we consider the case in which the exact function to optimize is not available or its evaluation is computationally demanding, but approximations of it are available at any prescribed accuracy level. The proposed method relies on a control of the accuracy level, and imposes an improvement of function approximations when the accuracy is detected to be too low to proceed with the optimization process. We prove global and local convergence and complexity of our procedure and show encouraging numerical results on test problems arising in data assimilation and machine learning.

Introduction

Let us consider the following nonlinear least-squares problem

min x∈R n f (x) = 1 2 F (x) 2 (1.1)
where F : R n → R N with N ≥ n, continuously differentiable. Let J(x) ∈ R N ×n be the Jacobian matrix of F (x) and g(x) ∈ R n the gradient of f (x). Let x * be a solution of (1.1). We are interested in large scale nonlinear least-squares problems for which we do not have access to exact values for function F and for the Jacobian matrix or in problems for which an evaluation of f is computationally demanding and can be replaced by cheaper approximations. In both cases, to recover x * , we assume we can solely rely on approximations f δ to f . We are interested in the case in which the accuracy level of these approximations can be estimated and improved when judged to be too low to proceed successfully with the optimization process.

Typical problems that fit in this framework are those arising in the broad context of derivativefree optimization, where models of the objective function may result from a possibly random sampling procedure, cf. [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF][START_REF] Conn | Introduction to derivative-free optimization[END_REF]. An example is given by data-fitting problems like those arising in machine learning, cf. [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF][START_REF] Gratton | Preconditioning and globalizing conjugate gradients in dual space for quadratically penalized nonlinear-least squares problems[END_REF], in which a huge amount of data is available, so that N is really large and optimizing f is usually very expensive. Moreover, in this context there is often an approximate form of redundancy in the measurements, which means that a full evaluation of the function or the gradient may be unnecessary to make progress in solving (1.1), see [START_REF] Friedlander | Hybrid deterministic-stochastic methods for data fitting[END_REF]. This motivates the derivation of methods that approximate the function and/or the gradient and even the Hessian through a subsampling. This topic has been widely studied recently, see for example [7-9, 15, 20, 21, 23]. In these papers the data-fitting problem involves a sum, over a large number of measurements, of the misfits. In [START_REF] Bollapragada | Exact and inexact subsampled newton methods for optimization[END_REF][START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF][START_REF] Krejić | Nonmonotone line search methods with variable sample size[END_REF][START_REF] Roosta-Khorasani | Sub-sampled newton methods 1: globally convergent algorithms[END_REF] exact and inexact Newton methods and line-search methods based on approximations of the gradient and the Hessian obtained through subsampling are considered, in [START_REF] Krejić | Inexact restoration approach for minimization with inexact evaluation of the objective function[END_REF] the problem is reformulated in terms of constrained optimization and handled with an Inexact Restoration technique. In [START_REF] Friedlander | Hybrid deterministic-stochastic methods for data fitting[END_REF] the stochastic gradient method is applied to the approximated problems and conditions on the size of the subproblems are given to maintain the rate of convergence of the full gradient method. In [START_REF] Blanchet | Convergence rate analysis of a stochastic trust region method for nonconvex optimization[END_REF][START_REF] Chen | Stochastic optimization using a trust-region method and random models[END_REF] a variant of the traditional trustregion method for stochastic nonlinear optimization problems is studied. A theoretical analysis is carried out with the help of martingale theory and under the fully-linear model assumption.

Examples of nonlinear least-squares problems for which the exact gradient is not available and is replaced by a random model arise in variational modelling for meteorology, such as 3D-Var and 4D-Var which are the dominant data assimilation least-squares formulations used in numerical weather prediction centers worldwide, cf. [START_REF] Courtier | A strategy for operational implementation of 4d-var, using an incremental approach[END_REF][START_REF] Trémolet | Model-error estimation in 4d-var[END_REF]. In this context tri-dimensional fields are reconstructed combining the information arising from present and past physical observations of the state of the atmosphere with results from the mathematical model, cf. [START_REF] Gratton | Preconditioning and globalizing conjugate gradients in dual space for quadratically penalized nonlinear-least squares problems[END_REF][START_REF] Weaver | Three-and four-dimensional variational assimilation with a general circulation model of the tropical pacific ocean. part i: Formulation, internal diagnostics, and consistency checks[END_REF]. The result of this minimization procedure is the initial state of a dynamical system, which is then integrated forward in time to produce a weather forecast. This topic has been studied for example in [START_REF] Bergou | Levenberg-marquardt methods based on probabilistic gradient models and inexact subproblem solution, with application to data assimilation[END_REF][START_REF] Gratton | Preconditioning and globalizing conjugate gradients in dual space for quadratically penalized nonlinear-least squares problems[END_REF][START_REF] Gratton | Observation thinning in data assimilation computations[END_REF]. In [START_REF] Gratton | Preconditioning and globalizing conjugate gradients in dual space for quadratically penalized nonlinear-least squares problems[END_REF] conjugate-gradients methods for the solution of nonlinear least-squares problems regularized by a quadratic penalty term are investigated. In [START_REF] Gratton | Observation thinning in data assimilation computations[END_REF] an observation-thinning method for the efficient numerical solution of large-scale incremental four dimensional (4D-Var) data assimilation problems is proposed, built exploiting an adaptive hierarchy of the observations which are successively added based on a posteriori error estimate. In [START_REF] Bergou | Levenberg-marquardt methods based on probabilistic gradient models and inexact subproblem solution, with application to data assimilation[END_REF] a Levenberg-Marquardt approach is proposed to deal with random gradient and Jacobian models. It is assumed that an approximation to the gradient is provided but only accurate with a certain probability and the knowledge of the probability of the error between the exact and the approximated gradient is used in the update of the regularization parameter.

Problems in which inaccurate function values occur, and do not necessarily arise from a sampling procedure, are those where the objective function evaluation is the result of a computation whose accuracy can vary and must be specified in advance. For instance, the evaluation of the objective function may involve the solution of a nonlinear equation or an inversion process. These are performed through an iterative process that can be stopped when a certain accuracy level is reached. Such problems are considered for example in [START_REF] Conn | Trust region methods[END_REF]Section 10.6], where a trust-region approach is proposed to solve them, provided that a bound on the accuracy level is available.

In this paper we present a modification of the approach proposed in [START_REF] Bergou | Levenberg-marquardt methods based on probabilistic gradient models and inexact subproblem solution, with application to data assimilation[END_REF], to obtain a method able to take into account inaccuracy also in the objective function, while in [START_REF] Bergou | Levenberg-marquardt methods based on probabilistic gradient models and inexact subproblem solution, with application to data assimilation[END_REF] it is assumed to have at disposal the exact function values. In our procedure, following [START_REF] Conn | Trust region methods[END_REF] and deviating from [START_REF] Bergou | Levenberg-marquardt methods based on probabilistic gradient models and inexact subproblem solution, with application to data assimilation[END_REF], we replace the request made in [START_REF] Bergou | Levenberg-marquardt methods based on probabilistic gradient models and inexact subproblem solution, with application to data assimilation[END_REF] on first-order accuracy of the gradient up to a certain probability with a control on the accuracy level of the function values. Then, we propose a Levenberg-Marquardt approach that aims at finding a solution of problem (1.1) considering a sequence of approximations f δ k of known and increasing accuracy. Moreover, having in mind large scale problems, the linear algebra operations will be handled by an iterative Krylov solver and inexact solutions of the subproblems will be sought for.

Let us outline briefly our solution method. We start the optimization process with a given accuracy level δ = δ 0 . We rely during the iterative process on a control that allows us to judge whether the accuracy level is too low. In this case the accuracy is changed, making possible the use of more accurate approximations of function, gradient and Jacobian in further iterations. We assume that we have access to approximate function and gradient values at any accuracy level. In the following we define the approximation of f at iteration k as

f δ k (x) = 1 2 F δ k (x) 2 , (1.2)
where F δ k is the approximation of F at iteration k. Moreover we denote by J δ k (x) ∈ R n×n the approximation to the Jacobian matrix of F (x) and with

g δ k (x) = J δ k (x) T F δ k (x) ∈ R n the gradient approximation.
We assume that there exists δ k ≥ 0 such that at each iteration k:

max{ f δ k (x k + p LM k ) -f (x k + p LM k ) , |f δ k (x k ) -f (x k )|} ≤ δ k . (1.3)
As the quality of both the approximations of f and g at x k depends on the distance max{

F δ k (x k )- F (x k ) , J δ k (x k ) -J(x k ) }, as follows: |f δ k (x k ) -f (x k )| ≤ 1 2 F δ k (x k ) -F (x k ) N j=1 |F j (x k ) + (F δ k ) j (x k ))|, g(x k ) -g δ k (x k ) ≤ J δ k (x k ) -J(x k ) F (x k ) + J δ k (x k ) F δ k (x) -F (x) ,
we can also assume that there exists K ≥ 0 such that at each iteration k:

g δ k (x k ) -g(x k ) ≤ Kδ k . (1.4)
We will refer to δ k as the accuracy level and to f δ k , g δ k , J δ k as approximated function, gradient and Jacobian matrix. If not differently specified, when the term accuracy is used, we refer to accuracy of such approximations.

Our intention is to rely on less accurate (and hopefully cheaper quantities) whenever possible in earlier stages of the algorithm, increasing only gradually the demanded accuracy, so as to obtain a reduced computational time for the overall solution process. To this aim we build a non-decreasing sequence of regularization parameters. This is needed to prevent the sequence of solution approximations from being attracted by a solution of the problems with approximated objective functions, cf. [START_REF] Bellavia | On an adaptive regularization for ill-posed nonlinear systems and its trust-region implementation[END_REF][START_REF] Hanke | A regularizing levenberg-marquardt scheme, with applications to inverse groundwater filtration problems[END_REF][START_REF] Kaltenbacher | Iterative regularization methods for nonlinear ill-posed problems[END_REF], and to allow inexactness in function and gradient till the last stage of the procedure. The obtained approach is shown to be globally convergent to first-order critical points. Along the iterations, the step computed by our procedure tends to assume the direction of the approximated negative gradient, due to the choice of generating a non-decreasing bounded above sequence of regularization parameters. Then, eventually the method reduces to a perturbed steepest descent method with step-length and accuracy in the gradient inherited by the updating parameter and accuracy control strategies employed. Local convergence for such a perturbed steepest descent method is proved, too. We stress that overall the procedure benefits from the use of a genuine Levenberg-Marquardt method till the last stage of convergence, gaining a faster convergence rate compared to a pure steepest descent method. Moreover this can be gained at a modest cost, thanks to the use of Krylov methods to solve the arising linear systems. These methods take a reduced number of iterations when the regularization term is large: in the last stage of the procedure the cost per iteration is that of a first-order method.

We are not aware of Levenberg-Marquardt methods for both zero and non-zero residual nonlinear least-squares problems with approximated function and gradient, for which both local and global convergence is proved. Contributions on this topic are given by [START_REF] Bergou | Levenberg-marquardt methods based on probabilistic gradient models and inexact subproblem solution, with application to data assimilation[END_REF] where the inexactness is present only in the gradient and in the Jacobian and local convergence is not proved and by [START_REF] Bellavia | On an adaptive regularization for ill-posed nonlinear systems and its trust-region implementation[END_REF][START_REF] Bellavia | On non-stationary tikhonov procedures for ill-posed nonlinear least squares problems[END_REF] where the Jacobian is exact and only local convergence is considered.

Importantly enough, the method and the related theory also apply to the situation where the output space of F δ k has smaller dimension than that of F , i.e. F δ k : R n → R K k with K k ≤ N for some k. This is the case for example when approximations to f stem from a subsampling technique and F δ k is obtained by selecting randomly some components of F . We assume it is possible to get a better approximation to f by adding more observations to the considered subset, i.e. increasing K k . We denote accordingly by

J δ k (x) ∈ R K k ×n the Jacobian matrix of F δ k (x).
The paper is organized as follows. We describe the proposed Levenberg-Marquardt approach in Section 2, focusing on the strategy to control the accuracy level. We analyse also the asymptotic behaviour of the sequence of regularization parameters generated. In Section 3 global convergence of the procedure to first-order critical points is proved. In Section 4, we show that the step computed by our procedure tends to asymptotically assume the direction of the approximated negative gradient and motivated by this asymptotic result, we prove local convergence for the steepest descent method we reduce to. In Section 5 we provide a global complexity bound for the proposed method showing that it shares its complexity properties with the steepest descent and trust-region methods. Finally, in Section 6 we numerically illustrate the approach on two test problems arising in data assimilation (Section 6.1) and in machine learning (Section 6.2). We show that our procedure is able to handle the inaccuracy in function values and find a solution of problem (1.1), i.e. of the original problem with exact objective function. Moreover we show that when the exact function is available, but it is expensive to optimize, the use of our accuracy control strategy allows us to obtain large computational savings.

The method

A Levenberg-Marquardt approach is an iterative procedure that at each iteration computes a step as the solution of the following linearized least-squares subproblem:

min p∈R n m k (x k + p) = 1 2 F δ k (x k ) + J δ k (x k )p 2 + 1 2 λ k p 2 , (2.5)
where λ k > 0 is an appropriately chosen regularization parameter. As we deal with large scale problems, we do not solve (2.5) exactly, but we seek for an approximate solution. We say that p approximately minimizes m k if it achieves the sufficient Cauchy decrease, i.e. if it provides at least as much reduction in m k as that achieved by the so-called Cauchy point, which is the minimizer of m k along the negative gradient direction [START_REF] Wright | Numerical optimization[END_REF]:

m k (x k ) -m k (x k + p) ≥ θ 2 g δ k (x k ) 2 J δ k (x k ) 2 + λ k , θ > 0. (2.6)
Achieving the Cauchy decrease is a sufficient condition to get global convergence of the method, so one can rely on approximated solutions of problem (2.5), see [START_REF] Conn | Trust region methods[END_REF]. A solution of (2.5) can alternatively be found solving the optimality conditions

(J δ k (x k ) T J δ k (x k ) + λ k I)p = -g δ k (x k ). (2.7)
Then, we approximately solve (2.7), i.e. we compute a step p such that

(J δ k (x k ) T J δ k (x k ) + λ k I)p = -g δ k (x k ) + r k ,
where vector

r k = (J δ k (x k ) T J δ k (x k ) + λ k I)p + g δ k (x k
) is the residual of (2.7). If the norm of the residual vector is small enough, p achieves the Cauchy decrease, as stated in the next Lemma.

Lemma 1 The inexact Levenberg-Marquardt step p LM k computed as

(J δ k (x k ) T J δ k (x k ) + λ k I)p LM k = -g δ k (x k ) + r k (2.8) for a residual r k satisfying r k ≤ k g δ k , with 0 ≤ k ≤ θ 2 λ k J δ k (x k ) 2 + λ k , (2.9)
for some θ 2 ∈ 0, 1 2 , achieves the Cauchy decrease (2.6),

with θ = 2(1 -θ 2 ) ∈ [1, 2).
Proof For the proof see [START_REF] Bergou | Levenberg-marquardt methods based on probabilistic gradient models and inexact subproblem solution, with application to data assimilation[END_REF], Lemma 4.1.

Let us describe now in details the way the accuracy level is controlled along the iterations. In classical Levenberg-Marquardt methods, at each iteration, if the objective function is sufficiently decreased, the step is accepted and the iteration is considered successful. Otherwise the step is rejected, λ k is updated and problem (2.5) is solved again. Here we assume that the objective function is not evaluated exactly. In our approach, it is desirable to have an accuracy level high enough to ensure that the achieved decrease in function values, observed after a successful iteration, is not merely an effect of the inaccuracy in these values, but corresponds to a true decrease also in the exact objective function. In [START_REF] Conn | Trust region methods[END_REF]Section 10.6], it is proved that this is achieved if the accuracy level δ k is smaller than a multiple of the reduction in the model:

δ k ≤ η 0 [m k (x k ) -m k (x k + p LM k )],
with η 0 > 0. We will prove in (2.14) and numerically illustrate in Section 6, that for our approach

m k (x k ) -m k (x k + p LM k ) = O(λ k p LM k 2 ).
(2.10)

According to this and following [START_REF] Conn | Trust region methods[END_REF], we control the accuracy level asking that

δ k ≤ κ d λ α k p LM k 2 , (2.11)
for constants κ d > 0 and α ∈

Given x 0 , δ 0 , κ d ≥ 0, α ∈ 1 2 , 1 , β > 1, η 1 ∈ (0, 1), η 2 > 0, λmax ≥ λ 0 > 0, γ > 1. Compute f δ 0 (x 0 ) and set δ -1 = δ 0 . For k = 0,
k ≤ κ d λ α k p LM k 2 , compute f δ k (x k + p LM k
), and set δ k+1 = δ k . Else reduce δ k : δ k = δ k β and go back to 1.

Compute ρ

δ k k (p LM k ) = f δ k-1 (x k )-f δ k (x k +p LM k ) m k (x k )-m k (x k +p LM k ) . (a) If ρ δ k k (p LM k ) ≥ η 1 , then set x k+1 = x k + p LM k and λ k+1 = min{γλ k , λmax} if g δ k (x k ) < η 2 /λ k , λ k if g δ k (x k ) ≥ η 2 /λ k . (b) Otherwise set x k+1 = x k , λ k+1 = γλ k and δ k+1 = δ k-1 .
Our approach is sketched in Algorithm 2.1. At each iteration k a trial step p LM k is computed using the accuracy level of the previous successful iteration. The norm of the trial step is then used to check condition (2.11). In case it is not satisfied, the accuracy level is increased in the loop at steps 1-2 until (2.11) is met. On the other hand, when the condition is satisfied it is not necessary to estimate the accuracy again for next iteration. The value δ k obtained at the end of the loop is used to compute f δ k (x k + p LM k ). Then, the ratio between the actual and the predicted reduction

ρ δ k k (p LM k ) = f δ k-1 (x k ) -f δ k (x k + p LM k ) m k (x k ) -m k (x k + p LM k ) (2.12)
is computed to decide whether to accept the step or not. Practically, notice that if at iteration k the accuracy level is changed, i.e. δ k = δ k-1 , the function is not evaluated again in x k to compute ρ δ k k (p LM k ), and the ratio is evaluated computing the difference between f δ k-1 (x k ) (evaluated at the previous step), and the new computed value f δ k (x k + p LM k ). The step acceptance and the updating of the regularization parameter are based on this ratio. A successful step is taken if

ρ δ k k (p LM k ) ≥ η 1 .
In such case, deviating from classical Levenberg-Marquardt and following [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF][START_REF] Bergou | Levenberg-marquardt methods based on probabilistic gradient models and inexact subproblem solution, with application to data assimilation[END_REF], λ k is increased if the norm of the gradient model is of the order of the inverse of the regularization parameter (condition g δ k (x k ) < η 2 /λ k in Algorithm 2.1), otherwise it is left unchanged. In case the step is unsuccessful λ k is increased and reductions of δ k performed at steps 1-2 are not taken into account. That is, the subsequent iteration k + 1 is started with the same accuracy level of iteration k (see step 3b).

First we prove the well-definedness of Algorithm 2.1. Specifically in Lemma 2, we prove that the loop at steps 1-2 of Algorithm 2.1 terminates in a finite number of steps. To this aim we need the following assumption: Assumption 1 Let {x k } be the sequence generated by Algorithm 2.1. Then there exists a positive constant κ J such that, for all k ≥ 0 and all

x ∈ [x k , x k + p LM k ], J δ k (x) ≤ κ J .
In standard Levenberg-Marquardt method it is customary to assume the boundedness of the norm of the Jacobian matrix, cf. [START_REF] Conn | Trust region methods[END_REF]. Here, we need the boundedness assumption on the norm of the Jacobian's approximations. In the applications of our interest, that we will present in the numerical results section, this assumption is met and κ J ∼ 1.

Lemma 2 Let Assumption 1 hold and let p LM k be defined as in Lemma 1. If x k is not a stationary point of f , the loop at steps 1-2 of Algorithm 2.1 terminates in a finite number of steps.

Proof If δ k tends to zero, g δ k (x k ) tends to g(x k ) from (1.4). Equation (2.9) yields k ≤ √ θ 2 , and from (2.8) it follows

p LM k = (J δ k (x k ) T J δ k (x k ) + λ k I) -1 (-g δ k (x k ) + r k ) ≥ ≥ (1 -k ) g δ k (x k ) J δ k (x k ) 2 + λ k ≥ (1 - √ θ 2 ) g δ k (x k ) κ 2 J + λ k . (2.13)
Then, lim inf

δ k →0 p LM k ≥ (1 - √ θ 2 ) κ 2 J + λ k g(x k ) > 0
as g(x k ) = 0, so for δ k small enough (2.11) is satisfied.

As far as the sequence of regularization parameters is concerned, we notice that it is bounded from below, as λ min = λ 0 ≤ λ k for all k. Moreover un upper bound λ max is provided for successful iterations in step 3a, so that the procedure gives rise to a sequence of regularization parameters with different behaviour than the one generated in [START_REF] Bergou | Levenberg-marquardt methods based on probabilistic gradient models and inexact subproblem solution, with application to data assimilation[END_REF]. It is indeed possible to prove that the bound is reached and for k large enough λ k = λ max on the subsequence of successful iterations, while in [START_REF] Bergou | Levenberg-marquardt methods based on probabilistic gradient models and inexact subproblem solution, with application to data assimilation[END_REF] the sequence is shown to diverge. The result is proved in the following Lemma. Proof If the result is not true, there exists a bound 0 < B < λ max such that the number of times that λ k < B happens is infinite. Because of the way λ k is updated one must have an infinity of iterations for which λ k+1 = λ k , and for them one has

ρ δ k k (p LM k ) ≥ η 1 and g δ k (x k ) ≥ η 2 /B.
Thus, from Lemma 1 and relation (2.6)

f δ k-1 (x k ) -f δ k (x k + p LM k ) ≥ η 1 (m k (x k ) -m k (x k + p LM k )) ≥ η 1 2 θ g δ k (x k ) 2 J δ k (x k ) 2 + λ k ≥ η 1 2 θ κ 2 J + B η 2 B 2 .
Since f δ k is bounded below by zero and the sequence {f δ k (x k+1 )} is decreasing and hence convergent, the number of such iterations cannot be infinite, hence we derive a contradiction. Then, for an infinite number of iterations λ k+1 > λ k and for the subsequence of successful iterations it exists k large enough for which λ k = λ max for all k ≥ k.

From the updating rule of λ k in Algorithm 2.1, the generated sequence of regularization parameters is non-decreasing. This result enables us to prove (2.10) and to motivate condition (2.11). From the model definition and (2.8) it holds

m k (x k ) -m k (x k + p LM k ) = - 1 2 (p LM k ) T (J δ k (x k ) T J δ k (x k ) + λ k I)p LM k -(p LM k ) T g δ k (x k ) = 1 2 J δ k (x k )p LM k 2 + 1 2 λ k p LM k 2 -(p LM k ) T r k .
Considering that from (2.9) and (2.13)

(p LM k ) T r k ≤ k p LM k g δ k (x k ) ≤ √ θ 2 1 - √ θ 2 (κ 2 J + λ k ) p LM k 2 ,
and that parameters λ k form a non-decreasing sequence, we can conclude that

m k (x k ) -m k (x k + p LM k ) = O(λ k p LM k 2 ). (2.14)
In the following section, we will prove that the sequence generated by Algorithm 2.1 converges globally to a solution of (1.1).

Global convergence

In this section we prove the global convergence of the sequence generated by Algorithm 2.1. We assume to compute an inexact Levenberg-Marquardt step according to the following assumption:

Assumption 2 Let p LM k satisfy (J δ k (x k ) T J δ k (x k ) + λ k I)p LM k = -g δ k (x k ) + r k for a residual r k satisfying r k ≤ k g δ k , with 0 ≤ k ≤ min θ 1 λ α k , θ 2 λ k J δ k (x k ) 2 + λ k (3.15)
where

θ 1 > 0, θ 2 ∈ 0, 1 2 and α ∈ 1 2 , 1 is defined in (2.11).
As stated in Lemma 1, this step achieves the Cauchy decrease. Then, the idea is to solve the linear systems (2.7) with an iterative solver, stopping the iterative process as soon as requirement (3.15) on the residual of the linear equations is met. The first bound in (3.15) will be used in the convergence analysis. We point out that the allowed inexactness level in the solution of the linear systems decreases as λ k increases. However, an upper bound on λ k is enforced, so we do not expect extremely small values of 1 λ α k , especially if α = 0.5 is chosen. Also, we point out that if λ k is large, the matrix in the linear systems is close to a multiple of the identity matrix and fast convergence of the Krylov iterative solver is expected.

We now report a result relating the step length and the norm of the approximated gradient at each iteration, that is going to be useful in the following analysis.

Lemma 4 Let Assumptions 1 and 2 hold. Then

p LM k ≤ 2 g δ k (x k ) λ k . (3.16) Proof Taking into account that from Assumption 2 r k ≤ k g δ k ≤ g δ k , it follows p LM k = (J T δ k J δ k + λ k I) -1 (-g δ k (x k ) + r k ) ≤ 2 g δ k (x k ) λ k .
In the following Lemma we establish a relationship between the exact and the approximated gradient which holds for λ k large enough. This relation shows that it is possible to control the accuracy on the gradient through the regularization parameter. Specifically, large values of λ k yield a small relative error on g δ k (x k ) .

Lemma 5 Let Assumptions 1 and 2 hold. For λ k sufficiently large, i.e. for

λ k ≥ νλ * = ν(2 δ 0 κ d K) 2 2-α ν > 1, (3.17)
it exists c k ∈ (0, 1) such that the following relation between the exact and the approximated gradient holds:

g(x k ) (1 + c k ) ≤ g δ k (x k ) ≤ g(x k ) (1 -c k ) , with c k = 2 K√ δ 0 κ d λ 1-α/2 k = λ * λ k 1-α/2 . (3.18)
Proof From (1.4), (2.11) and (3.16) it follows

| g(x k ) -g δ k (x k ) | ≤ g(x k ) -g δ k (x k ) ≤ K δ 0 δ k ≤ K δ 0 κ d λ α k p LM k 2 = K δ 0 κ d λ α/2 k p LM k ≤ 2 K δ 0 κ d g δ k (x k ) λ 1-α/2 k = c k g δ k (x k )
where we have set

c k = 2 K√ δ0κ d λ 1-α/2 k . Then, g(x k ) -g δ (x k ) ≤ c k g δ k (x k ) , (3.19)
(1 -c k ) g δ k (x k ) ≤ g(x k ) ≤ (1 + c k ) g δ k (x k ) , (3.20)
and for λ k > λ * , the thesis follows.

From the updating rule of the accuracy level δ k in step 2 of Algorithm 2.1, if δ k-1 is the successful accuracy level at iteration k -1, the successful accuracy level at iteration k is

δ k = δ k-1 β n k (3.21)
where n k ≥ 0 counts the number of times the accuracy is increased (i.e. δ k is decreased) in the loop at steps 1-2, that is finite from Lemma 2. We can also prove that the sequence {β n k } is bounded from above. To this aim, we need the following Assumption, which is standard in Levenberg-Marquardt methods, cf. [START_REF] Conn | Trust region methods[END_REF]:

Assumption 3 Assume that function f has Lipschitz continuous gradient with corresponding constant L > 0: g(x) -g(y) ≤ L x -y for all x, y ∈ R n .

Lemma 6 Let Assumptions 1, 2, 3, hold and λ * be defined in (3.17). Then, if λ k ≥ νλ * for ν > 1, there exists a constant β > 0 such that β n k ≤ β.

Proof Let δ k-1 be the successful accuracy level at iteration k -1. Then, it holds

δ k-1 ≤ κ d λ α k-1 p LM k-1 2 .
If in (3.21) n k ≤ 1 there is nothing to prove, so let assume

n k > 1. If n k > 1 it holds βδ k > κ d λ α k p LM k 2 .
From the updating rule at step 3 of Algorithm 2.1 it follows

λ k-1 ≤ λ k ≤ γλ k-1 . (3.22)
Using the first inequality in (3.22) and (2.11) we get from (3.21) that

β n k -1 = δ k-1 βδ k < κ d λ α k-1 p LM k-1 2 κ d λ α k p LM k 2 ≤ p LM k-1 2 p LM k 2 .
Then, from Assumption 2, recalling (3.18) and that k < √ θ 2 from (3.15), we have

β n k -1 ≤ (J δ k-1 (x k-1 ) T J δ k-1 (x k-1 ) + λ k-1 I) -1 (-g δ k-1 (x k-1 ) + r k-1 ) 2 (J δ k (x k ) T J δ k (x k ) + λ k I) -1 (-g δ k (x k ) + r k ) 2 ≤ ≤ J δ k (x k ) T J δ k (x k ) + λ k I 2 (1 - √ θ 2 ) 2 g δ k (x k ) 2 4 g δ k-1 (x k-1 ) 2 λ 2 k-1 ≤ ≤ 4 (1 - √ θ 2 ) 2 κ 2 J + λ k λ k-1 2 (1 + c k ) 2 (1 -c k-1 ) 2 g(x k-1 ) 2 g(x k ) 2 . By (3.22) it follows λ k-1 ≥ λ k γ > λ * γ .
This and c k < 1 yield

β n k -1 ≤ 16 (1 - √ θ 2 ) 2 κ 2 J λ min + γ 2 1 1 -γ 1-α/2 2 g(x k-1 ) g(x k ) 2 .
Let us now consider the term g(x k-1 ) g(x k ) . By (3.16), (3.18) and the Lipschitz continuity of the gradient we get:

g(x k-1 ) g(x k ) ≤1 + g(x k-1 ) -g(x k ) g(x k ) ≤ 1 + L p LM k g(x k ) ≤1 + 2L g δ k (x k ) λ k g(x k ) ≤ 1 + 2L (1 -c k )λ k ≤1 + 2L (1 -ν α 2 -1 )λ min .
We can then conclude that sequence β n k is bounded from above by a constant for λ k sufficiently large.

The result in Lemma 6 can be employed in the following Lemma, to prove that for sufficiently large values of the parameter λ k the iterations are successful.

Lemma 7 Let Assumptions 1, 2 and 3 hold. Assume that

λ k > max{νλ * , λ} (3.23)
with λ * defined in (3.17) and

λ = ϕ 1 -η 1 1 1-α ϕ = κ 2 J /λ min + 1 θ 2θ 1 λ 2α-1 min + 2L λ α min + 4(3 + β)κ d + 8κ d ḡ λ min , (3.24)
with η 1 , β, θ 1 , θ, α, L defined respectively in Algorithm 2.1, Lemma 6, Assumption 2, (2.6), (2.11) and Assumption 3, and

ḡ = κ J 2f δ0 (x 0 ). If x k is not a critical point of f then ρ δ k k (p LM k ) ≥ η 1 . Proof We consider 1 - ρ δ k k (p LM k ) 2 = -(p LM k ) T (J δ k (x k ) T J δ k (x k ) + λ k I)p LM k -2(p LM k ) T g δ k (x k ) 2(m k (x k ) -m k (x k + p LM k )) (3.25) + 1 2 F δ k (x k + p LM k ) 2 -1 2 F δ k-1 (x k ) 2 2(m k (x k ) -m k (x k + p LM k )) . (3.26)
From the Taylor expansion of f and denoting with R the reminder, we obtain

f δ k (x k + p LM k ) = f δ k (x k ) + (p LM k ) T g δ k (x k ) + (f δ k (x k + p LM k ) -f (x k + p LM k )) + (f (x k ) -f δ k (x k )) + ((p LM k ) T g(x k ) -(p LM k ) T g δ k (x k )) + R.
Then, from conditions (1.3), (2.11) and the fact that if

λ k > λ * from Lemma 6 δ k-1 = β n k δ k ≤ βδ k , it follows f δ k (x k + p LM k ) -f δ k-1 (x k ) =f δ k (x k ) -f δ k-1 (x k ) + (p LM k ) T g δ k (x k ) + R ≤(1 + β)κ d λ α k p LM k 2 + (p LM k ) T g δ k (x k ) + R, where R = (f δ k (x k + p LM k ) -f (x k + p LM k )) + (f (x k ) -f δ k (x k )) + (p LM k ) T (g(x k ) -g δ k (x k )) + R.
Moreover, by (1.3), (1.4) and (2.11) we can conclude that

|R| ≤ (2 + p LM k )κ d λ α k + L 2 p LM k 2 .
Then, from Lemma 4, Assumption 2 it follows that the numerator in (3.25)-(3.26) can be bounded above by

-(p LM k ) T (-g δ k (x k ) + r k ) -(p LM k ) T g δ k (x k ) + R + (1 + β)κ d λ α k p LM k 2 ≤ ≤ p LM k r k + κ d λ α k 2 + p LM k + L 2 p LM k 2 + (1 + β)κ d λ α k p LM k 2 ≤ ≤ 2θ 1 λ 1+α k + 2L λ 2 k + 4(3 + β)κ d λ 2-α k + 8κ d ḡ λ 3-α k g δ k (x k ) 2 , with ḡ = κ J 2f δ0 (x 0 ). From (2.6) it follows 1 - ρ δ k k (p LM k ) 2 ≤ κ 2 J /λ min + 1 θ 2θ 1 λ α k + 2L λ k + 4(3 + β)κ d λ 1-α k + 8κ d ḡ λ 2-α k ≤ ϕ λ 1-α k , with ϕ defined in (3.24) and from (3.23) ρ δ k k (p LM k ) ≥ 2η 1 > η 1 .
We can now state the following result, which guarantees that eventually the iterations are successful, provided that λ max > max{νλ * , λ}.

(3.27)

Corollary 1 Let Assumptions 1, 2 and 3 hold. Assume that λ max is chosen to satisfy (3.27).

Then, there exists an iteration index k such that the iterations generated by Algorithm 2.1 are successful for k > k. Moreover,

λ k ≤ max γ max{νλ * , λ}, λ max k > 0. (3.28)
Proof Notice that by the updating rules at step 3 of Algorithm 2.1, λ k increases in case of unsuccessful iterations and it is never decreased. Therefore, after a finite number of unsuccessful iterations it reaches the value max{νλ * , λ}. Moreover, condition (3.27) and the Algorithm's updating rules guarantee that λ k > max{νλ * , λ} for all the subsequent iterations. Then, by Lemma 7 it follows that eventually the iterations are successful. Finally, the parameter updating rules yield (3.28).

We are now ready to state and prove the global convergence of Algorithm 2.1 under the following further assumption: Assumption 4 Assume that λ max is chosen large enough to satisfy

λ max > γ max{νλ * , λ}. (3.29)
Notice that, under this assumption λ k ≤ λ max for all k > 0. In practice the choice of this value is not crucial. If a rather large value is set for this quantity the stopping criterion is usually satisfied before that value is reached. Moreover, since both λ * , λ depend on known algorithm's parameters, on the gradient Lipschitz constant L and on K in (1.4), assuming to be able to estimate these two latter quantities, it is possible to choose a value of λ max satisfying (3.29).

Theorem 1 Let Assumptions 1, 2, 3 and 4 hold. The sequences {δ k } and {x k } generated by Algorithm 2.1 are such that

lim k→∞ δ k = 0, lim k→∞ g(x k ) = 0.
Proof From the updating rule of the accuracy level, {δ k } is a decreasing sequence and so it is converging to some value δ * . Denoting with k s the first successful iteration and summing up over all the infinite successful iterations, from Lemma 1 and Assumption 4 we obtain

f δ k s -1 (x k s ) -lim inf k→∞ f δ k (x k+1 ) ≥ ksucc (f δ k-1 (x k ) -f δ k (x k+1 )) ≥ η 1 2 θ κ 2 J + λ max ksucc g δ k (x k ) 2 , so ksucc g δ k (x k )
2 is a finite number and g δ k (x k ) → 0 on the subsequence of successful iterations, so that lim inf k→∞ g δ k (x k ) = lim k→∞ g δ k (x k ) = 0, taking into account that by Corollary 1 the number of unsuccessful iterations is finite. Finally from (2.11) and (3.16) we have that

δ k ≤ κ d λ α k p LM k 2 ≤ 4κ d g δ k (x k ) 2 λ 2-α min ,
so we can conclude that δ k converges to zero and by (1.4) it follows that lim k→∞ g(x k ) = 0.

Local convergence

In this section we report on the local convergence of the proposed method. To this aim, it is useful to study the asymptotic behaviour of the inexact step. We first establish that, if p LM k satisfies (2.8), then asymptotically p LM k tends to assume the direction of the negative approximated gradient -g δ k (x k ). Then, we study the local convergence of the gradient method with a perturbed gradient step, where the accuracy in the gradient is driven by the accuracy control strategy employed.

Lemma 8 Let Assumptions 1, 2, 3 and 4 hold. Then

lim k→∞ (p LM k ) i + θ κ 2 J + λ k (g δ k (x k )) i = 0 for i = 1, . . . , n,
where (•) i denotes the i-th vector component.

Proof From (2.6)

θ 2 g δ k (x k ) 2 κ 2 J + λ k ≤ m k (x k ) -m k (x k + p LM k ) = -(p LM k ) T g δ k (x k ) - 1 2 (p LM k ) T (J δ k (x k ) T J δ k (x k ) + λ k I)p LM k ≤ -(p LM k ) T g δ k (x k ) - 1 2 λ k p LM k 2 .
Therefore, as θ ∈ [1, 2) from Lemma 1, it follows

θ g δ k (x k ) 2 κ 2 J + λ k + 2(p LM k ) T g δ k (x k ) + λ k θ p LM k 2 < 0, θ κ 2 J + λ k g δ k (x k ) + κ 2 J + λ k θ p LM k 2 ≤ κ 2 J θ p LM k 2 .
Then, from Lemma 4

θ κ 2 J + λ k g δ k (x k ) + p LM k 2 ≤ κ 2 J κ 2 J + λ k p LM k 2 ≤ 4κ 2 J g δ k (x k ) 2 κ 2 J λ 2 min
and the thesis follows as the right-hand side goes to zero when k tends to infinity from Theorem 1.

From Lemma 8, if λ k is large enough, p LM k tends to assume the direction of g δ k (x k ) with step-length θ κ 2 J +λ k . Then, eventually the method reduces to a perturbed steepest descent method with step-length and accuracy in the gradient inherited by the updating parameter and accuracy control strategies employed.

In the following theorem we prove local convergence for the steepest descent step resulting from our procedure. The analysis is inspired by the one reported in [22, §1.2.3], which is extended to allow inaccuracy in gradient values. It relies on analogous assumptions.

Theorem 2 Let x * be a solution of problem (1.1). Let Assumptions 1 and 3 hold and let {x k } be a sequence such that

x k+1 = x k + p SD k , k = 0, 1, 2, . . . with p SD k = -h(λ k )g δ k (x k ), (4.30)
the perturbed steepest descent step with step-length h(λ k ) = θ κ 2 J +λ k . Assume that there exists r > 0 such that f is twice differentiable in B r (x * ) and let H be its Hessian matrix. Assume that H(x) -H(y) ≤ M x -y for all x, y ∈ B r (x * ) and let 0 < l ≤ L < ∞ be such that lI H(x * ) LI. Assume that there exists an index k for which xk -x * < r and

λ k > max θ( L + l) 2 , λ * 1 + 2L l 2/(2-α) , (4.31)
where λ * is defined in (3.17) and r = min{r, l M }. Then for all k ≥ k the error is decreasing, i.e. x k+1 -x * < x k -x * , and x k -x * tends to zero.

Proof We follow the lines of the proof of Theorem 1.2.4 in [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF] for an exact gradient step, taking into account that our step is computed using an approximated gradient. As g(x * ) = 0,

g(x k ) = g(x k ) -g(x * ) = 1 0 H(x * + τ (x k -x * ))(x k -x * ) dτ := G k (x k -x * ),
where we have defined (3.18)). If we let k ≥ k, (4.31) ensures λ k > λ * , and c k < 1. Then, from (4.32) and the Lipschitz continuity of g we get

G k = 1 0 H(x * + τ (x k -x * )) dτ . From (4.30), x k+1 -x * = x k -x * -h(λ k )g(x k ) + h(λ k )(g(x k ) -g δ k (x k )) = = (I -h(λ k )G k )(x k -x * ) + h(λ k )(g(x k ) -g δ k (x k )). From (3.19) g δ k (x k ) -g(x k ) ≤ c k g δ k (x k ) ≤ c k g δ k (x k ) -g(x k ) + c k g(x k ) . (4.32) Notice that c k = λ * λ k 1-α 2 (see
(1 -c k ) g δ k (x k ) -g(x k ) ≤ c k g(x k ) -g(x * ) ≤ Lc k x k -x * .
Then, as (4.31) also yields λ

1-α 2 k -(λ * ) 1-α 2 ≥ 2L l (λ * ) 1-α 2 , it follows g δ k (x k ) -g(x k ) ≤ Lc k 1 -c k x k -x * ≤ l 2 x k -x * .
Let us denote e k = x k -x * . Then it holds

e k+1 ≤ I -h(λ k )G k e k + h(λ k ) g(x k ) -g δ k (x k ) ≤ I -h(λ k )G k e k + h(λ k )l 2 e k . ( 4

.33)

From [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF], Corollary 1.2.1

H(x * ) -τ M e k I H(x * + τ (x k -x * )) H(x * ) + τ M e k I.
Then,

l - e k 2 M I G k L + e k 2 M I, 1 -h(λ k ) L + e k 2 M I I -h(λ k )G k 1 -h(λ k ) l - e k 2 M I.
If we denote with

a k (h(λ k )) = 1 -h(λ k ) l - e k 2 M , b k (h(λ k )) = 1 -h(λ k ) L + e k 2 M , we obtain a k (h(λ k )) > -b k (h(λ k )) as by (4.31) h(λ k ) < 2 l+ L . Then it follows I -h(λ k )G k ≤ max{a k (h(λ k )), -b k (h(λ k ))} = 1 -h(λ k )l + M h(λ k ) 2 e k .
From (4.33)

e k+1 ≤ 1 - h(λ k )l 2 + M h(λ k )e k 2 e k < e k if e k < r = l M . Let us estimate the rate of convergence. Let us define q k = lh(λ k ) 2 and m k = M h(λ k ) 2 = q k r . Notice that as e k < r < q k +1 m k = 2 M h(λ k ) + l M , then 1 -m k e k + q k > 0. So e k+1 ≤ (1 -q k )e k + me 2 k = e k 1 -(m k e k -q k ) 2 1 -(m k e k -q k ) ≤ e k 1 -m k e k + q k 1 e k+1 ≥ 1 + q k -m k e k e k = 1 + q k e k -m k = 1 + q k e k - q k r , 1 e k+1 - 1 r ≥ (1 + q k ) 1 e k - 1 r ≥ (1 + q M ) 1 e k - 1 r > 0,
with q M = lθ 2(κ 2 J +λmax) . Then, we can iterate the procedure obtaining

1 e k ≥ 1 e k - 1 r ≥ (1 + q M ) k-k 1 ek - 1 r , e k ≤ 1 1 + q M k- k rek r -ek ,
and the convergence of x k -x * to zero follows.

Note that if in Algorithm 2.1 we choose

λ max > max γλ * , γ λ, θ( L + l) 2 , λ * 1 + 2L l 2/(2-α)
we have that it exists k such that for k ≥ k, all the iterations are successful and (4.31) is satisfied. Then, Theorem 2 shows the local behaviour of our procedure, provided that an accumulation point x * exists, at which the Hessian satisfies the assumptions in Theorem 2.

We point out however that overall the procedure benefits from the use of a genuine Levenberg-Marquardt method till the last stage of convergence, despite the use of increasingly large regularization parameters. We will see in the numerical results section that our method gains an overall faster convergence rate compared to a pure steepest descent method. Moreover this can be gained at a modest cost, as we solve the linear systems inexactly by an iterative solver. The number of inner iterations is small, even if the required inexactness level decreases with λ k . In fact, when the regularization term is large

J δ k (x k ) T J δ k (x k ) + λ k I λ k I.

Complexity

In this section we will provide a global complexity bound for the procedure sketched in Algorithm 2.1. The analysis is inspired by that reported in [START_REF] Zhao | Global complexity bound of the Levenberg-Marquardt method[END_REF]. We will prove that the number of iterations required to obtain an -accurate solution is O( -2 ).

Let us observe that in our procedure the regularization parameter at the current iteration depends on the outcome of the previous iteration and consequently let us define the following sets

S 1 = {k + 1 : ρ δ k k (p LM k ) ≥ η 1 ; g δ k (x k ) < η 2 /λ k }, (5.34)
S 2 = {k + 1 : ρ δ k k (p LM k ) ≥ η 1 ; g δ k (x k ) ≥ η 2 /λ k }, (5.35)
S 3 = {k + 1 : ρ δ k k (p LM k ) < η 1 }. (5.36) Let N i = |S i | for i = 1, 2, 3
, so that the number of successful iterations is N 1 + N 2 and the number of unsuccessful iterations is N 3 . Moreover S 1 can be split into two subsets

S 1 = A ∪ B = {k + 1 ∈ S 1 : γλ k < λ max } ∪ {k + 1 ∈ S 1 : γλ k ≥ λ max },
taking into account that if k + 1 ∈ S 1 from the updating rule at step 3a either λ k+1 = γλ k (A), or λ k+1 = λ max (B). The analysis is made under the following Assumption:

Assumption 5 Let us assume that the procedure sketched in Algorithm 2.1 is stopped when

g δ k (x k ) ≤ .
In the following Lemma we provide un upper bound for the number of successful iterations.

Lemma 9 Let Assumptions 1, 2, 3, 4 and 5 hold. Let k s be the index of the first successful iteration.

1. The number N 1 of successful iterations belonging to set S 1 is bounded above by:

N 1 ≤ f δ ks-1 (x ks ) 2 η 1 κ 2 J + λ max θ 2 = O( -2 ).
2. The number N 2 of successful iterations belonging to set S 2 is bounded above by a constant independent of :

N 2 ≤ f δ ks-1 (x ks ) 2 η 1 κ 2 J + λ max θ λ max η 2 2 .
Proof From (2.6), as λ k ≤ λ max for all k, it follows

m k (x k ) -m k (x k + p LM k ) ≥ θ 2 g δ k (x k ) 2 κ 2 J + λ max .
Then, as the iteration is successful

f δ k-1 (x k ) -f δ k (x k + p LM k ) ≥ η 1 (m k (x k ) -m k (x k + p LM k )) ≥ η 1 2 θ g δ k (x k ) 2 κ 2 J + λ max .
For all k it holds g δ k (x k ) 2 ≥ 2 and in particular for k ∈ S 2

g δ k (x k ) 2 ≥ η 2 λ max 2 .
Then

f δ ks -1 (x ks ) ≥ j∈S1∪S2 (f δj-1 (x j ) -f δj (x j+1 )) = j∈S1 (f δj-1 (x j ) -f δj (x j+1 )) + j∈S2 (f δj-1 (x j ) -f δj (x j+1 )) ≥ η 1 N 1 2 θ κ 2 J + λ max 2 + η 1 N 2 2 θ κ 2 J + λ max η 2 λ max 2 ,
and the thesis follows.

In the following Lemma we provide un upper bound for the number of unsuccessful iterations.

Lemma 10 Let Assumptions 1, 2, 3, 4 and 5 hold. The number of unsuccessful iterations N 3 is bounded above by a constant independent of :

N 3 ≤ log λmax λ0 log γ .
Proof Notice that from Assumption 4 it is not possible to have an iteration index in B before the last unsuccessful iteration. Then, if we denote with N the last unsuccessful iteration index, if k < N is a successful iteration, it belongs to A. Denoting with N a the number of such iterations, it follows

λ N = γ Na+N3 λ 0 ≤ λ max .
Then

N 3 ≤ N a + N 3 ≤ log λmax λ0 log γ ,
and the thesis follows.

Then, taking into account the results proved in the previous Lemmas, we can state the following complexity result, that shows that the proposed method shares the known complexity properties of the steepest descent and trust-region procedures.

Corollary 2 Let Assumptions 1, 2, 3, 4 and 5 hold, and let N T be the total number of iterations performed. Then,

N T ≤ f δ ks-1 (x ks ) 2 η 1 κ 2 J + λ max θ 1 2 + λ max η 2 2 + log λmax λ0 log γ = O( -2 ).
(5.37)

We underline that λ max and therefore the constant multiplying -2 in (5.37) may be large if κ J is large. On the other hand, in the applications of our interest, that we present in next section, it holds κ J 1.

Numerical results

In this section we analyse the numerical behaviour of the Levenberg-Marquardt method described in Algorithm 2.1. We show the results of its application to two large scale nonlinear least-squares problems, arising respectively from data assimilation and machine learning. These problems can be written as

min x∈R n f (x) = 1 2 F (x) 2 + 1 2 x 2 = N j=1 F j (x) 2 + 1 2 x 2 , (6.38)
with F j : R n → R, for j = 1, . . . , N .

In both test problems the inaccuracy in the function and gradient arises from the use of a subsampling technique. Then, at each iteration the approximations F δ k to F are built by selecting randomly a subset of the samples indices X k ⊆ {1, . . . , N } such that |X k | = K k for each k. For this reason we will denote Algorithm 2.1 as subsampled Levenberg-Marquardt method (SSLM). Each time condition (2.11) is not verified we increase the size of the subsampled set to improve the accuracy of the approximations. This is done in a linear way by a factor K * , so that if the loop 1-2 is performed n k times it holds

|X k+1 | = K n k * |X k |. (6.39)
Notice that other updates, different from linear, could be used affecting the speed of convergence of the procedure, see for example [START_REF] Bollapragada | Exact and inexact subsampled newton methods for optimization[END_REF][START_REF] Friedlander | Hybrid deterministic-stochastic methods for data fitting[END_REF]. Moreover, the subsampling is performed in a random way. In some cases, like for data assimilation problems, it is possible to devise more efficient strategies taking into account the particular structure of the problem. This leads to a quicker improvement in the accuracy level, the number of samples being the same, see [START_REF] Gratton | Observation thinning in data assimilation computations[END_REF], but this is out of the scope of this paper.

The procedure was implemented in Matlab and run using Matlab 2015a on an Intel(R) Core(TM) i7-4510U 2.00GHz, 16 GB RAM; the machine precision is m ∼ 2 • 10 -16 . We run SSLM with η 1 = 0.25, η 2 = 1.e -3, γ = 1.001, α = 0.9, λ max = 1.e + 6, λ min = 1. We remind that from the update at step 3 of Algorithm 2.1 the generated sequence of regularization parameters is increasing. This is needed to make the accuracy control (2.11) work. However a too quick growth would lead to a slow procedure. Then, the choice of γ in Algorithm 2.1 is relevant, as it determines the rate of growth of the parameters. In practice, it is advisable to choose a value of γ not too big. The chosen value of γ is shown to be effective in controlling the accuracy level without impacting negatively on the rate of convergence. On the other hand, the choice of λ max is not critical as, if a high value is chosen, the stopping criterion is satisfied before that value is reached.

In order to solve the linear subproblems (2.8) we employed the Matlab function cgls available at [START_REF] Saunders | Systems Optimization Laboratory[END_REF], that implements conjugate gradient (CG) method for least-squares. We set the stopping tolerance according to (3.15), where we have chosen θ 1 = θ 2 = 1.e -1. In both problems this choice corresponds to the tolerance k 1.e -1 along all the optimization process. We set to 20 the maximum number of iterations CG is allowed to perform. We will see in the numerical tests that the average number of CG iterations per nonlinear iteration is really low, and this maximum number is never reached.

In the following we are going to compare the performance of the proposed SSLM to that of three inexact methods, all of them used to solve the exact problem (6.38):

-Full Levenberg-Marquardt method (FLM), i.e. the procedure described in Algorithm 2.1, but run using all the available samples, so that K k = N for all k and δ k is zero along all the optimization process. -SLM, an inexact Levenberg-Marquardt method based on a standard update of the regularization parameters:

λ k+1 =      γ 1 λ k if ρ k (p LM k ) > η 2 , λ k if ρ k (p LM k ) ∈ [η 1 , η 2 ], γ 0 λ k if ρ k (p LM k ) < η 1 .
with λ 0 = 0.1, γ 0 = 2, γ 1 = 0.5, η 1 = 0.25, η 2 = 0.75. -An inexact Gauss-Newton method GN with λ k = 0 for all k.

For the three methods the linear algebra operations are handled as for the SSLM, i.e. the linear systems are solved with cgls with the same stopping criterion and maximum number of iterations set for SSLM. The aim of this comparison is to evaluate the effectiveness of the strategy we proposed to control the accuracy of function approximations. Our goal is to show that approximating the solution of (6.38) with the proposed strategy is more convenient in term of computational costs than solving the problem directly, and it does not affect the quality of the approximate solution found. First, we compare SSLM to FLM to show the savings arising from the use of approximations to the objective function, when the same rule for updating the regularization parameters is used. Then, we extend the comparison also to SLM and GN. This is done as we have to take into account that the specific update of the regularization parameters at step 3 of Algorithm 2.1 is designed to be used in conjunction with the strategy to control the accuracy (2.11). To solve the exact problem directly, it may be more effective to use a more standard update, that we expect to result in a quicker procedure.

To evaluate the performance of SSLM and compare it to that of the other solvers, we use two different counters, one for the nonlinear function evaluations and one for matrix-vector products involving the Jacobian matrix (transposed or not), which includes also the products performed in the conjugate gradients iterations. Notice that the counters are cost counters, i.e. they do not count the number of function evaluations or of matrix-vector products, but each function evaluation and each product is weighted according to the size of the samples set. The cost of a function evaluation or of a product is considered unitary when the full samples set is considered, otherwise it is weighted as |X k | N . We use a reference solution the approximation x * provided by FLM with stopping criterion g(x * ) < 1.e-8. We compared the solution approximations computed by all the other procedures to x * and we measured the distance by the Root Mean Square Error (RMSE):

RM SE = n i=1 (x * (i) -x(i)) 2 n .
In the Tables the columns heads have the following meanings: it: nonlinear iterations counter, CG it : average number of CG iterations per nonlinear iteration, cost f : function evaluations cost counter, cost p : matrix-vector products cost counter, |X it |: value of the cardinality of the samples set at the end of the process, RMSE: root mean square error, save f , save p : savings gained by SSLM in function evaluations and matrix-vector products respectively, compared to the FLM.

A data assimilation problem

In this section we consider the data assimilation problem described in [START_REF] Gratton | Observation thinning in data assimilation computations[END_REF]. We consider a onedimensional wave equation system, whose dynamics governed by the following nonlinear wave equation:

∂ 2 u(z, t) ∂t 2 - ∂ 2 u(z, t) ∂z 2 + f (u) = 0, (6.40) u(0, t) = u(1, t) = 0, (6.41) u(z, 0) = u 0 (z), ∂u(z, 0) ∂t = 0, (6.42) 0 ≤ t ≤ T, 0 ≤ z ≤ 1, (6.43)
where f (u) = µe νu . (6.44)

The system is discretized using a mesh involving n = 360 grid points for the spatial discretization and N t = 64 for the temporal one. We look for the initial state u 0 (z), which is possible to recover solving the following data assimilation problem:

min x∈R n 1 2 x -x b 2 B -1 + 1 2 Nt j=0 H j (x(t j )) -y j 2 R -1 j (6.45)
where, given a symmetric positive definite matrix M , x 2 M denote x T M x. Here x b ∈ R n is the background vector, which is an a priori estimate, y j ∈ R mj is the vector of observations at time t j and H j is the operator modelling the observation process at the same time. The state vector x(t j ) is the solution of the discretization of the nonlinear model (6.40)-( 6.43) at time t j . Matrices B ∈ R n×n and R j ∈ R mj ×mj represent the background-error covariance and the observation-error covariance at time t j respectively. In our test we build the background vector and the observations from a chosen initial true state x T by adding a noise following the normal distributions N (0, σ 2 b ) and N (0, σ 2 o ) respectively. We have chosen σ b = 0.2, σ o = 0.05 and we assume the covariances matrices to be diagonal: B = σ 2 b I and R j = σ 2 o I mj for each j. For further details on the test problem see [START_REF] Gratton | Observation thinning in data assimilation computations[END_REF]. We can reformulate (6.45) as a least-squares problem (6.38), defining

F (x) =    1 σo (H 0 (x(t 0 )) -y 0 ) . . . 1 σo (H Nt (x(t Nt )) -y Nt )    ,
where (H j (x(t j )) -y j ) ∈ R mj for j = 1, . . . , N t . We assume to have an observation for each grid point, so that the total number of observations is N = n•N t = 23040. The full problem (6.38) is obtained when all the observations are considered, in this case m j = 360 for every j. The approximations F δ k are obtained selecting randomly K k observations among the available ones, so that vectors (H j (x(t j )) -y j ) have dimension m j ≤ 360 and it may be m i = m j if i = j.

We consider two different problems of the form (6.45), corresponding to two different values of µ in (6.44). We consider a mildly nonlinear problem, choosing µ = 0.01 because this is usually the case in practical data assimilation applications and then we increase µ to 0.5 to consider the effect of the nonlinearity on our procedure.

We remind that we denote with x * the solution approximation found by FLM, computed asking g(x * ) < 1.e -8. If we compare this approximation to the true state x T we obtain

n i=1 (x * (i)-x T (i)) 2 n
= 5.2e -3. Then, we study the effect of the presence of inaccuracy in the function arising from the use of subsampling techniques and we compare the solution found by SSLM to x * . Taking into account (1.3) the accuracy level δ k is approximated in the following way. At the beginning of the iterative process δ 0 is set to |f δ0 (x 0 ) -f (x 0 )|. Then, it is left constant and updated only when condition (2.11) is not satisfied as follows

δ k |f δ k (x k ) -f (x k )|. ( 6.46)
We stress that this computation is not expensive as it does not affect the matrix-vector product counter and marginally contributes to the function evaluations counter. In fact, the evaluation of the full function is required only when condition (2.11) is not met and is performed just once in the loop 1-2, as x k is fixed.

In general a very accurate solution is not needed in practical applications, so the optimization process is stopped as soon as the residuals are detected to be Gaussian. As a normality test we employ the Anderson-Darling test, see [START_REF] Anderson | A test of goodness of fit[END_REF], which tests the hypothesis that a sample has been drawn from a population with a specified continuous cumulative distribution function Φ, in this case the Gaussian distribution. Assuming to have a sample of n ordered data

{x 1 ≤ x 2 ≤ • • • ≤ x n-1 ≤ x n },
a statistic is built in the following way:

W 2 n = -n - 1 n n i=1 (ln(Φ(x i )) + ln(1 -Φ(x n-i+1 ))).
If W 2 n exceeds a given critical value, then the hypothesis of normality is rejected with some significance level. We used the critical value for significance level 0.01 which is 1.092, see [START_REF] Stephens | Edf statistics for goodness of fit and some comparisons[END_REF].

The performance of our procedure is affected by mainly three choices: the cardinality of the starting set of observations K 0 , factor K * in (6.39) and the parameter κ d in (2.11). The choice of K * determines how much the accuracy is increased at each loop at steps 1-2 of Algorithm 2.1. A too small value could lead to a too low increase, gaining an accuracy level δ k that still does not satisfy condition (2.11), so that the loop should be performed n k > 1 times. Each time the accuracy is increased, the computation of a trial step is required, through the solution of a linear system (2.8) of increasing size, so it is advisable to consider a reasonable increase in the subsets size. Again too small values of κ d generally lead to too frequently rise the accuracy. In this section we investigate the effect of parameter κ d combined with different values of K 0 , while in the next section we analyse the effect of the choice of K * .

We run the procedure choosing two different values of K 0 , combined with different values of κ d , while K * is kept fixed, K * = 1.5. Tables 1 and2 refer to problem µ = 0.5 for K 0 = 2000 and K 0 = 5000 respectively, while Table 3 refers to test problem µ = 0.01 for K 0 = 2000 and K 0 = 7000. We also solved the two problems using FLM, GN and SLM. In these tables we report just figures corresponding to runs of FLM. On these problems indeed, FLM converges quickly and the update of the regularization parameter does not play a key role in the convergence. As a consequence, the behaviour of GN and SLM is really similar to that of FLM. Then, in the first column we report the results of the optimization process performed by FLM and in the last two rows the savings gained by SSLM in function evaluations save f and matrix-vector products save p .

We notice that SSLM requires on average a higher number of CG iterations than FLM and this is due to the need of recomputing the step when (2.11) is not satisfied. This number is affected by the choice of parameter κ d , as generally it decreases with κ d . This is less evident for µ = 0.5, while it is crystal clear for µ = 0.01. Moreover, the value of κ d does not affect the number of nonlinear iterations performed by SSLM, while it has a deep impact on the procedure cost, as we can see from the significant variation of function evaluations and matrix-vector products counters. We notice that in all cases our procedure is much less expensive than FLM, and consistent savings are provided by higher values of κ d . In these cases indeed the accuracy is increased less frequently, as condition (2.11) is more likely satisfied, and as a result the overall process is performed with less observations and is less expensive, at the cost however of a less accurate solution. Indeed, if κ d is too large the accuracy control strategy is not effective, the accuracy level may be never increased and the sequence may approach a solution of a problem with inaccurate function, that can be a bad approximation to that of (1.1). In Figure 1 we compare solution approximations for µ = 0.5 provided by: FLM (up left), SSLM with K 0 = 5000 and κ d = 10 (up right), SSLM with K 0 = 2000 (bottom left) and K 0 = 5000 (bottom right) and κ d = 10000. In all the plots the solid line represents the true state x T and the dotted line the computed solution approximation. It is evident that in the bottom left plot, corresponding to the last column of Table 1, the solution found is less accurate. In fact, due to the high value of κ d the accuracy is never increased and the problem is solved considering just the samples in the initial subset, which are not sufficient to obtain the same error in the approximate solution gained by the FLM. Then κ d should not be chosen too high, especially if K 0 is small. On the other hand, if K 0 is large enough one can expect to gain a low error in the solution approximation even with a higher κ d . For example K 0 = 5000 is large enough to obtain a good solution approximation, so the best performance is obtained with In Figure 2 we relate the savings gained with the corresponding error in the solution. The solid lines refer to Table 1 while the dotted ones to Table 2. In the left plot we report the savings in function evaluations (lines marked by stars) and in matrix-vector products (lines marked by circles), while in the right plot the RMSE, versus κ d . If K 0 = 5000 the error is almost the same for all choices of κ d but the savings increase with κ d , while if K 0 = 2000 the most significant savings are obtained choosing large κ d , but at the expense of a higher error.

Notice also that in the tests the final value |X it | is always less than N , which confirms that it is not necessary to use all the available observations to obtain a good solution approximation.

In Figure 3 we report as an example for problem µ = 0.5 and K 0 = 2000, κ d = 10 the behaviour of the accuracy level (left plot) and that of the error (right plot) through iterations. We underline that condition (2.11) is not violated at each iteration, then the accuracy is kept fixed for some consecutive iterations, and the evaluation of the full function, by the computation of the remaining components, is only sporadically necessary.

Regarding the choice µ = 0.01 we report statistics in Table 3. The problem is almost linear, so it is solved in few iterations. Due to the really low number of iterations, it is advisable to start with a rather large initial set to avoid converging to a solution of a problem with approximated objective function and to gain the same solution accuracy as FLM. In this case the procedure is less sensitive to the choice of parameter κ d than in the other case and only significant changes in κ d affect its performance. Also for this problem the use of SSLM provides significant savings compared to FLM.

A machine learning problem

In this section we consider a binary classification problem. Suppose to have at disposal a set of pairs {(z i , y i )} with z i ∈ R n , y i ∈ {-1, +1} and i = 1, . . . , N . We consider as a training objective function the logistic loss with l 2 regularization, see [START_REF] Bollapragada | Exact and inexact subsampled newton methods for optimization[END_REF]:

f (x) = 1 2N N i=1 log(1 + exp(-y i x T z i )) + 1 2N x 2 . (6.47)
Since this is a convex nonlinear programming problem, it could potentially be solved also by a subsampled Newton method. Here, for sake of gaining more computational experience with our approach, we reformulate it as a least-square problem (6.38), scaled by N , where F : R n → R N is given by

F (x) =    log(1 + exp(-y 1 x T z 1 )) . . . log(1 + exp(-y N x T z N ))    .
We consider the CINA dataset available at [START_REF]Causality workbench team[END_REF], for which n = 132 and that is divided in a training set of size N = 16033 and a testing set of size Ñ = 10000. We build the approximations f δ k as:

f δ k (x) = 1 2K k i∈X k log(1 + exp(-y i x T z i )) + 1 2K k x 2 .
We start the optimization process with K 0 = 1000 and we stop the procedure when g δ k (x k ) ≤ 1.e -4. Parameter κ d is set to 100. The results provided in this section are obtained without computing the accuracy level δ k as outlined in the previous subsection (see (6.46)), in order to avoid the evaluation of the full function f (x) each time condition (2.11) is not satisfied. Indeed, we can spare these evaluations by estimating the accuracy level in the following way:

δ k 2(N -K k ) K k , with K k = |X k |. (6.48)
This approximation is based on the observation that if the components F i (x) of F (x) were Gaussian, i / ∈X k F i (x) 2 would follow a Chi-squared distribution with standard deviation 2(N -K k ). Even if the normality assumption does not hold, this estimation works well in practice, as we will see in the following.

We study the effect of the choice of K * in (6.39) on the procedure performance. Then, we fix K 0 = 132 and κ d = 100. Once the problem is solved, the computed solution x † is used to classify the samples in the testing set. The classification error e te is defined as 1 2 Ñ Ñ i=1 log(1 + exp(-ŷ i x †T z i )), which consists of f (x † ) omitting the regularization term 1 2 Ñ x 2 , cf. [START_REF] Bollapragada | Exact and inexact subsampled newton methods for optimization[END_REF], and employing the estimations ŷi for y i , computed for z i in the testing set as ŷi The results are reported in Table 4. Concerning the three reference methods, we can notice that, as expected, the convergence rate of SLM or GN is faster compared to that of the FLM (the number of outer iterations performed is half of those required by FLM). However, the average number of CG iterations per outer iteration is more than the double. Indeed, the linear systems to be solved are less regularized, as the regularization parameters are smaller, so that the linear solver requires more iterations to converge. As a result, the cost of SLM and GN in terms of function evaluations is lower than that of FLM, but the cost in terms of matrix-vector products is comparable.

= +1 if σ(z i ) ≥ 0.5 -1 if σ(z i ) < 0.
The results reported in Table 4 show that for every choice of K * SSLM provides significant savings compared to all the reference methods, in terms of function evaluations (except for K * = 3.5) and especially in terms of matrix-vector products. The savings in percentage form in the last rows (save f , save p ) are computed with respect to FLM method. The RMSE and the testing error show that the quality of the approximate solutions is not affected by the use of the subsampling technique. The counters anyway are affected by the choice of K * , both too small and too large values lead to a more expensive SSLM procedure. The effect of small parameter values is clearly shown in Figure 4. In each plot are reported the values of n k (dotted line) and the number of CG iterations (dashed line) for each nonlinear iteration k, for K * = 1.1 (up left), K * = 2 (up right) and K * = 3.5 (bottom). The values of n k indicate how many times loop 1-2 in Algorithm 2.1 is performed (n k = 1 means that the accuracy in function values is increased once at iteration k, n k = 0 means that the accuracy is kept constant for that iteration). We notice that the number of CG iterations performed in a nonlinear iteration in which the accuracy is increased is much higher than that required by iterations in which it is kept constant, as the linear system (2.8) is solved more than once. When K * is small, the accuracy is increased of a small amount when condition (2.11) is not satisfied, and this leads to the need of increasing the accuracy of the approximations more often and then to perform a higher number of CG iterations, as it is shown in Table 4 and in Figure 4. On the other hand, with large values of K * the accuracy is increased less often, but a too large choice leads K k to quickly reach the maximum value N , so that many expensive iterations are performed and so again the computational costs are higher. In the left plot of Figure 5 we report values of |X k | along the iterations for different values of K * . We notice that when K * is small the size is often increased of a small amount while for larger values it raises quickly.

In the right plot of Figure 5 we compare the cost of matrix-vector products at each iteration of SSLM for various K * and of the FLM (solid line). We can see that significant savings are obtained at the beginning of the optimization process, due to the reduced size of the subproblems, which compensate the greater costs in the final stage, when the samples subsets are of size close to N and additional CG iterations are required when condition (2.11) is not satisfied. Notice that in all the tests performed the average number of CG iterations is generally low and the maximum number of allowed iterations is never reached. This is due to the low accuracy we solve the linear systems with, which anyway is enough to get convergence. Despite the use of an increasing sequence of regularization parameters, our method still gains the benefits of a quicker convergence compared to first-order methods. this is obtained at no great expense, as the number of CG iterations is extremely low. We used the Matlab function steep implementing the steepest descend method and available at [START_REF] Kelley | Iterative Methods for Optimization: Matlab Codes[END_REF] to solve the problem with exact objective function and after 1000 iterations the desired accuracy was not yet reached and the norm of the gradient was of the order of 1.e -3.

In the left plot of Figure 6 we consider SSLM with K * = 2 and we compare the approximation of the accuracy level provided by (6.46) (solid line) with that estimated through (6.48) (dashed line). The estimate is good enough to ensure the procedure run with the estimated accuracy level (SSLM est ) to achieve the same performance as that run approximating it via (6.46) (SSLM appr ), as it is shown in Table 5. We highlight that the use of (6.48) does not affect the quality of the classification process. Moreover it produces a saving in terms of function evaluations, even if also the approximation of the accuracy level through (6.46) is affordable, as increasing the accuracy, and so the evaluation of the full function, is needed only sporadically. For example for K * = 2 it is needed just four times along all the optimization process, as it is evident from Figure 4 or Figure 5.

Finally, in the right plot of Figure 6 we compare the decrease in the model m k (x k ) -m k (x k + p LM k ) to that of the term 1 2 λ k p LM k 2 used to approximate it in the control (2.11). As we have claimed in Section 2, the approximation is good, showing that in practice the assumption we made is verified.

Conclusions

In this paper we proposed an inexact Levenberg-Marquardt approach to solve nonlinear leastsquares problems with inaccurate function and gradient, assuming to be able to control the accuracy of the approximations. We proved that the proposed approach guarantees global convergence to a solution of the problem with exact objective function and that asymptotically the step tends to the direction of the negative approximated gradient. Then, we performed a local analysis for the perturbed steepest descent method we reduce to. The procedure was tested on two problems arising in machine learning and data assimilation. The results show that the implemented procedure is able to find a first-order solution of the problem with exact objective function and that the proposed strategy to control the accuracy level allows significant savings both in function evaluations and in matrix-vector products, compared to the same procedure performed on the problem with exact objective function. The provided numerical results also show the efficiency of the procedure in terms of inner Krylov iterations. Indeed, a very rough accuracy in the solution of the arising linear systems is imposed and as a result the number of performed Krylov iterations is quite low. Overall the method gains a faster convergence rate compared to a pure steepest descent method.

Lemma 3

 3 Let Assumption 1 hold and let p LM k be defined as in Lemma 1. It exists k ≥ 0 such that the regularization parameters {λ k } generated by Algorithm 2.1 satisfy λ k = λ max for any successful iteration k, with k ≥ k.

Fig. 1

 1 Fig. 1 Problem µ = 0.5. Comparison of true state (solid line) and computed solution (dotted line) computed by FLM (up left), SSLM with K 0 = 5000 and κ d = 10 (up right), SSLM with K 0 = 2000 (bottom left) and K 0 = 5000 (bottom right) with κ d = 10000.

Fig. 2

 2 Fig. 2 Left plot: save f (lines marked by a star) and savep (lines marked by a circle) for tests in Tables 1 (solid lines), 2 (dotted lines). Right plot: corresponding accuracy in problem's solution.

Fig. 3

 3 Fig. 3 Problem µ = 0.5, K 0 = 2000, κ d = 10. Left: log plot of the accuracy level versus iterations. Right: log plot the RMSE versus iterations.

5 ,

 5 where σ(z) = 1 1+exp(-z T x † ) . Notice that in these runs all the available training samples are used during the optimization process, so that the final value |X it | always reaches the maximum value N .

Fig. 4

 4 Fig. 4 Values of n k (dotted line) and number of CG iterations (dashed line) per nonlinear iteration, for K * = 1.1 (up left), K * = 2 (up right) and K * = 3.5 (bottom).

Fig. 5

 5 Fig. 5 Left: values of |X k | versus iterations for different values of K * . Right: Comparison of matrix-vector cost counter for FLM (solid line) and SSLM with various choices of K * .

  x k )-m k (x k +p k ) 1/2 λ k ||p k || 2

Fig. 6 2 (

 62 Fig. 6 Left: comparison of approximated accuracy level (solid line) and estimated accuracy level (dashed line) during run of SSLM. Right: decrease in the model m k (x k )-m k (x k +p LM k ) versus iterations (dashed line) compared to decrease of 1 2 λ k p LM k 2 (solid line).

  depends on the accuracy level. If condition(2.11) is not satisfied at iteration k, the uncertainty in the function values is considered too high and the accuracy is increased, i.e. δ k is decreased. We will prove in Lemma 2 that after a finite number of reductions condition (2.11) is met.

	1 2 , 1 . Parameter α in (2.11) is introduced to guarantee global
	convergence of the procedure, as shown in Section 3. Notice also that (2.11) is an implicit relation,
	as p LM k
	Algorithm 2.1: Levenberg-Marquardt method for problem (1.1)

Table 1

 1 Performance of SSLM for test problem µ = 0.5 and K 0 = 2000.

		FLM	κ d = 1	κ d = 10	κ d = 100	κ d = 1000	κ d = 10000
	it	9	11	12	12	12	11
	CG it	2.4	5.4	4.9	4.2	4.2	3.9
	cost f	10	9.7	6.1	3.3	3.2	2.0
	costp	67	46.1	26.8	14.9	13.5	10.3
	|X it |	23040	15188	6750	3000	3000	2000
	RMSE 1.2e-2 3.0e-2	2.8e-2	3.8e-2	4.4e-2	7.8e-2
	save f		3%	39%	67%	68%	80%
	savep		31%	60%	78%	80%	85%

Table 2

 2 Performance of SSLM for test problem µ = 0.5 and K 0 = 5000.

					FLM		κ d = 1	κ d = 10	κ d = 100	κ d = 1000	κ d = 10000
			it		9			11		11		12	12	12
			CG it	2.4		4.1		3.9		4.0	4.1	3.7
			cost f		10			9.1		6.5		5.1	4.9	3.6
			costp	67			54.8		37.2	34.6	32.9	27.3
			|X it |		23040		16875	11250	7500	7500	5000
			RMSE 1.2e-2 2.7e-2	3.0e-2	2.1e-2	2.1e-2	2.7e-2
			save f				9%		35%	49%	51%	64%
			savep				18%		44%	48%	51%	59%
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
						z

Table 3

 3 Performance of SSLM for test problem µ = 0.01.

			K 0 = 2000			K 0 = 7000
		FLM	κ d = 1	κ d = 10, 100 κ d = 1000	κ d = 1, 10	κ d = 100, 1000
	it	3	3	4	3	3	3
	CG it	3.0	12.3	9.5	6.0	5.7	4.0
	cost f	4	2.9	3.5	1.3	3.1	1.9
	costp	27	12.6	10.8	3.9	15.3	10.0
	|X it |	23040	6750	4500	2000	10500	7000
	RMSE 6.8e-3 2.0e-2	1.1e-2	3.4e-2	1.5e-2	1.6e-2
	save f		27%	12%	67%	22%	52%
	savep		53%	60%	85%	43%	63%

Table 4

 4 Performance of SSLM for machine learning test problem for different values of K * .

	K *

Table 5

 5 Comparison of subsampled Levenberg-Marquardt method with estimated accuracy level (first row, SSLMest) and accuracy level approximated by (6.46) (second row, SSLMappr).

	Solver	it	CG it	cost f	costp	|X it |	err	e te
	SSLMest	38 7.5	15.9	316.7	16000 5.4e-2 0.187
	SSLMappr	37 7.4	17.7	318.1	16000 5.7e-2 0.186

Acknowledgements We thank the authors of [17] for providing us the Matlab code for the data assimilation test problem.