A New Adaptive Accelerated Levenberg–Marquardt Method for Solving Nonlinear Equations and Its Applications in Supply Chain Problems

Abstract

In this paper, a new adaptive Levenberg–Marquardt method is proposed to solve the nonlinear equations including supply chain optimization problems. We present a new adaptive update rule which is a segmented function on the ratio between the actual and predicted reductions of the objective function to accept a large number of unsuccessful iterations and avoid jumping in local areas. The global convergence and quadratic convergence of the proposed method are proved by using the trust region technique and local error bound condition, respectively. In addition, we use the proposed algorithm to test on the symmetric and asymmetric linear equations. Numerical results show that the proposed method has good numerical performance and development prospects. Furthermore, we apply the algorithm to solve the fresh agricultural products supply chain optimization problems.

1. Introduction

With the development of science and technology, more and more fields are involved in the solution of nonlinear equation problems, such as chemistry, mechanics, economics and product management [1,2,3,4]. For example, decentralized decision models in supply chain management and gas pressure volume models in physics can be converted into the following nonlinear equations

$$F\left(x\right)=0,$$

(1)

where $F\left(x\right):{\mathbb{R}}^n\to {\mathbb{R}}^m$ is a continuously differentiable function. In particular, symmetric nonlinear equations with the Jacobian matrix symmetry also have a wide range of applications, such as the gradient mapping of unconstrained optimization problem, the Karush–Kuhn–Tucker (KKT) of equality constrained optimization problem, and other fields [5,6].

The steepest descent method, Newton method, quasi-Newton method, Gauss–Newton (GN) method are commonly used iterative methods for solving (1) [7,8,9,10]. The GN method is one of the most famous methods, when the Jacobian matrix is Lipschitz continuous and nonsingular at the solution of (1), the GN method has quadratic convergence. However, when the Jacobian matrix is singular or nearly singular, the GN method may not be well defined. In order to overcome this difficulty, the Levenberg–Marquardt (LM) method [11,12] for solving (1) was proposed. At the k-th iteration, the trial step is

$$d_k=-{(J_k^{\mathrm{T}}{J}_k+{\lambda }_{k}I)}^{-1}\left(J_k^{\mathrm{T}}{F}_k\right),$$

(2)

where $F_k=F\left(x_k\right)$, $J_k=J\left(x_k\right)$ is a Jacobian matrix of $F\left(x\right)$ at $x_k$, which may be a symmetric matrix or non-symmetric matrix, I is an identity matrix and the LM parameter ${\lambda }_k>0$.

The LM method ensures the uniqueness of solution of (1), and it also has quadratic convergence if $J_k$ is Lipschitz continuous, nonsingular at the solution, and ${\lambda }_k$ is selected appropriately. In this sense, the update of the LM parameter has a great impact on the performance and efficiency of algorithm, many effective LM parameters have been proposed. Yamashita and Fukushima [13] chose the LM parameter as ${\lambda }_k={\parallel F_k\parallel }^2$, and proved that the LM method had quadratic convergence under the local error bound condition and $J_k$ is Lipschitz continuous at the solution. However, when $\left\{x_k\right\}$ is far away from the solution set, ${\lambda }_k$ may be very large, which makes $d_k$ very small and reduces the efficiency of algorithm; when $\left\{x_k\right\}$ is sufficiently close to the solution set, ${\lambda }_k$ may be smaller than the machine epsilon and lose its role.

Based on these observations, Fan and Yuan [14] generalized the LM parameter in [13], and proved that the numerical results for choosing ${\lambda }_k=\parallel F_k\parallel$ is better than choosing ${\lambda }_k={\parallel F_k\parallel }^2$. Fan [15] first introduced the regularization factor ${\mu }_k$ into the LM method and chose ${\lambda }_k={\mu }_k\parallel F_k\parallel$, with numerical results showing that this choice of ${\lambda }_k$ provides the best performance. However, when $\left\{x_k\right\}$ is far away from the solution set, the choice of both LM parameters does not provide good results. Therefore, to avoid this situation, Fan and Pan [16] chose the LM parameter as ${\lambda }_k={\mu }_k\rho \left(x_k\right)$, in which ${\mu }_k$ is updated by a trust region technique. They defined $\rho \left(x_k\right)$ as a positive function of ${\mathbb{R}}^n\to {\mathbb{R}}_{+}$, i.e.,

$$\begin{array}{c}\hfill \rho \left(x_k\right)=\left\{\begin{array}{cc}\tilde{\rho }\left(x_k\right),\hfill & \mathrm{if}\tilde{\rho }\left(x_k\right)\le 1,\hfill \\ 1,\hfill & \mathrm{otherwise},\hfill \end{array}\right.\end{array}$$

where $\tilde{\rho }\left(x_k\right)=O(\parallel F_k{\parallel }^{\delta })$. This update strategy can obtain larger LM trial steps, so that the iterative sequence can quickly converge to the solution set when $\left\{x_k\right\}$ is far away from the solution set. Amini et al. [17] chose the LM parameter as

$$\begin{array}{c}\hfill {\lambda }_k=\frac{{\mu }_k\parallel F_k\parallel }{1+\parallel F_k\parallel }.\end{array}$$

It is clear that when $\left\{x_k\right\}$ is far away from the solution set and $\parallel F_k\parallel$ is very large, $\frac{\parallel F_k\parallel }{1+\parallel F_k\parallel }$ is close to 1, so ${\lambda }_k$ is close to ${\mu }_k$. The choice of ${\lambda }_k$ speeds up the efficiency of the algorithm more than previous LM parameters.

In addition to the above different choices of LM parameters, the introduction of adaptive technology also has a great impact on the LM method. As we all know, the ratio $r_k$ between the actual and predicted reductions of the objective function reflects the degree to which the approximate quadratic model approaches the value function. To make more use of information about the ratio, Fan and Yuan [18] proposed an adaptive LM method by selecting ${\lambda }_{k+1}={\mu }_{k+1}{\parallel F_{k+1}\parallel }^{\delta }$, ${\mu }_{k+1}={\mu }_{k}q\left(r_k\right)$, $q\left(r_k\right)$ is a continuous non-negative function about $r_k$, and $\sigma \in (0,2]$. The introduction of $q\left(r_k\right)$ avoids discontinuities when crossing the threshold $\frac{{\mu }_{k+1}}{{\mu }_k}$ of the ratio, and better numerical results can be obtained.

In fact, similar adaptive techniques have been proposed in the trust region algorithms. If $r_k$ is sufficiently greater than 1, the iteration is too successful at this time, then we can reduce ${\mu }_k$ to a very small value. Then, the algorithm will continue to perform a large number of consecutive unsuccessful iterations. On the other hand, if $r_k\to -\infty$, $d_k$ is a far-from-satisfactory trial step, then we can increase ${\mu }_k$ greatly. At this moment, the successive iteration points will be close to each other and the algorithm will converge slowly. Therefore, Hei [19] proposed an R-function by using an adaptive update strategy to update the trust region radius ${\Delta }_k$, i.e., ${\Delta }_{k+1}=R\left(r_k\right){\Delta }_k$. Furthermore, Walmag and Delhez [20] proposed a $\mathsf{\Lambda }$-function to update the trust region radius, i.e., ${\Delta }_{k+1}=\mathsf{\Lambda }\left(r_k\right){\Delta }_k$, where $\mathsf{\Lambda }$ is a non-negative and bounded function about $r_k$. On this basis, Lu et al. [21] argued that the consistency between the model and the objective function is not good enough in too-successful iterations, so an L-function was proposed to update the trust region radius. They showed that the L-function contains some favorable features of the R-function and the $\mathsf{\Lambda }$-function, and the method is more efficient in too-successful iterations. In this paper, we want to learn from the presentation of the L-function and provide a new adaptive strategy to update the LM parameter. Our innovations mainly include the following:

◊ A new adaptive accelerated LM method is proposed, which can improve the consistency between the model and the objective function in too-successful iterations by using the ratio information of the actual reduction to the predicted reduction;

◊ The new algorithm can solve the situation in which the iterative sequence is far away from the optimal solution set, accept a large number of unsuccessful iterations and avoid jumping in local areas, thus improving the efficiency and stability of the algorithm;

◊ The new adaptive accelerated LM method has global convergence and quadratic convergence under local error bound.

The rest of this paper is organized as follows. In Section 2, we describe in detail a new adaptive accelerated LM method which makes full use of the ratio information. Furthermore, we demonstrate that the new algorithm has global convergence under the appropriate conditions and maintains quadratic convergence under local error bound condition. In Section 3, numerical results are given, indicating that the new algorithm is efficient. The conclusion is given in the last section.

2. Methodology

2.1. The Adaptive Accelerated Levenberg–Marquardt Method

In this section, our main aim is to discuss how to update the LM parameter to propose a new adaptive accelerated LM method. It is easy to see from (2) that $d_k$ is the solution to the optimization problem

$$\underset{d\in {\mathbb{R}}^n}{min}\parallel F_k+J_k{d\parallel }^2+{\lambda }_k{\parallel d\parallel }^2\doteq {\psi }_k\left(d\right).$$

(3)

If

$${\Delta }_k=\parallel -{(J_k^{\mathrm{T}}{J}_k+{\lambda }_{k}I)}^{-1}\left(J_k^{\mathrm{T}}{F}_k\right)\parallel ,$$

(4)

then $d_k$ is also the solution of the subproblem

$$\begin{array}{ccc}\hfill & & \underset{d\in {\mathbb{R}}^n}{min}{\parallel F_k+J_{k}d\parallel }^2\doteq {\phi }_k\left(d\right),\hfill \\ & & \mathrm{s}.\mathrm{t}.\parallel d\parallel \le {\Delta }_k.\hfill \end{array}$$

(5)

Therefore, the LM method can be regarded as a trust region method, which implicitly modifies the trust region radius ${\Delta }_k$. The difference between the general trust region method and the LM method is that the LM method does not directly update the trust region radius, but updates the regularization factor ${\mu }_k$.

We define the actual reduction and predicted reduction of the merit function $\parallel F_k{\parallel }^2$ at the k-th iteration as

$$Are{d}_k=\parallel F_k{\parallel }^2-{\parallel F(x_k+d_k)\parallel }^2$$

(6)

and

$$Pre{d}_k={\phi }_k\left(0\right)-{\phi }_k\left(d_k\right).$$

(7)

The ratio between the actual and predicted reductions of the objective function is defined by

$$r_k=\frac{Are{d}_k}{Pre{d}_k}.$$

(8)

This ratio determines whether the trial step $d_k$ is accepted. Here, we choose the LM parameter as

$$\begin{array}{c}\hfill {\lambda }_{k+1}=\frac{{\mu }_{k+1}\parallel F_{k+1}\parallel }{1+\parallel F_{k+1}\parallel }.\end{array}$$

(9)

The usual empirical rules [22,23,24,25] of ${\mu }_{k+1}$ can be usually summarized as follows

$${\mu }_{k+1}=\left\{\begin{array}{cc}4{\mu }_k,\hfill & \mathrm{if}{r}_k<p_1,\hfill \\ {\mu }_k,\hfill & \mathrm{if}{p}_1\le r_k\le p_2,\hfill \\ \max \{\frac{{\mu }_k}{4},m\},\hfill & \mathrm{if}{r}_k>p_2,\hfill \end{array}\right.$$

(10)

where $m>0$ and $0<p_1<p_2<1$ are constants.

Iterations with $r_k$ greater than $p_2$ are very successful iterations. In this case, it is usually assumed that the approximation of the model function to the objective function is accurate and ${\mu }_k$ should be reduced. However, at too-successful iterations, i.e., $r_k$ is sufficiently greater than 1, the consistency between the model and the objective function is not good enough. Thus, we use an adaptive strategy to update the factor ${\mu }_{k+1}$, i.e., ${\mu }_{k+1}=K\left(r_k\right){\mu }_k$, where $K\left(r_k\right)$ is a function about $r_k$.

We construct $K\left(r_k\right)$ as follows:

$$K\left(r_k\right)=\left\{\begin{array}{cc}{\beta }_1+({\beta }_2-{\beta }_1)\exp (-{(\frac{-r_k+p_1}{p_1})}^2),\hfill & \mathrm{if}{r}_k\le p_1,\hfill \\ {\beta }_2,\hfill & \mathrm{if}{p}_1<r_k<p_2,\hfill \\ \frac{1-{\beta }_3\exp \left(p_2\right)}{1-\exp \left(p_2\right)}-\frac{(1-{\beta }_3)\exp \left(p_2\right)}{1-\exp \left(p_2\right)}\exp (-r_k+p_2)-\frac{1}{2},\hfill & \mathrm{if}{r}_k\ge p_2,\hfill \end{array}\right.$$

(11)

where $0<{\beta }_2<1<{\beta }_1\le {\beta }_3$ and $0<p_1<p_2<1$ are constants. Here, $K\left(r_k\right)$ satisfies the following properties   

$\left(1\right)\underset{r_k\to -\infty }{lim}K\left(r_k\right)={\beta }_1;$

$\left(2\right)\underset{r_k\to p_1}{lim}K\left(r_k\right)={\beta }_2;$

$\left(3\right)\underset{r_k\to p_2}{lim}K\left(r_k\right)=\frac{1}{2};$

$\left(4\right)\underset{r_k\to +\infty }{lim}K\left(r_k\right)=\frac{1-{\beta }_3\exp \left(p_2\right)}{1-\exp \left(p_2\right)}-\frac{1}{2}.$   

If we obtain a satisfactory trial step $d_k$ and ratio $r_k$, then we accept trial step $d_k$ and reduce ${\mu }_k$; otherwise, we reject trial step $d_k$ and increase ${\mu }_k$. At too-successful iterations, the actual reduction of the objective function obtained at iteration k is obviously greater than the predicted reduction. Although the current iteration allows the algorithm to progress towards the optimum, the approximation of the model function to the objective function is bad. Therefore, to avoid reducing ${\mu }_k$ too quickly, we use the K-function to update ${\mu }_k$.

According to the properties of the K-function, the rate of ${\mu }_k$ reduction is the fastest when $r_k$ is close to 1, i.e., when the model function provides an accurate local approximation of the objective function. The new idea we propose is to allow ${\mu }_k$ to be updated at a variable rate according to $r_k$, which would improve the efficiency and stability of the algorithm.

Based on the above analysis, we state a description of the new adaptive accelerated LM method (Algorithm 1) as follows.

In Algorithm 1, m is a given lower bound of the parameter ${\mu }_k$. It is introduced to prevent the step from being too large when the sequence is near the solution.

Algorithm 1 NAALM.

0. Given $x_0\in {\mathbb{R}}^n$, ${\mu }_0>m>0$, $0\le p_0<p_1<p_2<1$, $0<{\beta }_2<1<{\beta }_1\le {\beta }_3$, $\epsilon >0$. Let $k:=0$. 1. Compute $F_k$ and $J_k$. If $\parallel J_k^{\mathrm{T}}{F}_k\parallel \le \epsilon$, stop. Otherwise, compute ${\lambda }_k$ by (9). 2. Solving the following system $$(J_k^{\mathrm{T}}{J}_k+{\lambda }_{k}I)d=-J_k^{\mathrm{T}}{F}_k$$ (12) to determine $d_k$. 3. Compute $Pre{d}_k$, $Are{d}_k$ and $r_k$ by (6)–(8), respectively. 4. Set $$x_{k+1}=\left\{\begin{array}{cc}{x}_k+d_k,& \mathrm{if}{r}_k\ge p_0,\\ x_k,& \mathrm{if}{r}_k<p_0.\end{array}\right.$$ (13) 5. Choose ${\mu }_{k+1}$ as $${\mu }_{k+1}=max\{m,K\left(r_k\right){\mu }_k\},$$ (14) where $K\left(r_k\right)$ is given by (11). Set $k:=k+1$ and go to Step 1.

2.2. The Global Convergence

In this section, to obtain the global convergence of NAALM algorithm, we make the following assumption.

Assumption 1.

$F\left(x\right)$ is continuously differentiable, $F\left(x\right)$ and the Jacobian matrix $J\left(x\right)$ are Lipschitz continuous, i.e., there exist positive constants $L_1$ and $L_2$ such that

$$\parallel J\left(y\right)-J\left(x\right)\parallel \le L_1\parallel y-x\parallel ,\forall x,y\in {\mathbb{R}}^n,$$

(15)

and

$$\parallel F\left(y\right)-F\left(x\right)\parallel \le L_2\parallel y-x\parallel ,\forall x,y\in {\mathbb{R}}^n.$$

(16)

Lemma 1.

Let $d_k$ be computed by (12), then the inequality

$$Pre{d}_k\ge \parallel J_k^{\mathrm{T}}{F}_k\parallel min\left\{\parallel d_k\parallel ,\frac{\parallel J_k^{\mathrm{T}}{F}_k\parallel }{\parallel J_k^{\mathrm{T}}{J}_k\parallel }\right\}$$

(17)

holds for all $k\ge 0$.

Proof. 

From (7), for $\alpha \in [0,1]$, we have

$$\begin{array}{ccc}\hfill \parallel Pre{d}_k\parallel & =& \parallel F_k{\parallel }^2-{\parallel F_k+J_k{d}_k\parallel }^2\hfill \\ & \ge & \parallel F_k{\parallel }^2-{∥F_k-J_k\frac{\alpha \parallel d_k\parallel }{\parallel J_k^{\mathrm{T}}{F}_k\parallel }{J}_k^{\mathrm{T}}{F}_k∥}^2\hfill \\ & \ge & 2\alpha \parallel d_k\parallel \parallel J_k^{\mathrm{T}}{F}_k\parallel -{\alpha }^2\parallel d_k{\parallel }^2\parallel J_k^{\mathrm{T}}{F}_k\parallel ,\hfill \end{array}$$

(18)

then

$$\begin{array}{ccc}\hfill \parallel Pre{d}_k\parallel & \ge & \underset{0\le \alpha \le 1}{max}\left(2\alpha \parallel d_k\parallel \parallel J_k^{\mathrm{T}}{F}_k\parallel -{\alpha }^2\parallel d_k{\parallel }^2\parallel J_k^{\mathrm{T}}{F}_k\parallel \right)\hfill \\ & \ge & \parallel J_k^{\mathrm{T}}{F}_k\parallel min\left\{\parallel d_k\parallel ,\frac{\parallel J_k^{\mathrm{T}}{F}_k\parallel }{\parallel J_k^{\mathrm{T}}{J}_k\parallel }\right\}.\hfill \end{array}$$

(19)

The proof is complete. □

Theorem 1.

Under the conditions of Assumption 1, the sequence $\left\{x_k\right\}$ generated by NAALM algorithm satisfies

$$\underset{k\to \infty }{lim}\parallel J_k^{\mathrm{T}}{F}_k\parallel =0.$$

(20)

Proof. 

If the theorem is not true, then there exist a positive $\tau$ and infinitely many k such that

$$\parallel J_k^{\mathrm{T}}{F}_k\parallel \ge \tau .$$

(21)

Let $T_1$, $T_2$ be the sets of all indices that satisfy

$$\begin{array}{c}\hfill T_1=\left\{k\right|\parallel J_k^{\mathrm{T}}{F}_k\parallel \ge \tau \}\end{array}$$

and

$$\begin{array}{c}\hfill T_2=\left\{k\right|\parallel J_k^{\mathrm{T}}{F}_k\parallel \ge \frac{\tau }{2}\mathrm{and}{x}_{k+1}\ne x_k\}.\end{array}$$

Then, $T_1$ is an infinite set. In the following, we will derive the contradictions regarding whether $T_2$ is finite or infinite.

Case (I) $T_2$ is finite.

It follows from the definition of $T_2$ that the set

$$\begin{array}{c}\hfill T_3=\left\{k\right|\parallel J_k^{\mathrm{T}}{F}_k\parallel \ge \tau \mathrm{and}{x}_{k+1}\ne x_k\}\end{array}$$

is also finite. Let $\tilde{k}$ be the largest index of $T_3$. Then, we know that $x_{k+1}=x_k$ holds for all $k\in \{k>\tilde{k}|k\in T_1\}$. Define the indices set

$$\begin{array}{c}\hfill T_4=\{k>\tilde{k}|\parallel J_k^{\mathrm{T}}{F}_k\parallel \ge \tau \mathrm{and}{x}_{k+1}=x_k\}.\end{array}$$

Suppose $k\in T_4$. It is easy to see that $\parallel J_{k+1}^{\mathrm{T}}{F}_{k+1}\parallel \ge \tau$. Moreover, we have $x_{k+2}=x_{k+1}$. Otherwise, if $x_{k+2}\ne x_{k+1}$, then $k+1\in T_3$, which contradicts the fact that $\tilde{k}$ is the largest index of $T_3$. Hence, we have $k+1\in T_4$. By induction, we know that $\parallel J_k^{\mathrm{T}}{F}_k\parallel \ge \tau$ and $x_{k+1}=x_k$ hold for all $k>\tilde{k}$.

It now follows from Step 3 of the NAALM Algorithm that $r_k<p_0$ for all $k>\tilde{k}$, which imply

$${\mu }_k\to +\infty \mathrm{and}{\lambda }_k\to +\infty ,$$

(22)

due to (12)–(14) and $x_{k+1}=x_k$ for all $k>\tilde{k}$. Hence, we have

$$\underset{k\to \infty }{lim}{d}_k=0.$$

(23)

Furthermore, it follows from (21), (23) and Lemma 1 that

$$\begin{array}{ccc}\hfill |r_k-1|& =& \left|\frac{Are{d}_k}{Pre{d}_k}-1\right|\hfill \\ & =& \left|\frac{\parallel F_k+J_k{d}_k{\parallel }^2-{\parallel F(x_k+d_k)\parallel }^2}{Pre{d}_k}\right|\hfill \\ & =& \frac{\parallel F_k+J_k{d}_k\parallel O(\parallel d_k{\parallel }^2)+O(\parallel d_k{\parallel }^4)}{Pre{d}_k}\hfill \\ & \le & \frac{\parallel F_k+J_k{d}_k\parallel O(\parallel d_k{\parallel }^2)+O(\parallel d_k{\parallel }^4)}{\parallel J_k^{\mathrm{T}}{F}_k\parallel min\left\{\parallel d_k\parallel ,\frac{\parallel J_k^{\mathrm{T}}{F}_k\parallel }{\parallel J_k^{\mathrm{T}}{J}_k\parallel }\right\}}\hfill \\ & \le & \frac{O(\parallel d_k{\parallel }^2)}{\parallel d_k\parallel }\to 0,\hfill \end{array}$$

(24)

that is, $r_k\to 1$. In view of the updating rule of ${\mu }_k$, we know that there exists a positive constant $\tilde{m}>m$ such that ${\mu }_k<\tilde{m}$ holds for all sufficiently large k, which is a contradiction to (22). Hence, the supposition (21) cannot be true while $T_2$ is finite.

Case (II) $T_2$ is infinite.

It follows from Lemma 1 that

$$\begin{array}{ccc}\hfill \parallel F_1\parallel & \ge & \sum _{k\in T_2}(\parallel F_k{\parallel }^2-\parallel F_{k+1}{\parallel }^2)\hfill \\ & \ge & \sum _{k\in T_2}{p}_{0}Pre{d}_k\hfill \\ & \ge & \sum _{k\in T_2}{p}_0\parallel J_k^{\mathrm{T}}{F}_k\parallel min\left\{\parallel d_k\parallel ,\frac{\parallel J_k^{\mathrm{T}}{F}_k\parallel }{\parallel J_k^{\mathrm{T}}{J}_k\parallel }\right\}\hfill \\ & \ge & \sum _{k\in T_2}\frac{p_0\tau }{2},\hfill \end{array}$$

(25)

which gives

$$\sum _{k\in T_2}\parallel d_k\parallel <+\infty .$$

(26)

The above inequality, together with the Lipschitz conditions (15) and (16), implies that

$$\begin{array}{c}\hfill \sum _{k\in T_2}\left|\parallel J_k^{\mathrm{T}}{F}_k\parallel -\parallel J_{k+1}^{\mathrm{T}}{F}_{k+1}\parallel \right|<+\infty .\end{array}$$

(27)

Relation (27) and the fact that (21) holds for infinitely many k indicate that there exists a $\widehat{k}$ with $\parallel J_{\widehat{k}}^{\mathrm{T}}{F}_{\widehat{k}}\parallel \ge \tau$ such that

$$\begin{array}{c}\hfill \sum _{k\in T_2,k\ge \widehat{k}}\left|\parallel J_k^{\mathrm{T}}{F}_k\parallel -\parallel J_{k+1}^{\mathrm{T}}{F}_{k+1}\parallel \right|<\frac{\tau }{2}.\end{array}$$

By induction, we obtain that $\parallel J_k^{\mathrm{T}}{F}_k\parallel \ge \frac{\tau }{2}$ for all $k\ge \widehat{k}$. This result and (26) mean that

$$\underset{k\to \infty }{lim}\parallel d_k\parallel =0.$$

(28)

It follows from (12) and (13) that ${\mu }_k\to +\infty$. By the same analysis as (24), we know that ${\mu }_k\to 1$. Hence, there exists a positive constant $\overline{m}>m$ such that ${\mu }_k<\overline{m}$ holds for all large k, which introduces a contradiction. Therefore, the supposition (21) cannot be true when $T_2$ is infinite. The proof is complete. □

2.3. Local Convergence

In this section, we will study the local convergence properties of the NAALM algorithm by using the singular value decomposition (SVD) technique. We assume that the sequence $\left\{x_k\right\}$ generated by the NAALM algorithm converges to the nonempty solution set $X^{*}$ and lies in some neighborhood of $x^{*}\in X^{*}$. Firstly, we present some assumptions which the local convergence theory required.

Definition 1.

Let $\mathbb{N}\subset {\mathbb{R}}^n$ such that $\mathbb{N}\cap X^{*}\ne \varphi$, we say that $\parallel F\left(x\right)\parallel$ provides a local error bound on $\mathbb{N}$ for (1) if there exists a positive constant $c>0$ such that

$$c\mathrm{dist}(x,X^{*})\le \parallel F_k\parallel ,x\in \mathbb{N},$$

(29)

where $\mathrm{dist}(x,X^{*})$ is the distance from x to $X^{*}$.

Assumption 2.

(i) $F\left(x\right)$ is continuously differentiable, and $J\left(x\right)$ is Lipschitz continuous on $N(x^{*},b_1)$ with $b_1<1$, i.e., there exists a positive constant $L_1$ such that

$$\parallel J\left(y\right)-J\left(x\right)\parallel \le L_1\parallel y-x\parallel ,\forall x,y\in N(x^{*},b_1)=\left\{x\right|\parallel x-x^{*}\parallel \le b_1\}.$$

(30)

(ii) $F\left(x\right)$ provides a local error bound on some neighborhood of $x^{*}\in X^{*}$, i.e., there exists a positive constant $c_1>0$ such that

$$\parallel F\left(x\right)\parallel \ge c_1\mathrm{dist}(x,X^{*}),\forall x\in N(x^{*},b_1).$$

(31)

By the Lipschitzness of the Jacobian matrix proposed by (30), we have

$$\begin{array}{ccc}\hfill \parallel F\left(y\right)-F\left(x\right)-J\left(x\right)(y-x)\parallel & =& ∥{\int }_0^{1}J(x+t(y-x))(y-x)dt-J\left(x\right)(y-x)∥\hfill \\ & \le & \parallel y-x\parallel {\int }_0^1∥J(x+t(y-x\left)\right)-J\left(x\right)∥dt\hfill \\ & \le & L_1{\parallel y-x\parallel }^2,\hfill \end{array}$$

(32)

and

$$\parallel F\left(y\right)-F\left(x\right)\parallel \le L_2\parallel y-x\parallel ,\forall x,y\in N(x^{*},b_1),$$

(33)

where $L_2$ is a positive constant.

In the following, we use ${\overline{x}}_k$ to denote the vector in $X^{*}$ that satisfies

$$\parallel {\overline{x}}_k-x_k\parallel =\mathrm{dist}(x_k,X^{*}),\forall x,y\in N(x^{*},b_1).$$

(34)

To obtain the local convergence rate of $x_k$, we present some lemmas.

Lemma 2.

Under the conditions of Assumption 2, for all sufficiently large k, there exists a constant $c_2>0$ such that

$$\parallel d_k\parallel \le c_2\parallel {\overline{x}}_k-x_k\parallel .$$

(35)

Proof. 

According to (34), we have

$$\parallel {\overline{x}}_k-x^{*}\parallel \le \parallel {\overline{x}}_k-x_k\parallel +\parallel x_k-x^{*}\parallel \le 2\parallel x_k-x^{*}\parallel \le b_1,$$

(36)

which means that $\overline{x}\in N(x^{*},b_1)$. Following from (13),

$$\begin{array}{ccc}\hfill {\lambda }_k=\frac{{\mu }_k\parallel F_k\parallel }{1+\parallel F_k\parallel }& =& {\mu }_k\left(1-\frac{1}{1+\parallel F_k\parallel }\right)\hfill \\ & \ge & m\left(1-\frac{1}{1+c_1\parallel {\overline{x}}_k-x_k\parallel }\right)\hfill \\ & =& \frac{m{c}_1\parallel {\overline{x}}_k-x_k\parallel }{1+c_1\parallel {\overline{x}}_k-x_k\parallel },\hfill \end{array}$$

(37)

and we have from (32) that

$$\parallel F_k+J_k({\overline{x}}_k-x_k){\parallel }^2=\parallel F\left({\overline{x}}_k\right)-F_k-J_k({\overline{x}}_k-x_k){\parallel }^2\le L_1^2{\parallel {\overline{x}}_k-x_k\parallel }^4.$$

(38)

As $d_k$ is a minimizer of ${\psi }_k\left(d\right)$, we have

$$\begin{array}{ccc}\hfill \parallel d_k{\parallel }^2& \le & \frac{1}{{\lambda }_k}{\phi }_k\left(d_k\right)\hfill \\ & \le & \frac{1}{{\lambda }_k}{\phi }_k({\overline{x}}_k-x_k)\hfill \\ & =& \frac{1}{{\lambda }_k}(\parallel F_k+J_k({\overline{x}}_k-x_k){\parallel }^2+{\lambda }_k\parallel {\overline{x}}_k-x_k{\parallel }^2)\hfill \\ & \le & \frac{1+c_1\parallel {\overline{x}}_k-x_k\parallel }{m{c}_1\parallel {\overline{x}}_k-x_k\parallel }(L_1^2\parallel {\overline{x}}_k-x_k{\parallel }^4)+\parallel {\overline{x}}_k-x_k{\parallel }^2\hfill \\ & =& O(\parallel {\overline{x}}_k-x_k{\parallel }^2),\hfill \end{array}$$

then there exists a constant $c_2>0$ such that $\parallel d_k\parallel \le c_2\parallel {\overline{x}}_k-x_k\parallel$. The proof is completed. □

Lemma 3.

Under the conditions of Assumption 2, for all sufficiently large k, there exists a positive constant $M>m$ such that

$${\mu }_k\le M.$$

(39)

Proof. 

First, we show that for sufficiently large k, the following inequality holds

$$Pre{d}_k=\parallel F_k{\parallel }^2-\parallel F_k+J_k{d}_k{\parallel }^2\ge min\left\{\frac{c_1}{2c_2},\frac{c_1}{2}\right\}\parallel F_k\parallel \parallel d_k\parallel .$$

(40)

We consider two cases. In one case, if $\parallel {\overline{x}}_k-x_k\parallel \le \parallel d_k\parallel$, then the definition of $d_k$ and Assumption 2 imply that

$$\begin{array}{ccc}\hfill \parallel F_k\parallel -\parallel F_k+J_k{d}_k\parallel & \ge & \parallel F_k\parallel -\parallel F_k+J_k({\overline{x}}_k-x_k)\parallel \hfill \\ & \ge & c_1\parallel {\overline{x}}_k-x_k\parallel -L_1{\parallel {\overline{x}}_k-x_k\parallel }^2\hfill \\ & \ge & \frac{c_1}{2c_2}\parallel d_k\parallel .\hfill \end{array}$$

(41)

In the other case, if $\parallel {\overline{x}}_k-x_k\parallel >\parallel d_k\parallel$, then we have

$$\begin{array}{ccc}\hfill \parallel F_k\parallel -\parallel F_k+J_k{d}_k\parallel & \ge & \parallel F_k\parallel -∥F_k+\frac{\parallel d_k\parallel }{\parallel {\overline{x}}_k-x_k\parallel }{J}_k({\overline{x}}_k-x_k)∥\hfill \\ & \ge & \frac{\parallel d_k\parallel }{\parallel {\overline{x}}_k-x_k\parallel }\left(\parallel F_k\parallel -\parallel F_k+J_k({\overline{x}}_k-x_k)\parallel \right)\hfill \\ & \ge & \frac{\parallel d_k\parallel }{\parallel {\overline{x}}_k-x_k\parallel }\left(c_1\parallel {\overline{x}}_k-x_k\parallel -L_1{\parallel {\overline{x}}_k-x_k\parallel }^2\right)\hfill \\ & \ge & \frac{c_1}{2}\parallel d_k\parallel .\hfill \end{array}$$

(42)

Inequalities (41) and (42), together with Lemma 2 show that

$$\begin{array}{ccc}\hfill Pre{d}_k& =& (\parallel F_k\parallel +\parallel F_k+J_k{d}_k\parallel \left)\right(\parallel F_k\parallel -\parallel F_k+J_k{d}_k\parallel )\hfill \\ & \ge & \parallel F_k\parallel (\parallel F_k\parallel -\parallel F_k+J_k{d}_k\parallel )\hfill \\ & \ge & min\left\{\frac{c_1}{2c_2},\frac{c_1}{2}\right\}\parallel F_k\parallel \parallel d_k\parallel ,\hfill \end{array}$$

(43)

which gives (40). Hence, it follows from (40), Assumption 2 and Lemma 2 that

$$\begin{array}{ccc}\hfill |r_k-1|& =& \left|\frac{Are{d}_k}{Pre{d}_k}-1\right|\hfill \\ & =& \frac{\parallel F_k+J_k{d}_k\parallel O(\parallel d_k{\parallel }^2)+O(\parallel d_k{\parallel }^4)}{Pre{d}_k}\hfill \\ & \le & \frac{\parallel F_k\parallel O(\parallel d_k{\parallel }^2)+O(\parallel d_k{\parallel }^4)}{O(\parallel F_k\parallel \parallel d_k\parallel )}\hfill \\ & =& O(\parallel d_k\parallel )\to 0.\hfill \end{array}$$

Therefore, we have $r_k\to 1$, thus, there exists a constant $M>m$ such that ${\mu }_k\le M$ for all large k. The proof is completed. □

Without generality, we assume rank $\left(J\left(x^{*}\right)\right)=r$ for all $\overline{x}\in N(x^{*},b_1)\cap X^{*}$. Suppose the SVD of $J\left(\overline{x}\right)$ is

$$J\left(\overline{x}\right)=[{\overline{U}}_1,{\overline{U}}_2]\left[\begin{array}{cc}{\overline{\mathsf{\Sigma }}}_1& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\overline{V}}_1^{\mathrm{T}}\\ {\overline{V}}_2^{\mathrm{T}}\end{array}\right]={\overline{U}}_1{\overline{\mathsf{\Sigma }}}_1{\overline{V}}_1^{\mathrm{T}},$$

(44)

where ${\overline{\mathsf{\Sigma }}}_1=\mathrm{diag}({\overline{\sigma }}_1,{\overline{\sigma }}_2,\dots ,{\overline{\sigma }}_r)$ with ${\overline{\sigma }}_1\ge {\overline{\sigma }}_2\ge \cdots \ge {\overline{\sigma }}_r>0$ and $\overline{U}=[{\overline{U}}_1,{\overline{U}}_2]$, $\overline{V}=[{\overline{V}}_1,{\overline{V}}_2]$ are orthogonal matrices. Correspondingly, we consider SVD of $J\left(x_k\right)$ by

$$J\left(x_k\right)=[U_1,U_2,U_3]\left[\begin{array}{ccc}{\mathsf{\Sigma }}_1& 0& 0\\ 0& {\mathsf{\Sigma }}_2& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{V}_1^{\mathrm{T}}\\ V_2^{\mathrm{T}}\\ V_3^{\mathrm{T}}\end{array}\right]=U_1{\mathsf{\Sigma }}_1{V}_1^{\mathrm{T}}+U_2{\mathsf{\Sigma }}_2{V}_2^{\mathrm{T}},$$

(45)

where $U=[U_1,U_2,U_3]$, $V=[V_1,V_2,V_3]$ are orthogonal matrixes, ${\mathsf{\Sigma }}_1=\mathrm{diag}({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_r)$ with ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_r>0$ and ${\mathsf{\Sigma }}_2=\mathrm{diag}({\sigma }_{r+1},{\sigma }_{r+2},\cdots ,{\sigma }_{r+q})$ with ${\sigma }_{r+1}\ge {\sigma }_{r+2}\ge \cdots \ge {\sigma }_{r+q}>0$.

Lemma 4.

Under the conditions of Assumption 2, for all sufficiently large k, we have

$\left(a\right)\parallel U_1{U}_1^{\mathrm{T}}{F}_k\parallel \le O(\parallel {\overline{x}}_k-x_k\parallel )$;

$\left(b\right)\parallel U_2{U}_2^{\mathrm{T}}{F}_k\parallel \le O(\parallel {\overline{x}}_k-x_k{\parallel }^2)$;

$\left(c\right)\parallel U_3{U}_3^{\mathrm{T}}{F}_k\parallel \le O(\parallel {\overline{x}}_k-x_k{\parallel }^2)$;

$\left(d\right)\parallel F_k+J_k{d}_k\parallel \le O(\parallel {\overline{x}}_k-x_k{\parallel }^2)$.

Proof. 

The result $\left(a\right)$ follows immediately from (16). By (15) and the theory of matrix perturbation [26], we have

$$\parallel \mathrm{diag}({\mathsf{\Sigma }}_1-{\overline{\mathsf{\Sigma }}}_1,{\mathsf{\Sigma }}_2,0)\parallel \le \parallel J_k-J\left({\overline{x}}_k\right)\parallel \le L_1\parallel {\overline{x}}_k-x_k\parallel ,$$

which implies that

$$\parallel {\mathsf{\Sigma }}_1-{\overline{\mathsf{\Sigma }}}_1\parallel \le L_1\parallel {\overline{x}}_k-x_k\parallel \mathrm{and}\parallel {\mathsf{\Sigma }}_2\parallel \le L_1\parallel {\overline{x}}_k-x_k\parallel .$$

(46)

Let $s_k=-J_k^{+}{F}_k$, where $J_k^{+}$ is the pseudo-inverse of $J_k$. It is easy to see that $s_k$ is the least-squares solution of $min\parallel F_k+J_{k}s\parallel$, so we obtain from (32) that

$$\parallel U_3{U}_3^{\mathrm{T}}{F}_k\parallel =\parallel F_k+J_k{s}_k\parallel \le \parallel F_k+J_k({\overline{x}}_k-x_k)\parallel \le O(\parallel {\overline{x}}_k-x_k{\parallel }^2).$$

Let ${\overline{J}}_k=U_1{\mathsf{\Sigma }}_1{V}_1^{\mathrm{T}}$ and ${\overline{s}}_k=-{\overline{J}}_k^{+}{F}_k$. Since ${\overline{s}}_k$ is the least-squares solution of $min\parallel F_k+{\overline{J}}_{k}s\parallel$, it follows from (32) that

$$\begin{array}{ccc}\hfill \parallel (U_2{U}_2^{\mathrm{T}}+U_3{U}_3^{\mathrm{T}})F_k\parallel & =& \parallel F_k+{\overline{J}}_k{\overline{s}}_k\parallel \hfill \\ & \le & \parallel F_k+{\overline{J}}_k({\overline{x}}_k-x_k)\parallel \hfill \\ & \le & \parallel F_k+J_k({\overline{x}}_k-x_k)\parallel +\parallel ({\overline{J}}_k-J_k)({\overline{x}}_k-x_k)\parallel \hfill \\ & \le & L_1\parallel {\overline{x}}_k-x_k{\parallel }^2+\parallel U_2{\mathsf{\Sigma }}_2{V}_2^{\mathrm{T}}({\overline{x}}_k-x_k)\parallel \hfill \\ & \le & L_1\parallel {\overline{x}}_k-x_k{\parallel }^2+L_1\parallel {\overline{x}}_k-x_k\parallel \parallel {\overline{x}}_k-x_k\parallel \hfill \\ & \le & O(\parallel {\overline{x}}_k-x_k{\parallel }^2).\hfill \end{array}$$

Due to the orthogonality of $U_2$ and $U_3$, we obtain the result $\left(b\right)$.

Using (12) and (45), we obtain

$$d_k=-V_1{({\mathsf{\Sigma }}_1^2+{\lambda }_{k}I)}^{-1}{\mathsf{\Sigma }}_1{U}_1^{\mathrm{T}}{F}_k-V_2{({\mathsf{\Sigma }}_2^2+{\lambda }_{k}I)}^{-1}{\mathsf{\Sigma }}_2{U}_2^{\mathrm{T}}{F}_k,$$

and

$$\begin{array}{ccc}\hfill F_k+J_k{d}_k& =& F_k-U_1{\mathsf{\Sigma }}_1{({\mathsf{\Sigma }}_1^2+{\lambda }_{k}I)}^{-1}{\mathsf{\Sigma }}_1{U}_1^{\mathrm{T}}{F}_k-U_2{\mathsf{\Sigma }}_2{({\mathsf{\Sigma }}_2^2+{\lambda }_{k}I)}^{-1}{\mathsf{\Sigma }}_2{U}_2^{\mathrm{T}}{F}_k\hfill \\ & =& {\lambda }_k{U}_1{({\mathsf{\Sigma }}_1^2+{\lambda }_{k}I)}^{-1}{U}_1^{\mathrm{T}}{F}_k+{\lambda }_k{U}_2{({\mathsf{\Sigma }}_2^2+{\lambda }_{k}I)}^{-1}{U}_2^{\mathrm{T}}{F}_k+U_3{U}_3^{\mathrm{T}}{F}_k.\hfill \end{array}$$

Following from (13) and (33), the LM parameter satisfies

$${\lambda }_k=\frac{{\mu }_k\parallel F_k\parallel }{1+\parallel F_k\parallel }\le {\mu }_k\parallel F_k\parallel \le M{L}_2\parallel {\overline{x}}_k-x_k\parallel .$$

Since $\left\{x_k\right\}$ converges to the solution set $X^{*}$, we assume that $L_1\parallel {\overline{x}}_k-x_k\parallel \le \frac{{\overline{\sigma }}_r}{2}$ holds for all sufficiently large k. Then, it follows from (46) that

$$\parallel {\mathsf{\Sigma }}_1^{-1}\parallel \le \frac{1}{{\overline{\sigma }}_r-L_1\parallel {\overline{x}}_k-x_k\parallel }\le \frac{2}{{\overline{\sigma }}_r}.$$

It then follows from Lemmas 3 and 4 that

$$\begin{array}{ccc}\hfill \parallel F_k+J_k{d}_k\parallel & \le & {\lambda }_k\parallel {\mathsf{\Sigma }}_1^{-2}\parallel \parallel U_1^{\mathrm{T}}{F}_k\parallel +\parallel U_2^{\mathrm{T}}{F}_k\parallel +\parallel U_3{U}_3^{\mathrm{T}}{F}_k\parallel \hfill \\ & \le & \frac{4L_{2}M{\parallel {\overline{x}}_k-x_k\parallel }^2}{{\overline{\sigma }}_r^2}+O(\parallel {\overline{x}}_k-x_k{\parallel }^2)+O(\parallel {\overline{x}}_k-x_k{\parallel }^2)\hfill \\ & =& O(\parallel {\overline{x}}_k-x_k{\parallel }^2).\hfill \end{array}$$

(47)

The proof is completed. □

We can state the quadratic convergence result of the NAALM algorithm.

Theorem 2.

Let the sequence $\left\{x_k\right\}$ be generated by the NAALM algorithm, under Assumption 2, the sequence $\left\{x_k\right\}$ converges quadratically to a solution of nonlinear Equation (1).

Proof. 

It follows from Assumption 2, Lemma 2 and (47) that

$$\begin{array}{ccc}\hfill c_1\parallel {\overline{x}}_{k+1}-x_{k+1}\parallel & \le & \parallel F\left(x_{k+1}\right)\parallel \hfill \\ & =& \parallel F(x_k+d_k)\parallel \hfill \\ & \le & \parallel F_k+J_k{d}_k\parallel +O(\parallel d_k{\parallel }^2)\hfill \\ & =& O(\parallel {\overline{x}}_k-x_k{\parallel }^2).\hfill \end{array}$$

(48)

On the other hand, it is clear that

$$\parallel {\overline{x}}_k-x_k\parallel =\mathrm{dist}(x_k,X^{*})\le \parallel {\overline{x}}_{k+1}-x_k\parallel \le \parallel {\overline{x}}_{k+1}-x_{k+1}\parallel +\parallel d_k\parallel .$$

It follows from Lemma 2 that, for any sufficiently large k, we have

$$\parallel {\overline{x}}_k-x_k\parallel \le 2\parallel d_k\parallel \le O(\parallel {\overline{x}}_k-x_k\parallel ).$$

Therefore, $\parallel d_k\parallel =O(\parallel {\overline{x}}_k-x_k\parallel )$. This, along with (48), indicates that

$$\parallel d_{k+1}\parallel \le O(\parallel d_k{\parallel }^2),$$

which implies that $\left\{x_k\right\}$ is quadratically convergent to a solution of set $X^{*}$. The proof is completed. □

3. Numerical Results

In this section, the numerical performance of NAALM algorithm will be listed. All codes were written in MATLAB R2016b on a PC with 1.19 GHz, 8.00 GB RAM, using Windows 11 operation system. In this section, we will expand on the following two aspects. On the one hand, the effectiveness of the NAALM algorithm is illustrated by comparing it with other algorithms on some test questions. On the other hand, it shows that the NAALM algorithm has good development prospects by applying the algorithm to a fresh agricultural products supply chain problem.

3.1. Some Singular Nonlinear Equations Problems

The test problems are constructed by modifying the nonsingular problems given by Moré et al. [27], which have the following form as [28]:

$$\widehat{F}\left(x\right)=F\left(x\right)-J\left(x^{*}\right)A{\left(A^{\mathrm{T}}A\right)}^{-1}{A}^{\mathrm{T}}(x-x^{*}),$$

where $F\left(x\right)$ is the standard test function, $A\in {\mathbb{R}}^{n\times k}$ has full column rank with $0\le k\le n$ and x is a solution of the equation $F\left(x\right)=0$. According to the definition of $\widehat{F}\left(x\right)$, we obtain

$$\widehat{J}\left(x^{*}\right)=J\left(x^{*}\right)(I-A{\left(A^{\mathrm{T}}A\right)}^{-1}{A}^{\mathrm{T}}),$$

where $\widehat{J}\left(x^{*}\right)$ is Jacobian matrix of $F\left(x\right)$ at $x^{*}$ with rank $n-k$ and $\widehat{F}\left(x^{*}\right)=0$. In our test problems, some of $\widehat{J}\left(x^{*}\right)$ are symmetric matrices and some are non-symmetric matrices. Note that some roots of $\widehat{F}\left(x\right)$ may not be roots of $F\left(x\right)$. Similar to [28], we construct two sets of singular problems while $\widehat{J}\left(x^{*}\right)$ have rank $n-1$ or $n-2$, by choosing

$$A={[1,1,...,1]}^{\mathrm{T}}\in {\mathbb{R}}^{n\times 1},$$

and

$$A={\left[\begin{array}{cccccc}1& 1& 1& 1& ...& 1\\ 1& -1& 1& -1& ...& \pm 1\end{array}\right]}^{\mathrm{T}}\in {\mathbb{R}}^{n\times 2}.$$

We test our NAALM algorithm on some singular nonlinear equations, and compare it with the self-adaptive Levenberg–Marquardt algorithm (SLM) proposed in [18]. The main differences between these two algorithms are in the updating rule of ${\mu }_k$.

We set $p_0={10}^{-4}$, $p_1=\frac{1}{4}$, $p_2=\frac{3}{4}$, ${\beta }_1=\frac{5}{4}$, ${\beta }_2=\frac{1}{3}$, ${\beta }_3=\frac{6}{5}$, $m={10}^{-8}$, ${\mu }_0={10}^{-2}$, for all the tests. All test methods are terminated when $\parallel J_k^{\mathrm{T}}{F}_k\parallel \le {10}^{-5}$. The algorithm is considered to fail when the number of iterations exceeds 500. Considering the global convergence of the algorithms, we run each test problem for five starting points, $-10x_0$, $-x_0$, $x_0$, $10x_0$ and $100x_0$, where $x_0$ is given by [28]. For n as a variable, we take $n=500$, $n=1000$, respectively.

The performance profile of two algorithms. including the number of iterations (NI), function evaluations (NF), gradient evaluations (NG) and CPU time (CPU), is analyzed using the profiles of Dolan and Mor$\acute{\mathrm{e}}$ [29]. Let Y and W be the set of methods and test problems, $n_y$, $n_w$ be the number of methods and test problems, respectively. The performance profile $\psi :\mathbb{R}\to [0,1]$ is for each $y\in Y$ and $w\in W$ defined that $a_{w,y}>0$ is NI or NF or NG or CPU required to solve problems w by method y. Furthermore, the performance profile is obtained by

$$\begin{array}{c}\hfill {\psi }_y\left(\tau \right)=\frac{1}{n_w}size\{w\in W:lo{g}_2{r}_{w,y}\le \tau \},\end{array}$$

where $\tau >0$, $size\{·\}$ is the number of the elements in a set, and $r_{w,y}$ is the performance ratio defined as

$$\begin{array}{c}\hfill r_{w,y}=\frac{a_{w,y}}{min\{a_{w,y}:w\in Wandy\in Y\}}.\end{array}$$

Generally, the method whose performance profile plot is on the top right will represent the best method.

As can be seen from Figure 1, the NAALM algorithm is better than the SLM algorithm in terms of the number of iterations, especially when $\tau >2$, the curve of NAALM algorithm becomes stable, which indicates that NAALM algorithm can solve the problem only with fewer iterations. In terms of function evaluations, as shown in Figure 2, the NAALM algorithm curve in $\tau >1.75$, it has reached a stable state, while SLM algorithm can reach a stable state only when the curve coincides with that of NAALM algorithm at $\tau >2.75$; Figure 3 shows the performance diagram of the SLM algorithm and the NAALM algorithm in the Jacobian matrix. It can be seen that the NAALM algorithm can successfully solve test problems up to 98%, while SLM can only reach 94%, which shows that the NAALM algorithm can reduce the calculation times of the Jacobian matrix and save the calculation amount. Figure 4 shows the CPU time performance of the NAALM algorithm and the SLM algorithm. It can be seen from the figure that when $\tau <4.5$, the curves of the NAALM algorithm and the SLM algorithm are similar, but when $\tau >4.5$, both the NAALM algorithm and the SLM algorithm tend to be stable and coincide. Therefore, Figure 1, Figure 2, Figure 3 and Figure 4 show that the accelerated version of the LM algorithm proposed in this paper can not only converge to the solution quickly, but also reduce the computation amount of the Jacobian matrix.

3.2. Supply Chain Optimization Problems

The security and stability of the supply chain has a great impact on promoting high-quality and sustainable development of the economy. Therefore, supply chain has been applied to many fields, such as low-carbon supply chain, manufacturing green supply chain, food trade supply chain. In recent years, with the improvement of living standards, the quality of fresh agricultural products has attracted widespread attention from consumers. In order to meet the demand of consumers for high quality and low price of fresh agricultural products, we use the NAALM algorithm to study how suppliers and retailers make decisions to maximize both their own profits and the total profit of the fresh agricultural products supply chain under the decentralized policy.

In this supply chain, as the leader of Stackelberg game, fresh agricultural product suppliers supply the same variety of ordinary fresh agricultural products (ofp) and green fresh agricultural products (gfp) to retailers as followers, while retailers sell them to consumers. Suppliers need to choose the optimal wholesale price strategy of two fresh agricultural products, and retailers need to choose the optimal retail price strategy of two fresh agricultural products and determine the order quantity of two fresh agricultural products by market demand.

Without considering the impact of emergencies, the market demand for fresh agricultural products is relatively stable, and it is only related to price and freshness. Due to the substitution of the same varieties of ofp and gfp, there is a competitive relationship in the demand market. Based on the demand function theory of alternative price competition, it is assumed that the demand function of two fresh agricultural products is as follows

$$q_i=a-b\frac{p_i}{\theta }+r\frac{p_j}{\theta },i=1,2,j=3-i,$$

(49)

where $q_1$, $q_2$ represent the market demand of gfp and ofp, respectively, a represents the total potential market capacity of fresh agricultural products, $p_1$, $p_2$ represent the retail price of gfp and ofp, respectively, b is the price sensitivity coefficient, r is the competitive substitution coefficient of the two products, and it satisfies $b>r>0$, $\theta (0\le \theta \le 1)$ is the freshness of fresh produce when it arrives at the retailer’s store.

Under the decentralized policy, we regard suppliers and retailers as independent entities, and both with the goal of maximizing their respective interests. Now, the profit function of fresh agricultural products retailer is as follows

$$\underset{p_1,p_2}{max}{\pi }_s=(p_1-w_1)(a-b\frac{p_1}{\theta }+r\frac{p_2}{\theta })+(p_2-w_2)(a-b\frac{p_2}{\theta }+r\frac{p_1}{\theta }),$$

(50)

where $w_1$, $w_2$ represent the supply price of gfp and ofp, respectively, and the profit function of fresh agricultural products suppliers is as follows

$$\underset{w_1,w_2}{max}{\pi }_R=(w_1-\frac{c_1}{1-\beta })(a-b\frac{p_1}{\theta }+r\frac{p_2}{\theta }+(w_2-\frac{c_2}{1-\beta })(a-b\frac{p_2}{\theta }+r\frac{p_1}{\theta }),$$

(51)

where $\beta$$(0<\beta <1)$ represents the quantity loss of fresh produce when it reaches the retailer’s store, $c_1$, $c_2$ represents the unit production cost of gfp and ofp, respectively. Obviously, $p_1>p_1>0$ and $c_1>c_2>0$. We record the total profit of fresh agricultural products supply chain as follows:

$${\pi }_T={\pi }_s+{\pi }_R.$$

(52)

With reference to the setting of the parameters in the relevant literature [30], we set: $a=50$, $b=2$, $c_1=4$, $c_2=2$, $r=1.5$, $\beta =0.2$, $\theta =0.85$. These values satisfy the theoretical proof in [30] and can guarantee that the optimal value has practical significance. Now, we transform the unconstrained optimization problem (51) into a nonlinear equation problem, and then choose different initial points and use the NAALM algorithm to solve the nonlinear equation problem.

As can be seen from

Table 1. The optimal solution corresponding to different initial points by NAALM.

Initial Point	p1	p2	w1	w2	q1	q2	πT
(1;1;1;1)	45.0266	43.7205	65.065	64.329	10.5818	13.2479	1.4587 × 103
(10;10;10;10)	44.9648	43.7460	64.991	64.358	10.5163	13.2607	1.4587 × 103
(30;30;30;30)	44.9920	43.7496	64.956	64.376	10.5825	13.3037	1.4587 × 103
(50;50;50;50)	45.0158	43.7525	65.062	64.372	10.5125	13.3036	1.4587 × 103
(100;100;100;100)	44.9751	43.7553	65.955	64.269	10.5222	13.2679	1.4587 × 103

, with certain parameters, the NAALM algorithm can be used to solve the optimization problem, so as to obtain the optimal pricing strategy with maximum profit in the supply chain led by suppliers under the decentralized policy. In addition, the global convergence and robustness of the NAALM algorithm are verified according to different initial values and the number of iterations.

4. Conclusions

We constructed a new function that makes full use of the ratio information to update LM parameters adaptively. Based on this new LM parameter, we presented an adaptive accelerated Levenberg–Marquardt method for solving nonlinear equations. Furthermore, we showed the global convergence analysis of the proposed algorithm. Furthermore, the quadratic convergence is also obtained under the local error bound condition. Numerical experiments demonstrated that our method has good numerical performance. In addition, the application of the NAALM algorithm to a supply chain problem showed that the new algorithm has a good application prospect. We further highlight that the proposed NAALM algorithm can be used in other fields, such as the symmetric system of nonlinear equations. It is vital to note that the method’s convergence analysis in Hölderian local error bound condition will be taken into account in our future work.