New versions of Newton method: step-size choice, convergence domain and under-determined equations

TL;DR

This paper introduces new Newton-like algorithms with adaptive step-size choices that improve convergence properties and applicability to under-determined systems, including explicit complexity analysis and global convergence under certain conditions.

Contribution

It proposes novel Newton-type methods with adaptive step-sizes, extending applicability to under-determined systems and providing convergence domain estimates and complexity results.

Findings

New algorithms with improved convergence domains.

Global convergence under specific assumptions.

Explicit complexity bounds for the proposed methods.

Abstract

Newton method is one of the most powerful methods for finding solutions of nonlinear equations and for proving their existence. In its "pure" form it has fast convergence near the solution, but small convergence domain. On the other hand damped Newton method has slower convergence rate, but weaker conditions on the initial point. We provide new versions of Newton-like algorithms, resulting in combinations of Newton and damped Newton method with special step-size choice, and estimate its convergence domain. Under some assumptions the convergence is global. Explicit complexity results are also addressed. The adaptive version of the algorithm (with no a priori constants knowledge) is presented. The method is applicable for under-determined equations (with $m<n$, $m$ being the number of equations and $n$ being the number of variables). The results are specified for systems of quadratic equations, for composite mappings and for one-dimensional equations and inequalities.

Equations

$$P(x)=0,$$

$$x^{k+1}=x^k-P^{\prime }(x^k{)}^{-1}P(x^k).$$

$$x^{k+1}=x^k-\alpha P^{\prime }(x^k{)}^{-1}P(x^k)$$

$$\begin{array}[]{l}x^{k+1}=x^{k}-\alpha_{k}z^{k},\;\;k=0,1,\ldots\\ z^{k}\in{\rm{Arg}}\,\min_{z}\{||z||:P^{\prime}(x^{k})z=P(x^{k})\}.\end{array}$$

$$x^{k+1}=x^k-P^{\prime }(x^k{)}^{†}P(x^k),$$

$$g(x)=y,$$

$$g_i(x)=\frac{1}{2}(A_{i}x,x)+(b_i,x),A_i=A_i^T,b_i\in {\mathbb{R}}^n.$$

$$\begin{array}[]{l}x^{k+1}=x^{k}-z^{k},\;\;k=0,1,\ldots\\ \displaystyle z^{k}=\arg\min_{z}\{||P(x^{k})-P^{\prime}(x^{k})z||+M||z||^{2}\},\end{array}$$

$$||g(x^a)-g(x^b)-A(x^a-x^b)||\le C||x^a-x^b||$$

$$\min \{||z|{|}_1:P^{\prime }(x^k)z=P(x^k)\}.$$

$$\min \{||z|{|}_{\infty }:P^{\prime }(x^k)z=P(x^k)\}.$$

$$z^k=P^{\prime }(x^k{)}^{†}P(x^k).$$

$$Az=b,b\in {\mathbb{R}}^m,z\in {\mathbb{R}}^n$$

$$\Vert A^{T}h{\Vert }_{\ast }\ge {\mu }_0\Vert h{\Vert }_{\ast },{\mu }_0>0,$$

$$\hat{z}\in \mathrm{Arg}\min \{\Vert z\Vert :Az=b\}$$

$$H_k(\delta )=\sum _{\ell =k}^{\infty }{\delta }^{(2^{\ell })}.$$

$$c_{1}=H_{0}\Big{(}\frac{1}{2}\Big{)}\approx 0.8164215,$$

$$c_2=\underset{0\le r\le \frac{1}{4}}{\max }2H_0\left(\frac{1}{2}-r\right)+5r-4r^2-2c_1\approx \mathrm{0.0036003.}$$

$$0\le H_k(\delta )\le \frac{{\delta }^{(2^k)}}{1-{\delta }^{(2^k)}}=\frac{1}{{\delta }^{-(2^k)}-1}.$$

$$\begin{array}[]{l}x^{k+1}=x^{k}-\alpha_{k}z^{k},\;\;k=0,1,\ldots\\ z^{k}\in{\rm{Arg}}\,\min_{z}\{||z||:g^{\prime}(x^{k})z=g(x^{k})-y\}.\end{array}$$

$$\Vert g^{\prime }(x^a)-g^{\prime }(x^b)\Vert \le L\Vert x^a-x^b\Vert .$$

$$\Vert g^{\prime }(x{)}^{T}h{\Vert }_{\ast }\ge \mu \Vert h{\Vert }_{\ast },\forall h\in {\mathbb{R}}^m.$$

$$g(x+z)=g(x)+{\int }_0^1{g}^{\prime }(x+tz)zdt$$

$$u_{k+1}\le |1-\alpha |u_k+\frac{L{\alpha }^2}{2}\Vert z^k{\Vert }^2.$$

$$u_{k+1}\le |1-\alpha |u_k+\frac{L{\alpha }^2{u}_k^2}{2{\mu }^2}.$$

$$\Vert x^{k+s}-x^k\Vert \le \sum _{i=k}^{k+s-1}\Vert x^{i+1}-x^i\Vert \le q^k\frac{u_0}{(1-q)\mu }\to 0.$$

$$\Vert P^{\prime }(x^a)-P^{\prime }(x^b)\Vert \le L\Vert x^a-x^b\Vert .$$

$$\Vert P^{\prime }(x{)}^{T}h{\Vert }_{\ast }\ge \mu \Vert h{\Vert }_{\ast },\forall h\in {\mathbb{R}}^m.$$

$$u_{k+1}\le |1-{\alpha }_k|u_k+\frac{L{\alpha }_k^2\Vert z^k{\Vert }^2}{2},$$

$$u_{k+1}\le |1-{\alpha }_k|u_k+\frac{L{\alpha }_k^2{u}_k^2}{2{\mu }^2},$$

Full text

Keywords: Nonlinear equations; Newton method; under-determined equations; global convergence; adaptive algorithms; metric regularity

1 Introduction

Consider nonlinear equation

[TABLE]

written via the vector function $P:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$. There exists the huge bunch of literature on solvability of such equations and numerical methods for their solution, see e.g. the classical monographs [22, 4]. One of the most powerful methods is Newton method, going back to such giants as Newton, Cauchy, Fourier. The general form of the method is due to Kantorovich [14]; on history and references see [15, 33, 26, 5]. The basic version of Newton method for (1) is applicable when $P(x)$ is differentiable and $P^{\prime}(x)$ is invertible (this implies $m=n$):

[TABLE]

The method converges under some natural conditions, moreover it can be used for obtaining existence theorems for the solution (see references cited above). Unfortunately Newton method converges only locally: it requires a good initial approximation $x^{0}$ (so called “hot start”). Convergence conditions can be relaxed for damped Newton method

[TABLE]

with $0<\alpha<1$.

The advanced, generalized writing of Newton method is

[TABLE]

This variant relies on the solvability of the linear equation only, and it also admits non-constant step-size. It is applicable to under-determined systems of equations ($m<n$) and to non-linear equations in Banach space. The latter are outside of the scope of this paper but its analysis is essentially the same.

If $m=n$ and $P^{\prime}(x^{k})^{-1}$ exists, the method (3) coincides with classical Newton method for $\alpha_{k}=1$ and damped Newton method for $\alpha_{k}=\alpha<1$. Starting at some initial point $x^{0}$ the latter method converges to $x^{*}$ under some additional constraints on residual $||P(x^{0})||$ and function-related constants $L$, $\rho$, $\mu$, $\alpha$ (see Theorems below for rigorous conditions).

In explicit form Newton method for $m\neq n$ has been written by Ben-Israel [1]:

[TABLE]

where $A^{\dagger}$ stands for Moore-Penrose pseudoinverse of $A$. However the results in [1] are mostly oriented on over-determined systems, and the assumptions of the theorems in [1] are hard to check. Other publications on under-determined equations include [11, 21, 23, 31, 19]. Moreover there exist numerous literature on more general settings: equalities plus inequalities [29, 30], optimization problems [3, 10] with more general algorithms, which can be applied to solving of equations as particular case. In next Section 2 we discuss the under-determined finite-dimensional case in more details.

There is a very similar problem statement to (1), made in terms of the equation with the variable right-hand side

[TABLE]

with $g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ having a known solution $\overline{x}:g(\overline{x})=0$ and the variable $y$ as a parameter.

The question is: for which right-hand side part $y$ the equation is still feasible and what are the solutions? This problem arises in finding image of a mapping $\{g(x):x\in\mathbb{R}^{n}\}$, checking robustness/sensitivity of a solution or exploration problem of the image, etc. In general, this problem is hardly solvable, but we can provide local sufficient conditions of feasibility, imposed on $y$.

Trivially, equation (5) can be written in the form (1) with $P_{y}(x)=g(x)-y=0$ and $\|P_{y}(\overline{x})\|=\|y\|$. Thus conditions on feasibility of the equation $P_{y}(x)=0$ are exactly conditions on feasibility (5) with respect to right-hand side $y$, and vice-versa.

Let us explain the connection between (1) and (5) deeper. There are few approaches to treating feasibility of an equation. One of them is to prove existence of the solution by providing semi-local existence theorems. These involve conditions in some point and/or around it, and prove that if such conditions hold, then a solution exists. It is not necessary to provide tools for finding this solution.

Another way is to explore a constructive, algorithmic way of solving the equation, starting at a point $x^{0}$, resulting in a sequence $\{x^{k}\}$, and prove convergence to a solution $x^{k}\rightarrow x^{*}$ (e.g. fixed-point theorems). The convergence conditions are typically tied to the sequence, including starting point $x^{0}$. The conditions do not necessarily coincide with the conditions of semi-local existence theorems. Moreover, the conditions ensuring faster convergence of the algorithm are typically more strict, than the conditions of the semi-local existence theorems.

In Newton method theory these approaches are closely connected, and semi-local theorems are often proved via convergence of variants of Newton method. We also show it below in Theorem 3.2. This relation becomes very clear in comparison of equations (1) and (5). Naturally feasibility of (1) being solved by a Newton-like method (i.e. convergence of the Newton method started at $x^{0}$!) is stated in terms of norm of initial residual, say $\|P(x^{0})\|\leq s$. In terms of (5) the very same Newton-like method, being applied to the constructed $P_{y}(\cdot)$ and being started at the same point $\overline{x}=x^{0}$, converges for any fixed $y:\|y\|\leq s$. This analysis claims feasibility of (5) for all such $y$. Therefore finding the largest set of possible $y$ is essentially the same as the problem of finding the broadest residual range $P(x^{0})$. Without loss of generality, through the paper we assume $\overline{x}=0$, and switch between problems (1) and (5) as equivalent ones. The difference is clear from context.

We also examine some special cases of the nonlinear equations. One of them is the quadratic case, when all components of $g$ are quadratic functions:

[TABLE]

In this case we try to specify above results and design the algorithms to check feasibility of a vector $y\in R^{m}$.

The first goal of the present paper is to give explicit expressions of the method (3) for various norms and to provide simple, easily checkable conditions for convergence of the method. This also provides existence theorems: what is a feasible set $Y$ such that $y\in Y$ implies solvability of (5).

The second goal is to develop constructive algorithms for choosing step-sizes $\alpha_{k}$ to achieve fast and as global as possible convergence. We suggest different strategies for constructing algorithms, study their properties and provide explicit convergence conditions for the method and demonstrate its potentially global convergence.

The main contributions of the paper are threefold.

We propose the novel procedure for adjusting step-size $\alpha_{k}$. The strategy guarantees wide range of convergence (in some cases the algorithm converges globally) and fast rate of convergence (local quadratic convergence, typical for pure Newton method). Moreover explicit formula for method’s complexity is provided. 2. 2.

The choice of norms in the algorithm can be different, thus we arrive to different versions of the algorithm. For instance, Euclidean norms imply explicit form of desired direction $z^{k}$ (the same as in (4)) while $\ell_{1}$ norm provides sparse approximations etc. 3. 3.

We consider numerous applications, including under-determined cases, quadratic equations, one equation with $n$ variables.

Few words on comparison with known results. In the paper [23] results on solvability of nonlinear equations in Banach spaces and on application of Newton-like methods have been formulated in semi-local form. One of the results from [23] adopted to our notation and finite-dimensional case claims that if $P^{\prime}(x)$ exists and is Lipschitz on a ball $B$ centered in $x^{0}$ and estimate $||P^{\prime}(x)^{T}h||_{*}\geq\mu||h||_{*},\;\mu>0,\;\forall h$ holds on $B$, then equation (1) has a solution $x^{*}$ provided $||P(x^{0})||<\frac{\rho}{\mu}$. Another result deals with convergence of Newton method; however the method is not provided in explicit form. The condition on the derivative has extension to non-differentiable functions and multi-valued mapping an is known as metric regularity, used for proving existence theorems in different cases.

The paper, which contains the closest results to ours, is [21]. Nesterov addresses the same problem (1) and his method (in our notation) has the form

[TABLE]

where $M$ is the scalar parameter to be adjusted at each iteration. Nesterov’s assumptions are close to ours and his results on solvability of equations and on convergence of the method are similar. The main difference is the method itself; it is not clear how to solve the auxiliary optimization problem in Nesterov’s method, while finding $z^{k}$ in our method can be implemented in explicit form. Other papers on under-determined equations mentioned above either do not specify the technique for solving the linearized auxiliary equation, or restrict analysis with Euclidean norm and/or pure Newton step-size $\alpha_{k}=1$, see e.g. [23, 28, 31, 16].

The rest of the paper is organized as follows. Next section is introductory to the case of under-determined systems. In Section 3 we remind few notions and results. Next we prove simple solvability conditions for (1). In main Section 4 we propose few variants of general Newton algorithm (3), including adaptive ones and estimate their convergence rate. Some particular cases (scalar equations and inequalities, quadratic equations, problems with special structure) are treated in Section 5. Results of numerical simulation are exhibited in Section 6. Conclusion part finalizes the paper (Section 7).

2 Under-determined systems of equations

Under-determined equations attracted our attention by specific norm-dependency property. In case of $m<n$ norms in $\mathbb{R}^{n}$ (in the optimization sub-problem) and $\mathbb{R}^{m}$ (for residual) can be chosen arbitrarily, and they imply principally different forms and results of Newton method (3). Conditions on solvability and convergence look similar, but the results differ strongly.

Historically the case of under-determined equations ($m<n$) attracted much less attention than equations with the same number of equations and variables. The pioneering result is due to Graves [9] in more general setting of Banach spaces, for problem (5). Graves’ theorem for finite-dimensional case claims, that if condition

[TABLE]

holds in the ball of radius $\rho$, centered at zero, for a matrix $A$ with minimal singular value $\mu>C>0$, then a solution of the equation (5) exists provided $||y||$ is small enough, namely $||y||\leq\rho(\mu-C)$. The solution can be found via a version of modified Newton method, where next iteration requires solution of the linear equation with matrix $A$, see [6, 18] for details. The condition above gives rise to the mentioned metric regularity property. However, finding the matrix $A$ is still a problem itself.

First of all let us specify the subproblem of finding a vector $z^{k}$ in (3) for different norms of $x\in\mathbb{R}^{n}$. We skip simple verifications of the statements from convex analysis.

For $||x||=||x||_{1}$ vector $z^{k}$ is a solution of the problem

[TABLE] 2. 2.

For $||x||=||x||_{\infty}$ vector $z^{k}$ is a solution of the problem

[TABLE]

Both problems above can be easily reduced to linear programming. 3. 3.

For $||x||=||x||_{2}$ vector $z^{k}$ can be written explicitly

[TABLE]

In this case Newton method (3) coincide with (4). For $m\leq n$ Moore-Penrose pseudo-inverse of a matrix $A$ is written as $A^{\dagger}=A^{T}(AA^{T})^{-1}$, if $A$ has full row rank.

Thus in these (most important) cases algorithm (3) can be implemented effectively. Also the solution of the first two problems may be non-unique.

An important case is the scalar one, i.e. $m=1$. We specify general results for scalar equations and inequalities; the arising algorithms have much in common with unconstrained minimization methods. Finally we discuss nonlinear equations having some special structure. Then convergence results can be strongly enhanced.

3 Preliminaries and feasibility (existence) theorems

Key component in Newton method is the auxiliary convex optimization sub-problem, involving the linear constraint. Note that the constraint

[TABLE]

describes either a linear subspace, or the empty set. The classical result below (which goes back to Banach, see [14, 18, 21]) guarantees solvability of the linear equation (7) and gives an estimate of its solution. We prefer to provide the direct proof of the result because it is highly clear and short in finite-dimensional case. Suppose that spaces $\mathbb{R}^{n},\mathbb{R}^{m}$ are equipped with some norms, the dual norms are denoted $||\cdot||_{*}$ (for a linear functional $c$, associated with the vector of the same dimension, $\|c\|_{*}=\sup_{x:\|x\|=1}(c,x)$). Operator norm is subordinate with the vector norms, e.g. for $A:X\rightarrow Y$ we have $\|Ax\|_{Y}\leq\|A\|_{X,Y}\|x\|_{X}$. In most cases we do not specify vector norms; dual norms are obvious from the context. The adjoint operator $A^{*}$ is identified with matrix $A^{T}$.

Lemma 3.1**.**

If $A\in\mathbb{R}^{m\times n}$ satisfies condition

[TABLE]

for all $h\in\mathbb{R}^{m}$, then equation (7) has a solution for all $b\in\mathbb{R}^{m}$, and all solutions of optimization problem

[TABLE]

have bounded norms $\|\widehat{z}\|\leq\frac{1}{\mu_{0}}{\|b\|}$.

Proof.

Fix $b\in\mathbb{R}^{m}$ and denote $K=\{x\in\mathbb{R}^{n}:||x||\leq\frac{\|b\|}{\mu_{0}}\}$. This is a convex closed bounded set, and its linear image $Q=\{y\in\mathbb{R}^{m}:y=Ax,x\in K\}$ is convex closed bounded set as well. Suppose $b\notin Q$, then it can be strictly separated from $Q$: there exists $c\in\mathbb{R}^{m}:\max_{y\in Q}(c,y)<(c,b)$. But $\max_{y\in Q}(c,y)=\max_{x\in K}(c,Ax)=\max_{x\in K}(A^{T}c,x)=\frac{\|b\|}{\mu_{0}}\|A^{T}c\|_{*}\geq\|b\|\cdot\|c\|_{*}$, thus we get the contradiction: $\|b\|\cdot\|c\|_{*}<(c,b)$. Hence $b\in Q$, i.e. there exists $x\in\mathbb{R}^{n}:Ax=b,\|x\|\leq\frac{||b||}{\mu_{0}}\}$. A solution with the least norm obeys the same inequality.

∎

The Lemma is claiming that the matrix $A$ has full row rank equal to $m$ provided (8) holds. It is another way to say that the mapping $A:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ is onto mapping, i.e. covering all image space. In the case of Euclidean norms, parameter $\mu_{0}$ is the smallest singular value of the matrix $\mu_{0}=\sigma_{m}(A)$ (we denote singular values of a matrix in $\mathbb{R}^{m\times n}$ in decreasing order as $\sigma_{1}\geq\sigma_{2}\geq\ldots\geq\sigma_{m}$). In general case the conjugate operator $A^{*}$ is used instead of $A^{T}$, and $\mu$ is the metric regularity constant.

Some results below will exploit the sum of double exponentials functions $H_{k}\colon[0,1)\rightarrow\mathbb{R}_{+}$, cf. [23]:

[TABLE]

All functions $H_{k}(\cdot)$ are monotonically increasing and strictly convex. We also use two specific constants

[TABLE]

and

[TABLE]

Trivial approximations

[TABLE]

may be used for polynomial lower and upper bounds of $H_{0}(\delta)=\delta+\delta^{2}+\delta^{4}+\ldots+\delta^{(2^{k-1})}+H_{k}(\delta)$ with arbitrary precision. We also use property $H_{k}(\delta^{2})=H_{k+1}(\delta)$.

Below the problem of solvability of equation (5) is addressed. We apply algorithm (3) in the form

[TABLE]

with small $\alpha$ and prove that the iterations converge while the limit point is a solution. This techniques follows the idea from [23]. Remind that $\mathbb{R}^{n},\mathbb{R}^{m}$ are equipped with some norms, the dual norms are denoted $||\cdot||_{*}$.

Assumptions with respect to (5).

$\mathbf{A}.$ $g(0)=0$ (i.e $\overline{x}=0$), $g(x)$ is differentiable in the ball $B=\{x\in\mathbb{R}^{n}:\|x\|\leq\rho\}$, and its derivative $g^{\prime}(x)$ satisfies Lipschitz condition in $B$:

[TABLE]

$\mathbf{B}.$ The following inequality holds for all $x\in B$ and some fixed $\mu>0$:

[TABLE]

$\mathbf{C}.$ $\|y\|<\mu\rho$.

Theorem 3.2**.**

If conditions $\mathbf{A},\mathbf{B},\mathbf{C}$ hold then there exists a solution $x^{*}$ of (5), and $\|x^{*}\|\leq\frac{\|y\|}{\mu}$.

Proof.

We apply algorithm (11) with $\alpha>0$ small enough and $x^{0}=0$. The algorithm is well defined — condition $\mathbf{B}$ and Lemma 3.1 imply existence of solutions $z^{k}$ provided that $x^{k}\in B$; this is true for $k=0$ and will be validated recurrently. Standard formula

[TABLE]

combined with condition $\mathbf{A}$ provides for $x=x^{k},\;z=-\alpha z^{k}$ and $u_{k}=\|g(x^{k})-y\|$ recurrent relation

[TABLE]

Now condition $\mathbf{B}$ and Lemma 3.1 transform this estimate into

[TABLE]

Choose $\alpha=\varepsilon\frac{2\mu^{2}}{Lu_{0}}(1-\frac{u_{0}}{\mu\rho})$ with small $\varepsilon<1$ satisfying $0<\alpha<1$; it is possible due to condition $\mathbf{C}$. From the above inequality we get $u_{k+1}\leq u_{k}(1-\alpha+\alpha\varepsilon\frac{u_{k}}{u_{0}}(1-\frac{u_{0}}{\mu\rho}))$. For $k=0$ this implies $u_{1}<u_{0}$ and recurrently $u_{k+1}<u_{k}$. We also get $u_{k+1}\leq qu_{k},\;q=(1-\alpha+\alpha\varepsilon(1-\frac{u_{0}}{\mu\rho}))<1$. Thus $u_{k}\leq q^{k}u_{0}$ and $u_{k}\rightarrow 0$ for $k\rightarrow\infty$.

On the other hand we have $\|x^{k+1}-x^{k}\|=\alpha\|z^{k}\|\leq\frac{\|g(x^{k})-y\|}{\mu}=\frac{u_{k}}{\mu}\leq q^{k}\frac{u_{0}}{\mu}$. Hence for any $k,s$ and for $k\rightarrow\infty$

[TABLE]

It means that $x^{k}$ is a Cauchy sequence and converges to some point $x^{*}(\varepsilon)$. We had $g(x^{k})\rightarrow y$, thus continuity reasons imply $g(x^{*}(\varepsilon))=y$. Now, for all iterations we got $\|x^{k}-x^{0}\|=\|x^{k}\|\leq\sum_{j=0}^{k-1}\|x^{j+1}-x^{j}\|\leq\alpha\frac{u_{0}}{\mu(1-q)}=\frac{u_{0}}{\mu}\frac{1}{1-\varepsilon(1-\frac{u_{0}}{\mu\rho})}<\frac{u_{0}}{\mu}\frac{1}{1-(1-\frac{u_{0}}{\mu\rho})}=\rho$. Hence all iterations $x^{k}$ remain in the ball $B$ and our reasoning was correct. Finally under $\|x^{k}\|\leq\frac{u_{0}}{\mu}\frac{1}{1-\varepsilon(1-\frac{u_{0}}{\mu\rho})}$ we take $\varepsilon\rightarrow 0$, leading to $\|x^{k}\|\leq\frac{u_{0}}{\mu}$, so its limit point $x^{*}(\varepsilon)|_{\varepsilon\rightarrow 0}$. The limit point $x^{*}(\varepsilon)|_{\varepsilon\rightarrow 0}=x^{*}$ is a solution as well and $\|x^{*}\|\leq\frac{u_{0}}{\mu}$. ∎

Corollary 3.3**.**

If $\rho=\infty$ (that is conditions $\mathbf{A},\mathbf{B}$ hold on the entire space $\mathbb{R}^{n}$) then equation (5) has a solution for an arbitrary right-hand side $y$.

It is worth noting that if we apply pure Newton method (i.e. take $\alpha_{k}\equiv 1$), the conditions of its convergence are more restrictive: we need $\|y\|\leq\frac{2\mu^{2}}{L}$, that is we guarantee only local convergence even for $\rho=\infty$. This is a corollary of Newton-Mysovskikh theorem [14], which proof is valid for under-determined case as well, cf. also [23].

Corollary 3.4**.**

If $m=n$ and Condition $\mathbf{B}$ is replaced with $\|g^{\prime}(x)^{-1}\|\leq\frac{1}{\mu},x\in B$, then the statement of Theorem 3.2 holds true.

In this case our method (11) reduces to classical Newton method (2).

There exist numerous results on solvability of (5). Some of them are stronger than Theorem 3.2 and are based on general notion of metric regularity [7, 12]. We provided the proof based on our technique to exhibit its applicability to existence theorems.

4 Main algorithms

Here it is more convenient to use the main equation in form (1). In previous Section we proved solvability of equation by use of the algorithm with constant $\alpha_{k}\equiv\alpha>0$; choosing $\alpha$ smaller we obtained larger solvability domain.

However in this Section our goal is different — to reach the fastest convergence to a solution. For this purpose different strategies for design of step-sizes are needed. The basic policy is as follows. First, we rewrite assumptions in new notation. We remark that the assumptions in context of equation $P(x)=0$ are very much the same as $\mathbf{A},\mathbf{B}$.

$\mathbf{A^{\prime}}.$ $P(x)$ is differentiable in the ball $B=\{x\in\mathbb{R}^{n}:\|x-x^{0}\|\leq\rho\}$, and its derivative $P^{\prime}(x)$ satisfies Lipschitz condition in $B$:

[TABLE]

$\mathbf{B^{\prime}}.$ The following inequality holds for all $x\in B$ and some $\mu>0$:

[TABLE]

If $\mathbf{A^{\prime}},\mathbf{B^{\prime}}$ hold true, we have the same recurrent inequalities for $u_{k}=\|P({x^{k}})\|$:

[TABLE]

[TABLE]

the second one being just continuation of the first one based on the estimate $\|z^{k}\|\leq\frac{u_{k}}{\mu}$, compare with the calculations in the proof of Theorem 3.2. Now we can minimize right-hand sides of these inequalities over $\alpha_{k}$; it is natural to expect that such choice of step-size imply the fastest convergence of $u_{k}$ to zero and thus the fastest convergence of iterations $x_{k}$ to the solution. One of the main contributions of this paper is careful analysis of the resulting method.

If one applies such policy based on inequality (13), optimal $\alpha$ depends on $\mu,L$ (Algorithm 1 below). The values are hard to estimate in most applications, thus the method would be hard for implementation. Fortunately, we can modify the algorithm using parameter adjustment (Algorithm 2). On the other hand the same policy based on (12) requires just the value $L$, which is commonly available (Algorithm 3).

Thus we arrive to an algorithm which we call Newton method while in fact it is blended pure Newton with damped Newton with special rule for damping. In some relation it reminds Newton method for minimization of self-concordant functions [20]. Despite its simplicity, the idea of minimizing upper bound seems to be unexplored (or long forgotten) in Newton method theory with respect to equations. Authors found similar choice of step-size in [5], in different conditions and without explicit convergence bounds.

4.1 Newton method with known constants

If both constants $L$ and $\mu$ are known, then the step-size is taken as the minimizer of right-hand side of (13):

[TABLE]

Algorithm 1 (Basic Newton method)

[TABLE]

[TABLE]

The algorithm is well-defined, as soon $\|P(x^{k})\|=0$ means that a solution is already found (formally $z^{k}=0,\;\alpha_{k}=0$ thereafter). We remind that in calculation of $z^{k}$ any vector norm in $\mathbb{R}^{n}$ can be used, also any vector norm in $\mathbb{R}^{m}$ can be used for $\|P(x^{k})\|$, and constants $L,\mu$ must comply with these norms.

The update step in (15) can be written in less compact but more illustrative form:

[TABLE]

The latter case is a pure Newton step while the primal one is a damped Newton step. Direction $z^{k}$ calculation is the same in both stages. The result on convergence and rate of convergence is given below. We use upper ($\lceil\cdot\rceil$) and lower ($\lfloor\cdot\rfloor$) rounding to integer; constant $c_{1}\approx 0.8164$ was introduced in the end of Section 3. The theorem is followed by the corollary with simpler statements.

Theorem 4.1**.**

Suppose that Assumptions $\mathbf{A^{\prime}},\mathbf{B^{\prime}}$ hold and

[TABLE]

where $F^{{\rm{inv}}}_{1}(\cdot)$ is the inverse function for the continuous strictly increasing function $F_{1}(w)$, given by

[TABLE]

Then the sequence $\{x^{k}\}$ generated by Algorithm 1 converges to a solution $x^{*}:P(x^{*})=0$. The values $\|P(x^{k})\|$ are monotonically decreasing, and there are not more than

[TABLE]

iterations at Stage 1, then followed by Stage 2 steps. At $k$-th step the following estimates for the rate of convergence hold:

[TABLE]

where \overline{w}=\frac{L}{\mu^{2}}\|P(x^{0})\|-\frac{k_{\max}}{2}=\min\big{\{}\frac{L}{\mu^{2}}\|P(x^{0})\|, 1-\frac{1}{2}\big{\lceil}\frac{2L}{\mu^{2}}\|P(x^{0})\|\big{\rceil}+\frac{L}{\mu^{2}}\|P(x^{0})\|\big{\}}\in[0,1).

The Theorem’s statement may look quite involved, but both functions $H_{0}(\delta)$ and $F_{1}^{{\rm{inv}}}(p)$ are easily calculated in practice. In the interval of interest $\delta\in[0,\frac{1}{2}]$, the former function has rational approximation (10). The latter function can be evaluated with needed accuracy via binary search, as soon $F_{1}(w)$ is monotonically increasing on $R_{+}$.

Proof.

Assume that all $x^{k}\in B,\;k\geq 0$. Below we state condition enabling this assumption. Using $w_{k}=\frac{L}{\mu^{2}}\|P(x^{k})\|$ as the objective function, we rewrite (13) with generic step-size $\alpha$ as

[TABLE]

Its optimum over $\alpha$ is at $\alpha_{k}=\frac{1}{w_{k}}<1$, if $w_{k}>1$; and $\alpha_{k}=1$ otherwise; it is exactly (14).

During Stage 1 of damped Newton steps ($\alpha_{k}<1$), the objective function monotonically decreases as

[TABLE]

There are at most $k_{\max}=\max\{0,\;\lceil 2w_{0}\rceil-2\}$ iterations in the phase, say $\overline{k}$ ones, resulting in $w_{\overline{k}}\leq 1$. As soon $w_{k}$ reaches this unit threshold, the algorithm switches to Stage 2, pure Newton steps. Then recurrent relation (20) becomes

[TABLE]

so we can write

[TABLE]

For the second phase $\|x^{i+1}-x^{i}\|=\|z^{i}\|\leq\frac{1}{\mu}\|P(x^{i})\|=\frac{\mu}{L}w_{i}$ due Lemma 3.1, and for $\ell_{2}\geq\ell_{1}\geq 0$ holds

[TABLE]

The sequence $\{x^{k}\}$ is a Cauchy sequence because $H_{j}(\frac{w_{\overline{k}}}{2})\leq H_{j}(\frac{1}{2})\rightarrow_{j\rightarrow\infty}0$. It converges to a point $x^{*}:\|P(x^{*})\|=\lim_{k\rightarrow\infty}\|P(x^{k})\|=0$ due to continuity of $P$, with

[TABLE]

Next we are to estimate distance from points $x^{k}$ in Stage 1 to the limit solution point $x^{*}$. One-step distance for $k<\overline{k}$ is bounded by a constant: $\|x^{k+1}-x^{k}\|=\alpha_{k}\|z^{k}\|\leq\alpha_{k}{\frac{\mu}{L}}w_{k}=\frac{\mu}{L},$ and altogether

[TABLE]

Note that the formula also coincides with the upper bound (24) at $k=\overline{k}$. Exact number $\overline{k}$ of the steps in the first phase is not known, but we can replace it with the upper bound $k_{\max}$ in all estimates (21)–(25), due to monotonic decrease of $\{w_{k}\}$. We also have an upper bound for $w_{k_{\max}}\leq\overline{w}=w_{0}-\frac{1}{2}k_{\max}=w_{0}-\frac{\max\{0,\lceil 2w_{0}\rceil-2\}}{2}\in[0,1]$. Substituting $w_{k}=\frac{L}{\mu^{2}}\|P(x^{k})\|$ back we arrive to Theorem 4.1 bounds (19).

Finally we are to check our primal assumption of the algorithm-generated points $x^{k}$ being within $B$. This is guaranteed by one of two conditions, depending on whether the Algorithm starts from Stage 1 step or Stage 2 step.

In the first case $w_{0}>1$, and $\|x^{0}-x^{k}\|$ is bounded similarly to (25) as

[TABLE]

Here we also used $k_{\max}=\lceil 2w_{0}\rceil-2>0$ and the upper bound $w_{k_{\max}}\leq\overline{w}$. In other words, given $w_{0}=\frac{L}{\mu^{2}}\|P(x^{0})\|>1$, for $\{x^{k}\}\in B$ it is sufficient to satisfy \frac{L}{\mu}\rho\geq F_{1}\big{(}w_{0}\big{)}. This corresponds to (17b) part.

In the second case we have $w_{0}\leq 1$, and the algorithm makes pure Newton steps with $\alpha_{k}\equiv 1$ from the beginning. Then $\overline{k}=0$, $w_{k_{\max}}=w_{0}$ and from (23) follows

[TABLE]

Therefore if $w_{0}\leq 1$, then inequality $\|x^{0}-x^{k}\|\leq\rho,\;k\geq 0$ is satisfied if \frac{L}{\mu}\rho\geq 2H_{0}\big{(}\frac{w_{0}}{2}\big{)}{=F_{1}(w_{0})}. This corresponds to (17a) part.

Gluing the cases $w_{0}\leq 1$ and $w_{0}>1$ we arrive to the sufficient condition \rho\geq\frac{\mu}{L}F_{1}\big{(}{\frac{L}{\mu^{2}}}\|P(x^{0})\|\big{)}, resulting in $x^{k}\in B$. Due to $F_{1}(w)$ being strictly increasing this condition is equivalent to (16). ∎

Result on the rate of convergence means, roughly speaking, that after no more than $k_{\max}$ iterations one has very fast (quadratic) convergence. For good initial approximations $k_{\max}=0$, and pure Newton method steps are performed from the very start.

Corollary 4.2**.**

If $\rho=\infty$ (that is conditions $\mathbf{A^{\prime}},\mathbf{B^{\prime}}$ hold on the entire space $\mathbb{R}^{n}$), then Algorithm 1 converges to a solution of (1) for any $x^{0}\in\mathbb{R}^{n}$.

The following corollary provides simpler tight relaxed condition for Theorem 4.1. The idea is to develop an upper bound for (17), resulting in a lower bound on (16).

Corollary 4.3**.**

Condition (16) can be replaced with piece-wise linear one

[TABLE]

Proof.

In order to find a lower bound for $F_{1}^{{\rm{inv}}}(\cdot)$, we are to prove an upper bound for the function $F_{1}(w)$, (17). Both bounds are continuous strictly increasing functions. First we notice, that both (17a) and (17b) coincide at $w\in[\frac{1}{2},1]$, and it can be rewritten through the different junction point $w=1/2$ instead of $w=1$

[TABLE]

Due to convexity on interval $[0,\frac{1}{2}]$, function $H_{0}(\delta)$ is bounded by a secant segment:

[TABLE]

Here we used property H_{0}\big{(}\frac{1}{4}\big{)}=H_{0}\big{(}\frac{1}{2}\big{)}-\frac{1}{2}, which follows from the identity $H_{k}(x)=x^{(2^{k})}+H_{k}(x^{2}),\;\;x{\in[0,1)}$; and constant $c_{1}=H_{0}(\frac{1}{2})$ is introduced in Section 3.

Next we show that $F_{1}(w)\leq 2(w+c_{1}-1)$ for $w\geq\frac{1}{2}$. Indeed, following continuous function is periodic on $w,\;w\geq\frac{1}{2}$. It can be also written through the variable $r=\frac{1+2w-\lceil 2w\rceil}{2}\in(0,\frac{1}{2}]$.

[TABLE]

In the second row identity $H_{0}(\frac{1}{2})=\frac{1}{2}+H_{0}(\frac{1}{4})$ is used. The last inequality is due to convexity of $H_{0}(\delta)$, which is under secant segment on corresponding interval $[\frac{1}{4},\frac{1}{2}]$.

Thus we have the monotonically increasing piece-wise linear upper bound for (17)

[TABLE]

Its inverse $F_{1}^{{\rm{inv}},{\rm{lower}}}$ with property $F_{1}^{{\rm{inv}},{\rm{lower}}}(F_{1}^{{{\rm{upper}}}}(w))\equiv w,w\geq 0$ is the piece-wise linear lower bound for $F_{1}^{{\rm{inv}}}$:

[TABLE]

Substituting back $p=\frac{L}{\mu}\rho$ and $w=\frac{L}{\mu^{2}}\|P(x^{0})\|$ results in the Corollary statement. ∎

On Figure 1 the bound and its residuals are plotted. From the proof it is clear the bounds are tight.

The linear upper (lower) bounds of $F_{1}(\cdot),F_{1}^{{\rm{inv}}}(\cdot)$ for the interval $w\in[0,\frac{1}{2}]$ were chosen for consistency with linear bounds on $w\in[\frac{1}{2},\infty)$. Bound residual on these two intervals are also in the same order, cf. Figure 1. In Corollary 4.8 of Theorem 4.5 we present refined, quadratic approximation for $F_{1}(w)=2H_{0}(\frac{w}{2}),\;w\leq\frac{1}{2}$.

Corollary 4.4**.**

The upper bounds (19) may be simplified as well, using $\overline{w}\leq 1$ and thus $H_{0}(\frac{\overline{w}}{2})\leq H_{0}(\frac{1}{2})=c_{1}$:

[TABLE]

4.2 Adaptive Newton method

Presented Algorithm 1 explicitly uses two constants $\mu$ and $L$ but both enter into the algorithm as one parameter $\beta=\frac{\mu^{2}}{L}$. There is a simple modification allowing adaptively change an estimate of the parameter.

Input of the algorithm is an initial point $x^{0}$, the approximation $\beta_{0}$ and the parameter $0<q<1$.

Algorithm 2 (Adaptive Newton method)

Calculate

[TABLE] 2. 2.

If either

[TABLE]

or

[TABLE]

holds, go to Step 4. Otherwise 3. 3.

apply update rule $\beta_{k}\leftarrow q\beta_{k}$ and return to Step 1 without increasing the counter. 4. 4.

Take

[TABLE]

set $\beta_{k+1}=\beta_{k}$, increase counter $k\leftarrow k+1$, and go to Step 1.

Properties of Algorithm 2 are similar to Algorithm 1. We omit the formal proof of convergence; it follows the lines of the proof of Theorem 4.1 with respect to the properties:

• •

Algorithm 2 runs real steps at Step 4 and some number of fictitious steps resulting in update rule Step 3;

• •

$\beta_{k}$ is non-increasing sequence;

• •

if $\beta_{k}<\beta$ (the actual constant of the objective function), then Step 3 won’t appear and $\beta_{k}$ won’t decrease anymore. It means that there are at most $\widehat{k}=\max\{0,\;\lceil\log_{1/q}(\frac{\beta_{0}}{\beta})\rceil\}$ check steps. Minimal possible value of $\beta_{k}$ is $\beta_{\min}=q^{\widehat{k}}{\beta_{0}}$, and the number of Stage 1 steps is limited by $\widehat{k}_{\max}=\max\{0,\;\left\lceil 2\frac{\|P(x^{0})\|}{\beta_{\min}}\right\rceil-2\}$ as well;

• •

if Step 4 is performed with $\beta_{k}>\beta$ due to validity of a condition in Step 2, then $\|P(x^{k+1})\|$ decreases more than at the corresponding step with “optimal” step-size $\alpha_{k}=\min\{1,\;\frac{\beta}{\|P(x^{k})\|}\}$ (calculated with “true” value $\beta$).

Let us mention two other versions of adaptive Newton method. The first one uses increasing updates (e.g. $\beta_{k+1}=q_{2}\beta_{k}$ with $q_{2}>1$) in the end of Step 4, thus adapting the constant to current $x^{k}$. Also other decrease policies can be applied for $\beta_{k}$ in Step 3.

The alternative to the Algorithm 2 is line-search or Armijo-like rules for choosing step-size $\alpha_{k}$ to minimize objective function $\|P(x^{k}-\alpha z^{k})\|$ directly. It is known that this approach eventually leads to the quadratic convergence rate with pure Newton steps as well, but without any estimates [2]. The difference between is the following: in the proposed Algorithm 2 parameter $\beta$ is being monotonically tuned to the global problem-specific value, rather than trial-and-error procedure is performed at every iteration in the Armijo-like approach. We compare the alternatives in Example 1.

4.3 Method for $L$ known

Constant $\mu$, used in Assumptions $\mathbf{B},\mathbf{B^{\prime}},$ is rarely accessible. As said in the beginning of the section, we can use more accurate upper bound (12) instead of (13) for step-size choice. It results in the algorithm, which uses the Lipschitz constant only. The optimal step-size in this case is

[TABLE]

Algorithm 3 ($L$-Newton method)

[TABLE]

[TABLE]

The algorithm is well-defined, as condition $\|z^{k}\|=0$ holds only if $P(x^{k})=0$, i.e. a solution was found at the previous step. Formally we put $z^{k}=0,\;\alpha_{k}=1$ and $x^{k+1}=x^{k}$ thereafter.

For the Algorithm we also present similar convergence theorem and set of corollaries. We emphasize that while constant $\mu$ is still used in the bounds, Algorithm 3 does not depend on it.

Theorem 4.5**.**

Suppose that Assumptions $\mathbf{A^{\prime}},\mathbf{B^{\prime}}$ hold and

[TABLE]

where $F^{{\rm{inv}}}_{2}(\cdot)$ is the inverse function for the continuous strictly increasing function $F_{1}(w)$ given by

[TABLE]

Then the sequence $\{x^{k}\}$ generated by Algorithm 3 converges to a solution $x^{*}:P(x^{*})=0$. The values $\|P(x^{k})\|$ are monotonically decreasing, at $k$-th step the following estimates for the rate of convergence hold:

[TABLE]

where \overline{w}=\frac{L}{\mu^{2}}\|P(x^{0})\|-\frac{k_{\max}}{2}=\min\big{\{}\frac{L}{\mu^{2}}\|P(x^{0})\|, 1-\frac{1}{2}\big{\lceil}\frac{2L}{\mu^{2}}\|P(x^{0})\|\big{\rceil}+\frac{L}{\mu^{2}}\|P(x^{0})\|\big{\}}\in[0,1).

For proving the theorem we need the simple proposition about real sequences.

Proposition 4.6**.**

Consider two non-negative real sequences $w_{k}\geq 0,v_{k}\geq 0,\;k\geq 0$, and functions $h_{k}(v),f(v),k\geq 0$. Let $f(\cdot)$ be monotonically increasing function, being also a majorant function for $h_{k}(\cdot)$ with respect to $\{w_{k}\}$. Namely, we require $h_{k}(w_{k})\leq f(w_{k})$. If $w_{0}\leq v_{0}$ and the sequences satisfy

[TABLE]

then $w_{k}\leq v_{k},\;k\geq 0$.

The proposition is trivially proved by induction step $w_{k+1}\leq h_{k}(w_{k})\leq f(w_{k})\leq f(v_{k})=v_{k+1}$.

The proof of Theorem 4.5 resembles the proof of Theorem 4.1, and it uses majorization idea. Main issue is due to different step-size, now there is no clear separation between damped and pure Newton steps.

Proof.

We compare two discrete processes, both starting with the same value $v_{0}=w_{0}$. The first sequence is generated by recurrent equality $v_{k+1}=f(v_{k})$, where

[TABLE]

is monotonically increasing function on $v\geq 0$. The second process is $w_{k}=\frac{L}{\mu^{2}}\|P(x^{k})\|$, with $\{x^{k}\}$ generated by Algorithm 3. Assume that all $x^{k}\in B$ and thus Assumptions $\mathbf{A^{\prime}},\mathbf{B^{\prime}}$ hold. Then due to main inequality (12) and step-size (29)

[TABLE]

Here we used $\|z^{k}\|$ in parametric part $a^{2}_{k}\geq 0$ within introduced function $h_{k}(w)=\min_{\alpha}(|1-\alpha|w+a_{k}^{2}\frac{\alpha^{2}}{2})$. From $a_{k}=\frac{L}{\mu}\|z^{k}\|\leq\frac{L}{\mu^{2}}\|P(x^{k})\|=w_{k}$ by Assumption $\mathbf{B^{\prime}}$, the functions $|1-\alpha|w+a^{2}_{k}\frac{\alpha^{2}}{2}$ under minimization in $h_{k}(\cdot)$ are majorized by corresponding functions $|1-\alpha|v+w_{k}^{2}\frac{\alpha^{2}}{2}\leq|1-\alpha|v+\frac{\alpha^{2}}{2}v^{2}$ for $v\geq w_{k}$. Minimums of the functions over $\alpha$ (implicitly dependent on $w_{k},v$) also satisfy $h_{k}(v)\leq f(v)$ whenever $v\geq w_{k}$. It follows that $h_{k}(w_{k})\leq f(w_{k})$.

Therefore sequences $\{w_{k}\}$ and $\{v_{k}\}$, alongside with functions $h_{k}(\cdot),f(\cdot)$ satisfy Proposition 4.6, given $v_{0}=w_{0}=\frac{L}{\mu^{2}}\|P(x^{0})\|$. Now we have upper bound on $w_{k}$ through $v_{k}$ for all $k\geq 0$. Next we closely follow the lines and calculations of the proof of Theorem 4.1.

Analysis of $\{v_{k}\}$ is the same as analysis of the upper bound in Theorem 4.1. First, $v_{k}$ decreases, and the number of steps until $v_{k}$ reaches $1$ is the same $k_{\max}$, given by (18). Explicit expressions on $v_{k}$ are

[TABLE]

where $v_{k_{\max}}=\overline{w}=v_{0}-\frac{k_{\max}}{2}=\min\{w_{0},1-\frac{\lceil 2w_{0}\rceil-2w_{0}}{2}\}$ (remind that $v_{0}=w_{0}$ by definition). These expressions result in the bounds (32a) and (32c) on $\|P(x^{k})\|$, which are the same as in Theorem 4.1.

For the late steps with $k\geq k_{\max}$ we have $w_{k}\leq v_{k}\leq 1$, and thus $\alpha^{*}_{k}=1$ (just because $\frac{\|P(x^{k})\|}{L\|z^{k}\|^{2}}\geq\frac{\mu^{2}}{L\|P(x^{k})\|}=\frac{1}{w_{k}}\geq 1$). Then points $x^{k},k\geq k_{\max}$ form a Cauchy sequence like (23), and obey similar to (24) bound:

[TABLE]

This is (32d) bound, by the way the same as (19d) of Theorem 4.1.

What is the main difference from Theorem 4.1 proof, is the distance counting until $k_{\max}$-th step. In this case we assume $k_{\max}=\lceil 2w_{0}\rceil-2>0$. Due to the Algorithm’s step-size choice now we cannot be sure, whether $\alpha^{*}_{k}$ is always less than $1$ or not. For $i<k_{\max}$ we have $\|x^{i+1}-x^{i}\|=\alpha^{*}_{i}\|z^{i}\|\leq\|z^{i}\|\leq\frac{1}{\mu}\|P(x^{i})\|=\frac{\mu}{L}w_{i}\leq\frac{\mu}{L}v_{i}=\frac{\mu}{L}(v_{0}-\frac{1}{2}i)$, and arrive to (32b):

[TABLE]

In (32b) we used explicit formula for $\overline{w}=1+w_{0}-\frac{\lceil 2w_{0}\rceil}{2}$ in case $k_{\max}>0$.

The last part of the proof is checking assumption $x^{k}\in B$, i.e. $\|x^{0}-x^{k}\|\leq\rho$. From the derivation of bound (32b) above we have

[TABLE]

in case $k_{\max}>0$, i.e. $w_{0}\geq 1$, and $\|x^{0}-x^{k}\|\leq 2\frac{\mu}{L}H_{0}(\frac{w_{0}}{2})$ otherwise (from derivation of bound (32d)). Thus sufficient condition for $x^{k}\in B$ is

[TABLE]

which is equivalent to (30). ∎

Like for Algorithm 1, we can state few corollaries: on global convergence and some simple bounds.

Corollary 4.7**.**

If $\rho=\infty$ (that is conditions $\mathbf{A^{\prime}},\mathbf{B^{\prime}}$ hold on the entire space $\mathbb{R}^{n}$) then Algorithm 3 converges to a solution of (1) for any $x^{0}\in\mathbb{R}^{n}$.

There is a simpler tight relaxed condition for Theorem 4.5.

Corollary 4.8**.**

Condition (30) can be replaced with

[TABLE]

where constant $c_{3}=2c_{1}+c_{2}-1\approx 0.66885$.

The idea of the proof is the same as of Corollary 4.3: we derive upper bound $F_{2}^{{\rm{upper}}}(w)$ for function $F_{2}(w)$ (31). Then the sufficient condition for points $x^{k}$ generated by Algorithm 3 being inside $B$ is $\rho\geq\frac{\mu}{L}F_{2}^{{\rm{upper}}}(\frac{L}{\mu^{2}}\|P(x^{0})\|)$. Its inverse function $F_{2}^{{\rm{inv}},{\rm{lower}}}:F_{2}^{{\rm{inv}},{\rm{lower}}}(F_{2}^{{\rm{upper}}}(w))\equiv w,w\geq 0$ is a lower bound for $F_{2}^{{\rm{inv}}}(\cdot)$ then.

Proof.

Notice that both components of $F_{2}(w)$ are the same in $w\in[\frac{1}{2},1]$, so

[TABLE]

Let’s begin with the case $w>\frac{1}{2}$. We introduce the auxiliary function $r(w)=\frac{2w-\lceil 2w\rceil}{4}\in(0,\frac{1}{4}]$, related with the fractional part. This function is periodic on $w$, and

[TABLE]

Here we used definition (9) of constant $c_{2}$, introduced in Section 3.

By the definition of $H(\cdot)$ one can select terms up to quadratic in $H_{0}(\delta)=\delta+\delta^{2}+H_{2}(\delta)$, thus $H_{2}(\frac{1}{4})=H_{0}(\frac{1}{4})-\frac{5}{16}=H_{0}(\frac{1}{2})-\frac{1}{2}-\frac{5}{16}=c_{1}-\frac{13}{16}$. From convexity we have the upper linear bound for $H_{2}(\delta)\leq(4c_{1}-\frac{13}{4})\delta,\;\;\delta\in[0,\frac{1}{4}]$, and consequently for $2H_{0}(\frac{w}{2})$:

[TABLE]

In the last inequality we manually added a small positive linear term $2c_{2}w$ for continuity of the resulting upper bound. Combining two parts, we arrive to the increasing continuous upper bound for $F_{2}(\cdot)$:

[TABLE]

Using definition of $c_{3}=F_{2}^{{\rm{upper}}}(\frac{1}{2})=2c_{1}+c_{2}-1$, inverse of this function is the lower bound for $F_{2}^{{\rm{inv}}}(\cdot)$

[TABLE]

∎

The presented bounds are sharp and quite exact. Visually paired graphics of $F_{2}(w)$ and $F_{2}^{{\rm{upper}}}(w)$, $F_{2}^{{\rm{inv}}}(p)$ and $F_{2}^{{\rm{inv}},{\rm{lower}}}(p)$ looks the same, and its residuals are plotted on Figure 2.

Surprisingly enough, in practice the Algorithm 3 (and its adaptive modification) sometimes converges faster than Algorithm 1, possibly because direction-wise (along $z^{k}$) Lipschitz constant is less or equal than uniform Lipschitz constant of Assumption $\mathbf{A^{\prime}}$, and the convergence rate can be better.

The idea of adaptive algorithm with estimates $L_{k}$ works as well for Algorithm 3; including its modifications with increasing $L_{k}$.

4.4 Pure Newton method

For comparison let us specify convergence conditions of pure Newton method ($\alpha_{k}=1$).

Theorem 4.9**.**

Let conditions $\mathbf{A^{\prime}},\mathbf{B^{\prime}}$ hold. If $\delta=\frac{L}{2\mu^{2}}\|P(x^{0})\|<1$ and $\frac{2\mu}{L}H_{0}(\delta)\leq\rho$, then pure Newton method converges to a solution $x^{*}$ of (1), and

[TABLE]

It coincides with Corollary 1 of [23], proven in the Banach space setup (a misprint in [23] is corrected here). For $m=n$ case the result is a minor extension of Mysovskikh’s theorem [14].

5 Special Cases

In the section we outline few important cases in more detail, namely solving equations with special structure, solving scalar equations or inequalities, solvability of quadratic equations.

5.1 Structured problems

The problem is to solve equation $g(x)=y$ where $g(x)_{i}=\varphi(c_{i}^{T}x),\;c_{i}\in R^{n},i=1,\dots m$. Here $\varphi(t)$ is a twice differentiable scalar function,

[TABLE]

It is not hard to see that Assumptions $\mathbf{A},\mathbf{B}$ hold on the entire space $\mathbb{R}^{n}$ and Algorithm 1 converges, with Theorem 4.1 and Corollary 4.3 providing rate of convergence. The rate of convergence depends on estimates for $\mu,L$, which can be written as functions of $\mu_{\varphi},L_{\varphi}$ and minimal and maximal singular values $\sigma_{\min},\sigma_{\max}$ of matrix $C$ with columns $c_{i}$ (we suppose that $C$ has full rank, thus $\sigma_{\min}>0$). Indeed after simple calculations (see expression for $g^{\prime}(x)$ with $C$ below) we get

[TABLE]

However the special structure of the problem allows to get much sharper results.

Let’s use notation $P(x)=g(x)-y$. Indeed $P^{\prime}(x)=g^{\prime}(x)=D(x)C^{T},\;D(x)=\,{\rm{diag}}\,(\varphi^{\prime}(c_{i}^{T}x))$ and repeating the proof of Theorem 3.2 we get the equality

[TABLE]

where $D_{t}=D(x^{k}-\alpha tz^{k}),D=D(x^{k})$. Thus (recall $u_{k}=\|P(x^{k})\|$)

[TABLE]

Identity between spectral norm of a diagonal matrix and Euclidean norm of vector on the diagonal is used in the last line, followed by element-wise Lipschitz property of $\varphi^{\prime}(\cdot)$: $\|D_{t}-D\|=\|[\varphi^{\prime}(c_{i}^{T}(x^{k}-\alpha tz^{k}))-\varphi^{\prime}(c_{i}^{T}x^{k})]\|$ \leq\|[\big{|}\varphi^{\prime}(c_{i}^{T}(x^{k}-\alpha tz^{k}))-\varphi^{\prime}(c_{i}^{T}x^{k})\big{|}]\| \leq\|[L_{\varphi}\alpha t\big{|}c_{i}^{T}z^{k}\big{|}]\|=L_{\varphi}\alpha t\|[c_{i}^{T}z^{k}]\|=L_{\varphi}\alpha t\|C^{T}z^{k}\| for $t\geq 0$. Thus $\int_{0}^{1}||D_{t}-D||dt\leq\frac{1}{2}L_{\varphi}\alpha\|C^{T}z^{k}\|$.

But $P^{\prime}(x^{k})z^{k}=P(x^{k})$, thus $DC^{T}z^{k}=P(x^{k}),\;C^{T}z^{k}=D^{-1}P(x^{k})$ and hence $||C^{T}z^{k}||\leq\frac{u_{k}}{\mu_{\varphi}}$. We arrive to the inequality, very similar to (13), but with different constant $\gamma$

[TABLE]

Hence $u_{k+1}\leq u_{k}-\frac{1}{2\gamma}$ at Stage 1, thus this inequality does not depend on $C$! As the result we get estimates for the rate of convergence which are the same for ill-conditioned and well-conditioned matrices $C$. Of course this estimate is much better than standard one with \gamma=\big{(}\frac{\sigma_{\max}}{\sigma_{\min}}\big{)}^{\!2}\frac{L_{\varphi}}{\mu_{\varphi}^{2}} which follows from (34).

This example is just an illustrating one (explicit solution of the problem can be found easily), but it emphasizes the role of special structure in equations to solve. Numerical experiments with such problems are provided below, in Section 6.3.

5.2 One-dimensional case

Suppose we solve one equation with $n$ variables:

[TABLE]

Here [math] is not a minimal value of $f$, thus it is not a minimization problem! Nevertheless our algorithms will remind some minimization methods. This case has some specific features compared with arbitrary $m$. For instance calculation of $z^{k}$ may be done explicitly. Norm in image space is absolute value $|\cdot|$, and $\ell_{p}$ norms in pre-image space $\mathbb{R}^{n},\;p\in\{1,2,\infty\}$ can be considered. Then

[TABLE]

where $e^{j}=(0,\ldots,0,1,0,\ldots,0)^{T}$ is $j$-th orth vector, and $\,{\rm{sign}}\,(\cdot)$ function is coordinate-wise sign function, $\,{\rm{sign}}\,:\mathbb{R}^{n}\rightarrow\{-1,1\}^{n}$.

Constant $\mu$ (and $\mu_{0}$) are also calculated explicitly via conjugate (dual) vector norm $\mu=\min_{x\in B}\|\nabla f(x)\|_{*}$, $\mu_{0}=\|\nabla f(x^{0})\|_{*}$. For any norms $\|z^{k}\|=|f(x^{k})|/\|\nabla f(x^{k})\|_{*}$, and in Algorithm 3 damped Newton step is performed iff $\|\nabla f(x^{k})\|^{2}_{*}<L|f(x^{k})|$, otherwise pure Newton step is made.

If we choose $\ell_{1}$ norm, the method becomes coordinate-wise one. Thus, if we start with $x^{0}=0$ and perform few steps (e.g. we are in the domain of attraction of pure Newton algorithm) we arrive to a sparse solution of the equation.

In Euclidean case a Stage 1 step (damped Newton) of Algorithm 3 is

[TABLE]

which is exactly gradient minimization step for function $|f(x^{k})|$. Stage 2 (pure Newton) step is

[TABLE]

This reminds well-known subgradient method for minimization of convex functions. However in our case we do not assume any convexity properties, and the direction may be either gradient or anti-gradient in contrast with minimization methods!

5.3 Quadratic equations

Proceed to a specific nonlinear equation, namely the quadratic one. Then the function $g(x)$ may be written componentwise as (6), with gradients

[TABLE]

Obviously $g(0)=0$, the question is solvability of $g(x)=y$. There are some results on construction of the entire set of feasible points $Y=\{y:g(x)=y\}=g(\mathbb{R}^{n})$, including its convexity, see e.g. [24]. We focus on local solvability, trying to derive the largest ball inscribed in $Y$.

The derivative matrix $g^{\prime}(x)$ is formed row-wise as

[TABLE]

One has $g^{\prime}(0)=H,\;H$ being $m\times n$ matrix with rows $b_{i}$. We suppose $H$ has rank $m$ (recall that $m\leq n$), then its smallest singular value $\sigma_{\min}(H)>0$ serves as $\mu_{0}$.

The derivative $g^{\prime}(x)$ is linear on $x$, thus it has uniform Lipschitz constant $L$ on $\mathbb{R}^{n}$, and assumption $\mathbf{A}$ holds everywhere. There are several estimates for the Lipschitz constants, for example (for $\ell_{2}$ norm)

[TABLE]

from [25], where $\lambda_{\max}$ is the maximal eigenvalue of a matrix. Other estimates can be obtained via elaborate convex semidefinite optimization problem (SDP), cf. [32] for details.

Quadratic equations play significant role in power system analysis, because power flow equations are quadratic, see [17]. It is of interest to compare our estimates with some known results on solvability of power flow equations [34].

5.4 Solving systems of inequalities

Below we address some tricks to convert systems of inequalities into systems of equations.

First, if one seeks a solution of a system of inequalities

[TABLE]

then by introducing slack variables the problem is reduced to solution of the under-determined system of equations

[TABLE]

Similarly finding a feasible point for linear inequalities $x\geq 0,Ax=b,\;\;x\in\mathbb{R}^{n},\;b\in\mathbb{R}^{m}$ can be transformed to the under-determined system

[TABLE]

The efficiency of such reductions is unclear a priori and should be checked by intensive numerical study.

6 Numerical tests

We have performed several experiments to check effectiveness of the proposed approach for solving equations (1) and to compare it with known ones.

The first two experiments relate to the classical case $n=m$, i.e. the number of variables equals the number of equations. Algorithm 2 with adaptive parameter estimation is compared with well-known “backstepping” Armijo-like techniques for the damped Newton method. Namely, the competitor is the step-size proposed in [2].

[TABLE]

with some shrinkage parameter $q\in(0,1)$ and relaxation parameter $c\in(0,1)$. This method in some sense is similar to our Algorithm 2, but there are differences in step-size choice.

Two other examples relate to under-determined case, i.e. $n>m$. Example 3 is the illustration how to employ structure of the data as explained in Subsection 5.1. We show that such approach strongly accelerates convergence. Final Example 4 is an optimal control problem. Exploiting $L_{1}$ norms we construct sparse controls for the minimal-fuel problem.

6.1 Example 1

We studied the Fletcher-Powell system of equations [8]:

[TABLE]

for various dimensions $n$. The data are generated as proposed in original paper [8]: matrices $A,B\in\mathbb{R}^{n\times n}$ are random, then for some $x^{*}\in\mathbb{R}^{n}$ right-hand sides $E\in\mathbb{R}^{n}$ are calculated. The arising system of equations may have many solutions, however we have the guarantee that it is solvable. Then a set of 1000 initial points $x^{0}$ were randomly sampled (multistart policy). From each of the initial points, Newton algorithms were run with a) Armijo-kind step-size (35), and b) $\beta$-adaptive algorithm (Algorithm 2). Parameters of the algorithms were chosen as $\beta_{0}=100,q=0.95,c=0.8$.

Each run has a “success” or “fail” result. “Success” means that accuracy $||P(x^{k})||<10^{-8}$ is achieved, while “failure” is marked when either step-size threshold is $\gamma_{k}<10^{-13}$ attained or maximal number of iterations ($N=10000$) is performed. For each dimension $n$ data of algorithms’ outcomes were aggregated as following. For a random sample of equations (there were $100$ for each $n$) “success ratio” $r$ was calculated among $1000$ initial points as ratio of success runs and total number of runs: $N_{success}/1000$, both for Armijo-like approach ($r_{Armijo}$), and for our Algorithm 2 ($r_{Alg.~{}2}$). To emphasize comparison, we used ratio of ratios $r_{Alg.~{}2}/r_{Armijo}$ as indicator. Then $r_{Alg.~{}2}/r_{Armijo}$ values were imaged as the box-and-whisker plot for all dimensions (Figure 3). The middle line in a (quartile) box is the median, whiskers’ lengths are set to $0.05$ and $0.95$ percentiles, outlier data is dot-plotted – as soon there are $100$ points (for each of the samples), there are exactly 5 upper and 5 lower outliers).

We see that with high probability the ratio is larger than 1 for large dimensions. As a conclusion, our method finds a solution more often than Armijo-like approach.

The similar analysis was done for the function calls. For each of the samples the numbers of function calls (these can be many in one Newton step) are averaged over initial conditions for all runs, resulting in values $N_{Alg.~{}2},N_{Armijo}$. The ratios of the averaged function calls for Armijo-like step-size (35) and Algorithm 2 ($N_{Armijo}/N_{Alg.~{}2}$) were aggregated on the box-and-whisker plot on Figure 4. Typically Algorithm 2 admits much less functions evaluations compared with Armijo-like algorithm.

6.2 Example 2

The original problem is equality-constrained optimization:

[TABLE]

with scalar differentiable functions $f,h$, and $x\in R^{n}$. By use of Lagrange multiplier rule it is reduced to the solution of equations

[TABLE]

with new variable $X=[x^{T},\nu]^{T}\in R^{n+1}$. Note that we are interested at any solution of these equations, that is we do not distinguish minimum points and stationary points. The derivative is a block matrix

[TABLE]

We address the simplest case: minimization of a quadratic function (with symmetric $A$) on Euclidean unit sphere:

[TABLE]

Let the constraint be defined by $h(x)=\frac{1}{2}(x,x)-\frac{1}{2}$, then

[TABLE]

The experiment were run for different dimensions ($n=20,25,30,35,40,45,50$). For each dimension $100$ problems were randomly generated, and $1000$ initial points were randomly chosen for each problem. Then Algorithm 2 and Armijo-like step-size algorithm (35) were run as in the first example. The parameters and stopping criteria were the same, except for the parameter $q=0.85$, and number of Newton steps were bounded by $1000$. Matrix $A$ of the quadratic objective function is a positive semidefinite matrix, formed as $A=\frac{1}{2}MM^{T}.$ The auxiliary matrix $M\in\mathbb{R}^{n\times n+4}$ has coefficients, uniformly distributed on $[-0.9,2.1]$. Coefficients of the linear term $b_{i}$ are picked up from the scaled Gaussian variables: $b_{i}\sim 0.1\mathcal{N}(0,1)$. The initial conditions were chosen for the extended variable $X=[x^{T},\nu]^{T}$ as following:

• •

the first $n$ components ($x^{0}$, corresponding to the original variable $x$) are sampled from the uniform distribution on Euclidean sphere with radius $1$,

• •

the last, $n+1$ component (Lagrange multiplier $\nu$) is chosen as the “best approximation”, i.e. as the minimizer of the residual $\|P(X)\|^{2}=\|Ax^{0}+b+\nu x^{0}\|^{2}$. The explicit solution depends on the first $n$ components ($x^{0}$),

[TABLE]

Same as in Example 1, success ratios $r_{Armijo},r_{Alg.~{}2}$ and number of function calls $N_{Armijo},N_{Alg.~{}2}$ were calculated (averaged over initial points). On Figure 5 the ratios of $r_{Alg.~{}2}/r_{Armijo}$, and on Figure 6 the ratios $N_{Armijo}/N_{Alg.~{}2}$ were gathered in box-and-whisker plot.

Here it appears, that the success ration of Algorithm 2 prevails over Armijo-like approach only at dimension $50$ (still the ratio is about $1$, meaning that both algorithms behave quite similar); however, the average of function evaluations of our algorithm still lower for all cases.

6.3 Example 3

The problem is described in Section 5.1; it is to solve $g(x)=y$ where $g(x)_{i}=\varphi(c_{i}^{T}x),\;c_{i}\in R^{n},i=1,\dots m$. Here $\varphi(t)$ is twice differentiable scalar function,

[TABLE]

It has been explained in Section 5.1 that the special structure of the problem allows to get much sharper results.

Here we restrict ourselves with the single example to demonstrate how the methods work for medium-size problems ($n=40,m=21$). The equations have special structure as in Section 5.1:

[TABLE]

where

[TABLE]

Matrix $C$ with rows $c_{i}$, vectors $b,y$ were generated randomly. For function $\varphi(t)$ we have $\mu_{\varphi}=\max_{t}\varphi^{\prime}(t)\geq 0.5,\;L_{\varphi}=\max_{t}|\varphi^{\prime\prime}(t)|\leq 2$ for all $t$. Thus if we do not pay attention to the special structure of the problem we have $\mu\geq 0.5\sigma_{\min}(C),\;L\leq 2\sigma_{\max}(C)$. On the other hand if we take into account the structure we can replace $\frac{\mu^{2}}{L}$ ($=0.0012$ in the example) in Algorithm 1 (see Subsection 4.3) with $\frac{\mu_{\varphi}^{2}}{L_{\varphi}}=0.125$.

The results of simulations are as follows. When we apply Algorithm 1 with values $L,\mu$, it requires 6000 iterations to achieve the accuracy $||P(x^{k})||<10^{-12}$, while the same algorithm with $L_{\varphi},\mu_{\varphi}$ requires just 70 iterations. The similar result holds for Algorithm 3: the version exploiting $L$ requires 30 iterations, exploiting ${L_{\varphi}}$ — just 5 iterations. All algorithms which are not based on information on these constants (pure Newton, adaptive Algorithm 2) also converge in 5 iterations.

These results demonstrate how sensitive can be the proposed algorithms to a priori data and to the special structure of equations.

6.4 Example 4.

This test is devoted to the underdetermined systems of equations and, specifically, sparsity property, arising in optimal control problems. The behavior of a pendulum with force control $u$ is given by the second-order differential equation

[TABLE]

Given some initial condition $\phi(0),\dot{\phi}(0)$, the goal is to drive the pendulum to the specified terminal position and angular speed $[\phi(T),\dot{\phi}(T)]^{T}=b\in\mathbb{R}^{2}$ for the fixed time $T$. The secondary goal is to have sparse control and the least control capacity $\int_{0}^{T}|u(t)|dt$.

The model was discretized on the interval $[0,T]$. The discretized control $U$ has dimension $N=T/h-1$, where $h$ is the discretization step. The problem is to solve two equations $[\phi_{d}(T,U),\psi_{d}(T,U)]^{T}=b$, where $\phi_{d},\psi_{d}$ are the discrete counterparts of $\phi,\dot{\phi}$, in $N$ dimensional variable $U$.

The problem was solved by exploiting Algorithm 2 with specific choice of norm, namely - $\ell_{1}$-norm. First, it represents a discretized control capacity ($\|u\|=\sum_{i=0}^{N}|U^{(i)}|$). Second, it is known for its property of finding sparse solution. The initial approximation was $U=0$. Newton method converges in 3 steps, resulting in 5 non-zero components of control (i.e. the control should be applied at 5 time instants only). Moreover, the first Newton step reveals 2 components (time instants), which are sufficient to get to the goal, see Figure 7. Thus 2-impulse control (with impulses at $t=98$ and $t=153$) solves the problem.

Details on the simulation can be found in [27].

7 Conclusions and future research

New solvability conditions for under-determined equations (with wider solvability set) are proposed. The algorithms for finding a solution are easy to implement, they combine weaker assumptions on initial approximations and fast convergence rate. No convexity assumptions are required. The algorithms have large flexibility in using prior information, various norms and problem structure. It is worth mentioning that we do not try to convert the problem into optimization one. Combination of damped/pure Newton method is a contribution for solving classic $n=m$ problems as well.

There are numerous directions for future research.

We suppose that the auxiliary optimization problem for finding direction $z^{k}$ is solved exactly. Of course an approximate solution of the sub-problem suffices. 2. 2.

The algorithms provide a solution of the initial problem which is not specified a priori. Sometimes we are interested in the solution closest to $x^{0}$, i.e. $\min_{P(x)=0}\|x-x^{0}\|$. An algorithm for this purpose is of interest. 3. 3.

More general theory of structured problems (Section 5.1) is needed. 4. 4.

It is not obvious how to introduce regularization techniques into the algorithms.

Acknowledgments

The authors thank Yuri Nesterov and Alexander Ioffe for helpful discussions and references; the comments of the anonymous reviewers are highly acknowledged.

Funding

This work was supported by the Russian Science Foundation under Grant 16-11-10015.