Finding Zeros of H\"{o}lder Metrically Subregular Mappings via Globally Convergent Levenberg-Marquardt Methods

TL;DR

This paper introduces two globally convergent Levenberg-Marquardt algorithms for finding zeros of H"{o}lder metrically subregular mappings, with proven convergence, complexity bounds, and successful application to biological data.

Contribution

The paper develops two novel Levenberg-Marquardt methods with global convergence and complexity analysis for non-isolated zeros of H"{o}lder metrically subregular mappings.

Findings

Methods achieve global convergence to stationary points.

Worst-case complexity is O(ε^{-2}) evaluations.

Numerical results successfully applied to biological data.

Abstract

We present two globally convergent Levenberg-Marquardt methods for finding zeros of H\"{o}lder metrically subregular mappings that may have non-isolated zeros. The first method unifies the Levenberg- Marquardt direction and an Armijo-type line search, while the second incorporates this direction with a nonmonotone trust-region technique. For both methods, we prove the global convergence to a first-order stationary point of the associated merit function. Furthermore, the worst-case global complexity of these methods are provided, indicating that an approximate stationary point can be computed in at most $\mathcal{O}(\varepsilon^{-2})$ function and gradient evaluations, for an accuracy parameter $\varepsilon>0$. We also study the conditions for the proposed methods to converge to a zero of the associated mappings. Computing a moiety conserved steady state for biochemical reaction networks can be cast as the problem of finding a zero of a H\"{o}lder metrically subregular mapping. We report encouraging numerical results for finding a zero of such mappings derived from real-world biological data, which supports our theoretical foundations.

Figures

Tables

Table 1. ?tablename? 1: The list of 21 biological models, where the stoichiometric matrix N 𝑁 N is m × n 𝑚 𝑛 m\times n and r ​ a ​ n ​ k 𝑟 𝑎 𝑛 𝑘 rank is the rank of the matrix N 𝑁 N .

Model	$m$	$n$	$rank$	Model	$m$	$n$	$rank$
1. Ecoli_core	72	73	61	12. iMB745	525	598	490
2. iAF692	462	493	430	13. iNJ661	651	764	604
3. iAF1260	1520	1931	1456	14. iRsp1095	966	1042	921
4. iBsu1103	993	1167	956	15. iSB619	462	508	435
5. iCB925	415	558	386	16. iTH366	583	606	529
6. iIT341	424	428	392	17. iTZ479_v2	435	476	415
7. iJN678	641	669	589	18. iYL1228	1350	1695	1280
8. iJN746	727	795	700	19. L_lactis_MG1363	483	491	429
9. iJO1366	1654	2102	1582	20. Sc_thermophilis_rBioNet	348	365	320
10. iJP815	524	595	501	21. T_Maritima	434	470	414
11. iJR904	597	757	564

Equations

$$h(x)=0,x\in {\mathbb{R}}^m,$$

$$\underset{x\in {\mathbb{R}}^m}{\min }\psi (x),\text{with }\psi :{\mathbb{R}}^m\to \mathbb{R}\text{ given by }\psi (x):=\frac{1}{2}\Vert h(x){\Vert }^2,$$

$$\beta \mathrm{dist}(x,\Omega )\le \Vert h(x)\Vert ,\forall x\in \mathbb{B}(x^{\ast },\mathtt{r})$$

$$\beta \mathrm{dist}(x,\Omega )\le \Vert h(x){\Vert }^{\delta },\forall x\in \mathbb{B}(x^{\ast },\mathtt{r}),$$

$$\underset{d\in {\mathbb{R}}^m}{\min }{\varphi }_k(d),\text{with }{\varphi }_k:{\mathbb{R}}^m\to \mathbb{R}\text{ given by }{\varphi }_k(d):=\Vert \nabla h(x_k{)}^{T}d+h(x_k){\Vert }^2+{\mu }_k\Vert d{\Vert }^2,$$

$$\left(\nabla h(x_k)\nabla h(x_k{)}^T+{\mu }_{k}I\right)d_k=-\nabla h(x_k)h(x_k),$$

$${\mu }_k={\xi }_k\Vert h(x_k){\Vert }^{\eta }+{\omega }_k\Vert \nabla h(x_k)h(x_k){\Vert }^{\eta },$$

$$\psi (x_k+{\alpha }_k{d}_k)\le \psi (x_k)+\sigma {\alpha }_k\nabla \psi (x_k{)}^T{d}_k,$$

$$q_k(d):=\frac{1}{2}\Vert \nabla h(x_k{)}^{T}d+h(x_k){\Vert }^2.$$

$$r_k:=\frac{\psi (x_k)-\psi (x_k+d_k)}{q_k(0)-q_k(d_k)},$$

$$\nabla \psi (x_k{)}^T{d}_k=-d_k^{T}H(x_k)d_k<0,$$

$$D_k=\{\begin{array}{ll}\psi (x_0)& \mathrm{if}k=0,\\ (1-{\theta }_{k-1})\psi (x_k)+{\theta }_{k-1}{D}_{k-1}& \mathrm{if}k\ge 1,\end{array}$$

$$\Vert \nabla h(x)-\nabla h(y)\Vert \le L\Vert x-y\Vert ,\forall x,y\in {\mathbb{R}}^m.$$

$${\alpha }_k\ge \frac{\rho (1-\sigma ){\mu }^2}{\vartheta (L_0^2+\mu )}:=\hat{\alpha },$$

$$\vartheta :=\frac{1}{2}{L}_0^2+\frac{1}{2}{\overline{\alpha }}^2{\rho }^{-2}{L}^2{\mu }^{-2}{L}_0^2\Vert h(x_0){\Vert }^2+(1+\overline{\alpha }{\rho }^{-1}{L}_0^2{\mu }^{-1})L\Vert h(x_0)\Vert ,$$

$$0\le {\ell }_k\le \frac{\log (\hat{\alpha })-\log (\overline{\alpha })}{\log (\rho )}.$$

$$D_k-\psi (x_k)={\theta }_{k-1}(D_{k-1}-\psi (x_k))\ge 0,D_k-D_{k-1}=(1-{\theta }_{k-1})(\psi (x_k)-D_{k-1})\le 0.$$

$$\psi (x_k)\le D_k.$$

$$\psi (x_k)\le D_{k-1}\le D_{k-2}\le \dots \le D_0=\psi (x_0),$$

$$\Vert \nabla h(x_k)\Vert \le L_0,$$

$$\Vert H(x_k)\Vert =\Vert \nabla h(x_k)\nabla h(x_k{)}^T+{\mu }_{k}I\Vert \le \Vert \nabla h(x_k){\Vert }^2+{\mu }_k\le L_0^2+{\mu }_k,$$

$${\lambda }_{\min }\left(H(x_k{)}^{-1}\right)=\frac{1}{{\lambda }_{\max }\left(H(x_k)\right)}=\frac{1}{\Vert H(x_k)\Vert }\ge \frac{1}{L_0^2+{\mu }_k}.$$

$$\Vert d_k\Vert =\Vert H(x_k{)}^{-1}\nabla \psi (x_k)\Vert \le \Vert H(x_k{)}^{-1}\Vert \Vert \nabla \psi (x_k)\Vert =\frac{1}{{\lambda }_{\min }(H(x_k))}\Vert \nabla \psi (x_k)\Vert \le \frac{1}{{\mu }_k}\Vert \nabla \psi (x_k)\Vert .$$

$${\mu }_k={\xi }_k\Vert h(x_k){\Vert }^{\eta }+{\omega }_k\Vert \nabla h(x_k)h(x_k){\Vert }^{\eta }>{\xi }_{\min }{\epsilon }^{\eta },\text{leavevmode}\text{leavevmode}\text{leavevmode}\forall k\ge 0.$$

$${\mu }_k>\mu ,\text{leavevmode}\text{leavevmode}\text{leavevmode}\forall k\ge 0.$$

$$\Vert d_k\Vert \le {\mu }_k^{-1}\Vert \nabla \psi (x_k)\Vert \le {\mu }_k^{-1}\Vert \nabla h(x_k)\Vert \Vert h(x_k)\Vert \le {\mu }_k^{-1}{L}_0\Vert h(x_0)\Vert \le {\mu }^{-1}{L}_0\Vert h(x_0)\Vert .$$

$$\Vert h(x_k)+\alpha \nabla h(x_k{)}^T{d}_k\Vert \le \Vert h(x_k)\Vert +\alpha \Vert \nabla h(x_k)\Vert \Vert d_k\Vert \le (1+\alpha L_0^2{\mu }^{-1})\Vert h(x_0)\Vert .$$

$$\Vert \nabla h(x_k+t\alpha d_k)-\nabla h(x_k)\Vert \le L\alpha \Vert d_k\Vert \le \alpha L{L}_0{\mu }^{-1}\Vert h(x_0)\Vert .$$

$$\frac{1}{2}\Vert h(x_k)+\alpha \nabla h(x_k{)}^T{d}_k{\Vert }^2$$

$$\le \frac{1}{2}\Vert h(x_k){\Vert }^2+\alpha h(x_k{)}^T\nabla h(x_k{)}^T{d}_k+\frac{1}{2}{\alpha }^2\Vert \nabla h(x_k){\Vert }^2\Vert d_k{\Vert }^2$$

$$\le \frac{1}{2}\Vert h(x_k){\Vert }^2+\alpha h(x_k{)}^T\nabla h(x_k{)}^T{d}_k+\frac{1}{2}{\alpha }^2{L}_0^2\Vert d_k{\Vert }^2$$

Full text

Keywords: Nonlinear equation, Hölder metric subregularity, Non-isolated solutions, Levenberg-Marquardt methods, Global convergence, Worst-case global complexity, Biochemical reaction network kinetics.

AMS subject classifications: 90C26, 68Q25, 65K05

1 Introduction

We consider the problem of finding zeros of the nonlinear mapping $h:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$, i.e.,

[TABLE]

where $h$ is continuously differentiable and satisfies the Hölder metric subregularity (see bellow). The set of zeros of such mappings is denoted by $\Omega$, which is assumed to be nonempty.

A classical approach for finding a solution of (1) is to search for a minimiser of the nonlinear least-squares problem

[TABLE]

where $\|\cdot\|$ denotes the Euclidean norm. In order to guarantee the quadratic or superlinear convergence of many Newton-type schemes for solving (2), the existence of some constant $\beta>0$ satisfying

[TABLE]

is assumed, where $\mathbb{B}(x^{*},\mathtt{r})$ stands for the closed ball centered at $x^{*}$ with radius $\mathtt{r}>0$, cf. [27, 48]. Such an inequality is referred as an error bound (Lipschitzian error bound or metric regularity) condition. The notion of error bound has been very popular during the last few decades to study the local convergence of optimisation methodologies; however, there are many important mappings where (3) is not satisfied , see, e.g., [4, 34]. This motivated the authors of [4] to propose a weaker condition so-called the Hölder metric subregularity (Hölderian error bound), i.e.,

[TABLE]

for $\delta\in{]0,1]}$ and $\mathtt{r}\in{]0,1[}$. There are many mappings satisfying this condition, see, e.g., [4, 34] and references therein. See also Section 5 for a real-world nonlinear system satisfying (4), but not (3).

The Levenberg–Marquardt method is a standard technique used to solve (1), where, in the current point $x_{k}$ and for a positive parameter $\mu_{k}$, the convex subproblem

[TABLE]

is solved to compute a direction $d_{k}$ in which $\nabla h(x_{k})\in\mathbb{R}^{m\times n}$ is the gradient of $h$ at $x_{k}$. This requires finding the unique solution of the linear system

[TABLE]

where $I\in\mathbb{\mathbb{R}}^{m\times m}$ denotes the identity matrix. Then, the next iteration is be generated by $x_{k+1}=x_{k}+d_{k}$ and this scheme is continued until a stationary point of (2) is found, which may correspond to a zero $h$, when certain conditions are satisfied.

The choice of the parameter $\mu_{k}$ has substantial impacts on the global convergence, the local convergence rate, and the computational efficiency of Levenberg–Marquardt methods, cf. [4, 31, 38, 47, 48]. Hence, several ways to specify and to adapt this parameter have been proposed; see, e.g., [18, 19, 48]. A recently proposed Levenberg–Marquardt method by the authors [4] suggests an adaptive parameter $\mu_{k}$ based on the order $\delta\in\ ]0,1]$ of the Hölder metric subregularity (4), i.e.,

[TABLE]

where $\eta\in{]0,4\delta[}$, $\xi_{k}\in[\xi_{\min},\xi_{\max}]$ and $\omega_{k}\in[\omega_{\min},\omega_{\max}]$ with $\xi_{\min}+\omega_{\min}>0$. In [4], this Levenberg–Marquardt method, with adaptive regularisation (LM-AR), was presented and its local convergence was studied for Hölder metrically subregular mappings.

If one assumes that the starting point $x_{0}$ is close enough to a solution $x^{*}$ of (2), then the Levenberg–Marquardt method is known to be quadratically convergent if $\nabla h(x^{*})$ is nonsingular, in which case it is clearly convergent to a solution to (1). In fact, the nonsingularity assumption implies that the solution of the minimisation problem (2) must be locally unique; see [10, 28, 48]. However, assuming local uniqueness of the solution might be restrictive for many applications since the underlying mappings might have non-isolated zeros. Therefore, much attention has been devoted to the study of local convergence of the Levenberg–Marquardt method under local error bounds, which enables the solution of mappings with non-isolated zeros; see, e.g., [10, 18, 19, 48]. In particular, the local convergence of the Levenberg–Marquardt method was studied in [4] under the Hölder metric subregularity condition (4).

As is the case in many applications, one cannot provide a sufficiently close starting point $x_{0}$ to a solution $x^{*}$, and therefore the convergence of the Levenberg–Marquardt method is not guaranteed, which decreases the chance of practical applicability. To overcome this shortcoming, two globalisation techniques have been proposed to be combined with the Levenberg–Marquardt direction, namely, line search and trust-region; see, e.g., [3, 31, 30, 44]. Generally, a line search method finds a descent direction $d_{k}$, specifies a step-size $\alpha_{k}$, generates the new iteration $x_{k+1}=x_{k}+\alpha_{k}d_{k}$, and repeats this scheme until a stopping criterion holds. The step-size is usually determined by an inexact line search such as Armijo, Wolfe, or Goldstein backtracking schemes; see [14, 43]. In particular, the Armijo line search usually finds $\alpha_{k}$ using a backtracking procedure, which ends up with a step-size satisfying

[TABLE]

where $\sigma\in{]0,1[}$. In order to provide an outline for trust-region methods, let us define, firstly, the quadratic function $q_{k}:\mathbb{R}^{m}\rightarrow\mathbb{R}$ with

[TABLE]

Then, a Levenberg–Marquardt trust-region method solves the quadratic subproblem (5) to find a direction $d_{k}$, computes the ratio of the actual reduction to the predicted reduction

[TABLE]

and updates the parameter $\mu_{k}$ using $r_{k}$. For line search and trust-region methods, the global convergence to a first-order stationary point of $\psi$ can be guaranteed, which results to a monotone sequence of function values, i.e., $\psi(x_{k+1})\leq\psi(x_{k})$.

Regardless of the fact that the monotonicity seems natural for the minimisation goal, it has some drawbacks. We address two of them here: (i) the monotone method may lose its efficiency if iterations are trapped at the bottom of a curved narrow valley, where the monotonicity forces the iterations to trace the valley floor (causing very short steps or even an undesired zigzagging); (ii) the Armijo-type line search can break down for very small step-sizes because of rounding errors, when $\psi(x_{k}+\alpha_{k}d_{k})\simeq\psi(x_{k})$. In this case, the point $x_{k}$ may still be far from a stationary point of $\psi$; however, the Armijo rule cannot be satisfied due to indistinguishability of $\psi(x_{k}+\alpha_{k}d_{k})$ from $\psi(x_{k})$ in floating-point arithmetic. To overcome such limitations, the seminal article by Grippo et al. [22] addressed a variant of the Armijo rule (8) by substituting $\psi(x_{k})$ with $\psi_{l(k)}=\max_{0\leq j\leq m(k)}\psi(x_{k-j})$, where $m(0)=0$ and $0\leq m(k)\leq\min\{m(k-1)+1,M\}$, for some nonnegative integer constant $M$ for all $k\geq 1$. This does not guarantee the monotonicity condition $\psi(x_{k+1})\leq\psi(x_{k})$ and therefore called nonmonotone. Nonmonotonicity has also been studied for trust-region methods by replacing $\psi(x_{k})$ with a nonmonotone term; cf. [1]. On the basis of many studies in this area, nonmonotone methods have been recognised to be globally convergent and computationally efficient, even for highly nonlinear problems; see [1, 5] and references therein.

1.1 Motivation and contribution

Our analysis was motivated by the problem of finding moiety conserved steady states of deterministic equations representing the dynamical evolution of molecular species abundance in biochemical reaction networks. This problem can be considered as an application of finding zeros of a mapping $h:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$, that may not satisfy the local error bound (3). It was previously established [4] that this mapping is Hölder metrically subregular and that the merit function is real analytic using standard biochemical assumptions; cf. [7]. Applying a novel Levenberg–Marquardt algorithm with adaptive regularisation (LM-AR), to this problem, we proved local convergence to a zero of $h$ for all such networks if the sufficiently closeness of a starting point to $x^{*}$ can be assumed [4]. However, providing a starting point close enough to $x^{*}$ remains as a limitation in practice, as is the case for all local optimisation methods; see, e.g., [14, 16].

The global convergence and complexity of iterative methods has been the subject of intense debate within the nonlinear optimisation community over the last few decades. While the global convergence guarantees the convergence of the iteration sequence generated by a method for any given starting point $x_{0}$, the worst-case complexity provides an upper bound on the number iterations or function evaluations needed to reach a stationary point of the underlying objective function. These two factors are more important if the convexity or structured nonconvexity of the objective function is assumed; see, e.g., [8, 11, 39, 40, 41, 49]. In the particular case of solving nonlinear least-squares problems by Levenberg-Marquardt methods, there are less results about their global convergence and complexity, compared with the volume of literature concerning Newton-type methods; cf. [45, 46]. This motivates our aim to study the global convergence properties and complexity of two Levenberg-Marquardt methods using line search and trust-region techniques.

We analyse the global convergence, and investigate the complexity of, two Levenberg–Marquardt methods using line search scheme and trust-region globalisation techniques. For the first method, we use $\mu_{k}$ defined in (7), solve the linear system (6) to specify $d_{k}$, and combine this direction with a nonmonotone Armijo-type line search. We also propose a modified version of the Levenberg–Marquardt parameter (7), which is lower bounded, and combines the associated direction $d_{k}$ with a trust-region technique to adapt the Levenberg–Marquardt parameter. A global convergence analysis is provided for both methods. Moreover, we demonstrate that, for both methods, a first-order stationary point is attained after at most $\mathcal{O}(\varepsilon^{-2})$ iterations or function evaluations. We also illustrate some special mappings $h$ where the proposed methods are convergent to a solution to the nonlinear system (1). Finally, we demonstrate that the application of these two methods mappings derived from real-world biochemical reaction networks, from a diverse set of biological species, shows encouraging numerical results in practice. To the best of our knowledge, these two Levenberg–Marquardt methods are the first methods, globally convergent to a stationary point, for finding zeros of the mapping $h$ arising in the study of biological networks. All algorithms are made available within the COBRA Toolbox v.03 [25], an open source software package for modelling biochemical reaction networks.

This paper has five sections, besides this introductory section. Section 2 describes a globally convergent Levenberg–Marquardt line search method. Section 3 addresses a globally convergent Levenberg–Marquardt trust-region method, where in both sections the global convergence and complexity of these methods are analysed. In Section 4, finding a zero of some specific mappings with the proposed methods is discussed. Section 5 reports encouraging numerical results for a mapping appearing in biochemical reaction networks. Finally, conclusions and area for further research are identified in Section 6.

2 Levenberg–Marquardt line search method

For the sake of simplicity, we define $H(x_{k}):=\nabla h(x_{k})\nabla h(x_{k})^{T}+\mu_{k}I$. If $x_{k}$ is not a stationary point of $\psi$, from positive definiteness of $H(x_{k})$, we obtain

[TABLE]

which guarantees the descent property of $d_{k}$ at $x_{k}$. This motivates us to develop a globally convergent Levenberg–Marquardt method using (7). More precisely, we shall combine the Levenberg–Marquardt direction with a nonmonotone Armijo-type line search using the nonmonotone term

[TABLE]

where $\theta_{k-1}\in[\theta_{\min},\theta_{\max}]$ and $0\leq\theta_{\min}\leq\theta_{\max}<1$, cf. [2, 23].

A combination of the direction $d_{k}$ (given by solving (6) using the parameter (7)) with a nonmonotone Armijo-type line search using (12) leads to Algorithm 1.

In order to prove the global convergence of the sequence $\{x_{k}\}$ generated by LMLS to a stationary point of $\psi$, we assume that the next assumptions hold:

(A1)

The mapping $h$ is continuously differentiable and Hölder metrically subregular of order $\delta\in{]0,1]}$ at $(x^{*},0)$; i.e., there exist some constants $\beta>0$ and $\mathtt{r}>0$ such that (4) holds;

(A2)

The lower level set $\mathcal{L}(x_{0}):=\left\{x\in\mathbb{R}^{m}\>|\>\psi(x)\leq\psi(x_{0})\right\}$ is bounded;

(A3)

$\nabla h$ is Lipschitz continuous, i.e.,

[TABLE]

In the subsequent proposition, we first derive a lower bound for the step-size $\alpha_{k}$ and give a bound on the total number of function evaluations needed until the line search (Line 5 of LMLS) is satisfied.

Proposition 2.1**.**

Let $\left\{x_{k}\right\}$ be an infinite sequence generated by LMLS. Then,

(i)

$x_{k}\in\mathcal{L}(x_{0})$;

(ii)

if LMLS does not terminate at $x_{k}$, then

[TABLE]

with

[TABLE]

for $\mu\in{]0,\varepsilon^{\eta}[}$. Moreover, the inner loop of LMLS is terminated in a finite number of steps, denoted by $\ell_{k}$, which satisfies

[TABLE]

?proofname?

.

We prove Assertions (i) and (ii) by induction at the same time. Let us assume $i=1$. Since $D_{0}=\psi(x_{0})$, by the traditional results about the monotone Armijo line search, we have $\psi(x_{1})\leq D_{0}=\psi(x_{0})$. This implies that $x_{1}\in\mathcal{L}(x_{0})$. The proof of Assertion (ii) is similar to $i=k$, i.e., we therefore omit it.

We now assume Assertions (i) and (ii) hold for $i=1,\ldots,k-1$ and prove them for $i=k$. Since $x_{k-1}$ satisfies the line search and $\nabla\psi(x_{k-1})^{T}d_{k-1}<0$, similar to Lemma 2.3 in [2], we can show $\psi(x_{k})\leq D_{k-1}$. This and

[TABLE]

imply $D_{k}\leq D_{k-1}$ and

[TABLE]

Therefore,

[TABLE]

leading to $x_{k}\in\mathcal{L}(x_{0})$, i.e., Assertion (i) holds for $i=k$.

From (A1) and $x_{k}\in\mathcal{L}(x_{0})$, there exists some constant $L_{0}>0$ such that

[TABLE]

which implies

[TABLE]

leading to

[TABLE]

From the definition of $d_{k}$, we obtain

[TABLE]

Since LMLS does not stop at $x_{k}$, it holds $\|h(x_{k})\|>\varepsilon$ and $\|\nabla h(x_{k})h(x_{k})\|>\varepsilon$, which imply

[TABLE]

Let us consider a constant $\mu\in{]0,\xi_{\min}\varepsilon^{\eta}[}$, i.e.,

[TABLE]

We first derive a lower bound on the step-size $\alpha_{k}$. By (20), (18)  and (21), we get

[TABLE]

Therefore, for all $\alpha>0$, we have

[TABLE]

Further, for all $t\in[0,1]$ and $\alpha>0$, (A3) and (22) yield

[TABLE]

It follows from (18) that

[TABLE]

By this inequality, the Taylor expansion of $h(x_{k}+\alpha d_{k})$ around $x_{k}$, and the Cauchy–Schwarz inequality, for any $\alpha>0$, we come to

[TABLE]

This inequality, (18), (23), (24), and (A3) suggest

[TABLE]

From (16), we come to

[TABLE]

For $\alpha=\alpha_{k}/\rho\leq\overline{\alpha}/\rho$, we have

[TABLE]

which yields

[TABLE]

For $\alpha=\alpha_{k}/\rho$, the Armijo-type line search (Line 5 of LMLS) does not hold, i.e.,

[TABLE]

This and the inequality (27) lead to

[TABLE]

Substituting $\alpha=\alpha_{k}/\rho$, we have thanks to (20) and (19) that

[TABLE]

It follows from (28) and (21) that (13) is valid. Using $\alpha_{k}=\rho^{\ell_{k}}\overline{\alpha}$ and (13), we end up to

[TABLE]

which proves (14). ∎

The first main result of this section demonstrates some properties of the sequence $D_{k}$ and shows that any accumulation point of the sequence $\left\{x_{k}\right\}$ generated by LMLS is either a solution of (1) or a stationary point of $\psi$.

Theorem 2.2**.**

Let $\left\{x_{k}\right\}$ be an infinite sequence generated by LMLS. Then, for all $k\geq 0$, the following assertions hold:

(i)

$\left\{D_{k}\right\}$* is convergent and*

[TABLE]

(ii)

$\nabla\psi(x_{k})^{T}d_{k}\leq-c_{1}\|\nabla\psi(x_{k})\|^{2}$;

(iii)

LMLS either stops at finite number of iterations, satisfying $\|h(x_{k})\|\leq\varepsilon$ or $\|\nabla\psi(x_{k})\|\leq\varepsilon$, or generates an infinite sequence $\left\{x_{k}\right\}$ such that any accumulation point of this sequence is a stationary point of the merit function $\psi$, i.e.,

[TABLE]

?proofname?.

From (15), we have $D_{k}\leq D_{k-1}$. This, Proposition 2.1 (i) and (A2) imply that $\left\{D_{k}\right\}$ is convergent. Further, since $\theta_{k-1}\in[\theta_{\min},\theta_{\max}]$, with $\theta_{\max}\in[\theta_{\min},1[$, we have $1-\theta_{k-1}\geq 1-\theta_{\max}>0$. Taking limits from both sides of $D_{k}-D_{k-1}=(1-\theta_{k-1})(\psi(x_{k})-D_{k-1})$ when $k$ goes to infinity, we deduce (29).

It follows from (19) and the definition of $d_{k}$ that

[TABLE]

By the definition of $\mu_{k}$, (17), and (18), we get

[TABLE]

for all $k\geq 0$. This and (31) yield

[TABLE]

that is, Assertion (i) holds with $c_{1}:=\left(L_{0}^{2}+\overline{\mu}\right)^{-1}>0$.

Let us now prove the assertion (iii). If the algorithm stops in a finite number of iterations by either $\|h(x_{k})\|\leq\varepsilon$ or $\nabla\psi(x_{k})\leq\varepsilon$, the result is valid. Let us assume that the algorithm generates an infinite sequence $\left\{x_{k}\right\}$. For a fixed iteration $x_{k}$, the stopping criteria of LMLS do not hold, i.e., $\|h(x_{k})\|>\varepsilon$ and $\|\nabla h(x_{k})h(x_{k})\|>\varepsilon$. Therefore, from (21), we have

[TABLE]

It can be deduced from Line 5 of LMLS and Assertion (ii) that

[TABLE]

This and the assertion (ii) yield

[TABLE]

i.e., any accumulation point of $\left\{x_{k}\right\}$ is a stationary point of $\psi$. ∎

We continue the analysis of LMLS by providing the worst-case global and evaluation complexities of LMLS, which are upper bounds on the number of iterations and merit function evaluations required to get an approximate stationary point of $\psi$ satisfying $\|\nabla\psi(x)\|\leq\varepsilon$, for the accuracy parameter $\varepsilon$, respectively. Let us denote by $N_{i}(\varepsilon)$ and $N_{f}(\varepsilon)$ the total number of iterations and merit function evaluations of LMLS required to find and an $\varepsilon$-stationary point of (2) .

Theorem 2.3**.**

Let $\left\{x_{k}\right\}$ be the sequence generated by LMLS and (A2) and (A3) hold. Then,

(i)

the total number of iterations to guarantee $\|\nabla\psi(x_{k})\|\leq\varepsilon$ is bounded above and

[TABLE]

with $c_{2}:=c_{1}\sigma\widehat{\alpha}(1-\theta_{\max})$;

(ii)

the total number of function evaluations to guarantee $\|\nabla\psi(x_{k})\|\leq\varepsilon$ is bounded above and

[TABLE]

?proofname?.

To prove Assertion (i), we define

[TABLE]

which suggests

[TABLE]

Let us assume by contradiction that $N_{i}(\varepsilon)>\widehat{k}$, which means that the algorithm does not stop in $\widehat{k}$ iterations. From Line 5 of LMLS, (15), and Theorem 2.2 (ii), we obtain

[TABLE]

leading to

[TABLE]

which contradicts to (36). Therefore, (34) is valid.

Considering the bound on the number of merit function evaluations in step $k$ ($\ell_{k}$, given in Proposition 2.1), the following upper bound on the total number of merit function evaluations can be provided by

[TABLE]

giving the results. ∎

Theorem 2.3 implies that the worst-case global and evaluation complexities of LMLS to attain the approximate stationary point of $\psi$ are of the order $\mathcal{O}(\varepsilon^{-2})$, which is the same as the gradient method; see, e.g., [39]. However, in practice Levenberg-Marquardt methods usually performs much better than the gradient method.

Let us compute here the second derivative of $\psi$ at $x$, i.e.,

[TABLE]

where $S(x):=\sum_{i=1}^{m}h_{i}(x)\nabla^{2}h_{i}(x)$. Three types of the problem (1) are recognised with respect to the magnitude of $\|h(x^{*})\|$: (i) if $h(x^{*})=0$, the problem is called zero residual; (ii) if $\|h(x^{*})||$ is small, the problem is called small residual; and if $\|h(x^{*})\|$ is large, the problem is called large residual; see, e.g., [16]. Under the nonsingularity assumption of $\nabla h(x)$ at the limit point $x^{*}$ of $\{x_{k}\}$ and using (37), we investigate the superlinear convergence of $\left\{x_{k}\right\}$ generated by LMLS for zero residual problems, which is the same as the convergence rate given for quasi-Newton methods; see [15].

Theorem 2.4**.**

Let $\psi:\mathbb{R}^{m}\rightarrow\mathbb{R}$ be twice continuously differentiable on $\mathcal{L}(x_{0})$, and $\left\{x_{k}\right\}$ be the sequence generated by LMLS and (A1)-(A3) hold. If $\left\{x_{k}\right\}$ converges to $x^{*}$ and $\nabla h(x^{*})$ has full rank, then

[TABLE]

there exists $\overline{k}\geq 0$ such that $\alpha_{k}=1$ for all $k\geq\overline{k}$, and $\left\{x_{k}\right\}$ converges to $x^{*}$ superlinearly.

?proofname?.

Since $\nabla h(x^{*})$ has full rank, (30) implies $h(x^{*})=0$. This and (37) yield that $\nabla^{2}\psi(x^{*})$ is positive definite. Hence, $h(x^{*})=0$ leads to

[TABLE]

From (37), we obtain

[TABLE]

which implies

[TABLE]

Since $\psi$ is twice continuously differentiable in the compact set $\mathcal{L}(x_{0})$, $\|\nabla^{2}h_{i}(x)\|$ ($i=1,\ldots,m$) is bounded. This, the last inequality, and (39) give

[TABLE]

giving (38).

From Theorem 6.4 in [15] and (38), we have that (8) is valid with $\alpha_{k}=1$, for all $k$ sufficiently large. Therefore, the superlinear convergence of $\{x_{k}\}$ follows from Theorem 3.1 in [15]. ∎

3 Levenberg–Marquardt trust-region method

Trust-region methods are known to be effective for nonconvex optimisation problems (see [14]). Therefore, this section concerns with the development of a globally convergent Levenberg–Marquardt method using a trust-region technique and the investigation on its convergence analysis and complexity.

Let us start with some details of a trust-region globalisation technique that will be coupled with the Levenberg–Marquardt direction. We first draw your attention to some literature, e.g., [1, 3] and references therein, about the efficiency of nonmonotone trust-region methods compared to monotone ones for either optimisation or nonlinear systems. This motivates us to develop a nonmonotone Levenberg–Marquardt trust-region method for solving systems of nonlinear equations. To do so, we take advantage of the quadratic function $q_{k}$ (9) and define the ratio

[TABLE]

where the nonmonotone term $D_{k}$ defined by (12). In this ratio, the nominator is called nonmonotone reduction and the denominator is called the predicted reduction. Further, let us introduce a new Levenberg–Marquardt parameter that is a modified version of (7), i.e.,

[TABLE]

where $\mu_{k}$ is given by (7) with $\eta\in{]0,4\delta[}$, $\mu_{min}>0$, $\xi_{k}\in[\xi_{\min},\xi_{\max}]$, $\omega_{k}\in[\omega_{\min},\omega_{\max}]$ with $\xi_{\min}+\omega_{\min}>0$, and $\lambda_{k}$ is updated by

[TABLE]

in which $0<\rho_{2}<1<\rho_{1}$ and $0<\upsilon_{1}<\upsilon_{2}<1$ are some constants. A simple comparison between (7) and (41) indicates that $\widehat{\mu}_{k}$ is lower bounded and $\lambda_{k}$ helps to have a better control on the Levenberg–Marquardt parameter, which shows its effect on numerical performance of the method (see Section 5 for more details).

In our Levenberg–Marquardt trust-region method, we first determine $\widehat{\mu}_{k}$ (41), specify the direction $d_{k}$ by solving the linear system (6), and compute the ratio $\widehat{r}_{k}$ (40). If $\widehat{r}_{k}\geq\upsilon_{1}$, the trial point $d_{k}$ is accepted, i.e., $x_{k+1}=x_{k}+d_{k}$; otherwise, the parameter $\lambda_{k}$ should be increased by setting $\lambda_{k}=\rho_{1}\lambda_{k}$. In the case that $\widehat{r}_{k}\geq\upsilon_{2}$, the parameter $\lambda_{k}$ is decreased by setting $\lambda_{k+1}=\rho_{2}\lambda_{k}$. The final step will be the evaluation of stopping criteria, which here is either $\|h(x_{k+1})\|\leq\varepsilon$ or $\|\nabla\psi(x_{k+1})\|\leq\varepsilon$. We summarise this scheme in Algorithm 2.

In LMTR, the loop starts from Line 5 to Line 7 is called the inner loop and the loop starts from Line 3 to Line 14 is called the outer loop.

The subsequent proposition points out that the inner loop of LMTR is terminated after a finite number of steps and provides upper bounds for $\widehat{\mu}_{k}$ and $p_{k}$.

Proposition 3.1**.**

Let $\left\{x_{k}\right\}$ be an infinite sequence generated by LMTR and (A1)-(A3) holds. Then,

(i)

$q_{k}(0)-q_{k}(d_{k})\geq\frac{1}{2(L_{0}^{2}+\mu_{min})}\|\nabla\psi(x_{k})\|^{2}$;

(ii)

$q_{k}(0)-q_{k}(d_{k})\leq(\frac{1}{2}L_{0}^{2}+\widehat{\mu}_{k})\|d_{k}\|^{2}$;

(iii)

the inner loop is terminated in a finite number of steps. Moreover, if LMTR does not terminate at $x_{k}$, then

[TABLE]

with

[TABLE]

and

[TABLE]

?proofname?.

By the definition of $q_{k}$ in (9) and (19), we get

[TABLE]

giving Assertion (i).

It follows from (9) that

[TABLE]

proving Assertion (ii).

For the first part of the assertion (iii), we show that the inner loop is terminated after a finite number of steps. From Assertion (i) and (20), we obtain

[TABLE]

By (A3) and (20), for $t\in[0,1]$, we get

[TABLE]

By the Taylor expansion of $\psi(x_{k}+d_{k})$ around $x_{k}$, we come to

[TABLE]

From this, (23), and (47), it consequently holds

[TABLE]

Since $\widehat{r}_{k-1}\geq\upsilon_{1}$, we have $\psi(x_{k})\leq D_{k-1}$. This and

[TABLE]

imply $D_{k}\leq D_{k-1}$ and $\psi(x_{k})\leq D_{k}$, leading to $\psi(x_{k})\leq D_{k}\leq D_{k-1}\leq\ldots\leq D_{0}=\psi(x_{0})$, i.e.,

[TABLE]

It can be deduced from this and (49) that

[TABLE]

where $\widetilde{c}_{0}:=\frac{1}{2}L^{2}L_{0}^{2}\|h(x_{0})\|^{2}$, $\widetilde{c}_{1}:=L_{0}^{2}\|h(x_{0})\|$, and $\widetilde{c}_{2}:=\|h(x_{0})\|$. For sufficiently large $p_{k}$, we have $\widehat{\mu}_{k}=\rho_{1}^{p_{k}}\lambda_{k}\mu_{k}$. This, (10), and (46) yield

[TABLE]

It can be deduced from this and $\psi(x_{k})\leq D_{k}$ that $\widehat{r}_{k}\geq r_{k}\geq\upsilon_{1}$, for sufficiently large $p_{k}$, proving the first part of Assertion (iii).

In the second part of Assertion (iii), we provide upper bounds for $\widehat{\mu}_{k}$ and $p_{k}$. Let us denote by $\overline{d}_{k}$ the solution of the system (6) corresponding to the parameter $\overline{\mu}_{k}:=\rho_{1}^{p_{k}-1}\lambda_{k}\mu_{k}$ and set $\overline{x}_{k+1}=x_{k}+\overline{d}_{k}$. By (20) and (18), we get

[TABLE]

It follows from this and the triangle inequality that

[TABLE]

For all $t\in[0,1]$, (A3) and (53) imply

[TABLE]

From (18), (2), (48), and (11), we obtain

[TABLE]

Following $\psi(x_{k})\leq D_{k}$ and $\nabla\psi(x_{k})^{T}\overline{d}_{k}=-\overline{d}_{k}^{T}H_{k}\overline{d}_{k}\leq-\overline{\mu}_{k}\|\overline{d}_{k}\|^{2}$, it can be deduced

[TABLE]

It follows from (45) and the definition $\overline{\mu}_{k}$ that $\widehat{r}_{k}<\upsilon_{1}$ and

[TABLE]

Combining this inequality with that in (56) suggest

[TABLE]

leading to

[TABLE]

giving (42). Since $\widehat{\mu}_{k}=\rho_{1}^{p_{k}}\lambda_{k}\mu_{k}$, taking the logarithm from both sides of

[TABLE]

implies (43), completing the proof. ∎

We now draw your attention to the global convergence of the sequence $\left\{x_{k}\right\}$ generated by LMTR to a first-order stationary point $x^{*}$ of $\psi$ satisfying $\nabla\psi(x^{*})=0$. Let us first recall the following result for local convergence of the Levenberg–Marquardt method given in [4].

Theorem 3.2**.**

Let $\left\{x_{k}\right\}$ be the sequence generated by LMTR and (A1)-(A3) hold. Then, $\left\{D_{k}\right\}$ is convergent and

[TABLE]

Further, the algorithm either stops at finite number of iterations, satisfying $\|h(x_{k})\|\leq\varepsilon$ or $\nabla\psi(x_{k})\leq\varepsilon$, or generates an infinite sequence $\left\{x_{k}\right\}$ such that any accumulation point of this sequence is a stationary point of the merit function $\psi$, i.e.,

[TABLE]

?proofname?.

From (50) and (51), we have $D_{k}\leq D_{k-1}$ and $x_{k}\in\mathcal{L}(x_{0})$. Hence, the sequence $\{D_{k}\}$ is decreasing and bounded below, i.e., it is convergent. From $\theta_{k}\in[\theta_{\min},\theta_{\max}]$, with $\theta_{\max}\in[\theta_{\min},1[$, we obtain $1-\theta_{k}\geq 1-\theta_{\max}>0$. Taking limits when $k$ goes to infinity from $D_{k}\leq\psi(x_{k})\leq D_{k+1}$ gives (57).

If the algorithm stops in a finite number of iterations by either $\|h(x_{k})\|\leq\varepsilon$ or $\nabla\psi(x_{k})\leq\varepsilon$, the result is valid. If the algorithm generates the infinite sequence $\left\{x_{k}\right\}$, Proposition 3.1 (i) yields

[TABLE]

From this and (57), we obtain

[TABLE]

i.e., any accumulation point of $\left\{x_{k}\right\}$ is a stationary point of $\psi$. ∎

Let us continue this section by providing global and evaluation complexities of the sequence $\left\{x_{k}\right\}$ generated by LMTR using the results presented in Proposition 3.1.

Theorem 3.3**.**

Let $\left\{x_{k}\right\}$ be the sequence generated by LMTR and (A1)-(A3) hold. Then,

(i)

the total number of iterations to guarantee $\|\nabla\psi(x_{k})\|\leq\varepsilon$ is bounded above by

[TABLE]

where $\overline{c}_{3}:=\upsilon_{1}(1-\eta_{\max})/(2(L_{0}^{2}+\mu_{min}))$;

(ii)

the total number of function evaluations to guarantee $\|\nabla\psi(x_{k})\|\leq\varepsilon$ is bounded above by

[TABLE]

?proofname?.

To prove Assertion (i), we first define

[TABLE]

which is equivalent to

[TABLE]

Let us assume by contradiction that $N_{i}(\varepsilon)>\widetilde{k}$, which means that LMTR does not stop in $\widetilde{k}$ iterations. For a successful iteration $k$ of LMTR, it follows from (15) and Proposition 3.1 (i) that

[TABLE]

leading to

[TABLE]

which contradicts to (61), proving Assertion (i).

From (43) and (59), we obtain

[TABLE]

giving (60). ∎

We conclude this section by providing the local convergence rate of LMTR if the corresponding sequence $\{x_{k}\}$ is convergent to a solution of (1) under the  Lojasiewicz gradient inequality (see [35, 36]). To this end, the presence of the subsequent two lemmas are necessary in our local analysis of LMTR.

Lemma 3.4**.**

[7, Lemma 1]** Let $\left\{s_{k}\right\}$ be a sequence in $\mathbb{R}_{+}$ and let $\zeta,\nu$ be some nonnegative constants. Suppose that $s_{k}\to 0$ and that the sequence satisfies

[TABLE]

for all $k$ sufficiently large. Then

1. (i)

if $\zeta=0$, the sequence $\left\{s_{k}\right\}$ converges to [math] in a finite number of steps; 2. (ii)

if $\zeta\in\left]0,1\right]$, the sequence $\left\{s_{k}\right\}$ converges linearly to [math] with rate $1-\frac{1}{\nu}$; 3. (iii)

if $\zeta>1$, there exists $\varsigma>0$ such that, for all $k$ sufficiently large,

[TABLE]

Lemma 3.5**.**

[29, Theorem 2.5 and Lemma 2.3]** The sequence $\{x_{k}\}$ generated by LMTR with $p=0$ satisfies

[TABLE]

and

[TABLE]

Let us describe now the Łojasiewicz gradient inequality in the following definition.

Definition 3.6**.**

Let $\psi:U\to\mathbb{R}$ be a function defined on an open set $U\subseteq\mathbb{R}^{m}$, and assume that the set of zeros $\Omega:=\left\{x\in\mathbb{R}^{m},\,\psi(x)=0\right\}$ is nonempty. The function $\psi$ is said to satisfy the Łojasiewicz gradient inequality if for any critical point $\overline{x}$, there exist constants $\kappa>0,\varepsilon>0$ and $\theta\in[0,1[$ such that

[TABLE]

This inequality is valid for a large class of functions such as analytic, subanalytic, and semialgebraic functions, cf. [35, 36, 33]. See Section 5 for a mapping with a real analytic merit function, where finding zeros of this mapping is the main motivation of this study. Here, we further assume that

(A4)

the merit function $\psi$ satisfies the Łojasiewicz gradient inequality (64).

The next theorem is the third main result of this section, which provides the convergence of the sequences $\left\{\mathrm{dist}(x_{k},\Omega)\right\}$ and $\left\{\psi(x_{k})\right\}$ to [math] if an accumulation point $x^{*}$ of $\{x_{k}\}$ is a solution of the nonlinear system (1).

Theorem 3.7**.**

Suppose that (A4) holds and assume that the sequence $\{x_{k}\}$ generated by LMTR is convergent to a solution $x^{*}$ of the nonlinear system (1). Then,

1. (i)

for sufficiently large $k$, it holds $x_{k+1}=x_{k}+d_{k}$; 2. (ii)

there exist constants $s>0$, $\overline{s}>0$, and $k^{\prime}\in\mathbb{N}$ such that, for $x_{k^{\prime}}\in\mathbb{B}(x^{*},s)$,

[TABLE] 3. (iii)

if $\theta=0$, the sequences $\{\psi(x_{k})\}$ and $\{\mathrm{dist}(x_{k},\Omega)\}$ converge to [math] in a finite number of steps; 4. (iv)

if $\theta\in\left]0,\frac{1}{2}\right]$, the sequences $\{\psi(x_{k})\}$ and $\{\mathrm{dist}(x_{k},\Omega)\}$ converge linearly to [math]; 5. (v)

if $\theta\in\left]\frac{1}{2},1\right[$, there exist some positive constants $\varsigma_{1}$ and $\varsigma_{2}$ such that, for all large $k$,

[TABLE]

?proofname?.

Since $x^{*}$ is an accumulation point of $\left\{x_{k}\right\}$ and a solution of the nonlinear system (1), it can be deduced

[TABLE]

From this, Proposition 3.1 (i), Proposition 2.1 (ii), (49), $\widehat{\mu}_{k}\geq\mu_{min}$, and (65), we obtain

[TABLE]

which implies that there exists a $\overline{k}_{0}\in\mathbb{N}$ such that $r_{k}\geq\upsilon_{1}$. Hence, for all $k\geq\overline{k}_{0}$, it follows from $D_{k}\geq\psi(x_{k})$ that

[TABLE]

which means that $x_{k+1}=x_{k}+d_{k}$ with $p_{k}=0$ that justifies Assertion (i).

We divide the proof of Assertion (ii) into three parts. First, we will provide the values of $s$ and $\overline{s}$. Let us set $\varepsilon>0$ and $\kappa>0$ such that (64) holds and let $\overline{s}:=\min\{\mathtt{r},\varepsilon\}>0$. By the definition of $\widehat{\mu}_{k}$, (A2), and (42), we get

[TABLE]

By making $\overline{s}$ smaller if needed, we can guarantee

[TABLE]

Lipschitz continuity of $h$ and $\overline{s}\leq\mathtt{r}<1$, for all $x\in\mathbb{B}(x^{*},\overline{s})$, lead to

[TABLE]

We now define

[TABLE]

From $\overline{s}<1$ and $\theta\in[0,1[$, we obtain $s\leq\overline{s}$.

For $k^{\prime}\geq\overline{k}_{0}$, let us choose any $x_{k^{\prime}}\in\mathbb{B}(x^{*},s)$. Lemma 3.5 and $d_{k}=-H_{k}^{-1}\nabla h(x_{k})h(x_{k})$ imply, for all $k\in\mathbb{N}$,

[TABLE]

Next, let us show by induction that, for $i\in\mathbb{N}$,

[TABLE]

It follows from $x_{k^{\prime}}\in\mathbb{B}(x^{*},\overline{s})$ and (67) that

[TABLE]

leading to

[TABLE]

Then, from (69), one can deduce

[TABLE]

From the convexity of the function $\varphi(t):=-t^{1-\theta}$ with $t>0$, we come to

[TABLE]

This and (71) suggest

[TABLE]

It follows from $x_{0}\in\mathbb{B}(x^{*},s)\subseteq\mathbb{B}(x^{*},\overline{s})$ and (66) that $\|H_{0}\|\leq(L_{0}^{2}+\tau)$. Hence, by the Łojasiewicz gradient inequality (64), we get

[TABLE]

This, (73), and (68) yield

[TABLE]

which, proves the second assertion in (70) for $i=1$. Then, we have

[TABLE]

implying $x_{k^{\prime}+1}\in\mathbb{B}(x^{*},\overline{s})$. Now, let us assume that (70) holds for all $i=1,\ldots,k$. From $x_{k}\in\mathbb{B}(x^{*},\overline{s})$ and (67), it can be deduced

[TABLE]

leading to

[TABLE]

It follows from this and (69) that

[TABLE]

A combination of this inequality and (72) leads to

[TABLE]

Further, from $x_{k^{\prime}+k}\in\mathbb{B}(x^{*},\overline{s})$, (64), and (66), we obtain

[TABLE]

By the latter inequality and (75), we come to

[TABLE]

proving the second assertion in (70) for $i=k+1$. Then, it follows from (68) that

[TABLE]

Hence, the first assertion in (70) is valid for $i=k+1$.

Finally, we are in a position to show that Assertions (ii) is true. As shown in (70), $x_{k}\in\mathbb{B}(x^{*},\overline{s})$ for all $k\geq k^{\prime}$. This and (66), implies that $\|H_{k}\|\leq(L_{0}^{2}+\tau)$ for all $k\geq k^{\prime}$. Hence, for $k\geq k^{\prime}$, we have

[TABLE]

Then, by (74), we get

[TABLE]

From this and (64), it can be deduced

[TABLE]

which implies that $\{\psi(x_{k})\}$ converges to [math]. This ans the Hölder metric subregularity validate the statement of the assertion (ii).

Applying Lemma 3.4 with $s_{k}:=\psi(x_{k})$, $\nu:=2\kappa^{2}(L_{0}^{2}+\tau)$ and $\zeta:=2\theta$, we have that the convergence rate are dependent to $\theta$ as claimed in Assertions (iii)-(v). Therefore, the Hölder metric subregularity of $h$ implies that $\{\mathrm{dist}\left(x_{k},\Omega\right)\}$ converges to [math] with the rate given in (iii)-(v). ∎

4 Convergence to a solution of nonlinear systems

Let us emphasis that the algorithms LMLS and LMTR only guarantee the convergence of the sequence $\{x_{k}\}$ to a stationary point $x^{*}$ of the merit function $\psi$, which can be a local non-global minimiser of (2), i.e.,

[TABLE]

Therefore, the remainder of this section concerns with considering more restrictions on the mapping $h$ such that the global convergence of $\{x_{k}\}$ to a solution of (1) is guaranteed.

The next theorem extracts some classical results for cases that $\nabla h(x^{*})$ is nonsingular, which implies that $x^{*}$ is a solution of (1). Moreover, the worst-case global and evaluation complexities to attain solution of (1) are provided under the nonsingularity of $\nabla h(x)\nabla h(x)^{T}$ for all $x\in\mathcal{L}(x_{0})$. Under the assumption that all accumulation points of $\left\{x_{k}\right\}$ are solutions of (1) and $\nabla h(x^{*})$ is nonsingular for the accumulation point $x^{*}$, it is proved that the whole sequence $\left\{x_{k}\right\}$ converges to the isolated solution $x^{*}$ of (1).

Theorem 4.1**.**

Let $\left\{x_{k}\right\}$ be the sequence generated by LMLS or LMTR and (A1)-(A3) hold. Then

(i)

if $\nabla h(x^{*})$ is nonsingular at any accumulation point $x^{*}$ of $\left\{x_{k}\right\}$, then $x^{*}$ is a solution of the nonlinear system (1).

(ii)

if the matrix $\nabla h(x)\nabla h(x)^{T}$ is nonsingular for all $x\in\mathcal{L}(x_{0})$, i.e., there exists $\lambda>0$ such that $\lambda_{\min}(\nabla h(x)\nabla h(x)^{T})>\lambda$, then, for LMLS,

[TABLE]

and

[TABLE]

and, for LMTR,

[TABLE]

and

[TABLE]

(iii)

*if all accumulation points of $\left\{x_{k}\right\}$ are solutions of the nonlinear system (1), *$x^{*}$is an accumulation point of $\left\{x_{k}\right\}$ such that $\nabla h(x^{*})$ is nonsingular, and

[TABLE]

then $\left\{x_{k}\right\}$ converges to $x^{*}$.

?proofname?.

For any accumulation point $x^{*}$of $\left\{x_{k}\right\}$, it follows from Theorem 3.2 that $\nabla h(x^{*})h(x^{*})=0$. This, along with the nonsingularity of $\nabla h(x^{*})$, implies Assertion (i).

To prove Assertion (ii), we note that

[TABLE]

i.e., $\|h(x_{k})\|\leq\lambda^{-1/2}\|\nabla\psi(x_{k})\|$. This and Proposition 2.3 (i)-(ii) give (76) and (77), respectively. Similarly, (78) and (79) follow from this inequality and Proposition 3.3 (i)-(ii).

In order to prove Assertion (iii), let us assume that all accumulation points of $\left\{x_{k}\right\}$ are solutions of (1), $x^{*}$is an accumulation point such that $\nabla h(x^{*})$ is nonsingular, and (80) holds. From the inverse function theorem and the nonsingularity of $\nabla h(x^{*})$, there exists a neighborhood around 0 such that $h$ is invertible. Therefore, there exists a neighborhood $\mathbb{B}(x^{*},\mathtt{r}_{1})$ for $\mathtt{r}_{1}>0$ such that

[TABLE]

implying

[TABLE]

Since $x^{*}$is an accumulation point of $\left\{x_{k}\right\}$, $\mathbb{B}(x^{*},r_{1})$ contains an infinite number of iteration points of $\left\{x_{k}\right\}$. It remains to show that there exists $k_{2}\in\mathbb{N}$ such that $x_{k}\in\mathbb{B}(x^{*},r_{1})$, for all $k\geq k_{2}$. Hence, for an arbitrary $\varepsilon\in(0,r_{1})$, the set $\mathbb{B}(x^{*},r_{1})-\mathbb{B}(x^{*},\varepsilon)$ involves only a finite number of iterations of $\left\{x_{k}\right\}$, i.e., there exists $k_{3}\in\mathbb{N}$ such that

[TABLE]

It follows from (80) that there exists $k_{4}\in\mathbb{N}$ such that

[TABLE]

Let us set $k_{2}:=\max\left\{k_{3},k_{4}\right\}$ leading to

[TABLE]

giving the result. ∎

The mapping $h$ is called strictly monotone on $\mathbb{R}^{m}$ if

[TABLE]

In addition, the mapping $h$ is called strictly duplomonotone with constant $\overline{\tau}>0$ if

[TABLE]

whenever $h(x)\neq 0$; see [6, 44]. In the next result, we will show that if the mapping $h$ or $-h$ is strictly monotone (duplomonotone), then the sequence $\left\{x_{k}\right\}$ generated by LMLS converges to the unique solution of the nonlinear system (1).

Theorem 4.2**.**

Let $\left\{x_{k}\right\}$ be the sequence generated by LMLS or LMTR and (A1)-(A3) hold.

(i)

If the mapping $h$ or $-h$ is strictly monotone, then $\left\{x_{k}\right\}$ converges to the unique solution of the nonlinear system (1).

(ii)

If the mapping $h$ or $-h$ is strictly duplomonotone, then $\left\{x_{k}\right\}$ converges to a solution of the nonlinear system (1).

?proofname?.

In order to prove Assertion (i), let $h$ or $-h$ be strictly monotone and $x^{*}$ be an accumulation point of $\left\{x_{k}\right\}$. If $h$ is strictly monotone, for the points $x$ and $x+th$ with $t>0$ and $h\in\mathbb{R}^{m}$, we can deduce

[TABLE]

If $-h$ is strictly monotone, then

[TABLE]

By setting $z=h(x^{*})$ and $x=x^{*}$ in the last two inequalities, we get

[TABLE]

This, (30), and (58) imply $h(x^{*})=0$.

To prove Assertion (ii), let $h$ be strictly duplomonotone, which leads to

[TABLE]

If $-h$ is strictly duplomonotone, then

[TABLE]

The result follows from the last two inequalities at $x=x^{*}$, (30), and (58). ∎

Note that the strict monotonicity of $h$ does not implies the positive definiteness of $\nabla h(x)$. Therefore the results of Theorem 4.2 (i) is not a trivial consequence of Theorem 4.1 (i).

5 Application to biochemical reaction networks

In this section, we use the following notation: $\mathbb{Z}_{+}^{m\times n}:=\left\{A\in\mathbb{Z}^{m\times n}\leavevmode\nobreak\ |\leavevmode\nobreak\ A_{ij}\geq 0,\leavevmode\nobreak\ i=1,\ldots,m,\leavevmode\nobreak\ j=1,\ldots,n\right\}$, $\mathbb{R}_{+}^{m}:=\left\{a\in\mathbb{R}^{m}\leavevmode\nobreak\ |\leavevmode\nobreak\ a_{i}\geq 0,\leavevmode\nobreak\ i=1,\ldots,m\right\}$, and $\mathbb{R}_{++}^{m}:=\left\{a\in\mathbb{R}^{m}\leavevmode\nobreak\ |\leavevmode\nobreak\ a_{i}>0,\leavevmode\nobreak\ i=1,\ldots,m\right\}$. Let us consider a biochemical reaction network with $m$ molecular species and $n$ reversible elementary reactions111An elementary reaction is a chemical reaction for which no intermediate molecular species need to be postulated in order to describe the chemical reaction on a molecular scale.. We define forward and * reverse stoichiometric matrices*, $F,R\in\mathbb{\mathbb{Z}}_{+}^{m\times n}$, respectively, where $F_{ij}$ denotes the stoichiometry222Reaction stoichiometry is a quantitative relationship between the relative quantities of molecular species involved in a single chemical reaction. of the $i^{th}$ molecular species in the $j^{th}$ forward reaction and $R_{ij}$ denotes the stoichiometry of the $i^{th}$ molecular species in the $j^{th}$ reverse reaction. We assume that every reaction conserves mass, i.e., there exists at least a positive vector $l\in\mathbb{R}_{++}^{m}$ such that $(R-F)^{T}l=0$; cf. [21]. The matrix $N:=R-F$ represents net reaction stoichiometry and may be viewed as an incidence matrix of a directed hypergraph; see [32]. In practice, there are less molecular species than net reactions ($m<n$). We assume the cardinality of each row of $F$ and $R$ is at least one, and the cardinality of each column of $R-F$ is at least two. The matrices $F$ and $R$ are sparse and the sparsity pattern depends on the particular biochemical reaction network being modeled. It is here assumed that $\text{rank}([F,R])=m$, which is a requirement for kinetic consistency; cf. [20].

Let $c\in\mathbb{R}_{++}^{m}$ be a vector of molecular species concentrations. For nonnegative elementary kinetic parameters $k_{f},k_{r}\in\mathbb{R}_{+}^{n}$, elementary reaction kinetics for forward and reverse elementary reaction rates as $s(k_{f},c):=\exp(\ln(k_{f})+F^{T}\ln(c))$ and $r(k_{r},c):=\exp(\ln(k_{r})+R^{T}\ln(c))$, respectively, where $\exp(\cdot)$ and $\ln(\cdot)$ denote the respective componentwise functions; see, e.g., [7, 20]. Then, the system of differential equations

[TABLE]

shows the deterministic dynamical equation for time evolution of molecular species concentration. A vector $c^{*}$ is called a steady state if and only if $f(c^{*})=0$. Hence, $c^{*}$ is a steady state of the biochemical system if and only if

[TABLE]

where $\mathcal{N}(N)$ stands for the null space of $N$. The set of steady states $\Omega_{1}=\left\{c\in\mathbb{R}_{++}^{m},\,f(c)=0\right\}$ will be unchanged if $N$ is replaced by a matrix $\bar{N}$ with the same kernel. Suppose that $\bar{N}\in\mathbb{Z}^{r\times n}$ is the submatrix of $N$ whose rows are linearly independent, then $\mathrm{rank}\left(\bar{N}\right)=\mathrm{rank}(N)\eqqcolon r.$ If one replaces $N$ by $\bar{N}$ and transforms (5) into logarithmic scale, by letting $x\coloneqq\ln(c)\in\mathbb{R}^{m}$, $k\coloneqq[\ln(k_{f})^{T},\,\ln(k_{r})^{T}]^{T}\in\mathbb{R}^{2n}$, then the right-hand side of (5) can be translated to

[TABLE]

where $\left[\,\cdot\thinspace,\cdot\,\right]$ stands for the horizontal concatenation operator.

Let $L\in\mathbb{R}^{m-r,m}$ be a basis for the left nullspace of $N$, i.e., $L^{T}N=0$, where $\mathrm{rank}(N)=r$ and $\mathrm{rank}(L)=m-r$. The system satisfies moiety conservation if for any initial concentration $c_{0}\in\mathbb{R}_{++}^{m}$, it holds

[TABLE]

where $l_{0}\in\mathbb{R}_{++}^{m}.$ It is possible to compute $L$ such that each corresponds to a structurally identifiable conserved moiety in a biochemical reaction network; cf. [24]. Therefore, finding the moiety conserved steady state of a biochemical reaction network is equivalent to finding a zero of the mapping

[TABLE]

It was shown by the authors in Section 4.1 of [4] that the merit function $\psi$ satisfies Łojasiewicz gradient inequality (with an exponent $\theta\in[0,1[$) and the mapping $h$ is Hölder metrically subregular at $(x^{*},0)$, i.e., the assumption (A1) holds.

5.1 Computational results

We find zeros of the mapping (83) with a set of real-world biological data using LMLS and LMTR. In details, we compare the performance of LMLS and LMTR with some state-of-the-art algorithms on a set of 21 biochemical reaction networks given in Table 1. In Section 4.2 of [4], it is computationally shown that $\nabla h$ is rank-deficient or ill-conditioned at zeros of the mapping $h$ (83) for these biological models. This clearly justifies the reason of unsuccessful performance of many existing algorithms (e.g., gradient descent, Gauss-Newton, and trust-region methods) and vindicates the development of the two adaptive Levenberg-Mardquart methods (LMLS and LMTR) for such difficult problems.

All codes are written in MATLAB and runs are performed on a Dell Precision Tower 7000 Series 7810 (Dual Intel Xeon Processor E5-2620 v4 with 32 GB RAM). We compare LMLS and LMTR with

• $\bullet$

LM-YF: a Levenberg–Marquard line search method with $\mu_{k}=\|h(x_{k})\|^{2}$, given by Yamashita and Fukushima [48];

• $\bullet$

LM-FY: a Levenberg–Marquardline search method with $\mu_{k}=\|h(x_{k})\|$, given by Fan and Yuan [18];

• $\bullet$

LevMar: a Levenberg–Marquard trust-region method with $\mu_{k}=\|\nabla h(x_{k})h(x_{k})\|$, given by Ipsen et al. [26].

The codes of LMLS and LMTR are publicly available as a part of the COBRA Toolbox v3.0 [25]. Users can pass the solver name to the parameter structure of the MATLAB function optimizeVKmodels.m. For both LMLS and LMTR, on the basis of our experiments with the mapping (83), we set $\omega_{k}:=1-\xi_{k}$ and

[TABLE]

implying $\xi_{k}\in[10^{-10},0.95]$. We here use the starting point $x_{0}=0$ and consider the stopping criterion

[TABLE]

cf. [9]. We stop the algorithms if either (85) holds or the maximum number of iterations (say 10,000 for tuning $\eta$ and 100,000 for the comparison) is reached. While LMLS uses the parameters

[TABLE]

LMTR employs the parameters

[TABLE]

In our comparison, $N_{i}$, $N_{f}$ and $T$ denote the total number of iterations, the total number of function evaluations, and the running time, respectively. To illustrate the results, we used the Dolan and Moré performance profile [17] with the performance measures $N_{f}$ and $T$. In this procedure, the performance of each algorithm is measured by the ratio of its computational outcome versus the best numerical outcome of all algorithms. This performance profile offers a tool to statistically compare the performance of algorithms. Let $\mathcal{S}$ be a set of all algorithms and $\mathcal{P}$ be a set of test problems. For each problem $p$ and algorithm $s$, $t_{p,s}$ denotes the computational outcome with respect to the performance index, which is used in the definition of the performance ratio

[TABLE]

If an algorithm $s$ fails to solve a problem $p$, the procedure sets $r_{p,s}:=r_{\text{failed}}$, where $r_{\text{failed}}$ should be strictly larger than any performance ratio (86). Let $n_{p}$ be the number of problems in the experiment. For any factor $\tau\in\mathbb{R}$, the overall performance of an algorithm $s$ is given by

[TABLE]

Here, $\rho_{s}(\tau)$ is the probability that a performance ratio $r_{p,s}$ of an algorithm $s\in\mathcal{S}$ is within a factor $\tau$ of the best possible ratio. The function $\rho_{s}(\tau)$ is a distribution function for the performance ratio. In particular, $\rho_{s}(1)$ gives the probability that an algorithm $s$ wins over all other considered algorithms, and $\lim_{\tau\rightarrow r_{\mathrm{failed}}}\rho_{s}(\tau)$ gives the probability that algorithm $s$ solves all considered problems. Therefore, this performance profile can be considered as a measure of efficiency among all considered algorithms. In Figures 1 and 3, the number $\tau$ is represented in the $x$-axis, while $P(r_{p,s}\leq\tau:1\leq s\leq n_{s})$ is shown in the $y$-axis.

First, let us tune the parameter $\eta$ to get the best performance of LMLS and LMTR. To do so, we consider several versions of these algorithms corresponding to several levels of the parameter $\eta$ ($\eta=0.6,\leavevmode\nobreak\ 0.8,\leavevmode\nobreak\ 1.0,\leavevmode\nobreak\ 1.2,\leavevmode\nobreak\ 1.4$) and compare the results in Figure 1. From this figure, it is clear that $\eta=1.2$ attains the best results for both LMLS and LMTR. Therefore, we use $\eta=1.2$ for finding a zero of the mapping $h$ defined in (83); however, to solve a different mappings, one may tune this parameter carefully before any practical usage.

Next, we report the results of a comparison among LM-YF, LM-FY, LevMar, LMLS, and LMTR for finding a zero of $h$ (83) with respect to the total number of iterations ($N_{i}$), the total number of function evaluations ($N_{f}$), the mixed measure $N_{f}+3N_{i}$, and the running time ($T$) in Figure 2. From this figure, it can be seen that LMLS and LMTR outperform the others substantially with respect to all considered measures. Moreover, LMTR solves the problems even faster than LMLS; however, the slope of curve of LMLS indicates that its performance is much better than LM-YF, LM-FY, and LevMar, and its performance is close to the performance of LMTR. Surprisingly, both LMLS and LMTR are convergent to a zero of the mapping $h$ (83) not to a stationary point of the merit function $\psi$ given by (2). This clearly show the potential of LMLS and LMTR for finding the moiety conserved steady state of biochemical reaction networks.

Finally, we conclude this section by displaying the evolution of the merit function values during run of the considered algorithms. To this end, we illustrate the function values $\psi$ versus iterations in Figure 3 for the mapping (83) with the biological models iBsu1103 and iSB619. Here, we limit the maximum number of iterations to 1,000. From Figure 3, it can be seen that LMLS and LMTR perform much better than the others; however, the best performance is attained by LMTR.

6 Conclusion and further research

We have employed two globalisation techniques for Levenberg-Marquardt methods for finding a zero of Hölder metrically subregular mappings. First, we combined the Levenberg-Marquardt direction with a nonmonotone Armijo-type line search. Then, we modified the Levenberg-Marquardt parameter and combined the corresponding direction with a nonmonotone trust-region technique. Next, we studied the global convergence and the worst-case global and evaluation complexities or both methods, which are of the order $\mathcal{O}(\varepsilon^{-2})$. The worst-case behavior of the proposed methods, up to a factor, are equivalent to that of the steepest descent method for unconstrained optimisation, cf. [12, 39], which is not the best-known global complexity for nonconvex problems, cf. [13, 42]; however, practical usage of these methods show much better performance than the worse-case complexity, giving scope for future establishement of tigher complexity bounds. Finally, we have studied some special mappings that satisfy certain conditions for a stationary point to corresponds to a zero of a mapping, when obtained with the proposed methods.

We also investigate finding zeros of Hölder metrically subregular mappings that appear in modelling of biochemical reaction networks. Our numerical experiments establish the suitability of the proposed methods for a range of medium- and large-scale biochemical network problems. Nevertheless, biochemical reaction networks on the order of tens of millions of dimensions already exist [37], and the projection is for even larger models in the future. Therefore, considerable scope exists for development of accelerated solution methods.

Acknowledgements

This work was supported by the U.S. Department of Energy, Offices of Advanced Scientific Computing Research and the Biological and Environmental Research as part of the Scientific Discovery Through Advanced Computing program, grant #DE-SC0010429, the OPEN program from the Fonds National de la Recherche Luxembourg (FNR/O16/11402054), and the University of Luxembourg (IRP/OptBioSys).

?refname?

• [1]

Ahookhosh, M., and Amini, K.

An efficient nonmonotone trust-region method for unconstrained optimization.

Numerical Algorithms 59, 4 (2012), 523–540.

• [2]

Ahookhosh, M., Amini, K., and Bahrami, S.

A class of nonmonotone armijo-type line search method for unconstrained optimization.

Optimization 61, 4 (2012), 387–404.

• [3]

Ahookhosh, M., Amini, K., and Kimiaei, M.

A globally convergent trust-region method for large-scale symmetric nonlinear systems.

Numerical Functional Analysis and Optimization 36, 7 (2015), 830–855.

• [4]

Ahookhosh, M., Aragón Artacho, F. J., Fleming, R., and Phan, V.

Local convergence of the levenberg–marquardt methods under hölder metric subregularity.

Submitted (2017).

• [5]

Amini, K., Ahookhosh, M., and Nosratipour, H.

An inexact line search approach using modified nonmonotone strategy for unconstrained optimization.

Numerical Algorithms 66, 1 (2014), 49–78.

• [6]

Aragón Artacho, F., and Fleming, R.

Globally convergent algorithms for finding zeros of duplomonotone mappings.

Optimization Letter (2014), 1–16.

• [7]

Aragón Artacho, F., Fleming, R., and Phan, V.

Accelerating the DC algorithm for smooth functions.

Mathematical Programming (2017).

• [8]

Attouch, H., Bolte, J., Redont, P., and Soubeyran, A.

Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the kurdyka-łojasiewicz inequality.

Mathematics of Operations Research 35, 2 (2010), 438–457.

• [9]

Bellavia, S., Cartis, C., Gould, N., Morini, B., and Toint, P. L.

Convergence of a regularized euclidean residual algorithm for nonlinear least squares.

SIAM Journal on Numerical Analysis 48, 1 (2010), 1–29.

• [10]

Bellavia, S., and Morini, B.

Strong local convergence properties of adaptive regularized methods for nonlinear least squares.

IMA Journal of Numerical Analysis 35, 2 (2015), 947–968.

• [11]

Bolte, J., Sabach, S., and Teboulle, M.

Proximal alternating linearized minimization for nonconvex and nonsmooth problems.

Mathematical Programming 146, 1 (2014), 459–494.

• [12]

Cartis, C., Gould, N., and Toint, P. L.

On the complexity of steepest descent, newton’s and regularized newton’s methods for nonconvex unconstrained optimization.

SIAM Journal on Optimization 20, 6 (2010), 2833–2852.

• [13]

Cartis, C., Gould, N. I., and Toint, P. L.

Adaptive cubic regularisation methods for unconstrained optimization. part ii: worst-case function-and derivative-evaluation complexity.

Mathematical programming 130, 2 (2011), 295–319.

• [14]

Conn, A. R., Gould, N. I., and Toint, P. L.

Trust-Region Methods.

SIAM, 2000.

• [15]

Dennis, J., E., J., and Moré, J. J.

Quasi-newton methods, motivation and theory.

SIAM review 19, 1 (1977), 46–89.

• [16]

Dennis J.r., J. E., and Schnabel, R. B.

Numerical methods for unconstrained optimization and nonlinear equations.

SIAM, 1996.

• [17]

Dolan, E. D., and Moré, J. J.

Benchmarking optimization software with performance profiles.

Mathematical Programming 91, 2 (2002), 201–213.

• [18]

Fan, J., and Yuan, Y.

On the quadratic convergence of the levenberg-marquardt method without nonsingularity assumption.

Computing 74, 1 (2005), 23–39.

• [19]

Fischer, A.

Local behavior of an iterative framework for generalized equations with nonisolated solutions.

Mathematical Programming 94, 1 (2002), 91–124.

• [20]

Fleming, R. M., Vlassis, N., Thiele, I., and Saunders, M. A.

Conditions for duality between fluxes and concentrations in biochemical networks.

Journal of Theoretical Biology 409 (2016), 1 – 10.

• [21]

Gevorgyan, A., Poolman, M., and Fell, D.

Detection of stoichiometric inconsistencies in biomolecular models.

Bioinformatics 24, 19 (2008), 2245–2251.

• [22]

Grippo, L., Lampariello, F., and Lucidi, S.

A nonmonotone line search technique for newton’s method.

SIAM Journal on Numerical Analysis 23, 4 (1986), 707–716.

• [23]

Gu, N., and Mo, J.

Incorporating nonmonotone strategies into the trust region method for unconstrained optimization.

Computers & Mathematics with Applications 55, 9 (2008), 2158 – 2172.

• [24]

Haraldsdóttir, H. S., and Fleming, R. M.

Identification of conserved moieties in metabolic networks by graph theoretical analysis of atom transition networks.

PLoS Comput Biol 12, 11 (2016), e1004,999.

• [25]

Heirendt, L., e. a.

Creation and analysis of biochemical constraint-based models: the cobra toolbox v3.0.

Nature Protocols (accepted) (2018).

• [26]

Ipsen, I., Kelley, C., and Pope, S.

Rank-deficient nonlinear least squares problems and subst selection.

SIAM Journal on Numerical Analysis 49, 3 (2011), 1244–1266.

• [27]

Izmailov, A. F., and Solodov, M. V.

Newton-type methods for optimization and variational problems.

Springer Series in Operations Research and Financial Engineering. Springer, Cham, 2014.

• [28]

Kanzow, C., Yamashita, N., and Fukushima, M.

Levenberg-marquardt methods with strong local convergence properties for solving nonlinear equations with convex constraints.

Journal of Computational and Applied Mathematics 172, 2 (2004), 375–397.

• [29]

Karas, E. W., Santos, S. A., and Svaiter, B. F.

Algebraic rules for computing the regularization parameter of the levenberg–marquardt method.

Computational Optimization and Applications (2015), 1–29.

• [30]

Kelley, C.

Iterative Methods for Linear and Nonlinear Equations.

Frontiers Appl. Math. 16, SIAM, Philadelphia, 1999.

• [31]

Kelley, C.

Iterative Methods for Optimization.

Frontiers Appl. Math. 18, SIAM, Philadelphia, 1999.

• [32]

Klamt, S., Haus, U.-U., and Theis, F.

Hypergraphs and cellular networks.

PLoS Comput Biol 5, 5 (2009), e1000385.

• [33]

Kurdyka, K.

On gradients of functions definable in o-minimal structures.

Annales de l’institut Fourier 48, 3 (1998), 769–783.

• [34]

Li, G., Mordukhovich, B. S., and Phạm, T. S.

New fractional error bounds for polynomial systems with applications to hölderian stability in optimization and spectral theory of tensors.

Mathematical Programming 153, 2 (Nov 2015), 333–362.

• [35]

Lojasiewicz, S.

Une propriété topologique des sou s-ensembles analytiques réels.

Les équations aux dérivées partielles 117 (1963), 87–89.

• [36]

Lojasiewicz, S.

Ensembles semi-analytiques.

Université de Gracovie, 1965.

• [37]

Magnúsdóttir, S., Heinken, A., Kutt, L., Ravcheev, D. A., Bauer, E., Noronha, A., Greenhalgh, K., Jäger, C., Baginska, J., Wilmes, P., Fleming, R. M. T., and Thiele, I.

Generation of genome-scale metabolic reconstructions for 773 members of the human gut microbiota.

Nature Biotechnology 35, 1 (nov 2016), 81–89.

• [38]

Moré, J. J.

The levenberg-marquardt algorithm: implementation and theory.

In Numerical analysis. Springer, 1978, pp. 105–116.

• [39]

Nesterov, Y.

Introductory Lectures on Convex Optimization: A Basic Course.

Kluwer, Dordrecht, 2004.

• [40]

Nesterov, Y.

Modified gauss-newton scheme with worst-case guarantees for global performance.

Optimization Methods and Software 22, 3 (2007), 521–539.

• [41]

Nesterov, Y.

Gradient methods for minimizing composite functions.

Mathematical Programming 140, 1 (2013), 125–161.

• [42]

Nesterov, Y., and Polyak, B. T.

Cubic regularization of newton method and its global performance.

Mathematical Programming 108, 1 (2006), 177–205.

• [43]

Nocedal, J., and Wright, S.

Numerical optimization.

Springer, New York, 2006.

• [44]

Ortega, J., and Rheinboldt, W.

Iterative Solution of Nonlinear Equations in Several Variables.

Society for Industrial and Applied Mathematics, Jan. 2000.

• [45]

Ueda, K., and Yamashita, N.

On a global complexity bound of the levenberg-marquardt method.

Journal of Optimization Theory and Applications 147, 3 (Dec 2010), 443–453.

• [46]

Ueda, K., and Yamashita, N.

Global complexity bound analysis of the levenberg–marquardt method for nonsmooth equations and its application to the nonlinear complementarity problem.

Journal of Optimization Theory and Applications 152, 2 (Feb 2012), 450–467.

• [47]

Wright, S., and Holt, J. N.

An inexact levenberg-marquardt method for large sparse nonlinear least squres.

The Journal of the Australian Mathematical Society. Series B. Applied Mathematics 26, 04 (1985), 387–403.

• [48]

Yamashita, N., and Fukushima, M.

On the rate of convergence of the levenberg-marquardt method.

In Topics in Numerical Analysis, G. Alefeld and X. Chen, Eds., vol. 15. Springer Vienna, Vienna, 2001, pp. 239–249.

• [49]

Yudin, D., and Nemirovskii, A.

Problem Complexity and Method Efficiency in Optimization.

John Wiley and Sons, 1983.