Local convergence of the Levenberg-Marquardt method under H\"{o}lder metric subregularity

TL;DR

This paper investigates the local convergence of the Levenberg-Marquardt method for nonlinear equations under H"older metric subregularity, proposing an adaptive parameter and analyzing convergence with additional ojasiewicz inequality, supported by biochemical network applications.

Contribution

It introduces an adaptive Levenberg-Marquardt parameter formula and analyzes convergence under H"older metric subregularity and ojasiewicz inequality, extending existing convergence theory.

Findings

Convergence is established under H"older metric subregularity.

Numerical results confirm theoretical convergence for biochemical problems.

The method effectively computes steady states in biochemical networks.

Abstract

We describe and analyse Levenberg-Marquardt methods for solving systems of nonlinear equations. More specifically, we propose an adaptive formula for the Levenberg-Marquardt parameter and analyse the local convergence of the method under H\"{o}lder metric subregularity of the function defining the equation and H\"older continuity of its gradient mapping. Further, we analyse the local convergence of the method under the additional assumption that the \L{}ojasiewicz gradient inequality holds. We finally report encouraging numerical results confirming the theoretical findings for the problem of computing moiety conserved steady states in biochemical reaction networks. This problem can be cast as finding a solution of a system of nonlinear equations, where the associated mapping satisfies the \L{}ojasiewicz gradient inequality assumption.

Tables

Table 1. Table 1: Summary of the results of tuning the parameter η 𝜂 \eta for 1 with parameters ( 47 ) and η ∈ { 0.6 , 0.7 , 0.8 , 0.9 , 0.99 , 0.999 , 1 } 𝜂 0.6 0.7 0.8 0.9 0.99 0.999 1 \eta\in\{0.6,0.7,0.8,0.9,0.99,0.999,1\} to solve ( 46 ) in 20 biological models. For each model, the lowest number of iterations ( N i subscript 𝑁 𝑖 N_{i} ) and the lowest running time ( T 𝑇 T ) are displayed in bold.

Model	$m$	$n$	$r$	$\eta =0.6$		$\eta =0.7$		$\eta =0.8$		$\eta =0.9$		$\eta =0.99$		$\eta =0.999$		$\eta =1$
				$N_i$	$T$	$N_i$	$T$	$N_i$	$T$	$N_i$	$T$	$N_i$	$T$	$N_i$	$T$	$N_i$	$T$
1. Ecoli_core	72	73	61	210	0.10	191	0.08	179	0.07	159	0.06	136	0.06	138	0.04	139	0.09
2. iAF692	462	493	430	492	6.89	421	5.70	328	4.54	291	3.98	274	3.72	256	3.44	253	3.42
3. iAF1260	1520	1931	1456	473	98.75	410	83.41	357	74.05	334	68.65	257	53.69	271	57.05	268	54.88
4. iBsu1103	993	1167	956	421	34.34	356	28.46	313	25.14	254	20.34	232	18.38	218	17.13	226	18.04
5. iCB925	415	558	386	467	6.41	477	6.47	332	4.55	296	4.30	318	4.39	248	3.33	350	4.71
6. iIT341	424	428	392	388	4.33	333	3.48	253	2.68	226	2.29	207	2.13	226	2.34	212	2.19
7. iJN678	641	669	589	362	9.63	356	8.99	307	7.76	258	6.72	220	5.76	231	5.89	218	5.60
8. iJN746	727	795	700	470	17.05	376	13.29	301	10.66	255	8.93	256	9.26	231	8.17	251	9.10
9. iJO1366	1654	2102	1582	417	107.37	372	95.21	314	79.85	273	70.62	225	57.36	219	55.65	244	62.37
10. iJR904	597	757	564	441	11.99	420	11.15	346	8.80	279	7.26	245	6.34	249	6.39	253	6.98
11. iMB745	525	598	490	363	7.13	346	6.48	298	5.55	211	3.85	204	3.74	218	4.02	199	3.66
12. iNJ661	651	764	604	549	16.21	436	12.68	366	11.01	283	8.30	222	6.45	257	7.53	254	7.57
13. iRsp1095	966	1042	921	655	54.48	1057	92.98	374	28.20	333	25.03	301	21.89	301	22.63	336	25.52
14. iSB619	462	508	435	373	5.13	344	4.60	295	3.95	243	3.24	203	2.73	221	2.97	215	2.85
15. iTH366	583	606	529	349	7.27	338	6.91	279	5.67	219	4.46	212	4.33	204	4.16	207	4.21
16. iTZ479_v2	435	476	415	375	4.54	344	4.00	297	3.45	228	2.68	204	2.39	214	2.48	227	2.63
17. iYL1228	1350	1695	1280	846	136.52	550	87.18	387	60.96	326	51.24	294	46.28	318	49.99	293	46.18
18. L_lactis_MG1363	483	491	429	463	6.54	427	5.89	351	4.80	314	4.36	279	3.85	242	3.35	282	3.87
19. Sc_thermophilis_rBioNet	348	365	320	413	2.85	377	2.57	344	2.31	291	1.97	252	1.66	256	1.69	241	1.63
20. T_Maritima	434	470	414	407	4.83	318	3.72	281	3.26	233	2.70	231	2.65	215	2.49	212	2.44
Average				447	27.12	412	24.16	315	17.36	265	15.05	239	12.85	237	13.04	244	13.4

Table 2. Table 2: Summary of the results of LM-YF, LM-FY, LM-F, and 1 with parameters ( 47 ) and η = 0.999 𝜂 0.999 \eta=0.999 for solving ( 46 ) in 20 biological models. For each model, the lowest number of iterations ( N i subscript 𝑁 𝑖 N_{i} ) and the lowest running time ( T 𝑇 T ) are displayed in bold. On the bottom, we show the average among the successfully solved instances.

Model	$m$	$n$	$r$	LM-YF		LM-FY		LM-F		1
				$N_i$	$T$	$N_i$	$T$	$N_i$	$T$	$N_i$	$T$
1. Ecoli_core	72	73	61	238	0.10	10000	4.06	257	0.08	153	0.06
2. iAF692	462	493	430	1685	23.43	10000	135.63	1358	18.15	271	3.61
3. iAF1260	1520	1931	1456	8233	1726.92	10000	2066.38	2036	413.04	283	57.27
4. iBsu1103	993	1167	956	2396	187.24	10000	780.72	761	59.39	193	15.09
5. iCB925	415	558	386	1005	14.15	10000	131.01	10000	131.28	278	3.82
6. iIT341	424	428	392	1407	14.77	10000	103.09	944	9.71	222	2.26
7. iJN678	641	669	589	2218	56.33	10000	253.58	1043	26.34	229	5.82
8. iJN746	727	795	700	3107	108.72	10000	349.43	1069	37.18	217	7.53
9. iJO1366	1654	2102	1582	7716	1946.48	10000	2524.01	1066	268.38	232	58.33
10. iJR904	597	757	564	2789	72.26	10000	258.50	1231	31.74	262	6.75
11. iMB745	525	598	490	790	14.60	10000	181.40	1247	22.50	208	3.76
12. iNJ661	651	764	604	2635	76.62	10000	290.56	1357	39.45	360	10.44
13. iRsp1095	966	1042	921	3832	266.87	10000	694.21	10000	696.93	235	16.38
14. iSB619	462	508	435	1581	21.42	10000	133.94	814	10.82	233	3.12
15. iTH366	583	606	529	1641	33.78	10000	205.34	817	16.79	211	4.30
16. iTZ479_v2	435	476	415	1148	13.57	10000	117.46	713	8.34	221	2.61
17. iYL1228	1350	1695	1280	6070	956.92	10000	1565.75	10000	1567.47	272	42.94
18. L_lactis_MG1363	483	491	429	2180	30.17	10000	137.99	1231	17.06	287	4.10
19. Sc_thermophilis_rBioNet	348	365	320	1753	11.85	10000	68.54	935	6.33	244	1.61
20. T_Maritima	434	470	414	1169	14.33	10000	118.65	717	8.31	209	2.42
Average of successful				2680	279.53	—	—	1035	58.45	241	12.61

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,$$

$$\underset{d\in {\mathbb{R}}^m}{\min }{\varphi }_k(d),$$

$${\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),$$

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

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

$$\beta \mathrm{dist}\left(x,h^{-1}(\overline{y})\right)\le \Vert \overline{y}-h(x){\Vert }^{\delta },\forall x\in \mathbb{B}(\overline{x},r).$$

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

$$\mathrm{dist}(x,\Omega {)}^{\gamma }\le \varrho |\psi (x)|,\forall x\in C.$$

$$|\psi (x)-\psi (x^{\ast }){|}^{\theta }\le \kappa \Vert \nabla \psi (x)\Vert ,\forall x\in \mathbb{B}(x^{\ast },\epsilon ).$$

$$\frac{1}{\varrho }\mathrm{dist}(x,\Omega {)}^{\gamma }\le |\psi (x)|\le {\kappa }^{1/\theta }\Vert \nabla \psi (x){\Vert }^{1/\theta }={\kappa }^{1/\theta }\Vert \nabla h(x)h(x){\Vert }^{1/\theta }.$$

$$h(x_1,x_2,x_3,x_4):=\left(x_1+10x_2,\sqrt{5}(x_3-x_4),(x_2-2x_3{)}^2,\sqrt{10}(x_1-x_4{)}^2\right),$$

$$\mathrm{dist}(x_k,\Omega )=\Vert x_k\Vert =\frac{\sqrt{2}}{k}=\mathcal{O}(k^{-1}).$$

$$\frac{1}{2}\Vert h(x){\Vert }^2=\psi (x)\ge \beta \Vert x{\Vert }^{(4-1{)}^4+1}=\beta \Vert x{\Vert }^{82},\forall x\in \mathbb{B}(0_4,r).$$

$${\psi }_2(z):=\nabla h(\overline{x}{)}^T+D^{2}Ph(\overline{x})(z,\cdot ),$$

$$T_2:=\left\{z\in {\mathbb{R}}^m\mid \nabla h(\overline{x}{)}^{T}z=0_n\text{ and }{D}^{2}Ph(\overline{x})(z,z)=0_n\right\},$$

$${\psi }_2(z)$$

$$\Vert g(x)-g(y)\Vert \le L\Vert x-y{\Vert }^{\upsilon },\forall x,y\in \mathbb{B}(x^{\ast },r).$$

$$\Vert h(y)-h($$

$$=\Vert {\int }_0^1\nabla h(x+t(y-x){)}^T(y-x)dt-\nabla h(x{)}^T(y-x)\Vert$$

$$\le \Vert y-x\Vert {\int }_0^1\Vert \nabla h(x+t(y-x))-\nabla h(x)\Vert dt$$

$$\le L\Vert y-x{\Vert }^{1+\upsilon }{\int }_0^1{t}^{\upsilon }dt=\frac{L}{1+\upsilon }\Vert y-x{\Vert }^{1+\upsilon }.$$

$$\Vert x-x^{\ast }\Vert =\mathrm{dist}(x,\Omega ),\forall x\in \mathbb{B}(x^{\ast },r).$$

$$w_k:=x^{\ast }+\frac{r}{k\Vert z\Vert }z,\text{with }k=1,2,\dots .$$

$$\nabla h(x^{\ast }{)}^T(w_k-x^{\ast })=\frac{r}{k\Vert z\Vert }\nabla h(x^{\ast }{)}^{T}z=0.$$

$$\beta \Vert w_k-x^{\ast }\Vert$$

$$=\Vert h(w_k)-h(x^{\ast })-\nabla h(x^{\ast })(w_k-x^{\ast }){\Vert }^{\delta }$$

$$\le \frac{L^{\delta }}{(1+\upsilon {)}^{\delta }}\Vert w_k-x^{\ast }{\Vert }^{(1+\upsilon )\delta }.$$

$$\Vert w_k-x^{\ast }{\Vert }^{(1+\upsilon )\delta -1}\ge \frac{\beta (1+\upsilon {)}^{\delta }}{L^{\delta }},$$

$$\Vert h(y)-h(x)-\nabla h(x{)}^T(y-x)\Vert \le \frac{L}{1+\upsilon }\Vert y-x{\Vert }^{1+\upsilon },\forall x,y\in \mathbb{B}(x^{\ast },r).$$

$$\widetilde{r}:=\left\{\begin{array}[]{ll}\frac{r}{2},&\text{if }\xi_{\min}>0,\\ \min\left\{\frac{r}{2},\left(\frac{\beta^{2}(1+\upsilon)^{2\delta}}{2^{\delta}L^{2\delta}}\right)^{\frac{1}{2\delta(1+\upsilon)-2}}\right\},&\text{otherwise},\end{array}\right.$$

$$\varpi:=\left\{\begin{array}[]{ll}1,&\text{if }\xi_{\min}>0,\\ 2-\delta,&\text{otherwise}.\end{array}\right.$$

Full text

1 Introduction

For a given continuously differentiable mapping $h:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$, we consider the problem of finding a solution of the system of nonlinear equations

[TABLE]

We denote by $\Omega$ the set of solutions of this problem, which is assumed to be nonempty. Systems of nonlinear equations of type (1) frequently appear in the mathematical modelling of many real-world applications in the fields of solid-state physics Eilenberger , quantum field theory, optics, plasma physics hasegawa_plasma_1975 , fluid mechanics Whitham , chemical kinetics artacho_globally_2014 ; artacho_accelerating_2015 , and applied mathematics including the discretisation of ordinary and partial differential equations ortega_iterative_2000 .

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. This is a well-studied topic and there are many iterative schemes with fast local convergence rates (e.g., superlinear or quadratic) such as Newton, quasi-Newton, Gauss–Newton, adaptive regularised methods, and the Levenberg–Marquardt method. When $m=n$, to guarantee fast local convergence, these methods require an initial point $x_{0}$ to be sufficiently close to a solution $x^{*}$, and the matrix gradient of $h$ at $x^{*}$ (i.e., the transpose of the Jacobian matrix), denoted by $\nabla h(x^{*})$, to be nonsingular (i.e., full rank), cf. BelCGMT ; fischer2015globally ; nocedal_numerical_2006-1 ; ortega_iterative_2000 ; Yuan2015 .

The Levenberg–Marquardt method is a standard technique used to solve the nonlinear system (1), which is a combination of the gradient descent and the Gauss–Newton methods. More precisely, in each step, for a positive parameter $\mu_{k}$, the convex subproblem

[TABLE]

with $\phi_{k}:\mathbb{R}^{m}\rightarrow\mathbb{\mathbb{R}}$ given by

[TABLE]

is solved to compute a direction $d_{k}$, which is the unique solution to the system of linear equations

[TABLE]

where $I\in\mathbb{\mathbb{R}}^{m\times m}$ denotes the identity matrix. By choosing a suitable parameter $\mu_{k}$, the Levenberg–Marquardt method acts like the gradient descent method whenever the current iteration is far from a solution $x^{*}$, and behaves similar to the Gauss–Newton method if the current iteration is close to $x^{*}$. The parameter $\mu_{k}$ helps to overcome problematic cases where $\nabla h(x_{k})\nabla h(x_{k})^{T}$ is singular, or nearly singular, and thus ensures the existence of a unique solution to (4), or avoids very large steps, respectively. For $m=n$, the Levenberg–Marquardt method is known to be quadratically convergent to a solution of (1) if $\nabla h(x^{*})$ is nonsingular. In fact, the nonsingularity assumption implies that the solution to the minimisation problem (2) must be locally unique, see bellavia_strong_2015 ; kanzow_levenbergmarquardt_2004 ; alefeld_rate_2001 . However, assuming local uniqueness of the solution might be restrictive for many applications.

The notion of (local) error bound usually plays a key role in establishing the rate of convergence of the sequence of iterations generated by a given algorithm. This condition guarantees that the distance from the current iteration $x_{k}$ to the solution set $\Omega$, denoted by $\mathrm{dist}(x_{k},\Omega)=\inf_{y\in\Omega}\|x_{k}-y\|$, is less than the value of a residual function $R:\mathbb{R}^{m}\to\mathbb{R_{+}}$ at that point ($R(x_{k})$). The earliest publication using error bounds for solving a linear inequality system is due to Hoffman Hoffman_1952 , which was followed by many other authors, especially in optimisation. For more information about error bounds, we recommend the nice survey Pang_1997 .

For the particular case of nonlinear systems of equations, Yamashita and Fukushima alefeld_rate_2001 proved the local quadratic convergence of the Levenberg–Marquardt method with $\mu_{k}=\|h(x_{k})\|^{2}$ assuming a local error bound condition. More precisely, they assumed metric subregularity of $h$ around $(x^{*},0)$, which entails the existence of some constants $\beta>0$ and $r>0$ such that

[TABLE]

where $\mathbb{B}(x^{*},r)$ denotes the closed ball centered at $x^{*}$ with radius $r>0$. In this case, the residual function is given by $R(x):=\frac{1}{\beta}\|h(x)\|$. In those situations where the value of $\beta$ is known, the condition $\|h(x)\|\leq\varepsilon$ can be used as a stopping criterion for an iterative scheme, as it entails that the iterations must be close to a solution of (1).

Let us emphasise that, for $m=n$, the nonsingularity of $\nabla h(x^{*})$ implies that $x^{*}$ is locally unique and that (5) holds. Indeed, by the Lyusternik–Graves theorem (see, e.g., (DontchevRockafellar2014, , Theorem 5D.5), (Mordukhovich2006, , Theorem 1.57), or (CDK18, , Proposition 1.2)), the nonsingularity of $\nabla h(x^{*})$ is equivalent to the strong metric regularity of $h$ at $(x^{*},0)$, which implies strong metric subregularity of $h$ at $(x^{*},0)$. However, the latter does not imply the nonsingularity assumption and allows the solutions to be locally nonunique. This means that metric subregularity is a weaker assumption than the nonsingularity. In fact, for $m$ possibly different than $n$, strong metric subregularity of $h$ at $(x^{*},0)$ is equivalent to surjectivity of $\nabla h(x^{*})$ (see, e.g., (CDK18, , Proposition 1.2 and Theorem 2.6)). The successful use of the local error bound has motivated many researchers to investigate, under assumption (5), the local convergence of trust-region methods Fan2006 , adaptive regularised methods bellavia_strong_2015 , and the Levenberg–Marquardt method behling_effect_2013 ; fan_modified_2012 ; fan_quadratic_2005 , among other iterative schemes.

The main motivation for this paper comes from a nonlinear system of equations, the solution of which corresponds to a steady state of a given biochemical reaction network, which plays a crucial role in the modeling of biochemical reaction systems. These problems are usually ill-conditioned and require the application of the Levenberg–Marquardt method. As we numerically show in Section 4, $\nabla h$ is usually rank deficient at the solutions of (1). During our study of the properties of this problem, we were not able to show that the metric subregularity condition (5) is satisfied. However, taking standard biochemical assumptions artacho_accelerating_2015 , we can show that the corresponding merit function is real analytic and thus satisfies the Łojasiewicz gradient inequality and is Hölder metrically subregular around the solutions.

The local convergence of a Levenberg–Marquardt method under Hölder metric subregularity has been recently studied in Guo2015 ; ZhuLin16 . Nonetheless, the standard rules for the regularisation parameter have a very poor performance when they are applied for solving the nonlinear equation arising from the biochemical reaction network systems, as we show in a numerical experiment in Section 4. This motivated our quest to further investigate an adaptive Levenberg–Marquart method under the assumption that the underlying mapping is Hölder metrically subregular.

From the definition of the Levenberg–Marquardt direction in (4), we observe that a key factor in the performance of the Levenberg–Marquardt method is the choice of the parameter $\mu_{k}$, cf. izmailov2014newton ; Kelley1 . Several parameters have been proposed to improve the efficiency of the method. For example, Yamashita and Fukushima alefeld_rate_2001 took $\mu_{k}=\|h(x_{k})\|^{2}$, Fischer fischer2002local used $\mu_{k}=\|\nabla h(x_{k})h(x_{k})\|$, while Fan and Yuan fan_quadratic_2005 proposed $\mu_{k}=\|h(x_{k})\|^{\eta}$ with $\eta\in[1,2]$. Ma and Jiang MJ07 proposed a convex combination of these two types of parameters, namely, $\mu_{k}=\theta\|h(x_{k})\|+(1-\theta)\|\nabla h(x_{k})h(x_{k})\|$ for some constant $\theta\in{[0,1]}$. In a subsequent work, Fan and Pan fan_pan proposed the more general choice $\mu_{k}=\xi_{k}\rho(x_{k})$, where $\xi_{k}$ is updated by a trust-region technique, $\rho(x_{k})=\min\left\{\widetilde{\rho}(x_{k}),1\right\}$ and $\widetilde{\rho}:\mathbb{R}^{m}\to\mathbb{R}_{+}$ is a positive function such that $\widetilde{\rho}(x_{k})=O\left(\|h(x_{k})\|^{\eta}\right)$, with $\eta\in{]0,2]}$. Inspired by these works, and assuming that the function $h$ is Hölder metrically subregular of order $\delta\in\ ]0,1]$ and its gradient $\nabla h$ is Hölder continuous of order $\upsilon\in{]0,1]}$, in this paper we consider an adaptive parameter of the form

[TABLE]

where $\eta>0$, $\xi_{k}\in[\xi_{\min},\xi_{\max}]$ and $\omega_{k}\in[\omega_{\min},\omega_{\max}]$, for some constants $0\leq\xi_{\min}\leq\xi_{\max}$ and $0\leq\omega_{\min}\leq\omega_{\max}$ such that $\xi_{\min}+\omega_{\min}>0$.

In our first main result, Theorem 3.1, we provide an interval depending on $\delta$ and $\upsilon$ where the parameter $\eta$ must be chosen to guarantee the superlinear convergence of the sequence generated by the Levenberg–Marquardt method with the adaptive parameter (6). In our second main result, Theorem 3.2, under the additional assumption that the merit function $\psi$ defined in (2) satisfies the Łojasiewicz gradient inequality with exponent $\theta\in{]0,1[}$, we prove local convergence for every parameter $\eta$ smaller than a constant depending on both $\upsilon$ and $\theta$. As a consequence, we can ensure local convergence of the Lebenverg–Marquardt algorithm to a solution of (1) for all the above-mentioned biochemical networks as long as the parameter $\eta$ is chosen sufficiently small. To the best of our knowledge, this is the first such algorithm able to reliably handle these nonlinear systems arising in the study of biological networks. We successfully apply the proposed algorithm to nonlinear systems derived from many real biological networks, which are representative of a diverse set of biological species.

The remainder of this paper is organised as follows. In the next section, we particularise the Hölder metric subregularity for nonlinear equations and recall the Łojasiewicz inequalities. We investigate the local convergence of the Levenberg–Marquardt method under these conditions in Section 3. In Section 4, we report encouraging numerical results where nonlinear systems, arising from biochemical reaction networks, were quickly solved. Finally, we deliver some conclusions in Section 5.

2 Hölder metric subregularity and Łojasiewicz inequalities

Let us begin this section by recalling the notion of Hölder metric subregularity, which can be also defined in a similar manner for set-valued mappings (see, e.g., kruger_error_2015 ; CDK18 ).

Definition 1

A mapping $h:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ is said to be Hölder metrically subregular of order $\delta>0$ around $(\overline{x},\overline{y})$ with $\overline{y}=h(\overline{x})$ if there exist some constants $r>0$ and $\beta>0$ such that

[TABLE]

For any solution $x^{*}\in\Omega$ of the system of nonlinear equations (1), the Hölder metric subregularity of $h$ around $(x^{*},0)$ reduces to

[TABLE]

Therefore, this property provides an upper bound for the distance from any point sufficiently close to the solution $x^{*}$ to the nearest zero of the function.

Hölder metric subregularity around $(x^{*},0)$ is also called Hölderian local error bound ngai_global_2015 ; vui_global_2013 . It is known that Hölder metric subregularity is closely related to the Łojasiewicz inequalities, which are defined as follows.

Definition 2

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 U,\,\psi(x)=0\right\}$ is nonempty.

• (i)

The function $\psi$ is said to satisfy the Łojasiewicz inequality if for every compact subset $C\subset U$, there exist positive constants $\varrho$ and $\gamma$ such that

[TABLE]

• (ii)

The function $\psi$ is said to satisfy the Łojasiewicz gradient inequality if for any critical point $x^{*}$, there exist constants $\kappa>0,\varepsilon>0$ and $\theta\in{]0,1[}$ such that

[TABLE]

Stanisław Łojasiewicz proved that every real analytic function satisfies these properties lojasiewicz1965ensembles . Recall that a function $\psi:\mathbb{R}^{m}\to\mathbb{R}$ is said to be real analytic if it can be represented by a convergent power series. Fortunately, real analytic functions frequently appear in real world application problems. A relevant example in biochemistry is presented in Section 4.

Fact 1** ((lojasiewicz1965ensembles, , pp. 62 and 67))**

Every real analytic function $\psi:\mathbb{R}^{m}\to\mathbb{R}$ satisfies both the Łojasiewicz inequality and the Łojasiewicz gradient inequality.

Clearly, if the merit function $\psi(\cdot)=\frac{1}{2}\|h(\cdot)\|^{2}$ satisfies the Łojasiewicz inequality (8), then the mapping $h$ satisfies (7) with $\beta:=(2/\varrho)^{1/\gamma}$ and $\delta:=2/\gamma$; i.e., $h$ is Hölder metrically subregular around $(x^{*},0)$ of order $2/\gamma$. In addition, if $\psi(\cdot)$ satisfies the Łojasiewicz gradient inequality (9), then for any $\overline{x}\in\Omega$ and $x\in\mathbb{B}(\overline{x},\varepsilon)$, it holds

[TABLE]

The Łojasiewicz gradient inequality has recently gained much attention because of its role for proving the convergence of various numerical methods (e.g., BDL07 ; AB09 ; ABS13 ; artacho_accelerating_2015 ). The connection between this property and metric regularity of the set-valued mapping $\Psi(x):=[\psi(x),\infty[$ on an adequate set was revealed in BDLM10 , where it was also applied to deduce strong convergence of the proximal algorithm.

In some cases, for example when $\psi$ is a polynomial with an isolated zero at the origin, an order of the Hölder metric subregularity is known gwozdziewicz_lojasiewicz_1999 ; kurdyka_separation_2014 ; li_holder_2012 .

Fact 2** ((gwozdziewicz_lojasiewicz_1999, , Theorem 1.5))**

Let $\psi:\mathbb{R}^{m}\to\mathbb{R}$ be a polynomial function with an isolated zero at the origin. Then $\psi$ is Hölder metrically subregular around $(0,0)$ of order $\left((\mathrm{deg}\,\psi-1)^{m}+1\right)^{-1}$, where $\mathrm{deg}\,\psi$ denotes the degree of the polynomial function $\psi$.

The next example shows that the Powell singular function, which is a classical test function for nonlinear systems of equations, is not metrically subregular around its unique solution but is Hölder metrically subregular there. In addition, it demonstrates that the order given by Fact 2 is, in general, far from being tight.

Example 1

The Powell singular function more_1981 , which is the function $h:\mathbb{R}^{4}\to\mathbb{R}^{4}$ given by

[TABLE]

is (strongly) Hölder metrically subregular around $(0_{4},0)$ but does not satisfy the metric subregularity condition (5). We have $\Omega=\left\{0_{4}\right\}$ and $\nabla h(0_{4})$ is singular; thus, $h$ is not metrically regular around $(0_{4},0)$. Further, to prove that (5) does not hold, consider the sequence $\left\{x_{k}\right\}$ defined by $x_{k}=\left(0,0,\frac{1}{k},\frac{1}{k}\right)$. We see that $\left\{x_{k}\right\}\to 0_{4}$ and

[TABLE]

Since $\|h(x_{k})\|=\frac{\sqrt{26}}{k^{2}}=\mathcal{O}(k^{-2}),$ we conclude that (5) does not hold.

Consider the polynomial function $\psi(x):=\frac{1}{2}\|h(x)\|^{2}$ of degree $4$, which satisfies $\psi^{-1}(0)=0_{4}$. It follows from Fact 2 that there exist some constants $\beta>0$ and $r>0$ such that

[TABLE]

This implies that $h$ is Hölder metrically subregular of order $\delta=\frac{1}{41}$ around $(0_{4},0)$. Nonetheless, the order $\frac{1}{41}$ given by Fact 2 can be improved by using the theory of $2$-regularity: the function $h$ turns out to be $2$-regular at $0_{4}$, which implies by (izmailov2001, , Theorem 4) that (7) holds with $\delta=\frac{1}{2}$ (see also (izmailov2001, , Remark 7)). Recall that a twice differentiable mapping $h:\mathbb{R}^{m}\to\mathbb{R}^{n}$ is said to be 2-regular at the point $\overline{x}$ if the range of $\psi_{2}(z)$ is $\mathbb{R}^{n}$ for all $z\in T_{2}\setminus\{0\}$, where $\psi_{2}:\mathbb{R}^{m}\to\mathbb{R}^{n\times m}$ is defined for $z\in\mathbb{R}^{m}$ by

[TABLE]

$P$ is the projector in $\mathbb{R}^{n}$ onto the complementary subspace to the range of $\nabla h(\overline{x})^{T}$, and $D^{2}$ stands for the second-order (Fréchet) derivative.

Indeed, for any $z\in\mathbb{R}^{4}$, one has $\nabla h(0_{4})^{T}z=(z_{1}+10z_{2},\sqrt{5}(z_{3}-z_{4}),0,0)^{T}$, so the range of $\nabla h(0_{4})^{T}$ is $Y_{1}=\mathbb{R}^{2}\times\{0_{2}\}$, whose complementary subspace is $Y_{2}=\{0_{2}\}\times\mathbb{R}^{2}$. Then, $T_{2}=\{(-10t,t,0,0)^{T},t\in\mathbb{R}\}$ and for each $z\in T_{2}\setminus\{0_{4}\}$, one has

[TABLE]

which is full-rank for all $t\neq 0$. Therefore, the range of $\psi_{2}(z)$ is equal to $\mathbb{R}^{4}$ for all $z\in T_{2}\setminus\{0_{4}\}$, and the function $h$ is $2$-regular at $0_{4}$.$\Diamond$

There are many examples of smooth functions that are Hölder metrically subregular of order $\delta$ around some zero of the function and whose gradient is not full row rank at that point, cf. izmailov2001 ; izmailov2002 . Nonetheless, the following result restricts the possible values of $\delta$: if $x^{*}$ is an isolated solution in $\Omega$ (i.e., the function is Hölder strongly metrically subregular at $x^{*}$, cf. Mordukhovich2015 ; CDK18 ), and $\nabla h$ is Lipschitz continuous around $x^{*}$ then one must have $\delta\in{]0,1/2]}$ if $\delta\neq 1$. In fact, only Hölder continuity of $\nabla h$ is needed. Recall that a function $g:\mathbb{R}^{m}\to\mathbb{R}^{n}$ is said to be Hölder continuous of order $\upsilon\in{]0,1]}$ with constant $L>0$ around some point $x^{*}\in\mathbb{R}^{m}$ whenever there exist a positive constant $r$ such that

[TABLE]

When $\upsilon=1$, $g$ is said to be Lipschitz continuous with constant $L$ around $x^{*}$.

Proposition 1

Let $h:\mathbb{R}^{m}\to\mathbb{R}^{n}$ be a continuously differentiable function which is Hölder metrically subregular of order $\delta$ around some isolated solution $x^{*}\in\Omega=\{x\in\mathbb{R}^{m}:h(x)=0\}$. Assume further that $\nabla h$ is Hölder continuous around $x^{*}$ of order $\upsilon\in{]0,1]}$ and that $\nabla h(x^{*})$ is not full row rank. Then, it holds that $\delta\in\left]0,\frac{1}{1+\upsilon}\right]$.

Proof

Because of the Hölder continuity assumption and the mean value theorem, there are some positive constants $L$ and $r$ such that, for all $x,y\in\mathbb{B}(x^{*},r)$, it holds

[TABLE]

By using the fact that $x^{*}$ is an isolated solution, it is possible to make $r$ smaller if needed so that (7) holds and

[TABLE]

Since $\nabla h(x^{*})$ is not full row rank, there exists some $z\neq 0$ such that $\nabla h(x^{*})^{T}z=0$. Consider now the points

[TABLE]

Observe that

[TABLE]

As $w_{k}\in\mathbb{B}(x^{*},r)$ for all $k$, we deduce

[TABLE]

Thus, we get

[TABLE]

which implies that $\delta\leq\frac{1}{1+\upsilon}$, since $w_{k}\to x^{*}$, as claimed.∎

The next example shows that the full rank assumption in Proposition 1 is not redundant, and that the upper bound $\delta\leq\frac{1}{1+\upsilon}$ can be attained.

Example 2

Consider the continuously differentiable functions $h,\widehat{h}:\mathbb{R}\to\mathbb{R}$ given for $x\in\mathbb{R}$ by $h(x):=\frac{3}{4}\sqrt[3]{x^{4}}$ and $\widehat{h}(x):=\frac{3}{4}\sqrt[3]{x^{4}}+x$, whose solution sets are $\Omega=\{0\}$ and $\widehat{\Omega}=\left\{-\frac{64}{27},0\right\}$, respectively. Let $x^{*}:=0\in\Omega\cap\widehat{\Omega}$. Then, $h^{\prime}(x)=\sqrt[3]{x}$ and $\widehat{h}^{\prime}(x)=\sqrt[3]{x}+1$, which are both Hölder continuous around $x^{*}$ of order $\upsilon=\widehat{\upsilon}=\frac{1}{3}$. Observe that $h^{\prime}(0)=0$ while $\widehat{h}^{\prime}(0)=1$. Hence, it follows that $\widehat{h}$ is (Hölder) metrically subregular around $x^{*}$ of order $\widehat{\delta}:=1>\frac{1}{1+\widehat{\upsilon}}$, while it is easy to check that $h$ is Hölder metrically subregular around $x^{*}$ of order $\delta:=\frac{3}{4}=\frac{1}{1+\upsilon}$.$\Diamond$

3 Local convergence of the Levenberg–Marquardt method

In this section, to solve a nonlinear system of the form (1), we consider an adaptive Levenberg–Marquardt method and investigate its local convergence near a solution. Specifically, we consider the following Levenberg–Marquardt algorithm.

In order to prove the local convergence of algorithm 1 to some solution $x^{*}\in\Omega$, we assume throughout the paper that the next two conditions hold:

(A1)

There exists some constants $r\in{]0,1[}$, $\lambda>0$, $\beta>0$ and $\delta\in{]0,1]}$ such that the function $h$ is continuously differentiable and Lipschitz continuous with constant $\lambda$ on $\mathbb{B}(x^{*},r)$, and is Hölder metrically subregular of order $\delta$ around $(x^{*},0)$; that is, (7) holds.

(A2)

$\nabla h$ is Hölder continuous of order $\upsilon\in{]0,1]}$ with constant $L>0$ on $\mathbb{B}(x^{*},r)$.

Note that from (A1)-(A2) and the mean value theorem, see (2), it holds

[TABLE]

Let us define the constants

[TABLE]

and

[TABLE]

We begin our study with an analysis inspired by alefeld_rate_2001 , fischer2002local and Guo2015 . The following result provides a bound for the norm of the direction $d_{k}$ based on the distance of the current iteration $x_{k}$ to the solution set $\Omega$. This will be useful later for deducing the rate of convergence of 1.

Proposition 2

If $\xi_{\min}=0$, assume that $\delta>\frac{1}{1+\upsilon}$. Let $x_{k}\not\in\Omega$ be an iteration generated by 1 with $\eta\in{]0,2\delta(1+\upsilon)/\varpi[}$. Then, if $x_{k}\in\mathbb{B}(x^{*},\widetilde{r})$, the direction $d_{k}$ given by (4) satisfies

[TABLE]

where $\delta_{1}:=\min\left\{1+\upsilon-\frac{\eta\varpi}{2\delta},\,1\right\}$ and

[TABLE]

Proof

For all $k$, we will denote by $\overline{x}_{k}$ a vector in $\Omega$ such that $\|x_{k}-\overline{x}_{k}\|=\mathrm{dist}(x_{k},\Omega)$. Since $x_{k}\in\mathbb{B}(x^{*},r/2),$ we have

[TABLE]

which implies $\overline{x}{}_{k}\in\mathbb{B}(x^{*},r)$. Further,

[TABLE]

Observe that $\phi_{k}$ is strongly convex and the global minimiser of $\phi_{k}$ is given by (4). Then, we have

[TABLE]

From the definition of $\phi_{k}$ in (3), by (11) and (14), we deduce

[TABLE]

Let us assume first that $\xi_{\min}>0$. It follows from the definition of $\mu_{k}$ in (6) and (7) that

[TABLE]

leading to

[TABLE]

and this completes the proof of (12) for the case $\xi_{\min}>0$.

Let us consider now the case where $\xi_{\min}=0$, assuming then $\delta>\frac{1}{1+\upsilon}$. By (11), (7) and the Cauchy–Schwarz inequality, we have

[TABLE]

Thus, since $x_{k}\not\in\Omega$, we deduce

[TABLE]

Since $\delta>\frac{1}{1+\upsilon}$, we have

[TABLE]

and therefore

[TABLE]

This, together with the definition of $\mu_{k}$ in (6), implies

[TABLE]

Using (15), we obtain

[TABLE]

which completes the proof.∎

Remark 1

If $\delta>\frac{1}{1+\upsilon}$, by (16), we have that $\nabla h(x_{k})h(x_{k})=0$ implies $x_{k}\in\Omega$ whenever $x_{k}$ is sufficiently close to $x^{*}$.

The next result provides an upper bound for the distance of $x_{k+1}$ to the solution set $\Omega$ based on the distance of $x_{k}$ to $\Omega$.

Proposition 3

If $\xi_{\min}=0$, assume that $\delta>\frac{1}{1+\upsilon}$. Let $x_{k}\not\in\Omega$ and $x_{k+1}$ be two consecutive iterations generated by 1 with $\eta\in{]0,2\delta(1+\upsilon)/\varpi[}$. Then, if $x_{k},x_{k+1}\in\mathbb{B}(x^{*},\widetilde{r})$, we have

[TABLE]

where $\beta_{2}$ is a positive constant and

[TABLE]

Proof

Let $\overline{x}_{k}\in\Omega$ be such that $\|x_{k}-\overline{x}_{k}\|=\mathrm{dist}(x_{k},\Omega)$. From the definition of $\phi_{k}$ in (3) and the reasoning in (15), we obtain

[TABLE]

It follows from (A1) that there exists some constant $\widehat{L}$ such that $\|\nabla h(x)\|\leq\widehat{L}$ for all $x\in\mathbb{B}(x^{*},r)$. Then, by the definition of $\mu_{k}$ in (6) and the Lipschitz continuity of $h$, we have that

[TABLE]

which implies, thanks to (13),

[TABLE]

where $\zeta:=\min\left\{2(1+\upsilon),\,2+\eta\right\}$. By (7), (11), the latter inequality and Proposition 2, we get

[TABLE]

where

[TABLE]

Therefore,

[TABLE]

with $\delta_{2}$ given by (19) and $\beta_{2}:=\frac{1}{\beta}\widehat{\beta}_{2}^{\delta}$, giving the result.∎

The following proposition gives a different value of the exponent in (18).

Proposition 4

Assume that $\delta>\frac{1}{1+\upsilon}$. Let $x_{k}\not\in\Omega$ and $x_{k+1}$ be two consecutive iterations generated by 1 with $\eta\in{]0,2\delta(1+\upsilon)/\varpi[}$ and such that $x_{k},x_{k+1}\in\mathbb{B}(x^{*},\widetilde{r})$. Then, there exists a positive constant $\beta_{3}$ such that

[TABLE]

where

[TABLE]

Proof

Let $\overline{x}_{k},\overline{x}_{k+1}\in\Omega$ be such that $\|x_{k}-\overline{x}_{k}\|=\mathrm{dist}(x_{k},\Omega)$ and $\|x_{k+1}-\overline{x}_{k+1}\|=\mathrm{dist}(x_{k+1},\Omega)$. Assume that $x_{k+1}\not\in\Omega$ (otherwise, the inequality trivially holds). By (11), we have

[TABLE]

Thus, by the Cauchy–Schwarz inequality and (7), we get

[TABLE]

that is,

[TABLE]

Now, by (4), we have

[TABLE]

By (A1) and Proposition 2, it holds,

[TABLE]

It follows from (A1) that there exists some constant $\widehat{L}$ such that $\|\nabla h(x)\|\leq\widehat{L}$ for all $x\in\mathbb{B}(x^{*},r)$. Then, by the definition of $\mu_{k}$ in (6) and (A1), we get (20). Hence, by (24) and Proposition 2, we deduce

[TABLE]

where $\widehat{\beta}_{3}:=L\lambda\beta_{1}^{\upsilon}(\beta_{1}+1)+\widehat{L}L\beta_{1}^{1+\upsilon}{(1+\upsilon)}^{-1}+\left(\xi_{\max}+\omega_{\max}\widehat{L}^{\eta}\right)\lambda^{\eta}\beta_{1}$ and $\widehat{\delta}_{3}:=\min\{\eta+\delta_{1},\delta_{1}(1+\upsilon)\}$. Therefore, by (23),

[TABLE]

Since $\delta>\frac{1}{1+\upsilon}$, we have by (17) that

[TABLE]

Finally, by (25), we deduce

[TABLE]

whence,

[TABLE]

where $\beta_{3}:=\frac{4\widehat{\beta}_{3}}{\beta^{\frac{2}{\delta}}}$ and $\delta_{3}:=\frac{\widehat{\delta}_{3}\delta}{2-\delta}$. Since the expression for $\delta_{3}$ coincides with (22), the proof is complete.∎

Remark 2

(i) The bounds given by (18) and (21) are usually employed to analyse the rate of convergence of the sequence $\left\{x_{k}\right\}$ generated by 1. Observe that the values of $\delta_{2}$ and $\delta_{3}$ when $\xi_{\min}>0$ are greater or equal than their respective values when $\xi_{\min}=0$. A larger value of $\delta_{2}$ or $\delta_{3}$ would serve us to derive a better rate of convergence. To deduce a convergence result from Proposition 3, one needs to have $\delta_{2}>1$. This holds if and only if $\delta>\frac{1}{1+\upsilon}$ and $\eta\in\left]\frac{2}{\delta}-2,\frac{1}{\varpi}\left(2\delta(1+\upsilon)-\frac{2}{1+\upsilon}\right)\right[$, which imposes an additional requirement on the value of $\delta$ (to have a nonempty interval). For instance, when $\upsilon=1$, one must have $\delta>\frac{-1+\sqrt{33}}{8}$ if $\xi_{\min}>0$ and $\delta>\frac{-5+\sqrt{57}}{4}$ if $\xi_{\min}=0$. On the other hand, to guarantee that $\delta_{3}>1$, a stronger requirement would be needed, namely, $\delta>\frac{2}{2+\upsilon}\geq\frac{2}{3}$ and $\eta\in\left]\frac{2}{\delta}-2,\frac{1}{\varpi}\left(2\delta(1+\upsilon)-\frac{4-2\delta}{1+\upsilon}\right)\right[$. Nonetheless, it is important to observe that if $\delta=1$ one has that $\delta_{3}=1+\upsilon$ when $\eta\in[\upsilon,2\upsilon/\varpi]$, while $\delta_{2}=1+\upsilon$ only if $\eta=2\upsilon$ and $\varpi=1$. Therefore, if $\upsilon=\delta=1$, we can derive from Proposition 4 the quadratic convergence of the sequence for $\eta\in[1,2]$, which can only be guaranteed for $\eta=2$ by Proposition 3. In Figure 1, we plot the values of $\delta_{2}$ in Proposition 3 and $\delta_{3}$ in Proposition 4 when $\upsilon=1$ and $\xi_{\min}>0$.

(ii) The values of $\delta_{2}$ and $\delta_{3}$ are maximised when $\eta=\frac{2\upsilon\delta(2+\upsilon)}{\delta+\varpi(1+\upsilon)}$ and $\eta\in\left[\upsilon,\frac{2\upsilon\delta}{\varpi}\right]$, respectively, in which case $\delta_{2}=\delta+{\frac{\upsilon\delta^{2}\left(2+\upsilon\right)}{\delta+\varpi(1+\upsilon)}}$ and $\delta_{3}=\frac{(1+\upsilon)\delta}{2-\delta}$, having then $\delta_{2}\leq\delta_{3}$.

Remark 3

In light of Proposition 1, the extent of the results that can be derived from Propositions 3 and 4 is rather reduced when $x^{*}$ is an isolated solution and $\nabla h(x^{*})$ is not full rank, since it imposes $\delta\leq\frac{1}{1+\upsilon}$. Note that the function $F_{S}$ given as an example in (Guo2015, , Section 5) is Hölder metrically subregular of order $\delta=\frac{5}{6}>0.5$, but $\nabla F_{S}$ is not Lipschitz continuous around any zero of the function, so it does not satisfy (A2) for $\upsilon=1$ (and, therefore, it does not satisfy (Guo2015, , Assumption 4.1) either). However, with the additional assumption that the Łojasiewicz gradient inequality (9) holds, we will obtain local convergence for all $\delta\in{]0,1]}$ (see Theorem 3.2).

Next, we proceed to derive the main result of this section from Propositions 3 and 4, where we provide a region from which the parameter $\eta$ must be chosen so that superlinear convergence is guaranteed. Recall that a sequence $\left\{z_{k}\right\}$ is said to converge superlinearly to $z^{*}$ with order $q>1$ if $z_{k}$ converges to $z^{*}$ and there exists $K>0$ such that $\|z_{k+1}-z^{*}\|\leq K\|z_{k}-z^{*}\|^{q}$ for all $k$ sufficiently large.

Theorem 3.1

Assume that $\delta>\frac{1}{1+\upsilon}$ and $\eta\in\left]\frac{2}{\delta}-2,\frac{1}{\varpi}\left(2\delta(1+\upsilon)-\frac{2}{1+\upsilon}\right)\right[$. Then, there exists some $\overline{r}>0$ such that, for every sequence $\left\{x_{k}\right\}$ generated by 1 with $x_{0}\in\mathbb{B}(x^{*},\overline{r})$, one has that $\left\{\mathrm{dist}(x_{k},\Omega)\right\}$ is superlinearly convergent to [math] with order $\delta_{2}$ given by (19). Further, the sequence $\left\{x_{k}\right\}$ converges to a solution $\overline{x}\in\Omega\cap\mathbb{B}(x^{*},\widetilde{r})$, and if $\eta\leq\frac{2\upsilon\delta}{\varpi}$, its rate of convergence is also superlinear with order $\delta_{2}$. Moreover, if $\delta>\frac{2}{2+\upsilon}$ and $\eta<\frac{1}{\varpi}\left(2(1+\upsilon)\delta-\frac{4-2\delta}{1+\upsilon}\right)$, all the latter holds with order $\delta_{3}$ given by (22).

Proof

We assume that $x_{k}\not\in\Omega$ for all $k$ (otherwise, the statement trivially holds). Let $\delta_{1}$, $\beta_{1}$ be defined as in Proposition 2 and $\delta_{2}$, $\beta_{2}$ be defined as in the proof of Proposition 3. Since $\delta_{2}>1$, we have that $\delta_{1}\delta_{2}^{i}>i$ for all $i$ sufficiently large. As $\sum_{i=1}^{\infty}\left(\frac{1}{2}\right)^{i}=1$, we deduce that

[TABLE]

Define

[TABLE]

Note that $\overline{r}\in{]0,\widetilde{r}[}$, because $\widetilde{r}\in{]0,1[}$ and $\delta_{1}\leq 1$.

Pick any $x_{0}\in\mathbb{B}(x^{*},\overline{r})$ and let $\left\{x_{k}\right\}$ be an infinite sequence generated by 1. First, we will show by induction that $x_{k}\in\mathbb{B}(x^{*},\widetilde{r})$. It follows from $\overline{r}<1$ and (12) that

[TABLE]

Let us assume now that $x_{i}\in\mathbb{B}(x^{*},\widetilde{r})$ for $i=1,2,\dots,k$. Then, from Proposition 3 and the definition of $\overline{r}$, we have

[TABLE]

which yields

[TABLE]

The latter inequality, together with (12), (26) and (27), implies

[TABLE]

which completes the induction. Thus, we have shown that $x_{k}\in\mathbb{B}(x^{*},\widetilde{r})$ for all $k$, as claimed.

From Proposition 3, we obtain that $\left\{\mathrm{dist}(x_{k},\Omega)\right\}$ is superlinearly convergent to [math]. Further, it follows from (12) and (28) that

[TABLE]

Denoting by $s_{k}:=\sum_{i=1}^{k}\|d_{i}\|$, we have that $\left\{s_{k}\right\}$ is a Cauchy sequence. Then, for any $k,p\in\mathbb{N}\cup\{0\}$, we have

[TABLE]

which implies that $\left\{x_{k}\right\}$ is also a Cauchy sequence. Thus, the sequence $\left\{x_{k}\right\}$ converges to some $\overline{x}$. Since $x_{k}\in\mathbb{B}(x^{*},\widetilde{r})$ for all $k$ and $\left\{\mathrm{dist}(x_{k},\Omega)\right\}$ converges to [math], we have $\overline{x}\in\Omega\cap\mathbb{B}(x^{*},\widetilde{r})$.

Further, if $\eta\leq\frac{2\upsilon\delta}{\varpi}$ we have $\delta_{1}=1$ in Proposition 2, and by letting $p\to\infty$ in (29), we deduce

[TABLE]

Since $\{\mathrm{dist}(x_{k},\Omega)\}$ is superlinearly convergent to zero, for all $k$ sufficiently large, it holds that $\mathrm{dist}(x_{k+1},\Omega)\leq\frac{1}{2}\mathrm{dist}(x_{k},\Omega)$. Therefore, for $k$ sufficiently large, we have

[TABLE]

which proves the superlinear convergence of $x_{k}$ to $\overline{x}$ with order $\delta_{2}$.

Finally, the last assertion follows by the same argumentation, using $\delta_{3}$, $\beta_{3}$ and Proposition 4 instead of $\delta_{2}$, $\beta_{2}$ and Proposition 3, respectively.∎

Remark 4

Our results above generalise the results in Guo2015 ; ZhuLin16 , not only because in these works they assume $\nabla h$ to be Lipschitz continuous (i.e., $\upsilon=1$), but also because the parameter $\mu_{k}$ considered by these authors is equal to $\xi\|h(x_{k})\|^{\eta}$. Furthermore, in their convergence results, cf. (Guo2015, , Theorem 4.1 and Theorem 4.2) and (ZhuLin16, , Theorem 2.1 and Theorem 2.2), the authors assume $\delta>\max\left\{\frac{2}{3},\frac{2+\eta}{5}\right\}$ and $\delta>\max\left\{\frac{\sqrt{8\eta+1}+4\eta+1}{16},\frac{2}{2+\eta},\frac{1}{2+\eta}+\frac{\eta}{4},\frac{\eta+1}{4}\right\}>\frac{\sqrt{5}-1}{2}$, respectively, which both entail $\delta>\frac{-1+\sqrt{33}}{8}$, so we have slightly improved the lower bound on $\delta$ for the superlinear convergence in Theorem 3.1.

As a direct consequence of Theorem 3.1, whenever $\delta=\upsilon=1$ and $\eta\in[1,2]$, we can derive quadratic convergence of the sequence generated by 1.

Corollary 1

Assume that $\delta=1$ and $\eta\in{]0,2\upsilon]}$. Then, there exists $\overline{r}>0$ such that for every sequence $\left\{x_{k}\right\}$ generated by 1 with $x_{0}\in\mathbb{B}(x^{*},\overline{r})$, one has that $\left\{\mathrm{dist}(x_{k},\Omega)\right\}$ is superlinearly convergent to [math] with order

[TABLE]

Moreover, the sequence $\left\{x_{k}\right\}$ converges superlinearly with order $\delta_{3}$ to a solution $\overline{x}\in\Omega\cap\mathbb{B}(x^{*},\widetilde{r})$. Therefore, when $\upsilon=1$ and $\eta\in[1,2]$, the sequence $\left\{\mathrm{dist}(x_{k},\Omega)\right\}$ is quadratically convergent to [math], and the sequence $\left\{x_{k}\right\}$ converges quadratically to a solution $\overline{x}\in\Omega\cap\mathbb{B}(x^{*},\widetilde{r})$.

Remark 5

In particular, Corollary 1 generalizes (MJ07, , Theorem 3.7), where the authors prove quadratic convergence of the sequence $\{x_{k}\}$ by assuming $\delta=\upsilon=1$, and where the parameters in (6) are chosen as $\eta=1$, $\xi_{k}=\theta\in{[0,1]}$ and $\omega_{k}=1-\theta$, for all $k$.

Example 3 (Example 2 revisited)

Let $h$ and $\widehat{h}$ be the functions defined in Example 2. The function $h$ does not satisfy the assumptions of Theorem 3.1, since $\delta=\frac{1}{1+\upsilon}$. On the other hand, if $\widehat{\eta}\in\left]0,\frac{7}{6}\right[$ and the starting point $x_{0}$ is chosen sufficiently close to [math], Theorem 3.1 proves for the function $\widehat{h}$ the superlinear convergence of the sequence generated by 1 to [math] with order

[TABLE]

Note that, since the solution is locally unique, the additional assumption $\widehat{\eta}\leq 2\widehat{\upsilon}\hskip 1.0pt\widehat{\delta}=\frac{2}{3}$ is not needed. The order of convergence $\delta_{3}$ is thus maximised when $\widehat{\eta}\in\left[\frac{1}{3},\frac{2}{3}\right]$.$\Diamond$

The question of whether the sequence $\left\{\mathrm{dist}(x_{k},\Omega)\right\}$ converges to [math] when $\delta$ does not satisfy the requirements commented in Remark 2(i) remains open. However, with the additional assumption that $\psi$ satisfies the Łojasiewicz gradient inequality (which holds for real analytic functions), we can prove that the sequences $\left\{\mathrm{dist}(x_{k},\Omega)\right\}$ and $\left\{\psi(x_{k})\right\}$ converge to [math] for all $\delta\in(0,1]$ as long as the parameter $\eta$ is sufficiently small, and we can also provide a rate of convergence that depends on the exponent of the Łojasiewicz gradient inequality. This is the subject of the next subsection.

3.1 Convergence analysis under the Łojasiewicz gradient inequality

To prove our convergence result, we make use of the following two lemmas.

Lemma 1

Let $\left\{s_{k}\right\}$ be a nonnegative real sequence and let $\alpha,\vartheta$ be some nonnegative constants. Suppose that $s_{k}\to 0$ and that the sequence satisfies

[TABLE]

Then

1. (i)

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

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

if $\alpha>1$, there exists $\varsigma>0$ such that

[TABLE]

Proof

See (artacho_accelerating_2015, , Lemma 1).∎

Lemma 2

The sequence $\{x_{k}\}$ generated by 1 satisfies

[TABLE]

and

[TABLE]

Proof

This result is a straightforward modification of (karas2015algebraic, , Theorem 2.5 and Lemma 2.3), using (2) instead of the Lipschitz continuity of $\nabla h$.∎

In our second main result of this paper, under the additional assumption that the Łojasiewicz gradient inequality holds, we prove the convergence to [math] of the sequences $\left\{\mathrm{dist}(x_{k},\Omega)\right\}$ and $\left\{\psi(x_{k})\right\}$.

Theorem 3.2

Suppose that $\psi$ satisfies the Łojasiewicz gradient inequality (9) with exponent $\theta\in{]0,1[}$. Let

[TABLE]

Then, if $\eta\in{\left]0,\min\left\{\frac{2\upsilon}{\chi(1+\upsilon)},\frac{2(1-\theta)}{\chi}\right\}\right[}$, there exist some positive constants $s$ and $\overline{s}$ such that, for every $x_{0}\in\mathbb{B}(x^{*},s)$ and every sequence $\{x_{k}\}$ generated by 1, one has $\{x_{k}\}\subset\mathbb{B}(x^{*},\overline{s})$ and the two sequences $\{\psi(x_{k})\}$ and $\{\mathrm{dist}(x_{k},\Omega)\}$ converge to [math]. Moreover, the following holds:

1. (i)

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

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]

Proof

The proof has three key parts.

In the first part of the proof, we will set the values of $s$ and $\overline{s}$. Let $\varepsilon>0$ and $\kappa>0$ be such that (9) holds. Thus, one has

[TABLE]

Let $\overline{s}:=\min\{r,\varepsilon\}>0$. Then, by Assumption (A1), there exists some positive constant $M$ such that

[TABLE]

Since $\eta\in{\left]0,\frac{2\upsilon}{\chi(1+\upsilon)}\right[}$, it is possible to make $\overline{s}$ smaller if needed to ensure, for all $x\in\mathbb{B}(x^{*},\overline{s})$, that

[TABLE]

For all $x\in\mathbb{B}(x^{*},\overline{s})$, one has by the Lipschitz continuity of $h$ that

[TABLE]

since $\overline{s}\leq r<1$. Let

[TABLE]

Then, since $\overline{s}<1$ and $\theta+\frac{\eta\chi}{2}\in{]0,1[}$, we have $s\leq\overline{s}$.

In the second part of the proof, we will prove by induction that

[TABLE]

for all $i=1,2,\ldots$. Pick any $x_{0}\in\mathbb{B}(x^{*},s)$ and let $\{x_{k}\}$ be the sequence generated by 1. It follows from Lemma 2 that

[TABLE]

for all $k$, where $H_{k}=\nabla h(x_{k})\nabla h(x_{k})^{T}+\mu_{k}I$, since $d_{k}=-H_{k}^{-1}\nabla h(x_{k})h(x_{k})$. Since $x_{0}\in\mathbb{B}(x^{*},\overline{s})$, we have by (31), the definition of $\chi$ in (30) and (33) that

[TABLE]

which implies

[TABLE]

Therefore, from (37), we get

[TABLE]

Observe that the convexity of the function $\varphi(t):=-t^{1-\theta-\frac{\eta\chi}{2}}$ with $t>0$ yields

[TABLE]

for all $x,y\in\mathbb{R}^{m}\setminus\Omega$. By combining (39) with (40), we deduce

[TABLE]

Since $x_{0}\in\mathbb{B}(x^{*},s)\subseteq\mathbb{B}(x^{*},\overline{s})$, we have by (32) that $\|H_{0}\|\leq M$. Further, by the Łojasiewicz gradient inequality (9), it holds

[TABLE]

From the last inequality, together with (41), the first inequality in (38) and then (34), we obtain

[TABLE]

which, in particular, proves (36) for $i=1$. Hence,

[TABLE]

Therefore, $x_{1}\in\mathbb{B}(x^{*},\overline{s})$. Assume now that (35)–(36) holds for all $i=1,\ldots,k$. Since $x_{k}\in\mathbb{B}(x^{*},\overline{s})$, by (33) and the same argumentation as in (38), we have

[TABLE]

which implies

[TABLE]

Therefore, by (37), we get

[TABLE]

Combining the latter inequality with (40), we deduce

[TABLE]

Further, since $x_{k}\in\mathbb{B}(x^{*},\overline{s})$, from the Łojasiewicz gradient inequality (9) and (32), it holds

[TABLE]

From the last inequality and (43), we deduce

[TABLE]

which proves (36) for $i=k+1$. Hence, by (34), we have

[TABLE]

which proves (35) for $i=k+1$. This completes the second part of the proof.

In the third part of the proof, we will finally show the assertions in the statement of the theorem. From the second part of the proof we know that $x_{k}\in\mathbb{B}(x^{*},\overline{s})$ for all $k$. This, together with (32), implies that $\|H_{k}\|\leq M$ for all $k$. Thus,

[TABLE]

Therefore, by (42), we have

[TABLE]

It follows from the Łojasiewicz gradient inequality (9) and the last inequality that

[TABLE]

This implies that $\{\psi(x_{k})\}$ converges to [math]. By applying Lemma 1 with $s_{k}:=\psi(x_{k})$, $\vartheta:=2\kappa^{2}M$ and $\alpha:=2\theta$, we conclude that the rate of convergence depends on $\theta$ as claimed in (i)-(ii). Finally, observe that $\{\mathrm{dist}\left(x_{k},\Omega\right)\}$ converges to [math] with the rate stated in (i)-(ii) thanks to the Hölder metric subregularity of the function $h$.∎

Remark 6

Observe that every real analytic function satisfies the assumptions of Theorem 3.2, thanks to Fact 1 and the discussion after it in Section 2. Therefore, local sublinear convergence of 1 is guaranteed for all $\eta$ sufficiently small (i.e., whenever $\eta<\min\left\{\chi^{-1},2(1-\theta)\chi^{-1}\right\}$). This is the best that we can get with these weak assumptions, as we show in the next example.

Example 4 (Example 2 revisited)

Let $h(x)=\frac{3}{4}\sqrt[3]{x^{4}}$ be the function considered in Example 2. The function $h$ does not satisfy the assumptions of Theorem 3.1, but it verifies the ones of Theorem 3.2. Indeed, it is straightforward to check that $\psi(x)=\frac{1}{2}|h(x)|^{2}$ satisfies the Łojasiewicz gradient inequality (9) with exponent $\theta=\frac{5}{8}$. Since $\theta>\frac{1}{2}$, we can only guarantee the sublinear convergence of the sequence $\{x_{k}\}$ generated by 1 to [math] when $\eta\in{\left]0,\frac{1}{2\chi}\right[}={\left]0,\min\left\{\frac{1}{2\chi},\frac{3}{4\chi}\right\}\right[}$. In fact, this is the best convergence rate that we can get. Indeed, a direct computation gives us

[TABLE]

On the one hand, when $\xi_{min}>0$ and $\eta\in{\left]0,\frac{1}{2}\right[}$, we have $\frac{4\eta}{3}<\frac{2}{3}$. Therefore, it follows from (44) and $\xi_{k}\geq\xi_{\min}>0$ that

[TABLE]

which means that $\{x_{k}\}$ is sublinearly convergent to [math]. This coincides with what Theorem 3.2 asserts, since $\left]0,\frac{1}{2\chi}\right[=\left]0,\frac{1}{2}\right[$. On the other hand, when $\xi_{\min}=0$ and $\eta\in\left]0,\frac{2}{5}\right[$, sublinear convergence is also obtained from (44), which is exactly what Theorem 3.2 guarantees for all $\eta\in\left]0,\frac{1}{2\chi}\right[=\left]0,\frac{2}{5}\right[$.

4 Application to biochemical reaction networks

In this section, we introduce first a class of nonlinear equations arising in the study of biochemistry, cf. Fleming2012 . After that, we compare the performance of 1 with various Levenberg–Marquardt algorithms for finding steady states of nonlinear systems of biochemical networks on 20 different real data biological models.

4.1 Nonlinear systems in biochemical reaction networks

Consider a biochemical 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, that is, there exists at least one positive vector $l\in\mathbb{R}_{++}^{m}$ satisfying $(R-F)^{T}l=0$, cf. gevorgyan2008detection . The matrix $N:=R-F$ represents net reaction stoichiometry and may be viewed as the incidence matrix of a directed hypergraph, see Klamt2009 . We assume that there are less molecular species than there are net reactions, that is $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 particular sparsity pattern depends on the particular biochemical network being modeled. Moreover, we also assume that $\text{rank}([F,R])=m$, which is a requirement for kinetic consistency, cf. Fleming20161 .

Let $c\in\mathbb{R}_{++}^{m}$ denote a variable vector of molecular species concentrations. Assuming constant nonnegative elementary kinetic parameters $k_{f},k_{r}\in\mathbb{R}_{+}^{n}$, we assume 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., artacho_accelerating_2015 ; Fleming20161 . Then, the deterministic dynamical equation for time evolution of molecular species concentration is given by

[TABLE]

A vector $c^{*}$ is a steady state if and only if it satisfies

[TABLE]

Note that a vector $c^{*}$ is a steady state of the biochemical system if and only if

[TABLE]

here $\mathcal{N}(N)$ denotes the null space of $N$. Therefore, the set of steady states $\Omega=\left\{c\in\mathbb{R}_{++}^{m},\,f(c)=0\right\}$ is unchanged if we replace the matrix $N$ by a matrix $\bar{N}$ with the same null space. 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 (4.1) to 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 (4.1) is equal to the function

[TABLE]

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

Let $L\in\mathbb{R}^{(m-r)\times m}$ denote a basis for the left null space of $N$, which implies $LN=0$. We have $\mathrm{rank}(L)=m-r$. We say that the system satisfies moiety conservation if for any initial concentration $c_{0}\in\mathbb{R}_{++}^{m}$, it holds

[TABLE]

along the trajectory of (4.1), given an initial starting point $l_{0}\in\mathbb{R}_{++}^{m}$. It is possible to compute $L$ such that each row corresponds to a structurally identifiable conserved moiety in a biochemical network, cf. Haraldsdttir2016 . The problem of finding the moiety conserved steady state of a biochemical reaction network is equivalent to solving the nonlinear equation (1) with

[TABLE]

By replacing $f$ by $\bar{f}$ we have improved the rank deficiency of $\nabla f$, and thus the one of $h$ in (46). Nonetheless, as we demonstrate in Figure 5, $\nabla h$ is usually still far from being full rank at the solutions.

Let us show that $h$ is real analytic. Let $A:=\left[\bar{N},-\bar{N}\right]$ and $B:=[F,\,R]^{T}$. Then we can write

[TABLE]

where $Q=A^{T}A.$ Since $B_{ij}$ are nonnegative integers for all $i$ and $j$, we conclude that the function $\psi$ is real analytic (see Proposition 2.2.2 and Proposition 2.2.8 in parks1992primer ). It follows from Remark 6 that $\psi$ satisfy the Łojasiewicz gradient inequality (with some unknown exponent $\theta\in{[0,1[}$) and the mapping $h$ is Hölder metrically subregular around $(x^{*},0)$. Therefore, the assumptions of Theorem 3.2 are satisfied as long as $\eta$ is sufficiently small, and local sublinear convergence of 1 is guaranteed.

4.2 Computational experiments

In this subsection, we compare 1 with various Levenberg–Marquardt methods for solving the nonlinear system (1) with $h$ defined by (46) on 20 different biological models. These codes are available in the COBRA Toolbox v3 COBRA . In our implementation, all codes were written in MATLAB and runs were performed on Intel Core i7-4770 CPU 3.40GHz with 12GB RAM, under Windows 10 (64-bits). The algorithms were stopped whenever

[TABLE]

is satisfied or the maximum number of iterations (say 10,000) is reached. On the basis of our experiments with the mapping (46), we set

[TABLE]

The initial point is set to $x_{0}=0$ in all the experiments.

To illustrate the results, we use the Dolan and Moré performance profile Dolan2002 with the performance measures $N_{i}$ and $T$, where $N_{i}$ and $T$ denote the total number of iterations and the running time. 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 (48). 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 our first experiment, we explore for which parameter $\eta$ the best performance of 1 is obtained. To this end, we apply seven versions of 1 associated to each of the parameters $\eta\in\{0.6,0.7,0.8,0.9,0.99,\allowbreak 0.999,1\}$ to the nonlinear system (46) defined by 20 biological models. The results of this comparison are summarised in Table 1 and Figure 2, from where it can be observed that 1 with $\eta=0.999$ outperforms the other values of the parameters. It is also apparent that smaller values of $\eta$ are less efficient, although 1 successfully found a solution for every model and every value of $\eta$ that was tested. It is important to recall here that local convergence is only guaranteed by Theorem 3.2 for sufficiently small values of $\eta$, since the value of $\theta$ is unknown. Also, note that the local convergence for the value $\eta=1$ is not covered by Theorem 3.2 for our choice of the parameters, because it requires $\eta<\min\{1,2-2\theta\}$, since $\omega_{\min}=0$ in (47).

We now set $\eta=0.999$ and compare 1 with parameters (47) with the following Levenberg–Marquardt methods:

• $\bullet$

LM-YF: with $\mu_{k}=\|h(x_{k})\|^{2}$, given by Yamashita and Fukushima alefeld_rate_2001 ;

• $\bullet$

LM-FY: with $\mu_{k}=\|h(x_{k})\|$, given by Fan and Yuan fan_quadratic_2005 ;

• $\bullet$

LM-F: with $\mu_{k}=\|\nabla h(x_{k})h(x_{k})\|$, given by Fischer fischer2002local .

It is clear that all of these three methods are special cases of 1 by selecting suitable parameters $\xi_{k}$, $\omega_{k}$, and $\eta$. The results of our experiments are summarised in Table 2 and Figure 3. In Figures 3(a) and 3(b), we see that 1 is clearly always the winner, both for the number of iterations and the running time. Moreover, LM-F outperforms both LM-YF and LM-FY. In fact, LM-FY was not able to solve any of the considered problems within the 10,000 iterations.

In order to see the evolution of the merit function, we illustrate its value with respect to the number of iterations in Figure 4 for the mapping (46) with the biological models iAF692 and iNJ661. We limit the maximum number of iterations to 1,000. Clearly, 1 attains the best results, followed by LM-F. Both methods seem to be more suited to biological problems than LM-YF and LM-FY. We also show in Figure 4 the evolution of the value of the step size $\|d_{k}\|$. Both 1 and LM-F show a rippling behaviour, while the value of $\|d_{k}\|$ is nearly constant along the 1,000 iterations for LM-YF and LM-FY. Probably, this rippling behaviour is letting the first two methods escape from a flat valley of the merit function, while the two last methods get trapped there. Observe also that, by Lemma 2, one has that $\|d_{k}\|\leq\frac{1}{2}$ for LM-YF and $\|d_{k}\|\leq\frac{1}{2}\|h(x_{k})\|^{\frac{1}{2}}$ for LM-FY, while this upper bound can be larger for both 1 and LM-F.

In our last experiment, we find $10$ solutions of the nonlinear system (1) with 1 using $10$ random starting points $x_{0}\in{\left]-\frac{1}{2},\frac{1}{2}\right[^{m}}$ for each of the $20$ biological models and compute the rank of $\nabla h$ at each of these solutions. The results are shown in Figure 5, where we plot the rank deficiency of $\nabla h$ at each of the solutions. For all the models, except for the Ecoli_core, we observe that $\nabla h$ at the solutions found is far from being full rank. For the Ecoli_core, although $\nabla h$ had full rank at every solution found, the smallest eigenvalue at these solutions had a value around $10^{-9}$, making also this problem ill-conditioned. This explains the difficulties that most of the algorithms had for solving the nonlinear system (1) with $h$ defined by (46). Therefore, since we are dealing with a difficult problem, it is more meritorious the successfulness of 1 with parameters (47) for finding a solution of each of the $20$ models in less than 400 iterations (in less than one minute), as shown in Table 2.

5 Conclusion and further research

We have presented an adaptive Levenberg–Marquardt method for solving systems of nonlinear equations with possible non-isolated solutions. We have analysed its local convergence under Hölder metric subregularity of the underlying function and Hölder continuity of its gradient. We have further analysed the local convergence under the additional assumption that the Łojasiewicz gradient inequality holds. These properties hold in many applied problems, as they are satisfied by any real analytic function. One of these applications is computing a solution to a system of nonlinear equations arising in biochemical reaction networks, a problem which is usually ill-conditioned. We showed that such systems satisfy both the Hölder metric subregularity and the Łojasiewicz gradient inequality assumptions. In our numerical experiments, we clearly obtained a superior performance of our regularisation parameter, compared to existing Levenberg–Marquardt methods, for 20 different biological networks.

Several extensions to the present study are possible, the most important of which would be to develop a globally convergent version of the proposed Levenberg–Marquardt method. One approach, which is currently being investigated, would be to combine the scheme with an Armijo-type line search and a trust-region technique. This will be reported in a separate article AFV19 . It would also be interesting to analyse a regularisation parameter where the value of $\eta$ is updated at each iteration. The analysis of the convergence with such a parameter would be much more involved, so we leave this for future work.

Acknowledgements

We would like to thank Mikhail Solodov for suggesting the use of Levenberg–Marquardt methods for solving the system of nonlinear equations arising in biochemical reaction networks. Thanks also go to Michael Saunders for his useful comments on the first version of this manuscript. We are grateful to two anonymous reviewers for their constructive comments, which helped us improving the paper.

Appendix

See Tables 1 and 2 for the summary results of the comparisons.