Water resources management versus the world

1.   Introduction

In this paper, we consider the following bound constrained nonlinear systems of equations:

where $F(x) = (F_1(x), F_2(x), \cdots, F_n(x))^T$, and $F_i:\mathbb{R}^n\rightarrow \mathbb{R}$ is a nonlinear continuously differentiable function whose gradient is available. We denote by $F'(x) = (\nabla F_1(x), \nabla F_2(x), \cdots, \nabla F_n(x))^T$ the Jacobian matrix of $F$ at a given point $x.$ The set $\Omega\subseteq \mathbb{R}^n$ is defined as

for some given lower and upper bounds satisfying $-\infty\leq l_i < u_i\leq +\infty$ for all $i = 1, 2, \cdots, n.$

The bound constrained nonlinear equation of the type (1.1) is an important problem in the practical problems. There are a couple of different mathematical programming problems like Karush-Kuhn-Tucker systems and complementarity problems can be reformulated as the problem (1.1), see ^{[1,2,3,4,5,6,7,8,9]}. On the other hand, in many cases, the function $F_i(x)$ is not always defined on the whole space $\mathbb{R}^n,$ and one usually puts some suitable bounds on some or all of the variables.

The Newton type method is one of the most important numerical methods for problem (1.1) and many researchers are interested in this method ^{[7,10,11,12,13,14]}. Given a current iterate $x^k\in \Omega,$ the Newton method considers the least-squares solutions $d^k$ of the following nonlinear constrained equation:

We set the next iterate to be $x^{k+1} = x^{k}+d^k$ and call (1.2) the constrained Gauss-Newton method.

Another natural possibility is to consider solving the basic unconstrained Newton equation:

Denote the solution of (1.3) by $d^k_N$ if it exists and then define $x^{k+1}$ as the projection of $x^k+d^k_N$ onto $\Omega.$ This scheme can be called the projected Newton method.

On the other hand, there are many versions of the Newton type method, such as the constrained Levenberg-Marquardt method ^[6,15,16] usually used to solve the following subproblems:

where $\sigma$ is a positive constant.

Along with constrained versions of the methods in question, one can also consider their projected variants. The projected Levenberg-Marquardt method has been proposed in ^[15] and its iteration consists of finding the solution $d^k_{LM}$ of the unconstrained subproblem

and then defines the next iterate $x^{k+1}$ as the projection of $x^k+d^k_{LM}$ onto $\Omega.$

The Newton iteration can be costly, since partial derivatives must be computed and the linear system (1.3) must be solved at every iteration. This fact motivates the development of quasi-Newton methods ^[10,14,17] which are defined as the generalizations of (1.3) given by

In quasi-Newton methods, the matrices $B_k$ are intended to be approximations of $F'(x_k)$ and be updated by some quasi-Newton formulas. Another well known algorithm is the trust region type algorithm, for example, ^{[3,4,9,18,19,20,21]}.

Whether Newton method or quasi-Newton, one has to solve a linear system with full dimension, which will be expensive for large scale problems. To overcome this drawback, the active set methods are developed by many authors ^[7,12,13,22]. Since only a reduced dimension linear system to be dealt with at each iteration, the active set Newton methods are more efficient than the full Newton method especially for large scale problems.

To prove global convergence of the method outlined above, one often assumes that the iteration sequence is contained in a bounded set. If $l_i$ and $u_i$ are bounded and the algorithm generates a feasible sequence, the assumption holds naturally. Otherwise, one often makes an assumption that the level set is bounded. For unconstrained nonlinear equations system, M.Solodov designed a Newton method with projection technique, the method can generate a bounded iteration sequence without additional assumption and the global convergence is obtained. Motivated by the idea of M.Solodov ^[23], in this paper, we extend the method to constrained equations (1.1). By using this active set strategy, we only need to solve a linear system with reduced dimension at each iteration. The algorithm generates a bounded sequence automatically even if $l_i$ and $u_i$ are infinite. We obtain the global convergence and give numerical tests to show the efficiency of the proposed algorithm.

The paper is organized as follows: In section 2, we describe our algorithm in detail. In section 3, we prove the global convergence of the proposed algorithm. Some numerical tests are shown in Section 4 and a conclusion is given in section 5. Throughout this paper, we use $\|\cdot\|$ to denote the $2-$norm and $E$ denotes the identity matrix.

2.   Algorithm

We now describe our active set quasi-Newton method with projection technique in detail. To describe our algorithm, we introduce the definition of projection operator which is defined as a mapping form $\mathbb{R}^n$ to a nonempty closed convex subset $\Omega:$

A well-known property of the operator is that it is nonexpensive, namely,

Given a current iterate $x^k,$ let

where $\delta$ and $c$ are positive constants such that

and define the index sets

The precise statement of our algorithm is as follows:

Algorithm 2.1: (Active Set-type Quasi-Newton Method)

(S.0) Choose a positive definite matrix $B_k,$ $x^{0}\in[l, u],$ choose parameters $\beta\in(0, 1),$ $\lambda\in(0, 1),$ $\delta > 0,$ $c > 0,$ $\varepsilon > 0,$ $\mu_k > 0,$ and $\rho_k\in [0, 1),$ and set $k: = 0$.

(S.1) If $\|F(x^{k})\|\leq\varepsilon,$ stop.

(S.2) Try to compute a vector $d^{k}\in \mathbb{R}^n$ in the following way:

For $i\in\mathcal{A}_{k}$, set

For $i\in \mathcal{I}_{k}$, set solve the linear system

where

(S.3) Find $z^{k} = x^{k}+\alpha_{k}d^{k}$, where $\alpha_{k} = \beta^{m_{k}}$ with $m_{k}$ being the smallest nonnegative integer $m$ such that

(S.4) Compute

(S.5) Update $B_{k+1}$, set $k: = k+1$, go to (S.1).

Just as mentioned in ^[13], throughout this paper, we assume that the parameter $\delta > 0$ is chosen sufficiently small such that

This implies that we cannot have $x^k_i-l_i\leq\delta_k$ and $u_i-x^k_i\leq \delta_k$ for the same index $i\in \mathcal{A}_k.$

Our algorithm is somewhat different from the traditional active set Newton method as described in ^[13], where the search step $d^k$ in (S.2) is computed in the following formulas:

For $i\in\mathcal{A}_{k}$, set

For $i\in\mathcal{I}_{k}$, solve the linear system

As described in ^[13], in order to understand the formula for the computation of the components $d_i^k$ for $i \in \mathcal{I}_k$, note that, after a possible permutation of the rows and columns, ^[13] rewrite the standard (unconstrained) Newton equation $F'(x^k)d = -F(x^k)$ as

Here we replace (2.7) by (2.3) and (2.8) by (2.4), the main proposal is to guarantee that the inequality (2.5) holds. On the other hand, we compute $d_{\mathcal{I}_k}$ by (2.4) instead of (2.8) which can be seen as an inexact Newton method.

The matrix $B_k$ is updated by the well known rank two secant type formula updated by the well known BFGS formula

where $y_k = F(x^{k+1})-F(x^k)$ and $s_k = x^{k+1}-x^k.$

3.   Convergence

In the section, we prove the global convergence of Algorithm 2.1, we make the following assumption.

Assumption:

(A1) The function $F(x)$ is Lipschitz continuous and monotone, i.e., there exists a positive constant $L$ such that

and

(A2) The sequence of matrices $\{B_k\}$ is positive definite and bounded, i.e., there exists a positive constant $\kappa$ such that $\|B_k\|\leq \kappa$ for all $k.$

We first show that the algorithm is feasible, i.e., there exists a positive $m$ such that (2.5) holds.

Lemma 3.1. The Algorithm 2.1 is well defined.

Proof. We prove that the inequality (2.5) will hold with a nonnegative integer $m.$ Suppose that for some index $k$ this is not the case, which means, for all integer $m,$ we have

We further get

Now we take the limit of both sides of (3.4) as $m\rightarrow\infty$, when (3.4) holds which implies that $\lambda\geq 1,$ which contradicts the choice of $\lambda\in (0, 1).$ Hence we have that the inequality holds for some integer $m,$ and the whole algorithm is well defined.

In what follows, we assume that the algorithm generates an infinite iteration sequence. The following result shows that the algorithm generates a bounded sequence automatically and the proof is similar to Lemma 3.2 in ^[24] and we omit it here.

Theorem 3.2. Suppose assumptions (A1) and (A2) hold, sequences $\{x^k\}$ and $\{z^k\}$ are generated by Algorithm 2.1, then $\{x^k\}$ and $\{z^k\}$ are both bounded. Furthermore, for any $\overline{x}$ such that $F(\overline{x}) = 0,$ it holds that

and

Now we give the global convergence result of the Algorithm 2.1.

Lemma 3.3. Let $\{x^k\}$ be generated by Algorithm 2.1, assume Assumption (A1) and (A2) hold, and there exists constants $0 < \underline{\rho} < \overline{\rho} < 1,$ and $\underline{\mu} < \overline{\mu}$ such that $\underline{\rho}\leq\rho_k\leq\overline{\rho},$ and $\underline{\mu}\leq \mu_k\leq\overline{\mu}.$ Then $\{x^k\}$ converges to some $x^*$ such that $F(x^*) = 0.$

Proof. By the inequality (2.5), we have

By the definition of $d^k$, we have that

and

Combining (3.9) and (3.10), we can assume that there exists a positive constant $c_1$ such that

On the other hand, the definition of $d^k$ also gives that

From (2.4) and Assumption (A2), we have

Combining (3.12) and (3.13), we can assume that there exists a positive constant $c_2$ such that

Now by (3.8), we obtain

By the continuity of $F(x)$, the bound of sequence $\{z^k\}$ and (3.6), we have

We consider the two possible cases:

In the first case, the continuity of $F$ and the boundness of $\{x^k\}$ imply that the sequence $\{x^k\}$ has some accumulation point $x^*$ such that $F(x^*) = 0.$ Since $\overline{x}$ was an arbitrary solution, we can choose $\overline{x} = x^*$ in (3.5). The sequence $\{\|x^k-x^*\|\}$ converges and since $x^*$ is an accumulation point of $\{x^k\},$ it must be the case that $\{x^k\}$ converges to $x^*.$

Now consider the second case. From (3.14), we have

Hence by (3.16), we have

(The following proof is very similar to the last part in Theorem 2.1 ^[23], for complement, we list it here.) By the step rule, we have the inequality (2.5) is not valid for the value $\beta^{m_k-1},$ i.e.,

Let $k\rightarrow \infty,$ we get

Here $x^*$, $d^*$, $\rho^*$, $\mu^*$ denote the limits of the corresponding sequence respectively. On the other hand, by (3.4), we get

that contradicts the choice for $\lambda\in (0, 1).$ Hence the case $\liminf_{k\rightarrow\infty}\|F(x^k)\|$ is impossible.

This completes the proof.

4.   Numerical experiments

In this section, we demonstrate the numerical performance of Algorithm 2.1 (AQN) and its computational advantage by comparing with the modified Kanzow ^[13] ACTN method (denoted as AKP) and the classical Quasi-Newton method with project (denoted as CQN). All presented codes are written in MATLAB2019 and run on a PC with 3.30GHz CPU processor, 4.0GB memory and Windows 8 operation system.

We consider ten problems with dimension n = 1000, 5000, 10000. We use six different starting points, that is:

After several parameter selection experiments, we select the initial parameters that can make the three algorithms have better performance :

Set the terminating criterion for the iteration process as $||F(x_k)||\leq 10^{-6}$. The problems are listed as follows.

Problem 1. ^[25]

where $\Omega = \mathbb{R}^{n}_{+}.$

Problem 2. ^[25]

where $\Omega = \mathbb{R}^{n}_{+}.$

Problem 3. ^[25]

where $\Omega = \mathbb{R}^{n}_{+}.$

Problem 4. ^[25]

where $\Omega = \mathbb{R}^{n}_{+}.$

Problem 5. ^[25]

where $\Omega = \mathbb{R}^{n}_{+}.$

Problem 6. ^[25]

where $h = \frac{1}{n+1}$ and $\Omega = \mathbb{R}^{n}_{+}.$

Problem 7. ^[25]

where $\Omega = \mathbb{R}^{n}_{+}.$

Problem 8. ^[26]

where $\Omega = \mathbb{R}^{n}_{+}.$

Problem 9. ^[26]

where $\Omega = \mathbb{R}^{n}_{+}.$

Problem 10. ^[24]

where $\Omega = \mathbb{R}^{n}_{+}.$

Comprehensive results of our numerical experiment are presented in Tables 1–10. The columns of the presented tables have the following definitions:

IP: the initial points.

DIM: the dimension of the problem.

NI: the iterative number.

NF: the iterative number of function evaluation.

CPU: the CPU time in seconds when the algorithm terminate.

NORM: the final norm equation.

We denote result by '$-$' whenever the number of iterations exceeds 500 or the terminating criterion has not been satisfied. Among these results, none of the three methods were able to solve Problem 9 when initial point is $x_3 = (2, 2, ..., 2)^{T}$. Therefore, Table 9 does not include the case when the initial point is $x_3$. Meanwhile, in the drawing process, when the result was denoted by '$-$', its NI, NF, CPU and NORM are counted as $\infty$.

The performance of the three methods was evaluated using the performance profile which is presented by Dolan and Moré ^[27]. We comparing three methods with the same problem, dimension and initial point in an experiment, and recoding information of interest such as NI, NF, CPU and NORM.

We denote the set of problems as $\mathcal{P}$ and the set of methods as $\mathcal{M}$. For example, for each problem $p$ and method $m$, we define

Compare the performance on problem $p$ by method $m$ with the best performance by any method on this problem, that is, we use the performance ratio

We assume that a parameter $R\geq r_{p, m}$ for all $p, m$ is chosen, and $r_{p, m} = R$ if and only if method $m$ does not solve problem $p$. If method $m$ can solve problem $p$ successfully, we obtain an overall assessment of the performance between these methods. It can be described as follows:

where $n_p$ represents the number of elements in set $\mathcal{P}$, then $\rho_m(\tau)$ is the probability for method $m\in \mathcal{M}$ that a performance ratio $r_{p, m}$ is within a factor $\tau \in \mathbb{R}$ of the best possible ratio. The function $\rho_m$ is the (cumulative) distribution function for the performance ratio.

The performance profile $\rho_m: \mathbb{R}\longmapsto [0, 1]$ for a method is a nondecreasing, piecewise constant function, continuous from the right at each breakpoint. We are interested in methods with a high probability of solve success, then we need only to compare the values of $\rho_m(\tau)$ for all of the methods and choose the method with the largest, there means that we need to find which method's function $\rho_m$ first rearch the line $\rho_m(\tau) = 1$. In the same way, we can obtain the performance profile with respect to NI, NF and NORM.

As can be seen from the information in the Table 1–10, AQN has more stable solving performance and can solve more problems, such as what AKP cannot solve: $x_3$ and $x_5$ of Problem 3 when $n = 10000$; $x_2$ of Problem 3 when $n = 5000, 10000$; $x_2$ of Problem 5 when $n = 1000, 5000, 10000$; $x_3$ of Problem 5 when $n = 5000, 10000$; $x_2$ of Problem 8 when $n = 5000, 10000$; $x_2$ of Problem 10 when $n = 5000, 10000$. Compared with CQN, from Figure 1 and Figure 2, we can see that AQN reaches the line that $\rho_m(1) = 1$ before CQN, which demonstrates AQN has a faster solution time(CPU) and the solution results of the final norm equation(NORM) are more accurate. Although from Figure 3 and Figure 4, there shows the iterative number of CQN is less than AQN, in practice, we pay more attention to the advantage of solution time. To sum up, AQN has a more stable and faster solving performance.

5.   Conclusions

In this paper, we propose an active set quasi-Newton method for the solution of optimization problem with bound constraints. The implementation of the method uses the quasi-Newton step as a trial step and the project step as the correction step. By using active set technique, we only need to solve a reduced dimension linear equation at each iteration to generate the search direction. We prove that the generated sequence is bounded automatically and obtain the global convergence of the proposed algorithm. Meanwhile, compared with other algorithms, our method has the most stable performance. There are some questions that need studying in the near future. Firstly, it is possible to get the global convergence of the proposed algorithm without the assumption of the positive definite of the matrix $B_k$. Secondly, how to get the local convergence of the proposed algorithm especially under some weak condition such as the local error bound condition needs further studying.

Conflict of interest

All authors declare no conflicts of interest in this paper.