# Survey of Methods for Solving Systems of Nonlinear Equations, Part I: Root-
finding Approaches

Ilias S. Kotsireas
<EMAIL_ADDRESS>
Wilfrid Laurier University
Canada  Panos M. Pardalos
<EMAIL_ADDRESS>
Toxeus Systems LLC
Orlando, Florida, USA  Alexander Semenov
<EMAIL_ADDRESS>
University of Florida
Gainesville, Florida, USA  William T. Trevena
<EMAIL_ADDRESS>
University of Florida
Gainesville, Florida, USA  Michael N. Vrahatis
<EMAIL_ADDRESS>
University of Patras
Patras, Greece

###### Abstract

This paper presents a comprehensive survey of methods which can be utilized to
search for solutions to systems of nonlinear equations (SNEs). Our objectives
with this survey are to synthesize pertinent literature in this field by
presenting a thorough description and analysis of the known methods capable of
finding one or many solutions to SNEs, and to assist interested readers
seeking to identify solution techniques which are well suited for solving the
various classes of SNEs which one may encounter in real world applications.

To accomplish these objectives, we present a multi-part survey. In part one,
we focus on root-finding approaches which can be used to search for solutions
to a SNE without transforming it into an optimization problem. In part two, we
will introduce the various transformations which have been utilized to
transform a SNE into an optimization problem, and we discuss optimization
algorithms which can then be used to search for solutions. In part three, we
will present a robust quantitative comparative analysis of methods capable of
searching for solutions to SNEs.

Keywords: systems of nonlinear equations, localization of zeros, computation
of roots, topological degree, total number of solutions and extrema, interval
methods, symbolic computation, tensor methods, homotopy methods

##  1 Introduction

This paper presents part one of a survey on methods for finding one or many
solutions to a system of nonlinear equations (SNE):

${F}_{m}(x)={\mathit{\Theta}}_{m}\equiv(0,0,\ldots,0)^{\top}\quad\Longleftrightarrow\quad\left\\{\begin{aligned}
&f_{1}(x_{1},x_{2},\ldots,x_{n})=0,\\\\[2.84544pt]
&f_{2}(x_{1},x_{2},\ldots,x_{n})=0,\\\ &\kern 41.25641pt\vdots\\\
&f_{m}(x_{1},x_{2},\ldots,x_{n})=0,\end{aligned}\right.$ (1)

where
${F}_{m}=(f_{1},f_{2},\ldots,f_{m}):{\mathcal{D}}_{n}\subset{\mathbb{R}}^{n}\to{\mathbb{R}}^{m}$,
where $f_{1},f_{2},\ldots,f_{m}$ are real-valued continuous or continuously
differentiable functions on the domain ${\mathcal{D}}_{n}$, and where at least
one of $f_{1},f_{2},\ldots,f_{m}$ is nonlinear. For example, consider the
system of transcendental equations

${F}_{2}(x)={\mathit{\Theta}}_{2}\equiv(0,0)^{\top}\quad\Longleftrightarrow\quad\left\\{\begin{aligned}
f_{1}(x_{1},x_{2})&=x_{1}-x_{1}\sin(x_{1}+5x_{2})-x_{2}\cos(5x_{1}-x_{2})=0,\\\\[2.84544pt]
f_{2}(x_{1},x_{2})&=x_{2}-x_{2}\sin(5x_{1}-3x_{2})+x_{1}\cos(3x_{1}+5x_{2})=0,\end{aligned}\right.$
(2)

which is comprised of two transcendental equations of two unknowns (See Figure
1).

Figure 1: An example of a SNE with two transcendental equations of two
unknowns as introduced by Eq. (2):
(Blue):
$f_{1}(x_{1},x_{2})=x_{1}-x_{1}\sin(x_{1}+5x_{2})-x_{2}\cos(5x_{1}-x_{2})=0$;
(Red):
$f_{2}(x_{1},x_{2})=x_{2}-x_{2}\sin(5x_{1}-3x_{2})+x_{1}\cos(3x_{1}+5x_{2})=0$.
Solutions to this SNE are defined as the points where the blue and red
contours intersect. Finding all of the points within a certain region which
satisfy both equations is a challenging task.

Finding one or more solutions to a SNE is a challenging and ubiquitous task
faced in many fields including chemistry [1, 2, 3], chemical engineering [4],
automotive steering [5], power flow [6, 7], large-scale integrated circuit
designs [8], climate modeling [9], materials engineering [10], robotics [11,
12, 13, 14], nuclear engineering [15], image restoration [16], protein
interaction networks [8], neurophysiology [17], economics [18], finance [19],
applied mathematics [20], physics [21], finding string vacua [22], machine
learning [23, 24], geometric constraint solving (used in computer aided
design) [25], and geodesy [26, 27] among others. The problem of solving even a
system of polynomial equations has been proven to be NP-hard [28].
Furthermore, it has also been proven [29] that no general algorithm exists for
determining whether an integer solution exists for a polynomial equation with
a finite number of unknowns and only integer coefficients. The latter has been
known as Hilbert’s 10th problem.

###  1.1 Notation / Scientific Style

Throughout this paper, we utilize
$x=(x_{1},x_{2},\ldots,x_{n})^{\top}\in\mathcal{D}\subset\mathbb{R}^{n}$ to
denote a real vector within the bounded domain $\mathcal{D}$. Furthermore, we
utilize
$x^{*}=(x_{1}^{*},x_{2}^{*},\ldots,x_{n}^{*})^{\top}\in\mathcal{D}\subset\mathbb{R}^{n}$
to denote a real solution to a SNE such that all equations in the SNE are
satisfied ($F_{m}(x^{*})=0$). In an iterative method, we utilize
$x^{k}=(x_{1}^{k},x_{2}^{k},\ldots,x_{i}^{k},\ldots
x_{n}^{k})^{\top}\in\mathcal{D}\subset\mathbb{R}^{n}$ for $k=0,1,\ldots$ to
denote the vector found during the $k-$th iteration of the iterative method.
Here, $x^{k}_{i}$ denotes the $i-$th coordinate of the vector $x^{k}$.

###  1.2 Terminology

Although we refer to Eq. (1) as a system of nonlinear equations (SNE), such
systems have been referred to in a variety of different ways in literature.
For example, articles [30, 31, 32, 33, 34, 35] utilize the abbreviation “SNLE”
to refer to a system of nonlinear equations, and article [36] uses the
abbreviation “SoNE”. Other papers refer to Eq. (1) as a nonlinear system of
equations, and use the abbreviations “NSE” [37, 38] and “NLS” [39]. Eq. (1)
has also been referred to as a nonlinear equation system (NES) [40, 41, 42,
43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54]. The survey in article [55]
uses the terminology Nonlinear Equations (NEs) to refer to a system of one or
more nonlinear equations.

When $m>n$, a SNE can be referred to as an overdetermined SNE, and when $n>m$,
a SNE can be referred to as an underdetermined SNE. When $m=n$, a SNE can be
referred to as a square SNE [56]. Furthermore, a SNE is considered to be
consistent if a solution exists which satisfies all equations [25].

For a square]/ SNE, a solution
$x^{*}=(x_{1}^{*},x_{2}^{*},\ldots,x_{n}^{*})^{\top}$ of the SNE
${F}_{n}(x)={\mathit{\Theta}}_{n}$ or equivalently a zero $x^{*}$ of the
function ${F}_{n}(x)$ or a root $x^{*}$ of the function ${F}_{n}(x)$ is called
simple if for the determinant of the corresponding Jacobian matrix:

$J_{F_{n}}(x)\equiv
F^{\prime}_{n}(x)_{ij}\equiv\left\\{\frac{\displaystyle\partial
f_{i}(x)}{\displaystyle\partial
x_{j}}\right\\}_{ij}\equiv\left[\begin{array}[]{cccc}\vspace*{0.1cm}{\frac{\displaystyle\partial
f_{1}(x)}{\displaystyle\partial x_{1}}}&{\frac{\displaystyle\partial
f_{1}(x)}{\displaystyle\partial x_{2}}}&\cdots&{\frac{\displaystyle\partial
f_{1}(x)}{\displaystyle\partial x_{n}}}\\\\[5.69046pt]
{\frac{\displaystyle\partial f_{2}(x)}{\displaystyle\partial
x_{1}}}&{\frac{\displaystyle\partial f_{2}(x)}{\displaystyle\partial
x_{2}}}&\cdots&{\frac{\displaystyle\partial f_{2}(x)}{\displaystyle\partial
x_{n}}}\\\ \vdots&\vdots&\ddots&\vdots\\\\[2.84544pt]
{\frac{\displaystyle\partial f_{n}(x)}{\displaystyle\partial
x_{1}}}&{\frac{\displaystyle\partial f_{n}(x)}{\displaystyle\partial
x_{2}}}&\cdots&{\frac{\displaystyle\partial f_{n}(x)}{\displaystyle\partial
x_{n}}}\end{array}\right],\vspace*{0.1cm}$ (3)

at $x^{*}$ it holds that ${\rm det}\,J_{F_{n}}(x^{*})\neq 0$, otherwise it is
called multiple. The problem of conservation and decomposition of a multiple
root into simple roots in the case of systems of homogeneous algebraic
equations has been tackled in [57]. This approach can be applied to high
dimensional CAD where it is sometimes required to compute the intersection of
several hypersurfaces that are a perturbation of a set of original unperturbed
hypersurfaces.

When $F_{n}$ satisfies the monotonicity condition:

$\bigl{(}F_{n}(x)-F_{n}(y)\bigr{)}^{\top}(x-y)\geqslant 0,\kern
8.5359pt\forall\,x,y\in\mathbb{R}^{n},$ (4)

the corresponding SNE can be referred to as a system of monotone nonlinear
equations [58]. Furthermore, $F_{n}$ is considered to be Lipschitz continuous
if there exists $L>0$ such that

$\|F_{n}(x)-F_{n}(y)\|_{2}\leqslant L\|x-y\|_{2},\kern
5.69046pt\forall\,x,y\in\mathbb{R}^{n}.$ (5)

Many of the root finding methods described in Section 3 are guaranteed to
converge to a solution when applied to SNEs that satisfy both the monotonicity
and Lipchitz continuity conditions (the hybrid spectral methods introduced in
[59] for example).

###  1.3 Comparison to other surveys

Other surveys discussing solution techniques for SNEs include [60] and [55].
We have decided to conduct this comprehensive literature review because many
new solution techniques for SNEs have been introduced since the publication of
[60] in 1994, and because the recent survey presented in [55] focuses mainly
on methods which first convert a SNE into an optimization problem, and then
search for multiple solutions to the optimization problem using Intelligent
Optimization Algorithms (IOAs). The IOAs discussed in article [55] are
primarily metaheuristics for global optimization. Although the survey in
article [55] provides a very nice discussion of IOAs for solving SNEs
reformulated as optimization problems, many of the IOAs they discuss are only
introduced at a very high level, only eight IOAs were tested in their
computational study, and the IOAs were evaluated on SNEs comprised of 20
equations or less. Also, article [55] only briefly mentions methods which can
be used to search for solutions to SNEs without transforming them into
optimization problems.

We would like to present a broader survey which covers in detail the large set
of methods which can be used to solve a SNE without transforming it into an
optimization problem (i.e. homotopy and symbolic computation methods). These
methods are our main focus in part one of this survey. In part two, we will
expand upon article [55] by introducing additional reformulation techniques
and optimization algorithms which have been used to solve SNEs, and by
discussing in much more detail many optimization algorithms which were only
briefly introduced in article [55]. This will allow us to appropriately set
the stage for the comprehensive empirical study we will present in part three
of this survey. Furthermore, we believe it is imperative to introduce the
reader to a technique for determining the number of solutions to a SNE that
exist within a bounded domain. Such techniques are of critical practical
importance for those interested in finding all solutions to a SNE that exist
within a domain of interest.

###  1.4 Organization of this survey

We begin this survey by introducing a method which can be used to determine
the total number of solutions to a SNE that exist within a given bounded
domain. By determining the number of solutions to a SNE which exist within a
bounded domain of interest, in the event that no solutions exist, one can
avoid spending time and resources searching for solutions all together.
Alternatively, if one knows that a specific number of solutions exist to a SNE
within a bounded domain of interest, one can continue to search for solutions
until the desired number of solutions are found. Next, we introduce root-
finding methods which have been utilized in literature to search for solutions
to a SNE without transforming it into an optimization problem. We conclude our
paper by introducing additional methods which have been used to attempt to
solve SNEs, and by highlighting promising areas for future research.

##  2 Determining the number of solutions to a SNE in a bounded domain

The knowledge of all the solutions of a system of nonlinear equations and/or
all the extrema of a function is of major importance in various fields. The
total number of the solutions of a system of nonlinear equations can be
obtained by computing the topological degree. Suppose that the function
$F_{n}=(f_{1},f_{2},\ldots,f_{n})\colon{\mathcal{D}}_{n}\subset$
${\mathbb{R}}^{n}$ $\to$ ${\mathbb{R}}^{n}$ is defined and is two times
continuously differentiable in a bounded domain ${\mathcal{D}}_{n}$ of
${\mathbb{R}}^{n}$ with boundary $b({\mathcal{D}}_{n})$. Suppose further that
the solutions of $F_{n}(x)={\mathit{\Theta}}_{n}$ are not located on
$b({\mathcal{D}}_{n})$, and that they are simple (that the determinant of the
Jacobian of $F_{n}$ at these solutions is non-zero). Then the topological
degree of $F_{n}$ at ${\mathit{\Theta}}_{n}$ relative to ${\mathcal{D}}_{n}$
is denoted by ${\rm deg}[F_{n},{\mathcal{D}}_{n},$ ${{\mathit{\Theta}}_{n}}]$
and can be defined by the following relation:

${\rm
deg}[F_{n},{\mathcal{D}}_{n},{{\mathit{\Theta}}_{n}}]=\displaystyle\sum_{{\scriptstyle
x\in F_{n}^{-1}({{\mathit{\Theta}}_{n}})}}{\rm sgn}\,{\rm det}\,J_{F_{n}}(x),$

where ${\rm det}\,J_{F_{n}}(x)$ denotes the determinant of the Jacobian matrix
and ${\rm sgn}$ defines the three-valued sign function. The above definition
can be generalized when $F_{n}$ is only continuous [61].

It is evident that, since ${\rm deg}[F_{n},{\mathcal{D}}_{n},$
${{\mathit{\Theta}}_{n}}]$ is equal to the number of simple solutions of
$F_{n}(x)={\mathit{\Theta}}_{n}$ which give positive determinant of the
Jacobian matrix, minus the number of simple solutions which give negative
determinant of the Jacobian matrix, then the total number $N^{s}$ of simple
solutions of $F_{n}(x)={\mathit{\Theta}}_{n}$ can be obtained by the value of
${\rm deg}[F_{n},{\mathcal{D}}_{n},$ ${\mathit{\Theta}}_{n}]$ if all these
solutions have the same sign of the determinant of the Jacobian matrix. Thus,
Picard considered the following extensions of the function $F_{n}$ and the
domain ${\mathcal{D}}_{n}$ [62, 63]:

$F_{n+1}=(f_{1},f_{2},\ldots,f_{n},f_{n+1})\colon{\mathcal{D}}_{n+1}\subset{\mathbb{R}}^{n+1}\to{\mathbb{R}}^{n+1},$
(6)

where $f_{n+1}=y\,\,{\rm det}\,J_{F_{n}}$,
${\mathbb{R}}^{n+1}:x_{1},x_{2},\ldots,x_{n},y$, and ${\mathcal{D}}_{n+1}$ is
the direct product of the domain ${\mathcal{D}}_{n}$ with an arbitrary
interval of the real $y$-axis containing the point $y=0$. Then the solutions
of the following system of equations:

$\begin{array}[]{l}f_{i}(x_{1},x_{2},\ldots,x_{n})=0,\quad
i=1,2,\ldots,n,\\\\[2.84544pt] y\,\,{\rm
det}\,J_{F_{n}}(x_{1},x_{2},\ldots,x_{n})=0,\end{array}$

are the same simple solutions of $F_{n}(x)={\mathit{\Theta}}_{n}$ provided
that $y=0$. Obviously, the determinant of the Jacobian matrix obtained for the
function (6) is equal to $({\rm det}\,{J_{F_{n}}}(x))^{2}$ which is always
positive. Thus, the total number $N^{s}$ of the solutions of the system
$F_{n}(x)={\mathit{\Theta}}_{n}$ can be obtained by the following value of the
topological degree:

$N^{s}={\rm deg}[F_{n+1},{\mathcal{D}}_{n+1},{\mathit{\Theta}}_{n+1}].$

For example, in the one dimensional case, using the above Picard’s extensions
it is proved that the total number of simple solutions $N^{s}$ of the equation
$f(x)=0$, where $f\colon(a,b)\subset$ ${\mathbb{R}}$ $\to$ ${\mathbb{R}}$ is
twice continuously differentiable in a predetermined interval $(a,b)$, is
given by the following relation [62, 63]:

${N^{s}}=-\frac{1}{\pi}\left[\varepsilon\int_{a}^{b}\frac{f(x)\,f^{\prime\prime}(x)-{f^{\prime}}^{2}(x)}{f^{2}(x)+{\varepsilon}^{2}{f^{\prime}}^{2}(x)}\,dx+\arctan\left(\frac{\varepsilon
f^{\prime}(b)}{f(b)}\right)-\arctan\left(\frac{\varepsilon
f^{\prime}(a)}{f(a)}\right)\right],$ (7)

where $\varepsilon$ is a small positive constant. Note that ${N^{s}}$ was
shown to be independent of the value of $\varepsilon$. Also, the above
approach can be applied for computing the number of multiple solutions.
Obviously, the total number $N^{e}$ of the extrema of $f\in C^{3}$ i.e.
$x\in(a,b)$ such that $f^{\prime}(x)=0$ can be obtained using the above
formula (7). For details of the topological degree we refer the interested
reader to the books [64, 65, 66, 67, 61, 68]. Details of the computation of
the value of the topological degree and its usefulness as well as some
applications and issues related to the number of zeros can be found for
example in [69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83].

Article [84] also discusses utilizing the topological degree to determine the
existence of robust solutions to a SNE.

##  3 Root finding methods

This section discusses root finding methods which have been utilized to search
for solutions to SNEs. Specifically, Section 3.1 \- Section 3.10 discuss
methods which have been utilized to search for a single solution to a SNE, and
Section 3.11 discusses deflation techniques which can be used to compute
further solutions. Some of the root finding methods introduced in this section
are discussed in more detail than others, and after we introduce fundamental
root finding methods, we place a particular emphasis on promising methods
which do not appear to be widely discussed in recent literature. We also
briefly introduce some of the more well known root finding methods such as
Quasi-Newton methods, and we cite relevant literature for the interested
reader to refer to.

Although many of the root finding methods described in this section are
guaranteed to converge to a solution when applied to SNEs that satisfy certain
conditions such as monotonicity (Eq. (4)) and Lipschitz continuity (Eq. (5)),
none of the methods described in this section are guaranteed to converge to a
solution on all classes of SNEs. For example, the root finding method proposed
in article [58] was proven to converge on systems of monotone nonlinear
equations which satisfy additional regularity conditions, but the proposed
method is not guaranteed to converge on all SNEs.

Many of the root finding methods discussed in this section are specialized at
solving particular classes of SNEs. For example, Section 3.5 discusses a
tensor-free Chebyshev-Halley method designed to solve SNEs that are large,
sparse, and which may have ill-conditioned or singular Jacobian matrices at a
solution [85].

###  3.1 Classical Newton’s and Broyden’s methods for finding a single
solution

Newton’s method [86, 87, 61, 88] for finding roots is a well known approach
which can be utilized to find a single solution to a SNE. Starting at a “good”
initial guess or approximate solution
$x^{0}=(x_{1}^{0},x_{2}^{0},\ldots,x_{n}^{0})^{\top}$ of a system of $n$
nonlinear real equations in $n$ real unknowns, when ${F}_{n}(x)$ is Lipschitz
continuous and the Jacobian matrix $J_{{F}_{n}}(x)$ of ${F}_{n}(x)$ is
available, Newton’s method can be used to iteratively find elements of a
sequence converging towards a true solution
$x^{*}=(x_{1}^{*},x_{2}^{*},\ldots,x_{n}^{*})^{\top}$ with a quadratic rate of
convergence. Newton’s method attempts to calculate a new approximate solution
$x^{k+1}$ from $x^{k}$ via

$x^{k+1}=x^{k}-J_{{F}_{n}}\bigl{(}x^{k}\bigr{)}^{-1}\,{F}_{n}\bigl{(}x^{k}\bigr{)},\kern
14.22636ptk=0,1,2,\ldots$

Alternatively, instead of calculating the inverse of the Jacobian matrix, one
often can solve the system of linear equations:

$J_{{F}_{n}}\bigl{(}x^{k}\bigr{)}\,s^{k}=-{F}_{n}\bigl{(}x^{k}\bigr{)},\kern
14.22636ptk=0,1,2,\ldots$ (8)

to find the Newton step $s^{k}=x^{k+1}-x^{k}$. Thus, Newton’s method using a
initial guess $x^{0}$ for each iteration $k=0,1,2,\ldots$ performs the
following steps:

  * (a)

Solve: the system of linear equation
$J_{{F}_{n}}\bigl{(}x^{k}\bigr{)}s^{k}=-F_{n}\bigl{(}x^{k}\bigr{)}$ for
$s^{k}$.

  * (b)

Set: $x^{k+1}=x^{k}+s^{k}$.

In many cases, Step (b) is given by $x^{k+1}=x^{k}+\lambda s^{k}$ where
$\lambda$ is selected to guarantee decrease in $\|F_{n}\|$ (see [87]).

Many expansions upon Newton’s method for finding roots have been introduced in
literature, and most of these expansions are focused on accelerating the rate
of convergence. For example, articles [89, 90, 91, 92, 93, 94] present new
techniques based upon Newton’s method for finding roots and quadrature rules
to accelerate convergence. Newton’s method can also be augmented with higher
order information to accelerate convergence towards an approximation of a true
solution $x^{*}$. Overall, Newton’s method for finding roots seeks to
approximate the solution of a nonlinear system of equations by solving
iteratively a sequence of systems of linear equations. In cases where the
Jacobian matrix is singular, we refer the interested reader to articles [95,
96]. Convergence results for SNEs which have a Jacobian with a constant rank
in the vicinity of a solution are presented in articles [97, 98]. Dimension-
reducing modifications of Newton’s method that are able to tackle efficiently
and effectively, among others, almost linear systems can be found in [99, 100,
101, 102, 103].

In a similar approach, Broyden’s method [104, 86, 87, 61, 88] using an initial
guess $x^{0}$ and an initial matrix $B_{0}$ such that ${\rm det}\,B_{0}\neq
0$, for each iteration $k=0,1,2,\ldots$ performs the following steps:

  * (a)

Solve: the system of linear equations $B_{k}s^{k}=-F_{n}\bigl{(}x^{k}\bigr{)}$
for $s^{k}$.

  * (b)

Set: $x^{k+1}=x^{k}+s^{k}$.

  * (c)

Set: $y^{k}=F_{n}\bigl{(}x^{k+1}\bigr{)}-F_{n}\bigl{(}x^{k}\bigr{)}$.

  * (d)

Set: $\displaystyle
B_{k+1}=B_{k}+\displaystyle\frac{1}{\bigl{(}s^{k}\bigr{)}^{\top}s^{k}}\,\bigl{(}y^{k}-B_{k}s^{k}\bigr{)}\bigl{(}s^{k}\bigr{)}^{\top}.$

A good “choice” of $B_{0}$ is $B_{0}=J_{{F}_{n}}\bigl{(}x^{0}\bigr{)}$.
Obviously, by avoiding this choice Broyden’s method does not require the
computation of the Jacobian matrix. The effectiveness of Broyden’s method and
Newton’s method on different classes of SNEs has been studied, for example,
article [105] compares the effectiveness of Newton’s method, Broyden’s method,
and other methods at solving sparse SNEs.

There are many techniques that can be utilized to solve the system of linear
equations produced by each iteration of Newton’s or Broyden’s method for
finding roots. Two widely used methods for solving a system of linear
equations are described in the following sections, and we refer the interested
reader to [61] for more methods and details regarding these methods. For a
visualization of the basins of convergence for Newton’s and Broyden’s methods,
we refer the interested reader to article [106].

####  3.1.1 Generalized Minimum Residual method (GMRES)

Consider a system of linear equations of the form $Ax=b$. The GMRES method
[107] is an iterative method which at the $k-$th iteration uses the Arnoldi
iteration [108] to find an approximate solution $x^{k}$ in the $k-$th Krylov
subspace $\mathcal{K}_{k}$ which minimizes the residual $\|b-Ax^{k}\|_{2}$.
The $k$-th Krylov subspace of the GMRES method is

$\displaystyle\mathcal{K}_{k}=\mathcal{K}_{k}(A,r^{0})$ $\displaystyle={\rm
span}\\{r^{0},Ar^{0},A^{2}r^{0},\ldots,A^{k-1}r^{0}\\},$

where $r^{0}=b-Ax^{0}$ is the residual at the initial point $x^{0}$.

GMRES assumes that the matrix $A$ is invertible (and hence square) and that
$b$ is normalized such that $\|b\|_{2}=1$. It is worth mentioning that because
$A$ is assumed to be invertible, that this technique in its original form can
not be utilized with Newton’s method to find an approximate solution to a SNE
that has more equations than unknowns $(m>n)$.

####  3.1.2 Successive Overrelaxation method (SOR)

The SOR method is another iterative method for solving a system of linear
equations of the form $Ax=b$. If we assume that the diagonal elements $a_{ii}$
of $A$ are all non-zero, at the $k-$th iteration we can utilize the SOR
iteration to find a new approximate solution $x^{k+1}$ to $Ax=b$ by the
following update:

$x^{k+1}=x^{k}-\omega(D-\omega L)^{-1}(Ax^{k}-b),$

where $\omega>1$ is a relaxation parameter, and $A=D-L-U$ where $D,L$, and $U$
are diagonal, strictly lower triangular, and strictly upper triangular
matrices respectively. Here, $L$ and $U$ have zero diagonal elements, and the
assumption that the diagonal elements of $A$ are nonzero ensures that
$(D-L)^{-1}$ exists [61]. The Gauss-Seidel iteration is a special case of the
SOR where $\omega=1$.

Though the SOR method and the corresponding Jacobi method [61] can be utilized
to solve the series of linear systems produced by Newton’s method, the SOR and
Jacobi methods can also be extended to find a single solution to a SNE.

###  3.2 Solving a SNE using generalizations of iterative methods for linear
systems

The well known and widely used Gauss-Seidel iterative method for solving a
linear system of equations of the form $Ax=b$ can be generalized for solving
SNEs. Thus, if
${F}_{n}=(f_{1},f_{2},\ldots,f_{n}):{\mathcal{D}}\subset{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$,
then the basic step of the nonlinear Gauss-Seidel (NGS) iteration is to solve
(in analogy to linear case), the $i-$th nonlinear one-dimensional equation:

$f_{i}(x_{1}^{k+1},\ldots,x_{i-1}^{k+1},x_{i},x_{i+1}^{k},\ldots,x_{n}^{k})=0,\quad
k=0,1,\ldots,\quad i=1,2,\ldots,n.$ (9)

for $x_{i}$, and to set $x_{i}^{k+1}=x_{i}$. To obtain $x^{k+1}$ from $x^{k}$,
we solve successively the $n$ one-dimensional nonlinear equations (9) for
$i=1,2,\ldots,n$. If relaxation parameters $\omega_{k}$ are used, we may set
$x_{i}^{k+1}=x_{i}^{k}+\omega_{k}(x_{i}-x_{i}^{k})$ and the corresponding
method is called Nonlinear Successive Overrelaxation (NSOR) method. In an
analogous way, the $k-$th step of the Nonlinear Jacobi (NJ) iterative scheme
consists of solving the $i-$th nonlinear one-dimensional equation:

$f_{i}(x_{1}^{k},\ldots,x_{i-1}^{k},x_{i},x_{i+1}^{k},\ldots,x_{n}^{k})=0,\quad
k=0,1,\ldots,\quad i=1,2,\ldots,n.$ (10)

for $x_{i}$ and setting $x_{i}^{k+1}=x_{i}$, for $i=1,2,\ldots,n$. The main
characteristic of the nonlinear Jacobi approach is that it can be easily
parallelized.

In contrast to the linear case, in general, the analytic solutions of Eqs. (9)
and (10) are not available and an one-dimensional rootfinding method must be
applied that terminates after a suitable number of steps. Any type of one-
dimensional methods can be used leading to a large variety of combined
methods. On the other hand if many steps of these one-dimensional methods are
applied the whole procedure becomes cumbersome and thus in practice, in many
cases, a few steps or at least one step of these methods is applied. In this
case, for example, if we apply the one-dimensional Newton’s method, the
derivative-free one-dimensional secant method or the derivative-free one-
dimensional Steffensen method to Eqs. (9) and (10), we obtain respectively the
following methods for solving SNEs [61]:

(a) the one-step SOR-Newton method:

$x_{i}^{k+1}=x_{i}^{k}-\omega_{k}\,{\frac{f_{i}(x^{k,i})}{\partial_{i}f_{i}(x^{k,i})}},\quad
k=0,1,\ldots,\quad i=1,2,\ldots,n.$

where
$x^{k,i}=(x_{1}^{k+1},x_{2}^{k+1},\ldots,x_{i-1}^{k+1},x_{i}^{k},x_{i+1}^{k},\ldots,x_{n}^{k})^{\top}$
and where $\partial_{i}f_{i}(x)$ denotes the partial derivative of the
function $f_{i}(x)$ with respect to the variable $x_{i}$.

(b) the one-step Jacobi Newton method:

$x_{i}^{k+1}=x_{i}^{k}-\omega_{k}\,{\frac{f_{i}(x^{k})}{\partial_{i}f_{i}(x^{k})}},\quad
k=0,1,\ldots,\quad i=1,2,\ldots,n.$

(c) the derivative-free one-step SOR secant method:

$x_{i}^{k+1}=x_{i}^{k}-\omega_{k}{\frac{x_{i}^{k}-x_{i}^{k-1}}{f_{i}(x^{k,i})-f_{i}\bigl{(}x^{k,i}+(x_{i}^{k-1}-x_{i}^{k})\>\\!e^{i}\,\bigr{)}}}f_{i}\bigl{(}x^{k,i}\bigr{)},\quad
k=0,1,\ldots,\quad i=1,2,\ldots,n,$

where $e^{i}$ indicates the $i$-th column of the identity matrix $I_{n}$.

(d) the derivative-free one-step Jacobi secant method:

$x_{i}^{k+1}=x_{i}^{k}-\omega_{k}{\frac{x_{i}^{k}-x_{i}^{k-1}}{f_{i}(x^{k})-f_{i}\bigl{(}x^{k}+(x_{i}^{k-1}-x_{i}^{k})\>\\!e^{i}\,\bigr{)}}}f_{i}(x^{k}),\quad
k=0,1,\ldots,\quad i=1,2,\ldots,n.$

(e) the derivative-free one-step SOR Steffensen method:

$x_{i}^{k+1}=x_{i}^{k}-\omega_{k}{\frac{f_{i}\bigl{(}x^{k,i}\bigr{)}^{2}}{f_{i}(x^{k,i})-f_{i}\bigl{(}x^{k,i}-f_{i}(x^{k,i})\>\\!e^{i}\,\bigr{)}}},\quad
k=0,1,\ldots,\quad i=1,2,\ldots,n.$

(f) the derivative-free one-step Jacobi Steffensen method:

$x_{i}^{k+1}=x_{i}^{k}-\omega_{k}{\frac{f_{i}\bigl{(}x^{k}\bigr{)}^{2}}{f_{i}(x^{k})-f_{i}\bigl{(}x^{k}-f_{i}(x^{k})\>\\!e^{i}\,\bigr{)}}},\quad
k=0,1,\ldots,\quad i=1,2,\ldots,n.$ (11)

Furthermore, since in many problems of practical interest the functions values
are known only imprecisely, the traditional and widely applied one-dimensional
bisection method can be used for a suitable number of steps for computing
solutions of Eqs. (9) and (10). Specifically, for the computation of a zero of
a continuous function $f\colon[a,b]\subset$ ${\mathbb{R}}$ $\to$
${\mathbb{R}}$ the one-dimensional bisection method has be given in [80] by
the following sequence:

$x^{k+1}=x^{k}+{\rm sgn}f(x^{0})\,{\rm sgn}f(x^{k})\,(b-a)/2^{k+1},\qquad
x^{0}=a,\qquad k=0,1,\ldots\,.$ (12)

Similarly, instead of the above sequence we can also use the following one:

$x^{k+1}=x^{k}-{\rm sgn}f(x^{0})\,{\rm sgn}f(x^{k})\,(b-a)/2^{k+1},\qquad
x^{0}=b,\qquad k=0,1,\ldots\,.$ (13)

Obviously, the sequences (12) and (13) converge to a root $r\in(a,b)$ if for
some $x^{k},\,\,k=1,2,\ldots\,$ it holds that ${\rm sgn}f(x^{0})\,{\rm
sgn}f(x^{k})=-1.$ Also, the number of iterations $\nu$, that are required in
obtaining an approximate root $r^{*}$ such that
$|r-r^{*}|\leqslant\varepsilon$ for some $\varepsilon\in(0,1)$ is given by:

$\nu=\bigl{\lceil}{\log}_{2}\bigl{(}(b-a)\,{\varepsilon}^{-1}\bigr{)}\bigr{\rceil}.$
(14)

The one-dimensional bisection method always converges within the given
interval $(a,b)$ and is a globally convergent method. Furthermore, it has a
great advantage since it is worst-case optimal; i.e., it possesses
asymptotically the best possible rate of convergence in the worst case [109,
68]. Thus, it is guaranteed to converge to an approximate root with a
predetermined accuracy using the predefined number of iterations given by Eq.
(14). Notice that no other method has this characteristic. This method
actually requires only the signs of function values (and the gradient values
for the optimization case) to be correct. Thus, it requires only one bit of
information, namely the sign of a function value and consequently it can be
applied to problems with imprecise function values. Also, this method can be
generalized to tackle SNEs and optimization problems. In [80] a
straightforward generalization of the bisection method, named characteristic
bisection method has been presented. For a few details of the bisection and
the generalized bisection methods and for some applications we refer the
interested reader to [110, 111, 71, 74, 112, 79, 113, 80, 114, 82, 115, 116].

The main characteristic of all the above briefly described methods in this
section is that they can tackle $n$-dimensional SNEs using only one-
dimensional rootfinding methods. This issue is very important in cases where
the dimension $n$ is large. Also, the corresponding algorithms are very simple
to implement with a few lines of code. Notice that, the convergence properties
of all the above methods are well studied and analyzed (see for example [61]).
On the other hand, although the nonlinear iterative rootfinding methods have
been extensively studied, the optimization case has not been thoroughly
studied and analyzed. Details and some applications of the above described
methods can be found for example in [117, 118, 119, 114, 120].

###  3.3 Chebyshev-Halley methods

Although Newton’s method only has a quadratic rate of convergence, Halley’s
method possesses a cubic rate of convergence. Halley’s method, also known as
the method of tangent hyperbolas, utilizes information from the tensor of
second derivatives of ${F}_{n}(x)$ to accelerate convergence towards a local
minimizer $x^{*}$ of ${F}_{n}(x)$ from an initial guess $x^{0}$. Halley’s
method attempts to calculate a new approximate solution $x^{k+1}$ to a SNE
from $x^{k}$ using the following scheme:

$x^{k+1}=x^{k}-\Bigl{\\{}I+\frac{1}{2}L_{{F}_{n}}(x^{k})\,\bigl{[}I-\alpha
L_{{F}_{n}}(x^{k})\bigr{]}^{-1}\Bigr{\\}}\,{F}_{n}^{\prime}(x^{k})^{-1}{F}_{n}(x^{k}),$

where $\alpha$ is a real parameter, and
$L_{{F}_{n}}(x)={F}_{n}^{\prime}(x)^{-1}{F}_{n}^{\prime\prime}(x){F}_{n}^{\prime}(x)^{-1}{F}_{n}(x)$
is the degree of logarithmic convexity [121, 122, 123]. Therefore, Halley’s
method is applicable to cases where ${F}_{n}^{\prime\prime}(x)$ can be
computed.

Chebyshev’s method refers to the case where $\alpha=0$, the classical Halley’s
method refers to the case where $\alpha=\frac{1}{2}$, and the super-Halley
method refers to the case where $\alpha=1$. Details about Chebyshev-like
methods for solving SNEs can be found in [124].

###  3.4 Tensor methods utilizing higher order derivatives

There exists a class of tensor methods which utilizes information from higher
order derivatives to accelerate convergence to a single solution of a SNE. In
[125], Steihaug and Suleiman utilize the model:

$\|M^{k}(d)\|_{2}=\Bigl{\|}{F}_{n}(x^{k})+J_{{F}_{n}}(x^{k})d+\frac{1}{2}T^{k}dd\,\Bigr{\|}_{2}\leqslant\eta^{k}\|{F}_{n}(x^{k})\|_{2}.$
(15)

which seeks to determine a step $d^{k}$ at each iteration. In this model,
$T^{k}={F}_{n}^{\prime\prime}(x^{k})$ is the tensor of second derivatives of
${F}_{n}(x^{k})$, and $\eta^{k}\in[0,1)$. Steihaug and Suleiman proved that
any method that can be used to find $d^{k}$ in Eq. (15) while satisfying
$\|d^{k}\|_{2}=O(\|{F}_{n}(x^{k})\|_{2})$ is locally convergent. Furthermore,
Steihaug and Suleiman showed that the rate of convergence is at least:

  * (a)

Q-super-linear when $\eta^{k}\rightarrow 0$.

  * (b)

Q-quadratic when $\eta^{k}=O(\|{F}_{n}(x^{k})\|_{2})$.

  * (c)

Q-cubic when $\eta^{k}=O(\|{F}_{n}(x^{k})\|_{2}^{2})$.

  * (d)

Q-order $\min\bigl{\\{}\hat{p},\,3\bigr{\\}}$ when
$\eta^{k}=O(\|{F}_{n}(x^{k})\|_{2}^{\hat{p}-1}),\ \ 1<\hat{p}$.

Furthermore, in order to solve Eq. (15), Steihaug and Suleiman in [125]
introduce a class of inexact Chebyshev-Halley methods which under some
assumptions are locally convergent satisfying Eq. (15) and
$\|d^{k}\|_{2}=O(\|{F}_{n}(x^{k})\|_{2})$.

###  3.5 Tensor-free Chebyshev-Halley method

Instead of calculating the tensor term $T^{k}dd$ as performed in Eq. (15),
Eustaquio et al. in [85] avoid calculating the tensor term altogether by
introducing a general tensor-free Chebyshev-Halley method for solving SNEs.
Therefore, though Eq. (15) requires calculating ${F}_{n}^{\prime\prime}(x)$,
the method of Eustaquio et al. can be utilized in cases where
${F}_{n}^{\prime\prime}(x)$ does not exist. Also, this method does not require
the inequality in Eq. (15) to be satisfied. The general framework for the
inexact tensor-free Chebyshev-Halley class is the following:

  1. 1.

Given the SNE ${F}_{n}(x)={\mathit{\Theta}}_{n}$, the machine precision
$\varepsilon_{M}$, and $k=0$, select $\alpha\in{\mathbb{R}}$, a step length
$h>0$, an initial point $x^{0}\in{\mathbb{R}}^{n}$, a forcing term limit
$\tilde{\eta}\in(0,1)$, and a mapping
$C:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}^{n\times n}$.

  2. 2.

while  ${F}_{n}(x^{k})\neq{\mathit{\Theta}}_{n}$ do

    1. (a)

Select forcing term tolerances $\eta_{1}^{k}\in[0,\tilde{\eta})$ and
$\eta_{2}^{k}\in[0,\tilde{\eta})$.

    2. (b)

Compute $d_{1}$ such that
$\bigl{\|}J_{{F}_{n}}(x^{k})d_{1}+{F}_{n}(x^{k})\bigr{\|}_{2}\leqslant\eta_{1}^{k}\bigl{\|}{F}_{n}(x^{k})\bigr{\|}_{2}$.

    3. (c)

Compute $d_{2}$ such that
$\bigl{\|}\bigl{(}J_{{F}_{n}}(x^{k})+\alpha\,C(x^{k})\bigr{)}d_{2}+\frac{1}{2}C(x^{k})d_{1}\bigr{\|}_{2}\leqslant\eta_{2}^{k}\bigl{\|}\frac{1}{2}C(x^{k})d_{1}\bigr{\|}_{2}$.

    4. (d)

Set $x^{k+1}=x^{k}+d_{1}+d_{2}$.

    5. (e)

Set $k=k+1$.

  3. 3.

end while

where the authors utilized the mapping
$C(x^{k})=\frac{1}{h}\bigl{(}J_{{F}_{n}}(x^{k}+hd_{1})-J_{{F}_{n}}(x^{k})\bigr{)}$
with step length
$h=\sqrt{\varepsilon_{M}}\,\|x^{k}\|_{2}\big{/}\|d_{1}\|_{2}\in(\varepsilon_{M},0.5)$
and the forcing terms
$\eta_{1}^{k}=\min\bigl{\\{}10^{-8},\,\|{F}_{n}(x^{k})\|_{2}^{2}\bigr{\\}}\in(2\varepsilon_{M},10^{-8}]$
and
$\eta_{2}^{k}=\min\bigl{\\{}10^{-8},\,\|{F}_{n}(x^{k})\|_{2}\bigr{\\}}\in(2\varepsilon_{M},10^{-8}]$
to control the level of accuracy of the approximate solution to the SNE. In
comparison to Eq. (15), the mapping $C(x^{k})$ replaces the tensor term
$T^{k}d$ while still preserving a third-order rate of convergence. As shown in
the algorithm above, each iteration only requires approximately solving two
linear systems. Eustaquio et al. utilized Saad and Schultz’s Generalized
Minimum Residual method (GMRES) [107] to solve the two linear systems and
obtain the inexact Newton steps $d_{1}^{k}$ and $d_{2}^{k}$.

The authors proved that any method belonging to their class of inexact tensor-
free Chebyshev-Halley methods is locally convergent. Note that the system
solved in step 2(b) is analogous to system (8) solved in Newton’s method where
$d_{1}=x^{k+1}-x^{k}$.

Eustaquio et al. conducted a rigorous quantitative comparative analysis of
their method against the Inexact Newton method and the tensor methods
utilizing higher order derivatives presented by Steihaug and Suleiman [125].
Eustaquio et al. illustrated that the class of tensor-free Chebyshev-Halley
methods are capable of efficiently handling problems that have ill-conditioned
or singular Jacobian matrices at the solution.

Furthermore, Jarrat’s method [126] is defined as a specific case of the
tensor-free Chebyshev-Halley method where $\alpha=1$, $\eta_{1}^{k}=0$,
$\eta_{2}^{k}=0$, $h={2}/{3}$, and the mapping
$C(x^{k})={h^{-1}}\bigl{(}J_{{F}_{n}}(x^{k}+h\,d_{1})-J_{{F}_{n}}(x^{k})\bigr{)}$.
Jarrat’s method has a convergence rate of four. Articles using Jarratt-like
methods to solve SNEs can be found for example in [30, 127, 128].

###  3.6 Tensor methods not utilizing higher order derivatives

In addition to the class of tensor methods which utilize a tensor of higher
order derivatives, other tensor methods exist and have demonstrated
effectiveness at solving large, sparse, and ill-formed SNEs with singular
Jacobian matrices [129, 130, 131, 132, 133, 134, 135, 136, 137]. In [129]
Frank and Schnabel introduced tensor strategies for solving SNEs based around
solving the quadratic model:

$\min_{d\in{\mathbb{R}}^{n}}\|M^{k}(d)\|_{2}=\min_{d\in{\mathbb{R}}^{n}}\Bigl{\|}{F}_{n}(x^{k})+J_{{F}_{n}}(x^{k})d+\frac{1}{2}T^{k}dd\,\Bigr{\|}_{2},$
(16)

where $J_{{F}_{n}}(x)$ is the Jacobian matrix of ${F}_{n}(x)$, $d$ is the step
size, and $T^{k}\in{\mathbb{R}}^{n\times n\times n}$ is a carefully chosen
tensor.

Expanding upon this formulation, in [137], Bader presents three Krylov-based
methods for iteratively solving Eq, (16) to a specified tolerance. Bader’s
method selects $T^{k}$ such that the model interpolates $p\leqslant\sqrt{n}$
function values from its most recent history of iterates. By selecting $p=1$,
Bader reduces the tensor model about $x^{k}$ to

$\min_{d\in\mathcal{K}_{m}}\|M^{k}(d)\|_{2}=\min_{d\in\mathcal{K}_{m}}\Bigl{\|}{F}_{n}(x^{k})+J_{{F}_{n}}(x^{k})d+\frac{1}{2}a^{k}\bigl{(}(s^{k})^{\top}d\bigr{)}^{2}\Bigr{\|}_{2},$
(17)

where ${\mathcal{K}}_{m}$ is a $m$-dimensional Krylov subspace, and where
$a^{k}$, $s^{k}\in{\mathbb{R}}^{n}$ are given as follows:

$\displaystyle a^{k}$
$\displaystyle=\frac{2\bigl{(}{F}_{n}(x^{k-1})-{F}_{n}(x^{k})-J_{{F}_{n}}(x^{k})s^{k})}{((s^{k})^{\top}s^{k}\bigr{)}^{2}},$
$\displaystyle s^{k}$ $\displaystyle=x^{k-1}-x^{k}.$

The linear Krylov subspace method finds an approximate solution $x_{m}$ to the
linear system $Ax=b$ from an $m$-dimensional affine subspace
$x^{0}+\mathcal{K}_{m}$ where

$\displaystyle\mathcal{K}_{m}(A,r^{0})$ $\displaystyle={\rm
span}\\{r^{0},Ar^{0},A^{2}r^{0},\ldots,A^{m-1}r^{0}\\},$ $\displaystyle r^{0}$
$\displaystyle=b-Ax^{0},$

and $r^{0}$ is the residual at the initial point $x^{0}$.

The three techniques presented by Bader in [137] select $\mathcal{K}_{m}$ in
different ways as a subroutine within the following algorithm:

  1. 1.

Given the SNE ${F}_{n}(x)$, select an initial point $x_{0}$ and a maximum
number of iterations $k_{max}$.

  2. 2.

For $k=0,1,2,\ldots,k_{max}$, do:

    1. (a)

Select a forcing term tolerance $\eta^{k}\in[0,1)$.

    2. (b)

If $k=0$:

      1. i.

Calculate the Newton-GMRES [107] step $d_{N}$ based on the tolerance
$\eta^{k}$.

      2. ii.

Proceed to step 2e.

    3. (c)

Form the local tensor model (16).

    4. (d)

Calculate the approximate tensor step $d_{T}$ according to $\eta^{k}$ by
solving one of the three methods presented for selecting $\mathcal{K}_{m}$.

    5. (e)

Set $x^{k+1}=x^{k}+\alpha d$ where a linesearch strategy using the directions
$d_{T}$ and / or $d_{N}$ is used to select $d$ and $\alpha$.

    6. (f)

If $x^{k+1}$ is an acceptable approximate root of ${F}_{n}(x)$:

      1. i.

Stop.

Bader’s methods demonstrated effectiveness at solving large-scale SNEs,
especially those which have ill-conditioned or singular Jacobians at the
solution.

###  3.7 Quasi-Newton method for finding roots

Although Newton’s method for finding roots requires the Jacobian
$J_{{F}_{n}}(x)$, if $J_{{F}_{n}}(x)$ is not available, Quasi-Newton methods
can be utilized. A Quasi-Newton method is any method that utilizes an
approximation of $J_{{F}_{n}}(x)$ instead of the exact $J_{{F}_{n}}(x)$. A
classification of Quasi-Newton methods for solving SNEs is presented in
article [138]. Examples of Quasi-Newton methods include the methods discussed
in articles [139, 140].

One of the most popular Quasi-Newton methods is the
Broyden–Fletcher–Goldfarb–Shanno (BFGS) method [141, 142, 143, 144]. Popular
extension of BFGS method is limited memory BFGS (L-BFGS) [145]. Article [58]
combines a memoryless variant of the BFGS method with a projection technique
for solving SNEs, and the authors demonstrated the efficiency of the method by
solving large SNEs with up to 5000 dimensions in the order of milliseconds or
seconds, depending on the problem instance. Other Quasi-Newton methods include
the Symmetric Rank 1 (SR1) method and the Davidon–Fletcher–Powell (DFP) method
among others. Quasi-Newton methods were utilized to find solutions to SNEs in
articles [146, 147, 148, 149, 150, 151, 152, 153, 154].

###  3.8 Spectral methods

Spectral methods are another class of derivative-free methods which have been
used to solve SNEs, and they are particularly well suited for solving large
SNEs because they have a relatively low storage requirement [59]. Spectral
methods often utilize spectral parameters and coefficients to help determine
the search direction $d^{k}$ in the update step
$x^{k+1}=x^{k}+\alpha^{k}d^{k}$. For example, article [59] presents two
derivative-free hybrid spectral methods for solving SNEs, and proves that the
proposed methods will converge to a solution when applied to SNEs that satisfy
the conditions of monotonicity (Eq. (4)) and Lipschitz continuity (Eq. (5)).
Other articles which propose applying spectral methods for solving SNEs
include [155, 156].

###  3.9 Hermitian and skew-Hermitian splitting (HSS) based methods

First introduced in [157] for solving linear systems of equations, Hermitian
and skew-Hermitian splitting (HSS) based methods have been shown to be
effective at solving sparse SNEs [158] and SNEs which can be decomposed into a
linear part $Ax$ and a nonlinear part $\zeta(x)$ such that
$F_{m}(x)=Ax-\zeta(x)$ [159]. When the linear part $Ax$ is dominant over the
nonlinear part $\zeta(x)$, $F_{m}(x)$ can be referred to as a weakly nonlinear
system. By separating the linear and nonlinear elements of $F_{m}(x)$, certain
specialized techniques can potentially be utilized to find solutions to
$F_{m}(x)$ more quickly than by utilizing traditional root finding methods.
Some of these techniques are derivative-free, including the derivative-free
HSS based method introduced in [160] which is guaranteed to converge to a
solution on SNEs that satisfy certain conditions.

###  3.10 Levenberg-Marquardt method

The Levenberg-Marquardt algorithm [161, 162] was designed to solve a nonlinear
least squares problem that can be in general expressed as
$\varphi(x)=\frac{1}{2}\sum_{i=1}^{m}[f_{i}(x)]^{2}$ where
$\varphi:\mathbb{R}^{n}\rightarrow\mathbb{R}$. The Levenberg-Marquardt
algorithm is an iterative procedure which at each iteration calculates an
updated solution using the rule $x^{k+1}=x^{k}+\lambda d$, where the search
direction $d$ is found by solving equations of the form

$\bigl{(}{J}^{\top}_{{F}_{m}}{J}_{{F}_{m}}+\mu\,I_{n}\bigr{)}d=-{J}^{\top}_{{F}_{m}}{F}_{m},$

where ${J}_{{F}_{m}}\in\mathbb{R}^{m\times n}$ is the Jacobian matrix of
${F}_{m}$, $\mu\in\mathbb{R}$, $I_{n}$ is the identity matrix, and
${F}_{m}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$. In article [163], Fletcher
proposed the following modification to make the solution scale invariant:

$\bigl{(}{J}^{\top}_{{F}_{m}}{J}_{{F}_{m}}+\mu\,{\rm
diag}\\{{J}^{\top}_{{F}_{m}}{J}_{{F}_{m}}\\}\bigr{)}d=-{J}^{\top}_{{F}_{m}}{F}_{m},$

where the identity matrix $I_{n}$ is replaced with the diagonal matrix
consisting of the diagonal elements of ${J}^{\top}_{{F}_{m}}{J}_{{F}_{m}}$.
Article [164] proposed to use the Levenberg-Marquardt method for solving a SNE
modeling a port mooring structure. The Levenberg-Marquardt algorithm is widely
used in many applications, including artificial neural network training (for
example, see article [165]).

###  3.11 Deflation techniques for the computation of further solutions

As illustrated above, there is a plethora of methods for obtaining a single
solution of a system of $n$ nonlinear real equations in $n$ real unknowns.
Brown and Gearhart in [166] proposed deflation techniques for the computation
of further solutions of a system of nonlinear equations. Specifically, these
techniques proceed as follows: “Once a solution of a system of nonlinear
equations has been obtained a modified system is formed in such a way that it
retains those solutions of the original system which remain to be computed
except the solutions that has been already computed”. This procedure may be
applied sequentially until all solutions of the original system are obtained.

For example, assume a system of nonlinear equations
${F}_{n}(x)={\mathit{\Theta}}_{n}$, where
${F}_{n}=(f_{1},f_{2},\ldots,f_{n}):{\mathcal{D}}_{n}\subset{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$,
then to deflate out the $p$ already computed roots $r_{1}$,
$r_{2},\ldots,r_{p}$ in order to compute additional roots, the following norm
deflated function [166] is applied:

$\hat{f}_{i}(x)=\frac{1}{\prod_{j=1}^{p}\,\|x-r_{j}\|}\,{f}_{i}(x),\kern
14.22636pti=1,2,\ldots,n,$

or alternatively the following inner product deflated function [166] can be
used:

$\tilde{f}_{i}(x)=\frac{1}{\prod_{j=1}^{p}\,\bigl{\langle}\nabla
f_{i}(r_{j}),\,(x-r_{j})\bigr{\rangle}}\,{f}_{i}(x),\kern
14.22636pti=1,2,\ldots,n.$

The above inner product deflated function has proven to be useful in practice
when Newton’s method for root finding is used. The deflation techniques can
also be extended for the computation of multiple solutions. A study of the
deflation techniques for the one dimensional case can be found in article
[167]. Also, additional techniques for the optimization case related to the
above deflation techniques as well as some applications can be found for
example in the papers [168, 169, 170, 171, 172, 173].

##  4 Symbolic computation methods

Exact algorithmic methods for solving systems of nonlinear polynomial
equations have been developed in the realm of Symbolic Computation, also
called Computational Algebra, or Computer Algebra. Broadly speaking, this
particular subarea of Symbolic Computation can in fact be interpreted as a
constructive version of some parts of the well-established mathematical
theories of Commutative Algebra and Algebraic Geometry. We begin this section
by a brief presentation (via examples) of the main Symbolic Computation
methods for solving systems of nonlinear polynomial equations, namely
resultants and Gröbner bases. We end this section by citing a number of books
that discuss extensively all the developments in the area of Symbolic
Computation exact methods for solving systems of nonlinear polynomial
equations.

####  4.0.1 Resultants

Historically, the theory of resultants was developed with the aim to provide a
systematic means of elimination of variables/unknowns from a system of
nonlinear polynomial equations.

The resultant of two univariate polynomials of degrees $m$, $n$ respectively,

$p(x)=p_{m}x^{m}+\cdots+p_{1}+p_{0},\kern
22.76228ptq(x)=q_{n}x^{n}+\cdots+q_{1}+q_{0},$

is defined as the determinant of the $(n+m)\times(n+m)$ Sylvester matrix
associated to $p(x)$ and $q(x)$, namely:

${\rm res}_{x}\bigl{(}p(x),q(x)\bigr{)}={\rm
det}\left[\begin{array}[]{llll}a_{m}&a_{m-1}&&\\\ 0&a_{m}&a_{m-1}&\\\
0&0&a_{m}&a_{m-1}\\\ 0&0&0&a_{m}\\\ b_{n}&b_{n-1}&&\\\ 0&b_{n}&b_{n-1}&\\\
0&0&b_{n}&b_{n-1}\\\ 0&0&0&b_{n}\\\ \end{array}\right].$

One of the main properties of the resultant is that it is equal to zero, if
and only if the two polynomials have a common root.

The resultant of two bivariate polynomials $p=p(x,y)$,$q=q(x,y)$ with respect
to the variable $x$, or with respect to the variable $y$, is the determinant
of their associated Sylvester matrix, when $p$ and $q$ are considered as
polynomials in $x$, or in $y$ respectively.

Let us illustrate the concept of the resultant by a simple but instructive
example.

###### Example 1

Consider the system of two polynomials
$f_{1}(x_{1},x_{2})=0,f_{2}(x_{1},x_{1})=0$ in two variables $x_{1},x_{2}$,
given by:

$\begin{array}[]{l}f_{1}(x_{1},x_{2})=x_{1}^{2}+x_{2}^{2}-1,\\\\[2.84544pt]
f_{2}(x_{1},x_{2})=x_{1}-x_{2}-1.\\\ \end{array}$ (18)

The resultant of $f_{1},f_{2}$, with respect to $x_{1}$ is the determinant of
their $3\times 3$ Sylvester matrix when considered as polynomials in $x_{1}$,
namely, $f_{1}=1\cdot x_{1}^{2}+0\cdot x_{1}+(x_{2}^{2}-1)$, $f_{2}=1\cdot
x_{1}+(-x_{2}-1)$,

${\rm res}_{x_{1}}(f_{1},f_{2})=\det\left[\begin{array}[]{ccc}1&0&{{\it
x_{2}}}^{2}-1\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 1&-{\it x_{2}}-1&0\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&1&-{\it
x_{2}}-1\end{array}\right]=2x_{2}^{2}+2x_{2}.$

By factorizing the resultant and setting it equal to zero, we obtain $x_{2}=0$
and $x_{2}=-1$. By substituting each of these two values of $x_{2}$ back into
(18), we obtain the corresponding values of $x_{1}$ as: $x_{1}=1$ and
$x_{1}=0$. Therefore, we conclude that the system (18) possesses the two
solutions:

$(x_{1}=1,\,\,x_{2}=0),\kern 5.69046pt(x_{1}=0,\,\,x_{2}=-1),$

which can easily be verified to be correct. In terms of geometric
interpretation, the first equation $f_{1}(x_{1},x_{2})$ represent the unit
circle and the second equation $f_{1}(x_{1},x_{2})$ represents a straight
line, that intersects the unit circle. The two solutions of the system (18)
are the two points of intersection of the (red) unit circle and this (green)
straight line in Figure 2 below.

Figure 2: Geometric interpretation of System (18)

####  4.0.2 Gröbner bases

The theory of Gröbner bases was developed as a multivariate analogue of the
classical Gaussian elimination for systems of linear equations. Again we
illustrate the concept with an example and refer to the aforementioned books
for additional details and technicalities.

###### Example 2

Consider the system of three polynomials $f_{1}(x_{1},x_{2},x_{3})=0$,
$f_{2}(x_{1},x_{2},x_{3})=0$, $f_{3}(x_{1},x_{2},x_{3})=0$ in three variables
$x_{1},x_{2},x_{3}$, given by:

$\begin{array}[]{l}f_{1}(x_{1},x_{2},x_{3})=x_{1}^{2}+x_{2}x_{3}-2,\\\\[2.84544pt]
f_{2}(x_{1},x_{2},x_{3})=x_{1}x_{3}+x_{2}^{2}-3,\\\\[2.84544pt]
f_{3}(x_{1},x_{2},x_{3})=x_{1}x_{2}+x_{3}^{2}-5.\\\ \end{array}$ (19)

The lexicographical Gröbner basis of the ideal generated by these three
polynomials, with respect to the lexicographical ordering induced by
$x_{2}>x_{1}>x_{3}$ is given by the three polynomials:

$\begin{array}[]{l}8x_{2}^{8}-60x_{2}^{6}+142x_{2}^{4}-172x_{2}^{2}+1,\\\\[2.84544pt]
88x_{2}^{7}-680x_{2}^{5}+1674x_{2}^{3}-2081x_{2}+117x_{1},\\\\[2.84544pt]
152x_{2}^{7}-1132x_{2}^{5}+2700x_{2}^{3}-3403x_{2}+117x_{3}.\\\ \end{array}$
(20)

Now we notice that the first polynomial in (20) depends only on $x_{2}$ and is
of degree $8$. In addition, the second polynomial in (20) depends linearly on
$x_{1}$ and the third polynomial in (20) depends linearly on $x_{3}$. This
allows us to express $x_{1}$ and $x_{3}$ as polynomials in $x_{2}$. Therefore,
the solution process starts by finding the 8 roots of the first polynomial in
(20) and for each one of those roots, we find the unique values of $x_{1}$ and
$x_{3}$, given by the last two polynomials in (20). Here are the eight complex
roots (four real roots and two pairs of complex conjugate roots) of the first
polynomial in (20)

$\begin{array}[]{l}-2.16624183202470499\\\\[2.84544pt] \kern
6.25963pt2.16624183202470499\\\\[2.84544pt]
-0.07643337499454630\\\\[2.84544pt] \kern
6.25963pt0.07643337499454630\\\\[2.84544pt]
-1.329681781358029-0.606033421098925\,\mathrm{i}\\\\[2.84544pt]
-1.329681781358029+0.606033421098925\,\mathrm{i}\\\\[2.84544pt] \kern
6.25963pt1.329681781358029-0.606033421098925\,\mathrm{i}\\\\[2.84544pt] \kern
6.25963pt1.329681781358029+0.606033421098925\,\mathrm{i}\\\ \end{array}$

Note that the four real roots and the four complex roots come in pairs of the
form $(r,-r)$, a consequence of the fact that the first polynomial in (20)
contains only even powers of $x_{2}$. The presence of only even powers of the
indeterminate $x_{2}$ in the first polynomial in (20) is captured by the fact
that the order of the Galois group of this polynomial is smaller than the
order of the Galois group of a more “random/generic” polynomial of degree $8$.

A very readable, self-contained and pedagogical introduction to Gröbner bases
can be found in chapter 10 of [174].

####  4.0.3 Symbolic computation software systems

One of the major outcomes of the research area of Symbolic Computation
(Computer Algebra) in the past few decades, is the advent of the so-called
Computer Algebra Systems (CAS). A CAS is defined as a piece of mathematical
software that makes advanced functionalities available to the user, in a
transparent manner. Such advanced functionalities include univariate and
multivariate polynomial factorization, primality testing, integer
factorization, root finding, numerical and exact integration, visualization,
number theory, linear algebra, commutative algebra, commutative and non-
commutative settings, tools for optimization, graph theory, group theory,
coding theory, combinatorics, discrete mathematics and so forth. There is a
number of commercial and free (open source) CAS available today. Some of the
most well-known early CAS (such as Axiom and Macsyma) are largely deprecated
today. We focus our attention on those CAS that currently feature a
significant user base. We note that CAS are used by millions of research
mathematicians, physicists, chemists, astronomers, engineers, practitioners
and educators at academic, public, private and government institutions
worldwide. In what follows, we provide a brief summary of the main commercial
CAS in use today.

  * (a)

Maple is the flagship Canadian product in the area of mathematical software.
It is produced, maintained and distributed by the company Maplesoft, based in
Waterloo, Ontario, Canada https://www.maplesoft.com/. Maple features an easy-
to-learn underlying programming language as well as more than $150$ additional
packages that significantly expand its core functionalities.

  * (b)

Magma is the flagship Australian product in the area of mathematical software.
It is produced, maintained and distributed by the University of Sydney
http://magma.maths.usyd.edu.au/. It provides a mathematically rigorous
environment for defining and working with structures such as groups, rings,
fields, modules, algebras, schemes, curves, graphs, designs, codes and many
others. Magma also supports a number of databases designed to aid
computational research in those areas of mathematics which are algebraic in
nature.

  * (c)

Mathematica is a flagship American product in the area of mathematical
software. It is produced, maintained and distributed by the company Wolfram
Research, based in Champaign, Illinois, United States
https://www.wolfram.com/. Mathematica’s stated purposes include injecting
computational intelligence at every level, on every project by unifying
algorithms, data, notebooks, linguistics and deployment—enabling powerful
workflows across desktop, cloud, server and mobile.

  * (d)

Matlab is a flagship American product in the area of mathematical software.
https://www.mathworks.com. It is produced, maintained and distributed by the
company MathWorks, based in California, United States. MATLAB’s strong points
include matrix manipulations, plotting of functions and data, implementation
of algorithms, creation of user interfaces, and interfacing with programs
written in other languages.

Some of the most popular non-commercial mathematical software packages include
Sage, Singular/Plural, CoCoA, Macaulay 2 and so forth. Each one has their own
strengths and areas of particular focuses. For more information on CAS in
general, we refer the interested reader to the comprehensive book [175].

####  4.0.4 Symbolic computation books

In this section, we present a selection of books in Symbolic computation and
related areas. While there is inevitably some overlap among these books, there
are also different aspects of polynomial system solving and its applications
discussed in each one of them.

  1. 1.

The classic books [174, 176, 177] are general references for Computer Algebra.

  2. 2.

The series of books by Teo Mora [178, 179, 180, 181] is an encyclopedic
reference to the subject of polynomial system solving.

  3. 3.

The theory of Gröbner bases, see [182, 183, 184, 185], was initially developed
by Bruno Buchberger and further extended by several other researchers. The
Buchberger algorithm to compute Gröbner bases of polynomial systems is
implemented in every major Symbolic Computation software today. This includes
commercial software, such as Maple, Magma and Mathematica and open source
software, such as Singular [186], Plural [187], CoCoA [188, 189, 190, 191] and
the Macaulay2 software system [192].

  4. 4.

The book [193] discusses aspects of systems of polynomial equations from the
viewpoints of computational commutative algebra, discrete geometry,
elimination theory, real geometry, as well as their applications in various
domains such as partial differential equations, economics, probability, and
statistics.

  5. 5.

The book [194] and its second edition [195], describe useful algorithmic
aspects of systems with symmetries, treated using Invariant Theory. The book
[196] is concerned with polynomial systems with symmetries, that arise in the
context of dynamical systems. The french-language book [197] is devoted
exclusively to polynomial system solving methods. The french-language book
[198] contains a series of chapters on polynomial system solving.

  6. 6.

The books [199] and [200] describe the emerging area of Algebraic Statistics,
while the book [201] focuses on numerical aspects of polynomial system
solving. On the other hand, the book [202] focuses on theoretical aspects of
polynomial system solving. Furthermore, the books trilogy [203, 204, 205]
examine the foundation of polynomial system solving from the algebraic
geometry standpoint as well as the applications standpoint. The Galois group
of a polynomial is inextricably linked with the more general Galois Theory, we
mention the three books [206, 207, 208].

##  5 Homotopy / Continuation methods

Homotopy methods, also referred to as continuation methods, may be used for
finding solutions to polynomial equations and systems of polynomial equations
[209, 210, 211]. Polynomial equations often arise in kinematics and robotics
related problems [14, 12]; papers [6, 7] apply homotopy methods to solve
systems of power flow equations; large-scale integrated circuit designs and
protein-protein interaction equation are solved using homotopy method in [8],
[27] applies homotopy methods for solving SNEs arising in geodesy; [22]
applies homotopy for finding string vacua, [4] applies homotopy methods in
chemical engineering. Homotopy (or deformation) of a system of equations
${F}_{n}(x)={\mathit{\Theta}}_{n}\equiv(0,0,\ldots,0)^{\top}$ (cf. Eq. (1)) is
a function $H_{n}$ such that $H_{n}(x,1)=G_{n}(x)$, and
$H_{n}(x,0)={F}_{n}(x)$, where the roots of $G_{n}(x)$ are known. For example,
it is possible to choose a convex homotopy, $H_{n}(x,\lambda)=\lambda
G_{n}(x)+(1-\lambda){F}_{n}(x)$, where $\lambda\in[0,1]$ and trace the curve,
called the homotopy path, from a starting point $(x,1)$ to solution point
$(x,0)$. Thus, gradually deforming solutions of a starting system
$G_{n}(x)=\mathit{\Theta}_{n}$ into the solutions to the target system
${F}_{n}(x)=\mathit{\Theta}_{n}$. At each step of the process, a solution of
the current system $H_{n}(x,\tilde{\lambda})={\mathit{\Theta}}_{n}$ is used as
a starting solution to the next system
$H_{n}(x,\tilde{\lambda}+\Delta\lambda)={\mathit{\Theta}}_{n}$, which is
solved using Newton-type methods that require an invertible Jacobian. In the
case of polynomial systems it is trivial to find solutions of the starting
system $G_{n}(x)=\mathit{\Theta}_{n}$. Common approaches for defining
$G_{n}(x)$ include fixed point homotopy, where $G_{n}(x)=x-x^{0}$ (where
$x^{0}$ is a starting solution) and Newton homotopy, where
$G_{n}(x)=F_{n}(x)-F_{n}(x^{0})$. The latter is also referred to a global
homotopy, where $H_{n}(x,\lambda)={F}_{n}(x)-\lambda{F}_{n}(x^{0})$, where
$x^{0}$ is a starting solution [209]. An adaptive method for selecting the
steps of the homotopy path is presented in [212]. As per [213], advantages of
homotopy methods are 1) handling of singular solutions, 2) possibility to
obtain multiple solutions using one homotopy path, and 3) preserving Morse
indices for gradient systems. Paper [12] introduces a collision-based homotopy
continuation technique. Article [214] discusses the problem of divergent
homotopy paths and proposes an algorithm which performs projective path
tracking.

Article [215] describes the software package PHCpack for solving polynomial
systems using the homotopy method. Other software packages implementing
homotopy continuation methods include Bertini [216] and HomLab [217]. Paper
[218] describes a monodromy-based solver.

##  6 Interval methods

Interval arithmetic was introduced in article [219]. Books on interval methods
include [220] and [221]. A real interval $X$ is defined as a set of real
numbers between lower and upper bounds

$X=[a,b]=\\{x\in\mathbb{R}\,\,|\,\,a\leqslant x\leqslant b\\}.$

Interval arithmetic is a set of operations such as addition, subtraction,
multiplication, and division, defined on the intervals. Further, interval
functions, domain and range of which are the intervals, have been defined. In
addition to that, interval differentiation and integration have also been
proposed. The most important benefit of interval analysis is its
accountability for rounding errors due to limited machine precision. If a
value is represented as a single number, rounding errors occurring during the
computations may accumulate, thus leading to a wrong result. In interval
arithmetic a value is represented by a lower and upper bound which provide
reliable results during the computations. Further, interval arithmetic can be
used to model uncertainty, often arising in practical problems; for example,
uncertainty due to imprecise measurements [222]. Interval arithmetic was
standardized by the IEEE in 2015 [223]. Global optimization problems can be
solved using the interval branch-and-bound method which iteratively splits the
search space, and removes its parts that do not contain a global solution;
multiple splitting schemes have been proposed in the literature [224]. In this
case, application of interval arithmetic allows to guarantee, if a solution
exists within a region of interest. Paper [225] proposes to use interval
methods for inclusion and exclusion tests, where inclusion tests check
existence of the solution within an interval, and exclusion tests check its
non-existence. The interval Newton method [226] is a generalization of a
Newton method for interval arithmetic; it can be used to find zeros of a
function. Interval Newton method attempts to iteratively narrow down the new
interval solution $[x^{k+1}]$ from $[x^{k}]$ via

$[x^{k+1}]=[x^{k}]\cap\bigl{(}[x^{k}]-J_{{f}_{n}}\bigl{(}[x^{k}]\bigr{)}^{-1}\,{f}_{n}\bigl{(}[x^{k}]\bigr{)}\bigr{)},\kern
14.22636ptk=0,1,2,\ldots$

Paper [227] presents an interval branch-and-prune algorithm that is capable of
finding all solutions of a polynomial system, and the article evaluates the
proposed method on large systems with up to 320 variables. An interval method
was also effectively applied to solve large polynomial systems with up to 2500
variables in [228]. Article [229] provides a modification of the interval
branch-and-bound method and applies it to solving systems of equations. A
method of solving nonlinear equations by using interval arithmetic was
patented [230]. Article [231] presents details of implementation of a parallel
interval optimization algorithm and its application for solving systems of
nonlinear equations. A software package for solving equations using interval
methods is presented in [232] and [233]. Article [234] provides interval
extensions to Halley’s method (discussed in Section 3.3 of this paper), and
applies it to finding roots of a single nonlinear equation. Interval methods
are utilized to find all solutions of the kinematics SNEs in [235, 236].
Articles [237, 238] suggests to use linear programming (LP) for finding all
solutions of the SNE; the method in the paper is best suited to SNEs
consisting of linear equations with relatively few nonlinear terms. The method
is based on the surrounding nonlinear equations by rectangles, based on
interval arithmetic. Then, LP can be used to remove parts of the solution
space that do not contain the solutions. Paper [239] proposes to extend this
method and utilize LP narrowing, that is capable to solve large scale ($n$ =
50 000) separable SNEs. Article [240] proposes an extension of the interval
method for a SNE based on its transformation to separable form. Paper [241]
describes the UniCalc SNE solver software which is based on interval methods.
Another software package for interval optimization, that can also be used to
solve SNEs is RealPlayer [242].

##  7 Synopsis and concluding remarks

This article presents part one of a survey on methods for solving a system of
nonlinear equations (SNE). In part one we have presented a comprehensive
survey of methods which can be utilized to search for solutions to a SNE
without transforming a SNE into an optimization problem. Since many of the
SNEs that arise in real world applications are considered over a finite
bounded domain $\mathcal{D}$, we first introduced a technique which can be
utilized to determine the number of solutions to a SNE that exist within
$\mathcal{D}$. Then, we introduced a diverse set of root-finding methods which
can be used to search for solutions to a SNE. Next, we described additional
methods which have been used to search for solutions to SNEs including methods
from symbolic computation, homotopy / continuation methods, and interval
methods.

Analyzing this literature has led us to conclude that although there are a
variety of root-finding methods which are guaranteed to converge to a solution
when applied to SNEs that satisfy certain conditions (such as monotonicity
(Eq. 4) and Lipschitz continuity (Eq. 5)), for general SNEs, there is no
guarantee that any of the methods described in this paper will converge to a
solution in finite time. However, a new taxonomy of SNEs is needed to
facilitate the identification of new classes of tractable problems, and to
compare the performance of the methods that are most capable of solving them.
We are actively working on the development of such a taxonomy. Additionally,
we are actively exploring methods capable of solving systems of nonlinear
equations and inequalities.

In part one of this survey, we discussed methods for solving SNEs without
transforming them into optimization problems. In part two of this survey, we
will describe various transformations which can be utilized to transform a SNE
into an optimization problem, and we will discuss optimization algorithms
which can then be used to search for solutions. In part three of this survey,
we will present a robust quantitative comparative analysis of methods capable
of searching for solutions to SNEs.

## References

  * [1] Christodoulos A. Floudas.  Global optimization in design and control of chemical process systems.  Journal of Process Control, 10(2):125–134, 2000.
  * [2] Wagner F. Sacco and Nélio Henderson.  Finding all solutions of nonlinear systems using a hybrid metaheuristic with fuzzy clustering means.  Applied Soft Computing, 11(8):5424–5432, 2011.
  * [3] Astrid Holstad.  Numerical solution of nonlinear equations in chemical speciation calculations.  Computational Geosciences, 3(3):229–257, Dec 1999.
  * [4] Hugo Jiménez-Islas, Gloria M. Martínez-González, José L. Navarrete-Bolaños, José E. Botello-Álvarez, and J. Manuel Oliveros-Muñoz.  Nonlinear homotopic continuation methods: A chemical engineering perspective review.  Industrial & Engineering Chemistry Research, 52(42):14729–14742, Oct 2013.
  * [5] Nélio Henderson, Wagner F. Sacco, and Gustavo Mendes Platt.  Finding more than one root of nonlinear equations via a polarization technique: An application to double retrograde vaporization.  Chemical Engineering Research and Design, 88(5):551–561, 2010.
  * [6] Hsiao-Dong Chiang, Tian-Qi Zhao, Jiao-Jiao Deng, and Kaoru Koyanagi.  Homotopy-enhanced power flow methods for general distribution networks with distributed generators.  IEEE Transactions on Power Systems, 29(1):93–100, 2014.
  * [7] Dhagash Mehta, Hung Dinh Nguyen, and Konstantin Turitsyn.  Numerical polynomial homotopy continuation method to locate all the power flow solutions.  IET Generation, Transmission & Distribution, 10(12):2972–2980, 2016.
  * [8] Hsiao-Dong Chiang and Tao Wang.  Novel homotopy theory for nonlinear networks and systems and its applications to electrical grids.  IEEE Transactions on Control of Network Systems, 5(3):1051–1060, 2018.
  * [9] Chao Yang, Jianwen Cao, and Xiao-Chuan Cai.  A fully implicit domain decomposition algorithm for shallow water equations on the cubed-sphere.  SIAM Journal on Scientific Computing, 32(1):418–438, 2010.
  * [10] Matti Schneider, Daniel Wicht, and Thomas Böhlke.  On polarization-based schemes for the fft-based computational homogenization of inelastic materials.  Computational Mechanics, 64(4):1073–1095, Oct 2019.
  * [11] Yunong Zhang.  A set of nonlinear equations and inequalities arising in robotics and its online solution via a primal neural network.  Neurocomputing, 70(1):513–524, 2006.
  * [12] Amir Salimi Lafmejani, Ahmad Kalhor, and Mehdi Masouleh.  A new development of homotopy continuation method, applied in solving nonlinear kinematic system of equations of parallel mechanisms.  2015 3rd RSI International Conference on Robotics and Mechatronics (ICROM), 10 2015.
  * [13] David A Cox, John Little, and Donal O’Shea.  Robotics and automatic geometric theorem proving.  In Ideals, Varieties, and Algorithms, pages 291–343. Springer, 2015\.
  * [14] Ping Ji and Hongtao Wu.  Kinematics analysis of an offset 3-upu translational parallel robotic manipulator.  Robotics Auton. Syst., 42:117–123, 2003.
  * [15] Strachimir Cht. Mavrodiev and Maksym Deliyergiyev.  Modification of the nuclear landscape in the inverse problem framework using the generalized bethe-weizsäcker mass formula.  arXiv e-prints, pages arXiv–1602, 2016.
  * [16] Sani Aji, Poom Kumam, Punnarai Siricharoen, Auwal Bala Abubakar, and Mahmoud Muhammad Yahaya.  A modified conjugate descent projection method for monotone nonlinear equations and image restoration.  IEEE Access, 8:158656–158665, 2020.
  * [17] Jan Verschelde, Pierre Verlinden, and Ronald Cools.  Homotopies exploiting Newton polytopes for solving sparse polynomial systems.  SIAM Journal on Numerical Analysis, 31(3):915–930, 1994.
  * [18] Crina Grosan and Ajith Abraham.  A new approach for solving nonlinear equations systems.  IEEE Transactions on Systems, Man, and Cybernetics - Part A: Systems and Humans, 38(3):698–714, 2008.
  * [19] Ahmad Golbabai, Davood Ahmadian, and Mariyan Milev.  Radial basis functions with application to finance: American put option under jump diffusion.  Mathematical and Computer Modelling, 55(3):1354–1362, 2012.
  * [20] Guang Zhang and Liang Bai.  Existence of solutions for a nonlinear algebraic system.  Discrete dynamics in nature and society, 2009, 2009.
  * [21] Karol Kowalski and Karol Jankowski.  Towards complete solutions to systems of nonlinear equations of many-electron theories.  Phys. Rev. Lett., 81:1195–1198, Aug 1998.
  * [22] Dhagash Mehta.  Numerical polynomial homotopy continuation method and string vacua.  Advances in High Energy Physics, 2011:263937, Oct 2011.
  * [23] Yang Song, Chenlin Meng, Renjie Liao, and Stefano Ermon.  Nonlinear equation solving: A faster alternative to feedforward computation.  arXiv preprint arXiv:2002.03629, 2020.
  * [24] Yunfeng Cai and Ping Li.  Solving the robust matrix completion problem via a system of nonlinear equations.  In Silvia Chiappa and Roberto Calandra, editors, Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 4162–4172. PMLR, 26–28 Aug 2020.
  * [25] Stepan Yu. Gatilov.  Properties of nonlinear systems and convergence of the Newton-Raphson method in geometric constraint solving.  volume 32, pages 57––75. NCC, 2011.
  * [26] Béla Paláncz, Piroska Zaletnyik, Joseph L Awange, and Erik W Grafarend.  Dixon resultant’s solution of systems of geodetic polynomial equations.  Journal of Geodesy, 82(8):505–511, 2008.
  * [27] Béla Paláncz, Joseph L. Awange, Piroska Zaletnyik, and Robert H. Lewis.  Linear homotopy solution of nonlinear systems of equations in geodesy.  Journal of Geodesy, 84(1):79, Sep 2009.
  * [28] Christian Jansson.  An np-hardness result for nonlinear systems.  Reliable Computing, 4(4):345–350, Nov 1998.
  * [29] Yuri V. Matiyasevich.  Hilbert’s Tenth Problem.  MIT Press, Cambridge, MA, USA, 1993.
  * [30] Mustafa Q Khirallah and MA Hafiz.  Solving system of nonlinear equations using family of Jarratt methods.  International Journal of Differential Equations and Applications, 12(2), 2013.
  * [31] MQ Khirallah and MA Hafiz.  Novel three order methods for solving a system of nonlinear equations.  Bulletin of Mathematical Sciences & Applications, 1(2):01–14, 2012\.
  * [32] J Alikhani Koupaei and Seyed Mohammad Mehdi Hosseini.  A new hybrid algorithm based on chaotic maps for solving systems of nonlinear equations.  Chaos, Solitons & Fractals, 81:233–245, 2015.
  * [33] M. A. El-Shorbagy and Adel M. El-Refaey.  Hybridization of grasshopper optimization algorithm with genetic algorithm for solving system of non-linear equations.  IEEE Access, 8:220944–220961, 2020.
  * [34] Y Ramu Naidu and Akshay Kumar Ojha.  Solving multiobjective optimization problems using hybrid cooperative invasive weed optimization with multiple populations.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 48(6):821–832, 2016.
  * [35] Angel Kuri.  Solution of simultaneous non-linear equations using genetic algorithms.  WSEAS Transactions on Systems, 2, 01 2003.
  * [36] Lin Xiao, Zhijun Zhang, and Shuai Li.  Solving time-varying system of nonlinear equations by finite-time recurrent neural networks with application to motion tracking of robot manipulators.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 49(11):2210–2220, 2018.
  * [37] Huan-Tong Geng, Yi-Jie Sun, Qing-Xi Song, and Ting-Ting Wu.  Research of ranking method in evolution strategy for solving nonlinear system of equations.  In 2009 First International Conference on Information Science and Engineering, pages 348–351, 2009.
  * [38] Jun Pei, Zorica Dražić, Milan Dražić, Nenad Mladenović, and Panos M Pardalos.  Continuous variable neighborhood search (c-vns) for solving systems of nonlinear equations.  INFORMS Journal on Computing, 31(2):235–250, 2019.
  * [39] Abd allah Mousa and Islam Eldesoky.  Genls: Co-evolutionary algorithm for nonlinear system of equations.  Applied Mathematics and Computation, 197:633–642, 04 2008.
  * [40] Wenyin Gong, Yong Wang, Zhihua Cai, and Shengxiang Yang.  A weighted biobjective transformation technique for locating multiple optimal solutions of nonlinear equation systems.  IEEE Transactions on Evolutionary Computation, 21(5):697–713, 2017\.
  * [41] Wu Song, Yong Wang, Han-Xiong Li, and Zixing Cai.  Locating multiple optimal solutions of nonlinear equation systems based on multiobjective optimization.  IEEE Transactions on Evolutionary Computation, 19(3):414–431, 2015\.
  * [42] Sha Qin, Sanyou Zeng, Wei Dong, and Xi Li.  Nonlinear equation systems solved by many-objective hype.  In 2015 IEEE Congress on Evolutionary Computation (CEC), pages 2691–2696, 2015.
  * [43] Jinjin Guo, Binbin Qiu, and Yunong Zhang.  New-type dtz model for solving discrete time-dependent nonlinear equation system with robotic-arm application.  In 2020 10th International Conference on Information Science and Technology (ICIST), pages 153–162, 2020.
  * [44] Zuowen Liao, Wenyin Gong, and Ling Wang.  Memetic niching-based evolutionary algorithms for solving nonlinear equation system.  Expert Systems with Applications, 149:113261, 2020.
  * [45] Zuowen Liao, Wenyin Gong, Xuesong Yan, Ling Wang, and Chengyu Hu.  Solving nonlinear equations system with dynamic repulsion-based evolutionary algorithms.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 50(4):1590–1601, 2020.
  * [46] Weifeng Gao, Yuting Luo, Jingwei Xu, and Shengqi Zhu.  Evolutionary algorithm with multiobjective optimization technique for solving nonlinear equation systems.  Information Sciences, 541, 07 2020.
  * [47] Wenyin Gong, Yong Wang, Zhihua Cai, and Ling Wang.  Finding multiple roots of nonlinear equation systems via a repulsion-based adaptive differential evolution.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 50(4):1499–1513, 2020.
  * [48] Weifeng Gao, Genghui Li, Qingfu Zhang, Yuting Luo, and Zhenkun Wang.  Solving nonlinear equation systems by a two-phase evolutionary algorithm.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 51(9):5652–5663, 2021.
  * [49] Wei He, Wenyin Gong, Ling Wang, Xuesong Yan, and Chengyu Hu.  Fuzzy neighborhood-based differential evolution with orientation for nonlinear equation systems.  Knowledge-Based Systems, 182:104796, 2019.
  * [50] Jianye Wu, Wenyin Gong, and Ling Wang.  A clustering-based differential evolution with different crowding factors for nonlinear equations system.  Applied Soft Computing, 98:106733, 2021.
  * [51] Weifeng Gao and Yu Li.  Solving a new test set of nonlinear equation systems by evolutionary algorithm.  IEEE Transactions on Cybernetics, pages 1–10, 2021.
  * [52] Aijuan Song, Guohua Wu, and Witold Pedrycz.  Integrating variable reduction strategy with evolutionary algorithm for solving nonlinear equations systems.  CoRR, abs/2008.04223, 2020.
  * [53] Zuowen Liao, Wenyin Gong, Ling Wang, Xuesong Yan, and Chengyu Hu.  A decomposition-based differential evolution with reinitialization for nonlinear equations systems.  Knowledge-Based Systems, 191:105312, 2020.
  * [54] Shi Cheng, Junfeng Chen, Xiujuan Lei, and Yuhui Shi.  Locating multiple optima via brain storm optimization algorithms.  IEEE Access, 6:17039–17049, 2018.
  * [55] Wenyin Gong, Zuowen Liao, Xianyan Mi, Ling Wang, and Yuanyuan Guo.  Nonlinear equations solving with intelligent optimization algorithms: A survey.  Complex System Modeling and Simulation, 1(1):15–32, 2021.
  * [56] Masoud Ahookhosh, Keyvan Amini, and Somayeh Bahrami.  Two derivative-free projection approaches for systems of large-scale nonlinear monotone equations.  Numerical Algorithms, 64, 10 2013.
  * [57] Susumu Tanabé and Michael N. Vrahatis.  On perturbation of roots of homogeneous algebraic systems.  Mathematics of Computation, 255:1383–1402, 2006.
  * [58] Najib Ullah, Jamilu Sabi’u, and Abdullah Shah.  A derivative-free scaling memoryless broyden– fletcher–goldfarb–shanno method for solving a system of monotone nonlinear equations.  Numerical Linear Algebra with Applications, 28(5):e2374, 2021.
  * [59] Sani Aji, Poom Kumam, Aliyu Muhammed Awwal, Mahmoud Muhammad Yahaya, and Wiyada Kumam.  Two hybrid spectral methods with inertial effect for solving system of nonlinear monotone equations with application in robotics.  IEEE Access, 9:30918–30928, 2021.
  * [60] José Mario Martínez.  Algorithms for Solving Nonlinear Systems of Equations, pages 81–108.  Springer Netherlands, Dordrecht, 1994.
  * [61] James M. Ortega and Werner C. Rheinboldt.  Iterative Solution of Nonlinear Equations in Several Variables.  Society for Industrial and Applied Mathematics, January 2000.
  * [62] Émile Picard.  Sur le nombre des racines communes à plusieurs équations simultanées.  Journal de Mathématiques Pures et Appliquées, 8($4^{e}$ série):5–24, 1892.
  * [63] Émile Picard.  Traité d’Analyse.  Gauthier-Villars, Paris, 1922.
  * [64] Jane Cronin.  Fixed Points and Topological Degree in Nonlinear Analysis.  American Mathematical Society, Providence, Rhode Island, 1964.
  * [65] Noel G. Lloyd.  Degree Theory.  Oxford University Press, New York, 1978.
  * [66] Donal O’Regan, Yeol Je Cho, and Yu-Qing Chen.  Topological Degree Theory and Applications.  Taylor & Francis Group, Boca Raton, FL, 2006.
  * [67] Enrique Outerelo and Jesús M. Ruiz.  Mapping Degree Theory.  American Mathematical Society, Providence, Rhode Island, 2009.
  * [68] Krzysztof A. Sikorski.  Optimal Solution of Nonlinear Equations.  Oxford University Press, New York, 2001.
  * [69] Jeroen M. Bergamin, Tassos C. Bountis, and Michael N. Vrahatis.  Homoclinic orbits of invertible maps.  Nonlinearity, 15(5):1603–1619, 2002.
  * [70] Ioannis Z. Emiris, Bernard Mourrain, and Michael N. Vrahatis.  Sign methods for counting and computing real roots of algebraic systems.  Technical report, Rapport de Recherche No.3669, INRIA (Institut National de Recherche en Informatique et en Automatique), Sophia Antipolis, France, April 1999.
  * [71] Dimitris J. Kavvadias and Michael N. Vrahatis.  Locating and computing all the simple roots and extrema of a function.  SIAM Journal on Scientific Computing, 17(5):1232–1248, 19965.
  * [72] R. Baker Kearfott.  An efficient degree-computation method for a generalized method of bisection.  Numerische Mathematik, 32:109–127, 1979.
  * [73] Bernard Mourrain, Nicos G. Pavlidis, Dimitris K. Tasoulis, and Michael N. Vrahatis.  Determining the number of real roots of polynomials through neural networks.  Computers and Mathematics with Applications, 51(3-4):527–536, 2006\.
  * [74] Bernard Mourrain, Michael N. Vrahatis, and Yakoubsohn Jean-Claude.  On the complexity of isolating real roots and computing with certainty the topological degree.  Journal of Complexity, 18(2):612–640, 2002.
  * [75] Vassilis P. Plagianakos, Nicos K. Nousis, and Michael N. Vrahatis.  Locating and computing in parallel all the simple roots of special functions using pvm.  Journal of Computational and Applied Mathematics, 133(1-2):545–554, 2001.
  * [76] Chronis Polymilis, Graziano Servizi, Charalampos Skokos, Giorgio Turchetti, and Michael N. Vrahatis.  Topological degree theory and local analysis of area preserving maps.  Chaos, 13(1):94–104, 2003.
  * [77] Frank Stenger.  Computing the topological degree of a mapping in ${R}^{n}$.  Numerische Mathematik, 25:23–38, 1975.
  * [78] Martin Stynes.  An algorithm for numerical calculation of topological degree.  Applicable Analysis, 25:23–38, 1979.
  * [79] Michael N. Vrahatis.  An efficient method for locating and computing periodic orbits of nonlinear mappings.  Journal of Computational Physics, 119:105–119, 1995.
  * [80] Michael N. Vrahatis.  Solving systems of nonlinear equations using the nonzero value of the topological degree.  ACM Transactions on Mathematical Software, 14(4):312–329, Dec. 1988\.
  * [81] Michael N. Vrahatis.  A short proof and a generalization of Miranda’s existence theorem.  Proceedings of the American Mathematical Society, 107(3):701–703, 1989.
  * [82] Michael N. Vrahatis and Kosmas I. Iordanidis.  A rapid generalized method of bisection for solving systems of nonlinear equations.  Numerische Mathematik, 49(2-3):123–138, 1986.
  * [83] Michael N. Vrahatis, Graziano Servizi, Giorgio Turchetti, and Tassos C. Bountis.  A procedure to compute the fixed points and visualize the orbits of a 2D map.  Technical report, CERN SL/93-06 (AP), CERN-SL Division, European Organization for Nuclear Research (CERN), Geneva, Switzerland, February 1993.
  * [84] Peter Franek and Marek Krčál.  Robust satisfiability of systems of equations.  J. ACM, 62(4), sep 2015.
  * [85] Rodrigo G Eustaquio, Ademir A Ribeiro, and Miguel A Dumett.  A new class of root-finding methods in $r^{n}$: the inexact tensor-free chebyshev–halley class.  Computational and Applied Mathematics, 37(5):6654–6675, 2018.
  * [86] John E. Dennis and Robert B. Schnabel.  Numerical methods for unconstrained optimization and nonlinear equations.  SIAM, 1996.
  * [87] Carl T. Kelley.  Solving nonlinear equations with Newton’s method.  SIAM, 2003.
  * [88] Werner C. Rheinboldt.  Methods for solving systems of nonlinear equations.  SIAM, 1998.
  * [89] Marco Frontini and E Sormani.  Third-order methods from quadrature formulae for solving systems of nonlinear equations.  Applied Mathematics and Computation, 149(3):771–782, 2004.
  * [90] MT Darvishi and Ali Barati.  A fourth-order method from quadrature formulae to solve systems of nonlinear equations.  Applied Mathematics and Computation, 188(1):257–261, 2007.
  * [91] Hari M Srivastava, Javed Iqbal, Muhammad Arif, Alamgir Khan, Yusif S Gasimov, and Ronnason Chinram.  A new application of gauss quadrature method for solving systems of nonlinear equations.  Symmetry, 13(3):432, 2021.
  * [92] Muhammad Aslam Noor and Muhammad Waseem.  Some iterative methods for solving a system of nonlinear equations.  Computers & Mathematics with Applications, 57(1):101–106, 2009\.
  * [93] Adi Ben-Israel.  A modified newton-raphson method for the solution of systems of equations.  Israel Journal of Mathematics, 3(2):94–98, Jun 1965.
  * [94] N. E. Yudin.  Adaptive gauss–newton method for solving systems of nonlinear equations.  Doklady Mathematics, 104(2):293–296, Sep 2021.
  * [95] Adi Ben-Israel.  A newton-raphson method for the solution of systems of equations.  Journal of Mathematical Analysis and Applications, 15(2):243–252, 1966.
  * [96] José L. Hueso, Eulalia Martínez, and Juan R. Torregrosa.  Modified newton’s method for systems of nonlinear equations with singular jacobian.  Journal of Computational and Applied Mathematics, 224(1):77–83, 2009.
  * [97] Stepan Yu. Gatilov.  Using low-rank approximation of the Jacobian matrix in the Newton-Raphson method to solve certain singular equations.  J. Comput. Appl. Math., 272:8–24, 2014.
  * [98] Xiubin Xu and Chong Li.  Convergence of newton’s method for systems of equations with constant rank derivatives.  Journal of Computational Mathematics, 25(6):705–718, 2007.
  * [99] Theodoula N. Grapsa and Michael N. Vrahatis.  The implicit function theorem for solving systems of nonlinear equations in ${R}^{2}$.  International Journal of Computer Mathematics, 28(1-4):171–181, 1989.
  * [100] Theodoula N. Grapsa and Michael N. Vrahatis.  A dimension-reducing method for solving systems of nonlinear equations in ${R}^{n}$.  International Journal of Computer Mathematics, 32(3-4):205–216, 1990.
  * [101] Theodoula N. Grapsa and Michael N. Vrahatis.  A new dimension-reducing method for solving systems of nonlinear equations.  International Journal of Computer Mathematics, 55:235–244, 1995\.
  * [102] Theodoula N. Grapsa and Michael N. Vrahatis.  A dimension-reducing method for unconstrained optimization.  Journal of Computational and Applied Mathematics, 66:239–253, 1996\.
  * [103] Theodoula N. Grapsa, Michael N. Vrahatis, and Tassos C. Bountis.  Solving systems of nonlinear equations in ${R}^{n}$ using a rotating hyperplane in ${R}^{n+1}$.  International Journal of Computer Mathematics, 35(1-4):133–151, 1990.
  * [104] Charles G. Broyden.  A class of methods for solving nonlinear simultaneous equations.  Mathematics of Computation, 19(92):577–593, 1965.
  * [105] Márcia A. Gomes-Ruggiero, José Mario Martínez, and Antonio Carlos Moretti.  Comparing algorithms for solving sparse nonlinear systems of equations.  SIAM Journal on Scientific and Statistical Computing, 13(2):459–483, March 1992.
  * [106] Michael G. Epitropakis and Michael N. Vrahatis.  Studying the basin of convergence of methods for computing periodic orbits.  International Journal of Bifurcation and Chaos, 21(08):2079–2106, 2011.
  * [107] Youcef Saad and Martin H. Schultz.  Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems.  SIAM Journal on Scientific and Statistical Computing, 7(3):856–869, 1986.
  * [108] Walter E. Arnoldi.  The principle of minimized iterations in the solution of the matrix eigenvalue problem.  Quarterly of Applied Mathematics, 9:17–29, 1951.
  * [109] Krzysztof A Sikorski.  Bisection is optimal.  Numerische Mathematik, 40:111–117, 1982.
  * [110] Dimitris J. Kavvadias, Frosso S. Makri, and Michael N. Vrahatis.  Locating and computing arbitrarily distributed zeros.  SIAM Journal on Scientific Computing, 21(3):954–969, 1999.
  * [111] Dimitris J. Kavvadias, Frosso S. Makri, and Michael N. Vrahatis.  Efficiently computing many roots of a function.  SIAM Journal on Scientific Computing, 27(1):93–107, 2005.
  * [112] Michael N. Vrahatis.  An error estimation for the method of bisection in ${R}^{n}$.  Bulletin of the Greek Mathematical Society, 27:161–174, 1986.
  * [113] Michael N. Vrahatis.  CHABIS: A mathematical software package for locating and evaluating roots of systems of nonlinear equations.  ACM Transactions on Mathematical Software, 14(4):330–336, Dec. 1988\.
  * [114] Michael N. Vrahatis, George S. Androulakis, and George E. Manoussakis.  A new unconstrained optimization method for imprecise function and gradient values.  Journal of Mathematical Analysis and Applications, 197(2):586–607, 1996.
  * [115] Michael N. Vrahatis, Heinz Isliker, and Tassos C. Bountis.  Structure and breakdown of invariant tori in a 4-D mapping model of accelerator dynamics.  International Journal of Bifurcation and Chaos, 7(12):2707–2722, 1997.
  * [116] Dimitra-Nefeli A. Zottou, Dimitris J. Kavvadias, Frosso S. Makri, and Michael N. Vrahatis.  Manbis–a C++ mathematical software package for locating and computing efficiently many roots of a function: Theoretical issues.  ACM Transactions on Mathematical Software, 44(3):1–7, April 2018\.
  * [117] Stamatios-Aggelos N. Alexandropoulos, Panos M. Pardalos, and Michael N. Vrahatis.  Dynamic search trajectory methods for global optimization.  Annals of Mathematics and Artificial Intelligence, 88(1-3):3–37, 2020.
  * [118] Bassilis Boutsinas and Michael N. Vrahatis.  Artificial nonmonotonic neural networks.  Artificial Intelligence, 132(1):1–38, 2001.
  * [119] Michael N. Vrahatis, George S. Androulakis, John N. Lambrinos, and George D. Magoulas.  A class of gradient unconstrained minimization algorithms with adaptive stepsize.  Journal of Computational and Applied Mathematics, 114(2):367–386, 2000.
  * [120] Michael N. Vrahatis, George D. Magoulas, and Vassilis P Plagianakos.  Globally convergent modification of the quickprop method.  Neural Processing Letters, 12(2):159–169, 2000.
  * [121] M.A. Hernández and M.A. Salanova.  Indices of convexity and concavity. application to halley method.  Applied Mathematics and Computation, 103(1):27–49, 1999.
  * [122] MA Hernández and MA Salanova.  A family of chebyshev-halley type methods.  International Journal of Computer Mathematics, 47(1-2):59–63, 1993\.
  * [123] J.M. Gutiérrez and M.A. Hernández.  A family of chebyshev-halley type methods in banach spaces.  Bulletin of the Australian Mathematical Society, 55(1):113–130, 1997.
  * [124] DKR Babajee, MZ Dauhoo, MT Darvishi, A Karami, and Ali Barati.  Analysis of two chebyshev-like third order methods free from second derivatives for solving systems of nonlinear equations.  Journal of Computational and Applied Mathematics, 233(8):2002–2012, 2010.
  * [125] Trond Steihaug and Sara Suleiman.  Rate of convergence of higher order methods.  Applied Numerical Mathematics, 67:230–242, 2013.  NUMAN 2010.
  * [126] Ioannis K. Argyros, Dong Chen, and Qian Qingshan.  The jarratt method in banach space setting.  Journal of Computational and Applied Mathematics, 51(1):103–106, May 1994.
  * [127] Fayyaz Ahmad, Emran Tohidi, Malik Zaka Ullah, and Juan A Carrasco.  Higher order multi-step jarratt-like method for solving systems of nonlinear equations: Application to pdes and odes.  Computers & Mathematics with Applications, 70(4):624–636, 2015\.
  * [128] Thugal Zhanlav and Khuder Otgondorj.  Higher order jarratt-like iterations for solving systems of nonlinear equations.  Applied Mathematics and Computation, 395:125849, 2021.
  * [129] Paul David Frank.  Tensor methods for solving systems of nonlinear equations (numerical analysis, optimization).  PhD thesis, University of Colorado at Boulder, 1984.
  * [130] Robert B Schnabel and Paul D Frank.  Solving systems of nonlinear equations by tensor methods.  Technical report, Department of Computer Science, Colorado University at Boulder, 1986.
  * [131] Ali Bouaricha.  Solving large sparse systems of nonlinear equations and nonlinear least squares problems using tensor methods on sequential and parallel computers.  PhD thesis, University of Colorado at Boulder, 1992.
  * [132] Ali Bouaricha and Robert B Schnabel.  Tensor methods for large, sparse systems of nonlinear equations.  Preprint MCS-P473-1094, Mathematics and Computer Science Division, Argonne National Laboratory, 1994.
  * [133] Dan Feng and Thomas H Pulliam.  Tensor-gmres method for large systems of nonlinear equations.  SIAM Journal on Optimization, 7(3):757–779, 1997.
  * [134] Ali Bouaricha and Robert B Schnabel.  Algorithm 768: TENSOLVE: A software package for solving systems of nonlinear equations and nonlinear least-squares problems using tensor methods.  ACM Transactions on Mathematical Software (TOMS), 23(2):174–195, 1997.
  * [135] Ali Bouaricha and Robert B Schnabel.  Tensor methods for large sparse systems of nonlinear equations.  Mathematical programming, 82(3):377–400, 1998.
  * [136] Brett W Bader and Robert B Schnabel.  On the performance of tensor methods for solving ill-conditioned problems.  SIAM Journal on Scientific Computing, 29(6):2329–2351, 2004.
  * [137] Brett W. Bader.  Tensor-Krylov methods for solving large-scale systems of nonlinear equations.  SIAM Journal on Numerical Analysis, 43(3):1321–1347, 2005.
  * [138] Claude Brezinski.  A classification of quasi-Newton methods.  Numerical Algorithms, 33:123–135, 2004.
  * [139] Guangye Li.  Successive column correction algorithms for solving sparse nonlinear systems of equations.  Mathematical Programming, 43(1):187–207, Jan 1989.
  * [140] José Mario Martínez and Mário C. Zambaldi.  An inverse column-updating method for solving large–scale nonlinear systems of equations.  Optimization Methods and Software, 1(2):129–140, 1992.
  * [141] Charles G. Broyden.  The convergence of a class of double-rank minimization algorithms 1. General considerations.  IMA Journal of Applied Mathematics, 6(1):76–90, 03 1970.
  * [142] Roger Fletcher.  A new approach to variable metric algorithms.  The Computer Journal, 13(3):317–322, 01 1970.
  * [143] Donald Goldfarb.  A family of variable-metric methods derived by variational means.  Mathematics of Computation, 24(109):23–26, 1970.
  * [144] David F Shanno.  Conditioning of quasi-newton methods for function minimization.  Mathematics of computation, 24(111):647–656, 1970.
  * [145] Dong C. Liu and Jorge Nocedal.  On the limited memory bfgs method for large scale optimization.  Mathematical Programming, 45(1):503–528, Aug 1989.
  * [146] Ernesto G Birgin, Nataša Krejić, and José Mario Martínez.  Globally convergent inexact quasi-newton methods for solving nonlinear systems.  Numerical algorithms, 32(2):249–260, 2003.
  * [147] MK Dauda, AS Magaji, US Shehub, MA Usman, and MY Waziri.  A simple conjugate gradient type method for solving large-scale systems of nonlinear equations.  Malaysian Journal of Computing and Applied Mathematics, 3(2):25–35, 2020.
  * [148] David Ek.  Approaches to accelerate methods for solving systems of equations arising in nonlinear optimization.  PhD thesis, KTH Royal Institute of Technology, 2020.
  * [149] K Kamfa, MY Waziri, IM Sulaiman, M Mamat, and Hery Ibrahim.  A quasi-newton like method via modified rational approximation model for solving system of nonlinear equation.  Journal of Advanced Research in Dynamical and Control Systems, 2020\.
  * [150] Wah June Leong, Malik Abu Hassan, and Muhammad Waziri Yusuf.  A matrix-free quasi-newton method for solving large-scale nonlinear systems.  Computers & Mathematics with Applications, 62(5):2354–2363, 2011\.
  * [151] Ya-Zhong Luo, Guo-Jin Tang, and Li-Ni Zhou.  Hybrid approach for solving systems of nonlinear equations using chaos optimization and quasi-newton method.  Applied Soft Computing, 8(2):1068–1073, 2008.
  * [152] Mustafa Mamat, MK Dauda, MY Waziri, Fadhilah Ahmad, and Fatma Susilawati Mohamad.  Improved quasi-newton method via psb update for solving systems of nonlinear equations.  In AIP Conference Proceedings, volume 1782, page 030009. AIP Publishing LLC, 2016.
  * [153] José Mario Martínez.  Practical quasi-newton methods for solving nonlinear systems.  Journal of Computational and Applied Mathematics, 124(1):97–121, 2000.  Numerical Analysis 2000. Vol. IV: Optimization and Nonlinear Equations.
  * [154] Rosana Pérez and Véra Lucia Rocha Lopes.  Recent applications and numerical implementation of quasi-newton methods for solving nonlinear systems of equations.  Numerical Algorithms, 35(2):261–285, 2004.
  * [155] Alessandra Papini, Margherita Porcelli, and Cristina Sgattoni.  On the global convergence of a new spectral residual algorithm for nonlinear systems of equations.  Bollettino dell’Unione Matematica Italiana, pages 1–12, 2020.
  * [156] William La Cruz, José Mario Martínez, and Marcos Raydan.  Spectral residual method without gradient information for solving large-scale nonlinear systems of equations.  Mathematics of Computation, 75(255):1429–1449, April 2006.
  * [157] Zhong-Zhi Bai, Gene H. Golub, and Michael K. Ng.  Hermitian and skew-hermitian splitting methods for non-hermitian positive definite linear systems.  SIAM Journal on Matrix Analysis and Applications, 24(3):603–626, 2003.
  * [158] Qingbiao Wu and Minhong Chen.  Convergence analysis of modified newton-hss method for solving systems of nonlinear equations.  Numerical Algorithms, 64(4):659–683, 2013.
  * [159] Zhong-Zhi Bai and Xi Yang.  On hss-based iteration methods for weakly nonlinear systems.  Applied numerical mathematics, 59(12):2923–2936, 2009.
  * [160] Abdolreza Amiri, Alicia Cordero, Mohammad Taghi Darvishi, and Juan R Torregrosa.  A fast algorithm to solve systems of nonlinear equations.  Journal of Computational and Applied Mathematics, 354:242–258, 2019\.
  * [161] Kenneth Levenberg.  A method for the solution of certain non-linear problems in least squares.  Quarterly of Applied Mathematics, 2(2):164–168, 1944.
  * [162] Donald Marquardt.  An algorithm for least-squares estimation of nonlinear parameters.  SIAM Journal on Applied Mathematics, 11(2):431–441, 1963.
  * [163] Roger Fletcher.  A modified marquardt subroutine for nonlinear least squares.  Technical report, Atomic Energy Research Establishment, Harwell, England (United Kingdom), 1971.
  * [164] Valentina Kuzina and Alexander Koshev.  Modification of the levenberg–marquardt algorithm for solving complex computational construction problems.  In IOP Conference Series: Materials Science and Engineering, volume 960, page 032039. IOP Publishing, 2020.
  * [165] José de Jesús Rubio.  Stability analysis of the modified levenberg–marquardt algorithm for the artificial neural network training.  IEEE Transactions on Neural Networks and Learning Systems, 32(8):3510–3524, 2021.
  * [166] Kenneth M. Brown and William B. Gearhart.  Deflation techniques for the calculation of further solutions of a nonlinear system.  Numerische Mathematik, 16, 1971.
  * [167] James H. Wilkinson.  Rounding errors in algebraic processes.  Dover Publications, Inc., New York, 1994.
  * [168] Vassilis S. Kalantonis, Efstathios A. Perdios, Angela E. Perdiou, Omiros Ragos, and Michael N. Vrahatis.  Deflation techniques for the determination of periodic solutions of a certain period.  Astrophysics and Space Science, 288(2):591–599, 2003.
  * [169] Konstantinos E. Parsopoulos, Vasssilis P. Plagianakos, George D. Magoulas, and Michael N. Vrahatis.  Objective function “stretching” to alleviate convergence to local minima.  Nonlinear Analysis – Theory, Methods & Applications, 47(5):3419–3424, 2001.
  * [170] Konstantinos E. Parsopoulos and Michael N. Vrahatis.  Recent approaches to global optimization problems through particle swarm optimization.  Natural Computing, 1(2-3):235–306, 2002.
  * [171] Konstantinos E. Parsopoulos and Michael N. Vrahatis.  On the computation of all global minimizers through particle swarm optimization.  IEEE Transactions on Evolutionary Computation, 8(3):211–224, 2004\.
  * [172] Konstantinos E. Parsopoulos and Michael N. Vrahatis.  Particle swarm optimization and intelligence: Advances and applications.  Information Science Publishing (IGI Global), Hershey, PA, USA, 2010.
  * [173] Konstantinos E. Parsopoulos and Michael N. Vrahatis.  Deflection and stretching techniques for detection of multiple minimizers in multimodal optimization problems.  In Metaheuristics for Finding Multiple Solutions, M. Preuss, M.G. Epitropakis, J.E. Fieldsend and X. Li (eds.), Chapter 6, pages 129–144. Natural Computing Series, Springer Nature Switcherland AG, 2021.
  * [174] K. O. Geddes, S. R. Czapor, and G. Labahn.  Algorithms for computer algebra.  Kluwer Academic Publishers, Boston, MA, 1992.
  * [175] Michael J. Wester.  Computer Algebra Systems: A Practical Guide.  John Wiley & Sons, Chichester, United Kingdom, 1999.  https://math.unm.edu/~wester/cas/book/contents.html.
  * [176] J. H. Davenport, Y. Siret, and E. Tournier.  Computer algebra.  Academic Press, Ltd., London, second edition, 1993.  Systems and algorithms for algebraic computation, With a preface by Daniel Lazard, Translated from the French by A. Davenport and J. H. Davenport, With a foreword by Anthony C. Hearn.
  * [177] Joachim von zur Gathen and Jürgen Gerhard.  Modern computer algebra.  Cambridge University Press, Cambridge, third edition, 2013.
  * [178] Teo Mora.  Solving polynomial equation systems. I, volume 88 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, 2003.  The Kronecker-Duval philosophy.
  * [179] Teo Mora.  Solving polynomial equation systems. II, volume 99 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, 2005.  Macaulay’s paradigm and Gröbner technology.
  * [180] Teo Mora.  Solving polynomial equation systems. Vol. III. Algebraic solving, volume 157 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, 2015.
  * [181] Teo Mora.  Solving polynomial equation systems. Vol. IV. Buchberger theory and beyond, volume 158 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, 2016.
  * [182] William W. Adams and Philippe Loustaunau.  An introduction to Gröbner bases, volume 3 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 1994.
  * [183] Thomas Becker and Volker Weispfenning.  Gröbner bases, volume 141 of Graduate Texts in Mathematics.  Springer-Verlag, New York, 1993.  A computational approach to commutative algebra, In cooperation with Heinz Kredel.
  * [184] Ralf Fröberg.  An introduction to Gröbner bases.  Pure and Applied Mathematics (New York). John Wiley & Sons, Ltd., Chichester, 1997.
  * [185] Franz Winkler.  Polynomial algorithms in computer algebra.  Texts and Monographs in Symbolic Computation. Springer-Verlag, Vienna, 1996.
  * [186] Wolfram Decker and Christoph Lossen.  Computing in algebraic geometry, volume 16 of Algorithms and Computation in Mathematics.  Springer-Verlag, Berlin; Hindustan Book Agency, New Delhi, 2006.  A quick start using SINGULAR.
  * [187] Wolfram Decker and Gerhard Pfister.  A first course in computational algebraic geometry.  African Institute of Mathematics (AIMS) Library Series. Cambridge University Press, Cambridge, 2013.
  * [188] Martin Kreuzer and Lorenzo Robbiano.  Computational commutative algebra. 1.  Springer-Verlag, Berlin, 2000.
  * [189] Martin Kreuzer and Lorenzo Robbiano.  Computational commutative algebra. 2.  Springer-Verlag, Berlin, 2005.
  * [190] Martin Kreuzer and Lorenzo Robbiano.  Computational commutative algebra 1.  Springer-Verlag, Berlin, 2008.  Corrected reprint of the 2000 original.
  * [191] Martin Kreuzer and Lorenzo Robbiano.  Computational linear and commutative algebra.  Springer, Cham, 2016.
  * [192] Daniel R. Grayson and Michael E. Stillman.  Macaulay2, a software system for research in algebraic geometry.  Available at https://math.uiuc.edu/Macaulay2/.
  * [193] Bernd Sturmfels.  Solving systems of polynomial equations, volume 97 of CBMS Regional Conference Series in Mathematics.  Published for the Conference Board of the Mathematical Sciences, Washington, DC; by the American Mathematical Society, Providence, RI, 2002.
  * [194] Bernd Sturmfels.  Algorithms in invariant theory.  Texts and Monographs in Symbolic Computation. Springer-Verlag, Vienna, 1993.
  * [195] Bernd Sturmfels.  Algorithms in invariant theory.  Texts and Monographs in Symbolic Computation. Springer Wien, New York, Vienna, second edition, 2008.
  * [196] Karin Gatermann.  Computer algebra methods for equivariant dynamical systems, volume 1728 of Lecture Notes in Mathematics.  Springer-Verlag, Berlin, 2000.
  * [197] Mohamed Elkadi and Bernard Mourrain.  Introduction à la résolution des systèmes polynomiaux.  Mathématiques et Applications 59, 2007.
  * [198] Alin Bostan, Frédéric Chyzak, Marc Giusti, Romain Lebreton, Grégoire Lecerf, Bruno Salvy, and Éric Schost.  Algorithmes Efficaces en Calcul Formel.  HAL Archives Ouvertes, 2017.
  * [199] Giovanni Pistone, Eva Riccomagno, and Henry P. Wynn.  Algebraic statistics, volume 89 of Monographs on Statistics and Applied Probability.  Chapman & Hall/CRC, Boca Raton, FL, 2001.  Computational commutative algebra in statistics.
  * [200] Seth Sullivant.  Algebraic statistics, volume 194 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 2018.
  * [201] Hans J. Stetter.  Numerical polynomial algebra.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2004.
  * [202] Alicia Dickenstein and Ioannis Z. Emiris.  Solving Polynomial Equations. Foundations, Algorithms, and Applications, volume 14.  Algorithms and Computation in Mathematics, Springer, Berlin, Heidelberg, 2005.
  * [203] David A. Cox.  Applications of polynomial systems, volume 134 of CBMS Regional Conference Series in Mathematics.  American Mathematical Society, Providence, RI, 2020.
  * [204] David A. Cox, John Little, and Donal O’Shea.  Using algebraic geometry, volume 185 of Graduate Texts in Mathematics.  Springer, New York, second edition, 2005.
  * [205] David A. Cox, John Little, and Donal O’Shea.  Ideals, varieties, and algorithms.  Undergraduate Texts in Mathematics. Springer, Cham, fourth edition, 2015\.  An introduction to computational algebraic geometry and commutative algebra.
  * [206] David A. Cox.  Galois theory.  Pure and Applied Mathematics (Hoboken). John Wiley & Sons, Inc., Hoboken, NJ, second edition, 2012.
  * [207] Ian Stewart.  Galois theory.  CRC Press, Boca Raton, FL, fourth edition, 2015.
  * [208] Jean-Pierre Tignol.  Galois’ theory of algebraic equations.  World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, second edition, 2016.
  * [209] Eugene L. Allgower and Kurt Georg.  Introduction to Numerical Continuation Methods.  SIAM, 2003.
  * [210] Layne T. Watson.  Globally convergent homotopy algorithms for nonlinear systems of equations.  Nonlinear Dynamics, 1(2):143–191, Mar 1990.
  * [211] Alexander Morgan.  Solving Polynomial Systems Using Continuation for Engineering and Scientific Problems.  SIAM, 2009.
  * [212] Werner C Rheinboldt.  An adaptive continuation process for solving systems of nonlinear equations.  University of Maryland, 1975.
  * [213] Tianran Chen and Tien-Yien Li.  Homotopy continuation method for solving systems of nonlinear and polynomial equations.  Communications in Information and Systems, 15(2):119–307, 2015\.
  * [214] Tianran Chen and Tien-Yien Li.  Spherical projective path tracking for homotopy continuation methods.  Communications in Information and Systems, 12(3):195–220, 2012\.
  * [215] Jan Verschelde.  Algorithm 795: Phcpack: A general-purpose solver for polynomial systems by homotopy continuation.  ACM Trans. Math. Softw., 25(2):251–276, jun 1999.
  * [216] Daniel J. Bates, Jonathan D. Haunstein, Andrew J. Sommese, and Charles W. Wampler.  Numerically Solving Polynomial Systems with Bertini.  Society for Industrial and Applied Mathematics, USA, 2013.
  * [217] Andrew Sommese and Charles Wampler.  The Numerical Solution of Systems of Polynomials Arising in Engineering and Science.  03 2005.
  * [218] Timothy Duff, Cvetelina Hill, Anders Jensen, Kisun Lee, Anton Leykin, and Jeff Sommars.  Solving polynomial systems via homotopy continuation and monodromy.  IMA Journal of Numerical Analysis, 39(3):1421–1446, 04 2018.
  * [219] Ramon E Moore.  Interval analysis.  Prentice-Hall, 1966.
  * [220] Arnold Neumaier.  Interval Methods for Systems of Equations.  Encyclopedia of Mathematics and its Applications. Cambridge University Press, 1991.
  * [221] Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud.  Introduction to Interval Analysis.  Society for Industrial and Applied Mathematics, USA, 2009.
  * [222] Vladik Kreinovich.  Solving equations (and systems of equations) under uncertainty: how different practical problems lead to different mathematical and computational formulations.  Granular Computing, 1(3):171–179, Sep 2016.
  * [223] IEEE.  IEEE 1788-2015 – IEEE standard for interval arithmetic.  https://standards.ieee.org/standard/1788-2015.html, 2015.
  * [224] Helmut Ratschek and Jon G. Rokne.  Interval global optimization.  In Christodoulos A. Floudas and Panos M. Pardalos, editors, Encyclopedia of Optimization, Second Edition, pages 1739–1757. Springer, 2009\.
  * [225] Matthew Stuber, V. Kumar, and P. Barton.  Nonsmooth exclusion test for finding all solutions of nonlinear equations.  BIT Numerical Mathematics, 50:885–917, 12 2010.
  * [226] G. William Walster and Eldon R. Hansen.  Global optimization using interval analysis: Revised and expanded.  2007\.
  * [227] Pascal Van Hentenryck, David McAllester, and Deepak Kapur.  Solving polynomial systems using a branch and prune approach.  SIAM Journal on Numerical Analysis, 34(2):797–827, 1997.
  * [228] Kiyotaka Yamamura and Koki Suda.  An efficient algorithm for finding all solutions of separable systems of nonlinear equations.  BIT Numerical Mathematics, 47:681–691, 2007.
  * [229] Bartlomiej Jacek Kubica.  Excluding regions using sobol sequences in an interval branch-and-prune method for nonlinear systems.  Reliab. Comput., 19(4):385–397, 2013.
  * [230] G. William Walster and Eldon R. Hansen.  Us6915321b2: Method and apparatus for solving systems of nonlinear equations using interval arithmetic, 2002.
  * [231] Yves Papegay, David Daney, and Jean-Pierre Merlet.  Parallel implementation of interval analysis for equations solving.  In Jack Dongarra, Domenico Laforenza, and Salvatore Orlando, editors, Recent Advances in Parallel Virtual Machine and Message Passing Interface, pages 555–559, Berlin, Heidelberg, 2003. Springer Berlin Heidelberg.
  * [232] Max E. Jerrell.  Finding and verifying all solutions of a system of nonlinear equations.  IFAC Proceedings Volumes, 31(16):327–332, 1998.  IFAC Symposium on Computation in Economics, Finance and Engineering: Economic Systems, Cambridge, UK, 29 June - 1 July.
  * [233] Max E Jerrell and Wendy A Campione.  Finding and verifying all solutions of a system of nonlinear equations using public domain software.  2002\.
  * [234] Tahereh Eftekhari.  Interval extensions of the halley method and its modified method for finding enclosures of roots of nonlinear equations.  Computational Methods for Differential Equations, 8(2):222–235, 2020.
  * [235] Olivier Didrit, Michel Petitot, and Eric Walter.  Guaranteed solution of direct kinematic problems for general configurations of parallel manipulators.  IEEE Transactions on Robotics and Automation, 14(2):259–266, 1998\.
  * [236] A Castellet and F Thomas.  An algorithm for the solution of inverse kinematics problems based on an interval method.  In Advances in Robot Kinematics: Analysis and Control, pages 393–402. Springer, 1998.
  * [237] Kiyotaka Yamamura, Hitomi Kawata, and Ai Tokue.  Interval solution of nonlinear equations using linear programming.  BIT Numerical Mathematics, 38(1):186–199, Mar 1998.
  * [238] Kiyotaka Yamamura and Tsuyoshi Fujioka.  Finding all solutions of nonlinear equations using the dual simplex method.  Journal of Computational and Applied Mathematics, 152(1):587–595, 2003.  Proceedings of the International Conference on Recent Advances in Computational Mathematics.
  * [239] Kiyotaka Yamamura, Koki Suda, and Naoya Tamura.  Lp narrowing: A new strategy for finding all solutions of nonlinear equations.  Applied Mathematics and Computation, 215(1):405–413, 2009.
  * [240] Lubomir V Kolev.  An improved method for global solution of non-linear systems.  Reliable Computing, 5(2):103–111, 1999.
  * [241] A. Babichev, O. Kadyrova, T. Kashevarova, A. Leshchenko, and A. Semenov.  Unicalc, a novel approach to solving systems of algebraic equations.  Interval Computations, 1993, 01 1993.
  * [242] Laurent Granvilliers and Frédéric Benhamou.  Algorithm 852: Realpaver: An interval solver using constraint satisfaction techniques.  ACM Trans. Math. Softw., 32(1):138–156, mar 2006.