∎

11institutetext:  Xin-long Luo, Corresponding author 22institutetext: School
of Artificial Intelligence, Beijing University of Posts and
Telecommunications, P. O. Box 101, Xitucheng Road No. 10, Haidian District,
100876, Beijing China, 22email<EMAIL_ADDRESS>33institutetext: Hang
Xiao 44institutetext: School of Artificial Intelligence, Beijing University of
Posts and Telecommunications, P. O. Box 101, Xitucheng Road No. 10, Haidian
District, 100876, Beijing China, 44email<EMAIL_ADDRESS>
# Continuation Newton methods with deflation techniques and quasi-genetic
evolution for global optimization problems

Xin-long Luo$\ast$ Hang Xiao

(Received: date / Accepted: date)

###### Abstract

The global minimum point of an optimization problem is of interest in
engineering fields and it is difficult to be solved, especially for a
nonconvex large-scale optimization problem. In this article, we consider a new
memetic algorithm for this problem. That is to say, we use the determined
points (the stationary points of the function) as the initial seeds of the
evolutionary algorithm, other than the random initial seeds of the known
evolutionary algorithms. We compare it with the multi-start method (the built-
in subroutine GlobalSearch.m of the MATLAB R2020a environment), the branch-
and-bound method (Couenne of the state-of-the-art open-source solver for mixed
integer nonlinear programming problems), and two representative derivative-
free algorithms (CMA-ES and MCS), respectively. Numerical results show that
the proposed method performs well for the large-scale global optimization
problems, especially the problems of which are difficult to be solved by the
known global optimization methods.

###### Keywords:

continuation Newton method deflation technique nonconvex optimization global
optimization large-scale problem genetic evolution

###### MSC:

65K05 65L05 65L20

††journal: Journal of XXXX

##  1 Introduction

In this article, we are mainly concerned with the global minimum of the
unconstrained optimization problem

$\displaystyle\min_{x\in\Re^{n}}f(x),$ (1)

where $f:\;\Re^{n}\to\Re$ is a differentiable function. When problem (1) is
nonconvex and large-scale, it may have many local minimizers and it is
challenge to find its global minimum. For problem (1), there are some popular
global optimization methods, such as the multi-start methods Boender1987 ;
DSDW2016 ; RS2013 ; ULPG2007 , the branch-and-bound methods BLLMW2009 ;
Couenne2020 ; Sahinidis2021 ; TS2005 , the genetic evolution algorithms
Hansen2006 ; MTK1999 ; Mitchell1996 ; LW2001  and their memetic algorithms
Hart1994 ; MKCR2020 ; TRS2003 ; SGKZ2014 .

For multi-start methods, they use the local minimum method such as the trust-
region method DS1996 ; CGT2000 ; Yuan2015  or the line search method SY2006 ;
NW1999  to solve problem (1) and find its local minimum point $x^{\ast}$,
which satisfies the following first-order optimality condition

$\displaystyle F(x^{\ast})=\nabla f(x^{\ast})=0.$ (2)

Since the implementations of multi-start methods are simple and successful for
many real-world problems, many commercial global solvers are based on these
methods such as GlobalSearch.m of the MATLAB R2020a environment MATLAB ;
ULPG2007 . However, for the multi-start method, a starting point could easily
lie in the basin of attraction of a previously found local minimum point
RS2013 . Hence, the algorithm might find the same local minimum even if the
searches are initiated from points those are far apart.

The evolutionary algorithm Hansen2006 ; Moscato1989 ; SKM2018 ; SGKZ2014  is
another popular method and used by a broad community of engineers and
practitioners to solve real-world problems due to its simplicity and
attractive nature-inspired interpretations. However, for the evolutionary
algorithm, since there is not a guaranteed convergence and its computational
time is very heavy for the large-scale problem, it easily fails to find the
global minimum accurately.

In order to overcome the shortcomings of these two methods and make use of
their advantages, we revise the continuation Newton method AG2003 ; LXL2021 ;
LX2021  by using the deflation technique BG1971  to find the stationary points
of $f(x)$ as many as possible. Then, we select the first $L$ smallest
stationary points of $f(x)$ as the initial evolutionary seeds of the quasi-
genetic algorithm, other than the random initial seeds of the known
evolutionary algorithms MTK1999 ; Mitchell1996 ; Moscato1989 . After the
quasi-genetic algorithm evolves into several generations such as twenty
generations, we obtain a suboptimal point $x_{ag}^{\ast}$ near the global
minimum point. Finally, we use the continuation Newton method LXL2021 ; LX2021
with this suboptimal point $x_{ag}^{\ast}$ as the initial point to solve the
nonlinear system (2) of equations and obtain one stationary point
$x_{cn}^{\ast}$ of $f(x)$. Finally, we output
$x_{min}^{\ast}=\arg\min\\{f(x_{ag}^{\ast}),\;f(x_{cn}^{\ast})\\}$ as the
global minimum point of the algorithm.

The rest of this article is organized as follows. In the next section, we
convert the nonlinear system (2) into the continuation Newton flow. Then, we
give the continuation Newton method with the trust-region updating strategy to
follow the continuation Newton flow and obtain its steady-state solution, i.e.
a stationary point of $f(x)$. In section 3, we consider the deflation
technique to construct a new nonlinear function $G_{k}(\cdot)$ such that the
zeros of $G_{k}(x)$ are the stationary points of $f(x)$ and exclude the
previously found stationary points $x_{1}^{\ast},\,\ldots,\,x_{k}^{\ast}$ of
$f(x)$. In section 4, we give a quasi-genetic algorithm to obtain the
suboptimal point by using the evolutionary strategy of the convex combination
and using these stationary points $x_{1}^{\ast},\,\ldots,\,x_{K}^{\ast}$ as
the evolutionary seeds. Then, we refine suboptimal solution by using the
continuation Newton method from this suboptimal point to solve the system (2)
of nonlinear equations. In section 5, some promising numerical results of the
proposed method are reported, with comparison to the multi-start method (the
subroutine GlobalSearch.m of the MATLAB R2020a environment), the branch-and-
bound method (Couenne of the state-of-the-art open-source solver BLLMW2009 ;
Couenne2020 ), and two representative derivative-free algorithms RS2013  (the
covariance matrix adaptation evolution strategy (CMA-ES) Hansen2006 ;
Hansen2010  and the multilevel coordinate search (MCS) HN1999 ; Neumaier2000 )
for global optimization problems. Finally, some discussions and conclusions
are also given in section 6. $\|\cdot\|$ denotes the Euclidean vector norm or
its induced matrix norm throughout the paper.

##  2 Finding a stationary point with the continuation Newton method

For the completeness of the paper, we restate the continuation Newton method
with the trust-region updating strategy LXL2021 ; LX2021  for finding a
solution of the nonlinear system (2).

If we consider the damped Newton method for the nonlinear system (2)
Kelley2003 ; Kelley2018 ; OR2000 , we have

$\displaystyle x_{k+1}=x_{k}-\alpha_{k}J(x_{k})^{-1}F(x_{k}),$ (3)

where $J(x_{k})$ is the Jacobian matrix of $F(x_{k})$. We regard
$x_{k+1}=x(t_{k}+\alpha_{k})$, $x_{k}=x(t_{k})$ and let $\alpha_{k}\to 0$,
then we obtain the following continuous Newton flow AS2015 ; Branin1972 ;
Davidenko1953 ; LXL2021 ; Tanabe1979 :

$\displaystyle\frac{dx(t)}{dt}=-J(x)^{-1}F(x),\hskip 5.69054ptx(t_{0})=x_{0}.$
(4)

Actually, if we apply an iteration with the explicit Euler method AP1998  for
the continuous Newton flow (4), we obtain the damped Newton method (3). Since
the Jacobian matrix $J(x)=F^{\prime}(x)$ may be singular, we reformulate the
continuous Newton flow (4) as the following general formula:

$\displaystyle-J(x)\frac{dx(t)}{dt}=F(x),\hskip 5.69054ptx(t_{0})=x_{0}.$ (5)

The continuous Newton flow (5) is an old method and can be backtracked to
Davidenko’s work Davidenko1953  in 1953. After that, it was investigated by
Branin Branin1972 , Deuflhard et al DPR1975 , Tanabe Tanabe1979  and Abbott
Abbott1977  in 1970s, and applied to nonlinear boundary problems by Axelsson
and Sysala AS2015  recently. The continuous and even growing interest in this
method originates from its some nice properties. One of them is that the
solution $x(t)$ of the continuous Newton flow (5) converges to the steady-
state solution $x^{\ast}$ from any initial point $x_{0}$, as described by the
following property 1.

###### Property 1

(Branin Branin1972  and Tanabe Tanabe1979 ) Assume that $x(t)$ is the solution
of the continuous Newton flow (5), then $E(x(t))=\|F(x)\|^{2}$ converges to
zero when $t\to\infty$. That is to say, for every limit point $x^{\ast}$ of
$x(t)$, it is also a solution of the nonlinear system (2). Furthermore, every
element $F^{i}(x)$ of $F(x)$ has the same convergence rate $exp(-t)$ and
$x(t)$ can not converge to the solution $x^{\ast}$ of the nonlinear system (2)
on the finite interval when the initial point $x_{0}$ is not a solution of the
nonlinear system (2).

Proof. Assume that $x(t)$ is the solution of the continuous Newton flow (5),
then we have

$\displaystyle\frac{d}{dt}\left(e^{t}F(x)\right)=e^{t}J(x)\frac{dx(t)}{dt}+e^{t}F(x)=0.$

Consequently, we obtain

$\displaystyle F(x(t))=F(x_{0})e^{-t}.$ (6)

From equation (6), it is not difficult to know that every element $F^{i}(x)$
of $F(x)$ converges to zero with the convergence rate $e^{-t}$ when
$t\to\infty$. Thus, if the solution $x(t)$ of the continuous Newton flow (5)
belongs to a compact set, it has a limit point $x^{\ast}$ when $t\to\infty$,
and this limit point $x^{\ast}$ is also a solution of the nonlinear system
(2).

If we assume that the solution $x(t)$ of the continuous Newton flow (5)
converges to the solution $x^{\ast}$ of the nonlinear system (2) on the finite
interval $(0,\,T]$, from equation (6), we have

$\displaystyle F(x^{\ast})=F(x_{0})e^{-T}.$ (7)

Since $x^{\ast}$ is a solution of the nonlinear system (2), we have
$F(x^{\ast})=0$. By substituting it into equation (7), we obtain

$\displaystyle F(x_{0})=0.$

Thus, it contradicts the assumption that $x_{0}$ is not a solution of the
nonlinear system (2). Consequently, the solution $x(t)$ of the continuous
Newton flow (5) can not converge to the solution $x^{\ast}$ of the nonlinear
system (2) on the finite interval. ∎

###### Remark 1

The inverse $J(x)^{-1}$ of the Jacobian matrix $J(x)$ can be regarded as the
pre-conditioner of $F(x)$ such that the solution elements
$x^{i}(t)\,(i=1,\,2,\ldots,n)$ of the continuous Newton flow (4) have the
roughly same convergence rates and it mitigates the stiff property of the ODE
(4) (the definition of the stiff problem can be found in HW1996  and
references therein). This property is very useful since it makes us adopt the
explicit ODE method to follow the trajectory of the Newton flow.

Actually, if we consider $F(x)=Ax$, from the ODE (5), we have

$\displaystyle A\frac{dx}{dt}=-Ax,\;x(0)=x_{0}.$ (8)

By integrating the linear ODE (8), we obtain

$\displaystyle x(t)=e^{-t}x_{0}.$ (9)

From equation (9), we know that the solution $x(t)$ of the ODE (8) converges
to zero exponentially with the same rate $e^{-t}$ when $t$ tends to infinity.

When the Jacobian matrix $J(x)$ is singular or nearly singular, the ODE (5) is
the system of differential-algebraic equations (DAEs) and its trajectory can
not be efficiently followed by the general ODE method such as the backward
differentiation formulas (the built-in subroutine ode15s.m of the MATLAB
R2020a environment AP1998 ; HW1996 ; MATLAB ; SGT2003 ). Thus, we need to
construct the special method to handle this problem.

Since the continuous Newton flow (5) is intrinsically a nonlinear diminishing
system for the energy function $f(x(t))=\|F(x(t))\|^{2}$, it can be integrated
by the strong stability preserving methods GS1998 ; GST2001  and the steady-
state solution $x^{\ast}$ can be obtained after the long time integration.
Here, we consider another continuation approach based on the traditional
optimization methods for problem (5). We expect that the new method has the
global convergence as the homotopy continuation methods and the fast
convergence rate near the solution $x^{\ast}$ as the merit-function methods.
In order to achieve these two aims, we construct the special continuation
Newton method with the new step size $\alpha_{k}=\Delta t_{k}/(1+\Delta
t_{k})$ LXL2021 ; LY2021 ; LXLZ2021 ; LX2021 ; LX2021B  and the time step
$\Delta t_{k}$ is adaptively adjusted by the trust-region updating strategy
for problem (5).

Firstly, we apply the implicit Euler method to the continuous Newton flow (5)
AP1998 ; HW1996 , then we obtain

$\displaystyle J(x_{k+1})\left(x_{k+1}-x_{k}\right)=-{\Delta
t_{k}}F(x_{k+1}).$ (10)

The scheme (10) is an implicit formula and it needs to solve a system of
nonlinear equations at every iteration. To avoid solving the system of
nonlinear equations, we replace the Jacobian matrix $J(x_{k+1})$ with
$J(x_{k})$ and substitute $F(x_{k+1})$ with its linear approximation
$F(x_{k})+J(x_{k})(x_{k+1}-x_{k})$ into equation (10), respectively. Thus, we
obtain the following explicit continuation Newton method:

$\displaystyle J(x_{k})s_{k}^{N}=-F(x_{k}),\;s_{k}=\frac{\Delta
t_{k}}{1+\Delta t_{k}}s_{k}^{N},\;x_{k+1}=x_{k}+s_{k},$ (11)

where $J_{k}=J(x_{k}),\;F_{k}=F(x_{k})$.

###### Remark 2

The continuation Newton method (11) is similar to the damped Newton method (3)
if we let $\alpha_{k}=\Delta t_{k}/(1+\Delta t_{k})$ in equation (11).
However, from the view of the ODE method, they are different. The damped
Newton method (3) is obtained by the explicit Euler method applied to the
continuous Newton flow (5), and its time-stepping size $\alpha_{k}$ is
restricted by the numerical stability HW1996 . That is to say, for the linear
test equation $dx/dt=-\lambda x$, its time step size $\alpha_{k}$ is
restricted by the stable region $|1-\lambda\alpha_{k}|\leq 1$. Therefore, the
large time step $\alpha_{k}$ can not be adopted in the steady-state phase. The
continuation Newton method (11) is obtained by the implicit Euler method and
its linear approximation applied to the continuous Newton flow (5), and its
time step size $\Delta t_{k}$ is not restricted by the numerical stability for
the linear test equation. Therefore, the large time step can be adopted in the
steady-state phase for the continuation Newton method (11), and it mimics the
Newton method near the solution $x^{\ast}$ such that it has the fast local
convergence rate. The most of all, $\alpha_{k}=\Delta t_{k}/(\Delta t_{k}+1)$
in equation (11) is favourable to adopt the trust-region updating strategy for
adaptively adjusting the time step $\Delta t_{k}$ such that the continuation
Newton method (11) accurately traces the trajectory of the continuous Newton
flow in the transient-state phase and achieves the fast convergence rate near
the equilibrium point $x^{\ast}$.

Another issue is how to adaptively adjust the time step $\Delta t_{k}$ at
every iteration. There is a popular way to control the time step size based on
the trust-region updating strategy Deuflhard2004 ; Higham1999 ; LL2010 ;
Luo2010 ; LKLT2009 ; Luo2012 ; LLS2021 ; LY2021 . Its main idea is that the
time step size $\Delta t_{k+1}$ will be enlarged when the linear model
$F(x_{k})+J(x_{k})s_{k}$ approximates $F(x_{k}+s_{k})$ well, and $\Delta
t_{k+1}$ will be reduced when $F(x_{k})+J(x_{k})s_{k}$ approximates
$F(x_{k}+s_{k})$ badly. Thus, by using the relation (11), we enlarge or reduce
the time step size $\Delta t_{k}$ at every iteration according to the
following ratio:

$\displaystyle\rho_{k}=\frac{\|F(x_{k})\|-\|F(x_{k}+s_{k})\|}{\|F(x_{k})\|-\|F(x_{k})+J(x_{k})s_{k}\|}=\frac{\|F(x_{k})\|-\|F(x_{k}+s_{k})\|}{(\Delta
t_{k}/(1+\Delta t_{k}))\|F(x_{k})\|}.$ (12)

A particular adjustment strategy is given as follows:

$\displaystyle\Delta t_{k+1}=\begin{cases}2\Delta
t_{k},\;&\text{if}\;0\leq\left|1-\rho_{k}\right|\leq\eta_{1},\\\ \Delta
t_{k},\;&\text{else if}\;\eta_{1}<\left|1-\rho_{k}\right|<\eta_{2},\\\
\frac{1}{2}\Delta t_{k},\;&\text{others},\end{cases}$ (13)

where the constants are selected as $\eta_{1}=0.25,\;\eta_{2}=0.75$ according
to our numerical experiments. When $\rho_{k}\geq\eta_{a}$, we accept the trial
step and set

$\displaystyle s_{k}=\frac{\Delta t_{k}}{1+\Delta
t_{k}}s_{k}^{N},\;x_{k+1}=x_{k}+s_{k},$ (14)

where the Newton step $s_{k}^{N}$ is solved by equation (11) and $\eta_{a}$ is
a small positive number such as $\eta_{a}=1.0\times 10^{-6}$. Otherwise, we
discard the trial step and set

$\displaystyle x_{k+1}=x_{k},\;s_{k+1}^{N}=s_{k}^{N},\;F_{k+1}=F_{k}.$ (15)

###### Remark 3

This new time-stepping selection based on the trust-region updating strategy
has some advantages compared to the traditional line search strategy Luo2005 .
If we use the line search strategy and the damped Newton method (3) to track
the trajectory $x(t)$ of the continuous Newton flow (5), in order to achieve
the fast convergence rate in the steady-state phase, the time step
$\alpha_{k}$ of the damped Newton method is tried from 1 and reduced by the
half with many times at every iteration. Since the linear model
$F(x_{k})+J(x_{k})s_{k}^{N}$ may not approximate $F(x_{k}+s_{k}^{N})$ well in
the transient-state phase, the time step $\alpha_{k}$ will be small.
Consequently, the line search strategy consumes the unnecessary trial steps in
the transient-state phase. However, the time-stepping selection method based
on the trust-region updating strategy (12)-(13) can overcome this shortcoming.

For a system of nonlinear equations, the computational time of the Jacobian
matrix is heavy if we update the Jacobian matrix $J(x_{k})$ at every
iteration. In order to save the computational time of the Jacobian evaluation,
similarly to the switching preconditioned technique LXLZ2021 ; LX2021 , we set
$J_{k+1}=J_{k}$ when $F_{k}+J_{k}s_{k}$ approximates $F(x_{k}+s_{k})$ well.
Otherwise, we update $J_{k+1}=J(x_{k+1})$. An effective updating strategy is
give by

$\displaystyle
J_{k+1}=\begin{cases}J_{k},\;\text{if}\;|1-\rho_{k}|\leq\eta_{1}\;\text{and}\;\|s_{k}\|_{\infty}<1,\\\
J(x_{k+1}),\;\text{otherwise},\end{cases}$ (16)

where $\rho_{k}$ is defined by equation (12) and $\eta_{1}=0.25$. In practice,
in order to save the computational time of decomposing the matrix $J_{k+1}$
when $J_{k}$ performs well, i.e. $|1-\rho_{k}|\leq\eta_{1}$, according to the
updating formula (16), we set $J_{k+1}=J_{k}$ in (11).

For a real-world problem, the analytical Jacobian $J(x_{k})$ may not be
offered. Thus, in practice, we replace the Jacobian matrix $J(x_{k})$ with its
difference approximation as follows:

$\displaystyle J(x_{k})\approx\left[\frac{F(x_{k}+\epsilon
e_{1})-F(x_{k})}{\epsilon},\,\ldots,\,\frac{F(x_{k}+\epsilon
e_{n})-F(x_{k})}{\epsilon}\right],$ (17)

where $e_{i}$ represents the unit vector whose elements equal zeros except for
the $i$-th element which equals 1, and the parameter $\epsilon$ can be
selected as $10^{-6}$ according to our numerical experiments.

According to the above discussions, we give the detailed descriptions of the
continuation Newton method and the trust-region updating strategy for finding
a stationary point of $f(x)$ (i.e. a solution of the nonlinear system (2)) in
Algorithm 1. The global convergence and local suplinear convergence of
Algorithm 1 can be found in LXL2021 ; LX2021 .

Algorithm 1 Continuation Newton methods and the trust-region updating strategy
for finding a solution of nonlinear equations (The CNMTr method)

0: Function $F:\;\Re^{n}\to\Re^{n}$, the initial point $x_{0}$ (optional), the
Jacobian $J$ of $F$ (optional), the tolerance $\epsilon$ (optional).

0: An approximation solution $x^{\ast}$ of nonlinear equations.

1: Set the default $x_{0}=\text{ones}(n,\,1)$ and $\epsilon=10^{-6}$ when
$x_{0}$ or $\epsilon$ is not provided by the calling subroutine.

2: Initialize the parameters:
$\eta_{a}=10^{-6},\;\eta_{1}=0.25,\;\gamma_{1}=2,\;\eta_{2}=0.75,\;\gamma_{2}=0.5,\;\text{maxit}=400$.

3: Set $\Delta t_{0}=2$, flag_success_trialstep = 1, $\text{itc}=0,\;k=0$.

4: Evaluate $F_{k}=F(x_{k})$.

5: Set $\rho_{k-1}=0,\;s_{k-1}=0$.

6: while (itc $<$ maxit) do

7: if (flag_success_trialstep == 1) then

8: Set itc = itc + 1.

9: Compute $\text{Res}_{k}=\|F_{k}\|_{\infty}$.

10: if ($\text{Res}_{k}<\epsilon$) then

11: break;

12: end if

13: if
$\left((|1-\rho_{k-1}|>0.25)\;\text{or}\;(\|s_{k-1}\|_{\infty}>1)\right)$ then

14: Evaluate $J_{k}=J(x_{k})$.

15: else

16: Set $J_{k}=J_{k-1}$;

17: end if

18: By solving $J_{k}s_{k}^{N}=-F_{k}$; we obtain the Newton step $s_{k}^{N}$.

19: end if

20: Set $s_{k}={\Delta t_{k}}/{(1+\Delta
t_{k})}\,s_{k}^{N},\;x_{k+1}=x_{k}+s_{k}$.

21: Evaluate $F_{k+1}=F(x_{k+1})$.

22: if $\left(\|F_{k}\|<\|F_{k+1}\|\right)$ then

23: $\rho_{k}=-1$;

24: else

25: Compute the ratio $\rho_{k}$ from equation (12).

26: end if

27: Adjust the time step size $\Delta t_{k+1}$ according to the trust-region
updating strategy (13).

28: if $(\rho_{k}\geq\eta_{a})$ then

29: Accept the trial point $x_{k+1}$. Set flag_success_trialstep = 1.

30: else

31: Set $x_{k+1}=x_{k}$, $F_{k+1}=F_{k}$, $s_{k+1}^{N}=s_{k}^{N}$,
flag_success_trialstep = 0.

32: end if

33: Set $k\longleftarrow k+1$.

34: end while

##  3 The deflation technique

By using Algorithm 1 (i.e., the CNMTr method) in section 2, we can find a
stationary point $x_{1}^{\ast}$ of $f(x)$. Our strategy is to use Algorithm 1
repeatedly, until the stationary points $x_{k}^{\ast}\,(k=1,\,2,\,\ldots,\,K)$
of $f(x)$ are found as many as possible. Then, we select the first $L$
minimizers of $f(x)$ from the set
$\left\\{x_{1}^{\ast}\,\ldots,x_{K}^{\ast}\right\\}$ as the initial seeds of
the genetic algorithm in section 4.

If we directly use Algorithm 1 and the multi-start techniques, a starting
point could easily lie in the basin of attraction of a previously found local
minimum and algorithm 1 might find the same local minimum even if the searches
are initiated from points those are far apart. In order to overcome this
disadvantage, we use the deflation technique BG1971  described as follows.

Assume that $x_{k}^{\ast}$ is a zero point of $G_{k-1}(x)$. Then, we construct
another function $G_{k}(\cdot)$ via eliminating the zero point $x_{k}^{\ast}$
of $G_{k-1}(x)$ as follows:

$\displaystyle G_{k}(x)=\frac{1}{\|x-x_{k}^{\ast}\|}G_{k-1}(x),$ (18)

where $G_{0}(x)=F(x)=\nabla f(x)$. Thus, from equation (18), we obtain

$\displaystyle
G_{k}(x)=\frac{1}{\|x-x_{1}^{\ast}\|\cdots\|x-x_{k}^{\ast}\|}F(x),$ (19)

where $x_{i}^{\ast}\,(i=1:k)$ are the zero points of $F(x)$.

When $F(x_{i}^{\ast})=0$ and $J(x_{i}^{\ast})$ is nonsingular, Brown and
Gearhard BG1971  have proved

$\displaystyle\lim_{j\to\infty}\inf\|G_{k}(x_{j})\|>0$

for any sequence $x_{j}\to x_{i}^{\ast}$, where $G_{k}(x)$ is defined by
equation (19). Actually, we can obtain the following estimation of the lower
bound of $G_{k}(x)$ when $x$ belongs to the neighborhood of $x_{i}^{\ast}$.

###### Lemma 1

Assume that $x_{i}^{\ast}$ is a zero point of the continuously differential
function $F(\cdot)$ and its Jacobian matrix $J(x_{i}^{\ast})$ satisfies

$\displaystyle\|J(x_{i}^{\ast})y\|\geq c_{l}\|y\|,\;\forall\;y\in R^{n}.$ (20)

Then there exists a neighborhood $B_{\delta_{i}}(x_{i}^{\ast})$ of
$x_{i}^{\ast}$ such that

$\displaystyle\|G_{k}(x)\|\geq\frac{c_{l}}{2}\frac{1}{\|x-x_{1}^{\ast}\|\cdots\|x-x_{i-1}^{\ast}\|\|x-x_{i+1}^{\ast}\|\cdots\|x-x_{k}^{\ast}\|},\;\forall
x\in B_{\delta_{i}}(x_{i}^{\ast})$ (21)

holds for all $x\in B_{\delta_{i}}(x_{i}^{\ast})$, where $G_{k}(x)$ is defined
by equation (19). That is to say, $x_{i}^{\ast}$ is not the zero point of
$G_{k}(x)$.

Proof. Since $x_{i}^{\ast}$ is a zero point of $F(x)$, according to the first-
order Taylor expansion, we have

$\displaystyle F(x)$
$\displaystyle=F(x_{i}^{\ast})+\int_{0}^{1}J(x_{i}^{\ast}+t(x-x_{i}^{\ast}))(x-x_{i}^{\ast})dt$
$\displaystyle=\int_{0}^{1}J(x_{i}^{\ast}+t(x-x_{i}^{\ast}))(x-x_{i}^{\ast})dt,$
(22)

where $J(x)=F^{\prime}(x)$ is the Jacobian matrix of $F(x)$. According to the
nonsingular assumption (20) of $J(x_{i}^{\ast})$ and the continuity of
$J(\cdot)$, we know that there exists a neighborhood
$B_{\delta_{i}}(x_{i}^{\ast})=\\{x:\|x-x_{i}^{\ast}\|\leq\delta_{i}\\}$ of
$x_{i}^{\ast}$ such that

$\displaystyle\|J(x)-J(x_{i}^{\ast})\|\leq\frac{c_{l}}{2}$ (23)

holds for all $x\in B_{\delta_{i}}(x_{i}^{\ast})$. From equations (22)-(23),
we have

$\displaystyle\|F(x)-J(x_{i}^{\ast})(x-x_{i}^{\ast})\|=\left\|\int_{0}^{1}(J(x_{i}^{\ast}+t(x-x_{i}^{\ast}))-J(x_{i}^{\ast}))(x-x_{i}^{\ast})dt\right\|$
$\displaystyle\leq\int_{0}^{1}\|J(x_{i}^{\ast}+t(x-x_{i}^{\ast}))-J(x_{i}^{\ast})\|\|x-x_{i}^{\ast}\|dt\leq\frac{c_{l}}{2}\|x-x_{i}^{\ast}\|,\;\forall
x\in B_{\delta_{i}}(x_{i}^{\ast}).$ (24)

Furthermore, from the triangle inequality $\|x-y\|\geq\|x\|-\|y\|$ and the
assumption (20), we have

$\displaystyle\|F(x)-J(x_{i}^{\ast})(x-x_{i}^{\ast})\|\geq\|J(x_{i}^{\ast})(x-x_{i}^{\ast})\|-\|F(x)\|\geq
c_{l}\|x-x_{i}^{\ast}\|-\|F(x)\|.$ (25)

Thus, from equations (24)-(25), we obtain

$\displaystyle\|F(x)\|\geq\frac{c_{l}}{2}\|x-x_{i}^{\ast}\|,\;\forall x\in
B_{\delta_{i}}(x_{i}^{\ast}).$ (26)

By substituting inequality (26) into equation (19), we obtain

$\displaystyle\|G_{k}(x)\|\geq\frac{c_{l}}{2}\frac{1}{\|x-x_{1}^{\ast}\|\cdots\|x-x_{i-1}^{\ast}\|\|x-x_{i+1}^{\ast}\|\cdots\|x-x_{k}^{\ast}\|}>0,\;\forall
x\in B_{\delta_{i}}(x_{i}^{\ast}).$

That is to say, $x_{i}^{\ast}$ is not a zero point of $G_{k}(x)$. ∎

In practice, we need to avoid the overflow or the underflow when we compute
$G_{k}(x)$ and expect that the continuation Newton method (Algorithm 1) can
find the zero points of $F(x)$ as many as possible. Thus, according to our
numerical experiments, we modify $G_{k}(x)$ defined by equation (18) as

$\displaystyle G_{k}(x)=\frac{\alpha_{k}}{\|x-x_{k}^{\ast}\|_{1}}G_{k-1}(x),$
(27)

where $G_{0}(x)=F(x)$ and $\alpha_{k}$ is defined by

$\displaystyle\alpha_{k}=\begin{cases}n,\;\text{if}\;\|x_{k}^{\ast}\|_{1}\leq
10^{-6},\\\ \|x_{k}^{\ast}\|_{1},\;\text{otherwise}.\end{cases}$ (28)

Since $G_{k}(x)$ has the special structure and its sub-differential exists
except for finite points $x_{1}^{\ast},\,\ldots,\,x_{k}^{\ast}$, from
equations (27)-(28), we can compute its sub-differential $\partial G_{k}(x)$
of $G_{k}(x)$ as follows:

$\displaystyle\partial
G_{k}(x)=\frac{\alpha_{1}\cdots\alpha_{k}}{\|x-x_{1}^{\ast}\|_{1}\cdots\|x-x_{k}^{\ast}\|_{1}}\frac{\partial
F(x)}{\partial
x}+\frac{\alpha_{1}\cdots\alpha_{k}}{\|x-x_{1}^{\ast}\|_{1}\cdots\|x-x_{k}^{\ast}\|_{1}}F(x)p(x)^{T}$
$\displaystyle\hskip
5.69054pt=\frac{\alpha_{1}\cdots\alpha_{k}}{\|x-x_{1}^{\ast}\|_{1}\cdots\|x-x_{k}^{\ast}\|_{1}}\left(\frac{\partial
F(x)}{\partial x}+F(x)p(x)^{T}\right),$ (29)

where the column vector function $p(\cdot)$ is defined by

$\displaystyle
p(x)^{T}=-\left(\frac{\text{sgn}\left(x-x_{1}^{\ast}\right)}{\|x-x_{1}^{\ast}\|_{1}}+\cdots+\frac{\text{sgn}\left(x-x_{k}^{\ast}\right)}{\|x-x_{k}^{\ast}\|_{1}}\right),$
(30)

the sign function $\text{sgn}(\cdot)$ of which is defined by

$\displaystyle\text{sgn}(\beta)=\begin{cases}1,\;\text{if}\;\beta>0,\\\
0,\;\text{if}\;\beta=0,\\\ -1,\;\text{otherwise}.\end{cases}$ (31)

Now, by using the strategy (27) of eliminating zero points of $F(x)$
repeatedly, we can find the most of stationary points of $f(x)$. Then, we
select the minimizer from the stationary points of $f(x)$ as the approximation
of the global minimizer of problem (1). We describe its details in Algorithm
2.

###### Remark 4

From equation (27), we know that $F(x_{l})$ may converge to a positive number
when $\lim_{l\to\infty}G_{k}(x_{l})=0$ and
$\lim_{l\to\infty}\|x_{l}-x_{i}^{\ast}\|=\infty$. One of the reasons is that
the Newton method converges to the same solution of $F(x)$ and can not
converge to other solutions of $F(x)$ when it is initialized in the special
region Abbott1977 . Thus, the continuation Newton method with the deflation
technique may not find all the stationary points of $f(x)$ from the same
initial point $x_{0}$. In order to find the stationary point of $f(x)$ as many
as possible, we use Algorithm 2 to find its stationary points from several far
apart initial points sequentially such as $x_{0}=\text{ones}(n,1)$,
$x_{0}=[\text{ones}(\frac{n}{2},1);\,-\text{ones}(\frac{n}{2},1)]$,
$x_{0}=[-\text{ones}(\frac{n}{2},1);\,\text{ones}(\frac{n}{2},1)]$ and
$x_{0}=-\text{ones}(n,1)$. According to our numerical experiments, the
continuation Newton method with the deflation technique (i.e. Algorithm 2) can
find all the stationary points of $f(x)$ when the number of its stationary
points is less than thirty and Algorithm 2 is initialized from those four
points in the most cases.

Algorithm 2 Continuation Newton methods with the deflation technique for
multi-solutions of nonlinear equations (The CNMDT method)

0: The gradient $F(x)=\nabla f(x)$ of $f(x)$, the four initial points
$x_{0}=\text{ones}(n,1)$,
$x_{0}=[\text{ones}(n/2,1);\,-\text{ones}(n/2,\,1)]$,
$x_{0}=[-\text{ones}(n/2,1);\,\text{ones}(n/2,\,1)]$ and
$x_{0}=-\text{ones}(n,\,1)$, the tolerance $\epsilon$ (optional), the Jacobian
$J(x)$ of $F(x)$ (optional).

0: The $K$ approximation zero points $x_{1}^{\ast},\cdots,x_{K}^{\ast}$ of
$F(\cdot)$ and the minimizer $x_{ap}^{\ast}$ of $f(x)$.

1: Set Number_of_Stationary_Points = 0.

2: Call Algorithm 1 to find a zero of $F$ and we denote it as $x_{1}^{\ast}$.

3: Set Number_of_Stationary_Points = Number_of_Statonary_Points +1.

4: Set $G_{0}(x)=F(x)$.

5: Set Flag_of_Fail = 0 and $k=0$.

6: while (Flag_of_Fail == 0) do

7: Call Algorithm 1 to find a zero of $G_{k+1}(\cdot)$ and we denote it as
$x_{k+2}^{\ast}$, where $G_{k+1}(\cdot)$ is defined by equations (27)-(28).

8: if (“A zero of $G_{k+1}$ is obtained by Algorithm 1 successfully.”) then

9: Set Flag_of_Fail = 0.

10: Set Number_of_Stationary_Points = Number_of_Stationary_Points +1.

11: else

12: Set Flag_of_Fail = 1.

13: end if

14: Try again starting at step 6 with $k\longleftarrow k+1$.

15: end while

16: Set $f_{ap}=f(x_{1}^{\ast})$ and $x_{ap}^{\ast}=x_{1}^{\ast}$.

17: for k = 2 : Number_of_Stationary_Points do

18: if $f(x_{k}^{\ast})<f_{ap}$ then

19: Set $f_{a}=f(x_{k}^{\ast})$ and $x_{ap}^{\ast}=x_{k}^{\ast}$.

20: end if

21: end for

##  4 The quasi-genetic evolution and the local search method

By using the continuation Newton method with the deflation technique
(Algorithm 2) in section 3, we can find the most of stationary points of
$f(x)$ and the suboptimal minimizer $x_{ap}^{\ast}$ of $f(x)$ when it has
infinite stationary points. According to our numerical experiments, the global
minimizer is one of the found stationary points or in the convex region
generated by the found stationary points for the most problems. Therefore, we
use the idea of memetic algorithms Moscato1989 ; SGKZ2014  to approach the
global minimum further. Namely, we use the heuristic crossover of the genetic
algorithm to evolve for approaching the global minimum of $f(x)$ from its
stationary points, which are solved by Algorithm 2. We describe it in detail
as follows.

Firstly, we use Algorithm 2 to find the stationary points of $f(x)$ as many as
possible, and we denote these stationary points as
$x_{k}^{\ast}\,(k=1,\,2,\,\ldots,\,K)$. Then, we select the first $L$
minimizers of $f(x)$ from the set
$\left\\{x_{1}^{\ast},\,\ldots,\,x_{K}^{\ast}\right\\}$ and denote these $L$
stationary points as $x_{k}^{0}\,(k=1,\,2,\,\ldots,\,L)$. We set
$S_{0}=\left\\{x_{1}^{0},\,\ldots,\,x_{L}^{0}\right\\}$ and adopt them as the
initial seeds of the quasi-genetic algorithm. When the number of the found
stationary points is less than that of the initial seeds (i.e. $K<L$), we
supplement $L$ determined seeds such as
$\\{\text{zeros}(n,\,1),\;0.1*\bar{S}_{0},\;\bar{S}_{0},\;10*\bar{S}_{0},\;100*\bar{S}_{0},\;1000*\bar{S}_{0},\;\ldots\\}$,
where the set
$\bar{S}_{0}=\\{\text{ones}(n,\,1),\;-\text{ones}(n,\,1),\;[\text{ones}(n/2,\,1);-\text{ones}(n/2,\,1)],\;[-\text{ones}(n/2,\,1)$;
$\text{ones}(n/2,\,1)]\\}$.

Secondly, we select any two elements $x_{i}^{0}$ and $x_{j}^{0}$ from the seed
set $S_{0}$, and use the convex crossover to generate next offspring. Namely,
we set

$\displaystyle\bar{x}_{k}^{1}=\frac{1}{2}\left(x_{i}^{0}+x_{j}^{0}\right),\;i,\,j=1:L,\;k=1:\frac{1}{2}L(L-1).$
(32)

Similarly, we select the first $L$ minimizers of $f(x)$ from the set
$S_{0}\bigcup\left\\{\bar{x}_{1}^{1},\,\ldots,\,\bar{x}_{L(L-1)/2}^{1}\right\\}$
and denote these first $L$ minimizers as $x_{k}^{1}\,(k=1,\,\ldots,\,L)$. We
set $S_{1}=\left\\{x_{1}^{1},\,\ldots,\,x_{L}^{1}\right\\}$ and regard them as
the good population.

Repeatedly, after the $M$-th evolution, we generate the offspring
$\left\\{\bar{x}_{1}^{M},\,\ldots,\,\bar{x}_{L(L-1)/2}^{M}\right\\}$ and
select the first $L$ minimizers of $f(x)$ from the set
$S_{M-1}\bigcup\left\\{\bar{x}_{1}^{M},\,\ldots,\,\bar{x}_{L(L-1)/2}^{M}\right\\}$.
We denote these first $L$ minimizers as $x_{k}^{M}\,(k=1,\,\ldots,\,L)$, set
$S_{M}=\left\\{x_{1}^{M},\,\ldots,\,x_{L}^{M}\right\\}$ and stop evolving. We
denote the minimizer of $f(x)$ in the set $S_{M}$ as $x_{ag}^{\ast}$.

Finally, in order to improve the convergence of the quasi-genetic evolution,
we use the continuation Newton method (i.e. Algorithm 1) from the initial
point $x_{ag}^{\ast}$ to find the solution of $F(x)$ and denote this solution
as $x_{cn}^{\ast}$, where we set $F(x)=\nabla f(x)$ in Algorithm 1. We output
the minimum point $x_{min}^{\ast}$ of $f(x)$ between $x_{cn}^{\ast}$ and
$x_{ag}^{\ast}$, i.e.

$\displaystyle
x_{min}^{\ast}=\arg\min\left\\{f(x_{cn}^{\ast}),\;f(x_{ag}^{\ast})\right\\}.$

We can select the parameter $L=\min\\{20,\;K\\}$ and $M=10$ according to our
numerical experiments. Now, we give the detailed descriptions of the
continuation Newton method with the deflation technique and the quasi-genetic
algorithm for the global optimization problem in Algorithm 3.

Algorithm 3 Continuation Newton methods with the deflation technique and the
quasi-genetic evolution for the global optimization problem (the CNMGE method)

0: The gradient $F(x)=\nabla f(x)$ of $f(x)$, the initial point $x_{0}$
(optional), the tolerance $\epsilon$ (optional), the Jacobian $J(x)$ of $F(x)$
(optional).

0: An approximation minimizer $x_{min}^{\ast}$ of $f(x)$.

1: Call Algorithm 2 to find the stationary points of $f(x)$ as many as
possible and denote these stationary points as
$x_{k}^{\ast}\,(k=1,\,2,\,\ldots,\,K)$.

2: Set $L=20$, $l=0$ and the maximum evolution iterations $M=10$.

3: Select the first $L$ minimizers of $f(x)$ from the set
$\left\\{x_{1}^{\ast}\,\ldots,\,x_{K}^{\ast}\right\\}$ and denote these first
$L$ minimizers as $x_{k}^{0}\,(k=1,\,2,\,\ldots,\,L)$. Set
$S_{0}=\left\\{x_{1}^{0},\,\ldots,\,x_{L}^{0}\right\\}$. When $K<L$, we
supplement $L$ determined seeds such as
$\\{\text{zeros}(n,\,1),\;0.1*\bar{S}_{0},\;\bar{S}_{0},\;10*\bar{S}_{0},\;100*\bar{S}_{0},\;1000*\bar{S}_{0},\;\ldots\\}$,
where
$\bar{S}_{0}=\\{\text{ones}(n,\,1),\;-\text{ones}(n,\,1),\;[\text{ones}(n/2,\,1);-\text{ones}(n/2,\,1)],\;[-\text{ones}(n/2,\,1);\,\text{ones}(n/2,\,1)]\\}$.

4: while ($l<M$) do

5: Set $k=1$.

6: for $i=1:L$ do

7: for $j=i+1:L$ do

8: Set $\bar{x}_{k}^{l+1}=\frac{1}{2}\left(x_{i}^{l}+x_{j}^{l}\right)$;

9: Set $k=k+1$;

10: end for

11: end for

12: Select the first $L$ minimizers of $f(x)$ from the set
$S_{l}\bigcup\left\\{\bar{x}_{1}^{l+1},\,\ldots,\,\bar{x}_{L(L+1)/2}^{l+1}\right\\}$,
and denote these first $L$ minimizers as $x_{i}^{l+1}\,(i=1,\,\ldots,\,L)$.
Set $S_{l+1}=\left\\{x_{1}^{l+1},\,\ldots,\,x_{L}^{l+1}\right\\}$.

13: Set $l=l+1$.

14: end while

15: Select the smallest minimizer $x_{ag}^{\ast}$ of $f(x)$ from $S_{M}$.

16: Set $F(x)=\nabla f(x)$ and call Algorithm 1 with the initial point
$x_{ag}^{\ast}$ to find the solution $x_{cn}^{\ast}$ of $F(x)$.

17: Set
$x_{min}^{\ast}=\arg\min\left\\{f(x_{cn}^{\ast}),\;f(x_{ag}^{\ast})\right\\}$.

##  5 Numerical Experiments

In this section, some numerical experiments are conducted to test the
performance of the continuation Newton method with the deflation technique and
the quasi-genetic evolution (i.e. Algorithm 3, the CNMGE method) for the
global optimization problems. We compare CNMGE with the multi-start algorithm
(the built-in subroutine GlobalSearch.m of the MATLAB R2020a environment
MATLAB ; ULPG2007 ), the branch-and-bound method (Couenne of the state-of-the-
art open-source solver BLLMW2009 ; Couenne2020 ), and two representative
derivative-free algorithms RS2013  (CMA-ES Hansen2006 ; Hansen2010  and MCS
HN1999 ; Neumaier2000 ) for $68$ unconstrained global optimization problems
which can be found in Andrei2008 ; AD2005 ; EVNTBHA2019 ; LM2004 ; MGH1981 ;
SB2020 . The codes are executed by a HP notebook with the Intel quad-core CPU
and 8Gb memory except for Couenne, which is executed by the NEOS server
CMM1998 ; Dolan2001 ; GM1997 ; NEOS2021 .

GlobalSearch is a stable and efficient global optimization solver and is
widely used in engineering field. Couenne is a state-of-the-art open source
solver for global optimization problems. The covariance matrix adaptation
evolution strategy (CMA-ES, its MATLAB code can be downloaded from CMAES2012 )
and the multilevel coordinate search (MCS, its MATLAB code can be downloaded
from MCS2000 ) are two representative derivative-free algorithms. Therefore,
we select these four methods as the basis for comparison. The termination
condition of Algorithm 1 for finding a solution of the nonlinear system (2) is

$\displaystyle\|F(x_{k})\|_{\infty}\leq 10^{-6}.$ (33)

The scales of test problems vary from dimension $1$ to $1000$. The numerical
results are arranged in Tables 1-4. Since the global minima of some test
problems are unknown, we mark their global minima of these problems with the
notation “$\backslash$” in Tables 1-4. For these problems, we regard that the
method fails to find their global minima if its found solutions are greater
than the minimum solutions. The number of CMA-ES’s evolution generations is
set to $10^{5}$. Due to the randomness of CMA-ES, we repeatedly calculate
every problem ten times with CMA-ES, and select the minimum as its optimal
value in Tables 1-4, respectively. The computational time of CNMGE,
GlobalSearch, Couenne, CMA-ES and MCS is illustrated by Figures 1-2.

When the unconstrained optimization problem (1) has many local minima, the
numerical tests show that Algorithm 2 can find at least $30$ stationary points
of $f(x)$, which are useful to make Algorithm 3 find the global minimum of
$f(x)$ for the difficult problem. This result can be illustrated by Problem 1.
Problem 1 LM2004  is a molecular potential energy problem and its objective
function has the following form:

$\displaystyle f(\omega)=\sum_{i=1}^{n}\left(1\ +\ cos(3\omega_{i})\ +\
\frac{(-1)^{i}}{\sqrt{10.60099896-4.141720682(cos\omega_{i})}}\right),$ (34)

where $\omega_{i}$ is the torsion angle and $n+3$ is the number of atoms or
beads in the given system. The number of its local minimizers increases
exponentially with the size of the problem, which characterizes the principal
difficult in minimizing molecular potential energy functions for the
traditional global optimization methods. CNMGE can find 91 stationary points
of problem (34) with $n=1000$ and its global minimum
$f(\omega^{\ast})=-41.118303$ at
$\omega^{\ast}=(1.039195,\,3.141593,\,\ldots,\,1.039195,\,3.141593)$. However,
for problem 1 with 50 atoms, the traditional global optimization methods such
as the interval analysis method kearfott1996 ; LM2004  or the multi-start
method MATLAB ; ULPG2007  are difficult to find its global minimum. Couenne,
CMA-ES and MCS all fail to solve this problem.

For problem $44$, reference SB2020  gives its global minimum value (i.e.
$-9.596407E+02$) at $x=(512,\;404.2319)$, when it is evaluated on the interval
$x_{i}\in[-512,\;512]\;(i=1,\,2)$. Since these three methods do not consider
the constraints $x_{i}\in[-512,\;512]\;(i=1,\,2)$ and regard this problem as
an unconstrained optimization problem, CNMGE can find a smaller value of
$f(x)$ than that of reference SB2020 .

From Tables 1-4, we find that CNMGE performs well for unconstrained global
optimization problems and finds their global minima. From Tables 3-4, we find
that GlobalSearch, Couenne, CMA-ES and MCS all perform well for small-scale
problems. However, from Tables 1-2, we find that GlobalSearch, CMA-ES and MCS
all fail to find the global minima of the large-scale problems about 17/34.
Couenne can solve the most problems efficiently and it can not solve the
large-scale problems efficiently about 12/34 (Couenne fails to solve 3 large-
scale problems and 9 large-scale problems are solved by it about 8 hours).
Therefore, CNMGE is more robust than GlobalSearch, Couenne, CMA-ES and MCS to
find the global minimum of the unconstrained optimization problem.

From Tables 1-2 and Figure 1, we find that the computational time of CNMGE is
slightly more than that of GlobalSearch and Couenne for the large-scale
unconstrained optimization problem. From Tables 1-2 and Figure 1, we find that
the computational time of CNMGE and GlobalSearch is about $1/10$ of that of
CMA-ES and MCS for the large-scale problem. From Tables 3-4 and Figure 2, we
find that the computational time of CNMGE, CMA-ES and MCS is about $1/5$ of
that of GlobalSearch for the small-scale problem. Therefore, the computational
time of CNMGE is not serious and far less than that of CMA-ES and MCS for the
large-scale problem, since CNMGE only consumes several seconds to three
hundred seconds for large-scale problems.

Table 1: Numerical results of large-scale problems (no.1-17).

Problem | CNMGE | GlobalSearch | Couenne | CMA-ES | MCS | Global Minimum
---|---|---|---|---|---|---
min $f(x)$ | min $f(x)$ | min $f(x)$ | min $f(x)$ | min $f(x)$
| (CPU time (s)) | (CPU time (s)) | (CPU time (s)) | (CPU time (s)) | (CPU time (s)) |
1\. Molecular Energy | -4.1118E+01 | -2.8761E-01 | -4.1118E+01 | 1.0008E+01 | -3.6903E-01 | -4.1118E+01
Problem (n = 1000) LM2004  | (2.02E+01) | (1.82E+01) | (2.68E+04) | (2.98E+02) | (2.01E+02) |
|  | (failed) | (too long) | (failed) | (failed) |
2\. Ackley Function | 1.2690E-07 | 8.8817E-16 | 6.7078E-06 | 3.6253 | 0 | 0
(n = 1000) SB2020  | (7.13E+01) | (1.81E+01) | (1.77E+03) | (7.71E+01) | (1.28E+04) |
|  |  |  | (failed) | (too long) |
3\. Levy Function | 1.1477E-15 | 0.3651 | 1.1997E-15 | 1.4998E-32 | 2.2419E-09 | 0
(n = 1000) SB2020  | (2.08E+02) | (3.74E+00) | (1.43E-01) | (1.51E+02) | (1.10E+02) |
|  | (failed) |  |  |  |
4\. Schwefel Function | 1.2728E-02 | 4.7903E+03 | 3.1688E+03 | 1.7096E+05 | 1.1844E+05 | 1.2728E-02
(n = 1000) SB2020  | (3.43E+02) | (4.34E+01) | (2.42E+04) | (2.45E+02) | (2.60E+02) |
|  | (failed) | (failed) | (failed) | (failed) |
5\. Rastrigin Function | 0 | 0 | 8.5893E-26 | 1.0000E+03 | 0 | 0
(n = 1000) SB2020  | (5.14E+00) | (3.48E+00) | (2.82E+01) | (1.84E+02) | (1.25E+02) |
|  |  |  | (failed) |  |
6\. Styblinski-Tang | -3.9166E+04 | -2.5013e+04 | -3.9166E+04 | -8.1869E+03 | -3.9166E+04 | -3.9166E+04
Function (n = 1000) SB2020  | (6.58E+01) | (4.52E+00) | (2.55E+04) | (1.69E+02) | (3.17E+02) |
|  | (failed) | (too long) | (failed) |  |
7\. Trid Function | -1.6717E+08 | -2.4920e+04 | -1.6716E+08 | -9.9900E+02 | -4.2784E+05 | -1.6716E+08
(n = 1000) SB2020  | (4.86E+01) | (4.34E+00) | (2.64E+04) | (1.84E+02) | (5.45E+01) |
|  | (failed) | (too long) | (failed) | (failed) |
8\. Sum Squares Function | 0 | 0 | 3.7885E-44 | 2.3517E-16 | 0 | 0
(n = 1000) SB2020  | (5.86E+00) | (3.84E+00) | (4.14E-02) | (6.20E+01) | (1.26E+04) |
|  |  |  |  | (too long) |
9\. Sphere Function | 0 | 3.9800E-11 | 1.4569E-31 | 2.9222E-16 | 0 | 0
(n = 1000) SB2020  | (4.52E-01) | (3.22E+00) | (2.66E-03) | (5.91E+01) | (1.26E+04) |
|  |  |  |  | (too long) |
10\. Rotated Hyper-Ellipsoid | 0 | 0 | 2.1437E-45 | 3.7091E-16 | 0 | 0
Function (n = 1000) SB2020  | (9.57E+00) | (1.54E+01) | (3.26E-01) | (6.64E+02) | (1.43E+04) |
|  |  |  |  | (too long) |
11\. Zakharov Function | 0 | 0 | 1.3612E-14 | 1.5541E+04 | 1.7016E+04 | 0
(n = 1000) SB2020  | (1.92E+01) | (5.17E+00) | (3.20E+00) | (1.61E+02) | (1.44E+02) |
|  |  |  | (failed) | (failed) |
12\. Dixon-Price Function | 9.0315E-14 | 7.2931E-01 | 1.0905E-15 | 6.6667E-01 | 6.6667E-01 | 0
(n = 1000) SB2020  | (7.67E+00) | (4.03E+00) | (2.13E-01) | (9.77E+01) | (8.45E+02) |
|  | (failed) |  | (failed) | (failed) |
13\. Rosenbrock Function | 0 | 3.0630E+01 | 9.8899E+02 | 9.7702E+02 | 2.1796E+02 | 0
(n = 1000) MGH1981 ; SB2020  | (7.83E+00) | (3.82E+00) | (1.07E+02) | (1.77E+02) | (2.18E+02) |
|  | (failed) | (failed) | (failed) | (failed) |
14\. Powell Function | 0 | 0 | 5.6167E-10 | 4.9855E-03 | 2.1627E+01 | 0
(n = 1000) MGH1981 ; SB2020  | (5.60E+01) | (4.72E+00) | (9.05E-02) | (2.63E+02) | (2.05E+02) |
|  |  |  |  | (failed) |
15\. Quartic With Noise | 5.0000E-01 | 5.0000E-01 | 5.0000E-01 | 5.0000E-01 | 5.0000E-01 | 5.0000E-01
Function (n = 1000) AD2005  | (7.23E+00) | (3.89E+00) | (9.90E-02) | (5.75E+01) | (1.27E+04) |
|  |  |  |  | (too long) |
16\. Schubert Function | -1.2031E+04 | -9.4947E+03 | -1.2031E+04 | -4.2396E+03 | -2.6399E+03 | -1.2031E+04
(n = 1000) AD2005  | (2.13E+02) | (6.30E+01) | (2.62E+04) | (1.66E+02) | (3.34E+02) |
|  | (failed) | (too long) | (failed) | (failed) |
17\. Raydan 1 Function | 5.0050E+04 | 5.0050E+04 | 5.0050E+04 | 5.0050E+04 | 5.0050E+04 | $\backslash$
(n = 1000) Andrei2008  | (7.81E-01) | (1.03E+01) | (1.15E-01) | (2.83E+02) | (1.22E+04) | $\backslash$
|  |  |  |  | (too long) |

Table 2: Numerical results of large-scale problems (no.18-34).

Problem | CNMGE | GlobalSearch | Couenne | CMA-ES | MCS | Global Minimum
---|---|---|---|---|---|---
min $f(x)$ | min $f(x)$ | min $f(x)$ | min $f(x)$ | min $f(x)$
| (CPU time (s)) | (CPU time (s)) | (CPU time (s)) | (CPU time (s)) | (CPU time (s)) |
18\. Raydan 2 Function | 1.0000E+03 | 1.0000E+03 | 1.0000E+03 | 1.3083E+03 | 1.0000E+03 | $\backslash$
(n = 1000) Andrei2008  | (1.28E+00) | (7.60E+00) | (9.68E-02) | (2.93E+02) | (1.31E+04) | $\backslash$
|  |  |  | (failed) | (too long) |
19\. Extended Tridiagonal 1 | 3.0726E-07 | 2.9700E+02 | 9.9721E+02 | 1.0000E+03 | 4.0306E-01 | $\backslash$
Function(n = 1000) Andrei2008  | (1.19E+01) | (4.51E+00) | (5.61E+02) | (2.66E+02) | (1.59E+02) | $\backslash$
|  | (failed) | (failed) | (failed) | (failed) |
20\. Extended quadratic penalty | 3.9900E+03 | 3.9962E+03 | 3.9900E+03 | 3.9900E+03 | 3.9912E+03 | $\backslash$
QP1 Function (n = 1000) Andrei2008  | (9.56E+00) | (3.30E+00) | (2.58E+04) | (1.05E+02) | (7.88E+01) | $\backslash$
|  | (failed) | (too long) |  |  |
21\. Extended quadratic penalty | 6.5901E-21 | 1.0000E+04 | 3.6937E-22 | 1.3031E-01 | 6.6117E+02 | $\backslash$
QP2 Function (n = 1000) Andrei2008  | (1.24E+01) | (3.66E+00) | (2.48E+00) | (6.64E+01) | (6.79E+01) | $\backslash$
|  | (failed) |  | (failed) | (failed) |
22\. Quadratic QF2 Function | -1.0001E+00 | -1.0001E+00 | -1.0001E+00 | -1.0001E+00 | -1.0001E+00 | $\backslash$
(n = 1000) Andrei2008  | (1.24E+01) | (3.66E+00) | (4.10E+01) | (5.25E+01) | (7.79E+01) | $\backslash$
23\. Extended PSC1 Function | 3.8660E+02 | 5.0000E+02 | 3.8660E+02 | 3.8660E+02 | 4.0143E+02 | $\backslash$
(n = 1000) Andrei2008  | (5.72E+00) | (3.32E+00) | (2.82E+04) | (9.50E+01) | (9.22E+01) | $\backslash$
|  | (failed) | (too long) |  | (failed) |
24\. Extended BD1 Function | 0 | 0 | 5.0327E-23 | 1.0616E-15 | 1.2762E-23 | $\backslash$
(n = 1000) Andrei2008  | (2.45E+01) | (2.50E+00) | (3.49E-01) | (9.98E+01) | (9.52E+01) | $\backslash$
25\. Extended Cliff Function | 9.9893E+01 | 2.3779E+02 | 9.9893E+01 | 1.0179E+02 | 2.9165E+02 | $\backslash$
(n = 1000) Andrei2008  | (9.25E+00) | (5.44E+00) | (1.85E+01) | (2.44E+02) | (7.05E+01) | $\backslash$
|  | (failed) |  |  | (failed) |
26\. Perturbed quadratic diagonal | 0 | 0 | 1.1750E-42 | 1.3102E-14 | 0 | $\backslash$
Function (n = 1000) Andrei2008  | (6.53E+00) | (3.17E+00) | (1.23E+00) | (8.59E+01) | (1.27E+04) | $\backslash$
|  |  |  |  | (too long) |
27\. Extended Hiebert Function | 6.3646E-18 | 1.2500E+12 | 4.6675E-21 | 1.1715E+08 | 6.5817E+06 | $\backslash$
(n = 1000) Andrei2008  | (1.16E+01) | (2.21E+01) | (1.91E-01) | (1.22E+02) | (1.41E+02) | $\backslash$
|  | (failed) |  | (failed) | (failed) |
28\. Extended TET Function | 1.2796E+03 | 1.2796E+03 | 1.2796E+03 | 1.2796E+03 | 1.2797E+03 | $\backslash$
(n = 1000) Andrei2008  | (1.04E+01) | (3.65E+00) | (2.67E+04) | (1.75E+02) | (2.17E+02) | $\backslash$
|  |  | (too long) |  |  |
29\. Diagonal 1 Function | -2.7068E+06 | -2.2288E+06 | -2.7068E+06 | -2.7068E+06 | -2.6997E+06 | $\backslash$
(n = 1000) Andrei2008  | (4.31E+00) | (3.54E+00) | (5.33E+01) | (2.52E+02) | (7.52E+01) | $\backslash$
|  | (failed) |  |  | (failed) |
30\. Diagonal 3 Function | -5.0050E+05 | -4.8990E+05 | -5.0050E+05 | -4.8991E+05 | -4.9623E+05 | $\backslash$
(n = 1000) Andrei2008  | (2.59E+01) | (1.83E+01) | (7.34E-02) | (2.17E+02) | (2.79E+02) | $\backslash$
|  | (failed) |  | (failed) | (failed) |
31\. Diagonal 5 Function | 6.9315E+02 | 6.9315E+02 | 6.9315E+02 | 6.9315E+02 | 6.9315E+02 | $\backslash$
(n = 1000) Andrei2008  | (7.88E-01) | (1.16E+01) | (2.81E+04) | (1.65E+02) | (1.25E+04) | $\backslash$
|  |  | (too long) |  | (too long) |
32\. Extended Maratos Function | -5.0031E+02 | -5.0026E+02 | -5.0031E+02 | -1.5082E+02 | -5.0031E+02 | $\backslash$
(n = 1000) Andrei2008  | (7.71E+00) | (3.63E+00) | (1.60E+01) | (1.74E+02) | (8.24E+01) | $\backslash$
|  |  |  | (failed) |  |
33\. EG2 Function | -9.9950E+02 | -9.2965E+02 | -9.9950E+02 | -8.0038E+02 | -9.9917E+02 | $\backslash$
(n = 1000) Andrei2008  | (1.08E+02) | (2.71E+00) | (5.76E-02) | (6.47E+01) | (1.01E+02) | $\backslash$
|  | (failed) |  | (failed) |  |
34\. SINQUAD Function | 0 | 0 | 2.4405E-13 | 0 | 5.5349E-01 | $\backslash$
(n = 1000) Andrei2008  | (2.02E+01) | (3.30E+00) | (1.12E+00) | (2.05E+02) | (1.70E+02) | $\backslash$
|  |  |  |  | (failed) |

Table 3: Numerical results of small-scale problems (no.35-55).

Problem | CNMGE | GlobalSearch | Couenne | CMA-ES | MCS | Global Minimum
---|---|---|---|---|---|---
min $f(x)$ | min $f(x)$ | min $f(x)$ | min $f(x)$ | min $f(x)$
| (CPU time (s)) | (CPU time (s)) | (CPU time (s)) | (CPU time (s)) | (CPU time (s)) |
35\. Griewank Function | 1.9318E-14 | 1.7187E-12 | 1.5929E-26 | 6.6613E-16 | 0 | 0
(n = 10) MGH1981 ; SB2020  | (5.90E-01) | (2.66E+00) | (1.14E-03) | (1.33E+00) | (6.38E-01) |
36\. Levy Function N. 13 | 2.5205E-26 | 1.0981E-01 | 2.8774E-16 | 1.3498E-31 | 1.3498E-31 | 0
(n = 2) SB2020  | (5.67E-01) | (2.85E+00) | (8.48E-03) | (6.78E-01) | (7.01E-01) |
|  | (failed) |  |  |  |
37\. Hosaki Function | -2.3458E+00 | -1.1277E+00 | -2.3458E+00 | -2.3458E+00 | -2.3458E+00 | -2.3458E+00
(n = 2) AD2005  | (5.67E-01) | (2.62E+00) | (2.80E+04) | (5.70E-01) | (3.26E-01) |
|  | (failed) | (too long) |  |  |
38\. Beale Function | 4.7472E-14 | 0 | 5.0180E-23 | 1.3347E-17 | 1.5473E-22 | 0
(n = 2) MGH1981 ; SB2020  | (6.70E-01) | (3.80E+00) | (2.68E-03) | (6.14E-01) | (2.82E-01) |
39\. Easom Function | -1.0000E+00 | -1.0000E+00 | -1.0000E+00 | -3.0308E-05 | -1.0000E+00 | -1.0000E+00
(n = 2) SB2020  | (3.79E-01) | (2.64E+00) | (2.55E-03) | (3.71E-02) | (6.13E-01) |
|  |  |  | (failed) |  |
40\. Price Function | 0 | 1.1149E-12 | 1.1663E-37 | 5.5711E-13 | 0 | 0
(n = 2) AD2005  | (2.30E-01) | (5.21E+00) | (1.46E-03) | (9.81E-01) | (1.42E-01) |
41\. Branin Function | 3.9789E-01 | 3.9789E-01 | 3.9789E-01 | 3.9789E-01 | 3.9789E-01 | 3.9789E-01
(n = 2) SB2020  | (2.64E-01) | (2.21E+00) | (8.57E-03) | (6.35E-01) | (4.91E-01) |
42\. Trecanni Function | 0 | 1.8313E-15 | 6.0525E-50 | -3.5527E-15 | 0 | 0
(n = 2) AD2005  | (2.64E-01) | (2.43E+00) | (8.27E-01) | (5.41E-01) | (4.46E-01) |
43\. Booth Function | 1.3604E-15 | 4.1193E-15 | 0 | 1.8582E-18 | 0 | 0
(n = 2) SB2020  | (2.25E-01) | (2.45E+00) | (1.05E-03) | (5.49E-01) | (2.48E-01) |
44\. Matyas Function | 0 | 6.9378E-16 | 0 | 7.9654E-18 | 0 | 0
(n = 2) SB2020  | (2.25E-01) | (2.45E+00) | (1.05E-03) | (5.74E-01) | (9.18E-02) |
45\. McCormick Function | -1.9132E+00 | -1.9132E+00 | -1.9132E+00 | -1.9132E+00 | -1.9132E+00 | -1.9132E+00
(n = 2) SB2020  | (2.45E-01) | (2.99E+00) | (7.47E-03) | (5.37E-01) | (2.30E-01) |
46\. Power Sum Function | 3.0180E-22 | 2.0472E-07 | 4.2680E-22 | 5.4319E-12 | 1.4325E-06 | 0
(n = 4) SB2020  | (2.06E-01) | (4.19E+00) | (4.76E-03) | (4.46E+00) | (1.07) |
47\. Colville Function | 0 | 5.3436E-13 | 0 | 0 | 3.8343E-19 | 0
(n = 4) SB2020  | (3.49E-01) | (3.80E+00) | (2.39E-01) | (1.09E+00) | (5.95E-01) |
48\. Schaffer Function N. 2 | 1.6653E-16 | 0 | 0 | 2.9258E-01 | 0 | 0
(n = 2) SB2020  | (2.49E-01) | (4.17E+00) | (1.10E-03) | (1.32E+00) | (6.62E-01) |
49\. Bohachevsky Function | 0 | 0 | -5.5511E-17 | 0 | 0 | 0
(n = 2) SB2020  | (1.57E+00) | (2.43E+00) | (2.43E-02) | (6.54E-01) | (1.66E-01) |
50\. Three-Hump Camel | 0 | 3.2640E-16 | 1.4758E-18 | 7.5064E-18 | 0 | 0
Function (n = 2) SB2020  | (3.84E-01) | (2.32E+00) | (2.10E-01) | (5.50E-01) | (1.15E-01) |
51\. Six-Hump Camel | -1.0316E+00 | -1.0316E+00 | -1.0316E+00 | -1.0316E+00 | -1.0316E+00 | -1.0316E+00
Function (n = 2) SB2020  | (2.57E-01) | (2.29E+00) | (9.38E-02) | (6.07E-01) | (3.74E-01) |
52\. Drop-Wave Function | -1.0000E+00 | -1.0000E+00 | -1.0000E+00 | -9.8389E-01 | -1.0000E+00 | -1.0000E+00
(n = 2) SB2020  | (8.42E-01) | (3.12E+00) | (2.59E-02) | (7.34E-01) | (5.65E-01) |
53\. Perm Function 0, d, $\beta$ | 0 | 6.3025E+02 | 5.3868E-10 | 3.0872E-16 | 1.0130E+02 | 0
(n = 4) SB2020  | (8.72E-01) | (5.79E+00) | (2.10E-01) | (1.11E+00) | (9.28E-01) |
|  | (failed) |  |  | (failed) |
54\. Hartmann 3-D Function | -3.8628E+00 | -3.8628E+00 | -3.8628E+00 | -3.8628E+00 | -3.8628E+00 | -3.8628E+00
(n = 3) SB2020  | (3.75E-01) | (2.77E+00) | (4.50E-01) | (8.20E-01) | (3.93E-01) |
55\. Trefethen 4 Function | -3.3069E+00 | -2.8376E+00 | -3.3069E+00 | -3.2965E+00 | -2.7612E+00 | -3.3069E+00
(n = 2) AD2005  | (2.80E+00) | (5.91E+00) | (7.88E-01) | (8.23E-01) | (7.37E-01) |

Table 4: Numerical results of small-scale problems (no.56-68).

Problem | CNMGE | GlobalSearch | Couenne | CMA-ES | MCS | Global Minimum
---|---|---|---|---|---|---
min $f(x)$ | min $f(x)$ | min $f(x)$ | min $f(x)$ | min $f(x)$
| (CPU time (s)) | (CPU time (s)) | (CPU time (s)) | (CPU time (s)) | (CPU time (s)) |
56\. Zettl Function | -3.7912E-03 | -3.7912E-03 | -3.7912E-03 | -3.7912E-03 | -3.7912E-03 | -3.7912E-03
(n = 2) AD2005  | (2.85E-01) | (2.36E+00) | (1.80E-02) | (6.59E-01) | (2.57E-01) |
57\. Exp2 Function | 3.6945E-14 | 4.8793E-14 | 8.0737E-20 | 3.0512E-18 | 2.8282E-21 | 0
(n = 2) AD2005  | (2.51E-01) | (2.81E+00) | (7.61E-03) | (6.26E-01) | (2.44E-01) |
58\. Hansen Function | -1.7654E+02 | -1.7654E+02 | -1.7654E+02 | -1.7654E+02 | -1.4548E+02 | -1.7654E+02
(n = 2) AD2005  | (2.95E-01) | (2.98E+00) | (9.89E-01) | (8.64E-01) | (7.22E-01) |
|  |  |  |  | (failed) |
59\. Schaffer Function N. 4 | 2.9258E-01 | 2.9258E-01 | 2.9258E-01 | 2.9387E-01 | 4.9289E-01 | 2.9258E-01
(n = 2) SB2020  | (2.17E+00) | (4.60E+00) | (1.49E-01) | (9.21E-01) | (9.21E-01) |
60\. Holder Table Function | -1.9209E+01 | -1.9209E+01 | -1.9209E+01 | -1.9209E+01 | -1.9209E+01 | -1.9209E+01
(n = 2) SB2020  | (5.38E-01) | (3.36E+00) | (3.81E-01) | (6.05E-01) | (5.86E-01) |
61\. Gramacy $\&$ Lee (2012) | -8.6901E-01 | -8.6901E-01 | -8.6901E-01 | -8.6901E-01 | -6.6328E-01 | -8.6901E-01
Function (n = 1) SB2020  | (1.41E+00) | (5.49E-01) | (2.31E-03) | (5.67E-01) | (4.56E-01) |
62\. Eggholder Function | -2.0408E+03 | -9.3533E+02 | -9.5964E+02 | -9.3432E+02 | -5.5737E+02 | -9.5964E+02
(n = 2) AD2005  | (3.53E-01) | (3.53E+00) | (2.98E-02) | (8.48E-01) | (4.70E-01) |
|  | (failed) |  | (failed) | (failed) |
63\. Michalewicz Function | -1.8013E+00 | -1.8013E+00 | -1.8013E+00 | -1.8013E+00 | -1.8013E+00 | -1.8013E+00
(n = 2) SB2020  | (4.13E-01) | (4.32E+00) | (3.62E-01) | (5.97E-01) | (3.11E-01) |
64\. Box-Betts Exponential | 0 | 6.8000E-17 | 6.3148E-16 | 1.1671E-18 | 0 | 0
Quadratic Function (n = 3) AD2005  | (4.70E-01) | (2.76E+00) | (6.52E-03) | (8.48E-01) | (1.59E-01) |
65\. Cross-in-Tray Function | -2.0626E+00 | -2.0626E+00 | -2.0626E+00 | -2.0626E+00 | -1.8899E+00 | -2.0626E+00
(n = 2) SB2020  | (2.75E-01) | (3.88E+00) | (7.46E-04) | (5.69E-01) | (3.40E-01) |
66\. Himmeblau Function | 1.8778E-16 | 0 | 2.2376E-25 | 1.0669E-18 | 3.4597E-20 | 0
(n = 2) SB2020  | (1.65E+00) | (2.43E+00) | (2.32E-03) | (5.95E-01) | (4.00E+00) |
67\. Forrester et al. (2008) | -6.0207E+00 | -6.0207E+00 | -6.0207E+00 | -6.0207E+00 | -6.0207E+00 | -6.0207E+00
Function (n = 1) SB2020  | (2.70E-01) | (2.75E+00) | (6.78E-02) | (4.10E-01) | (2.65E-01) |
68\. Goldstein-Price | 3.0000E+00 | 3.0000E+00 | 3.0000E+00 | 3.0000E+00 | 3.0000E+00 | 3.0000E+00
Function (n = 2) SB2020  | (4.73E-01) | (2.71E+00) | (3.61E+00) | (7.32E-01) | (3.39E-01) |

Figure 1: The computational time of CNMGE, GlobalSearch, Couenne, CMA-ES and
MCS for the large-scale problems. Figure 2: The computational time of CNMGE,
GlobalSearch, Couenne, CMA-ES and MCS for the small-scale problems.

##  6 Conclusions

Based on the deflation technique and the genetic evolution, we give an
efficient continuation Newton for solving the global minimum of an
unconstrained optimization problem. The numerical results show that CNMGE
(Algorithm 3) performs well to solve the large-scale problems, especially the
problems of which are difficult to be solved by the known global optimization
methods. Therefore, CNMGE can be regarded as an alternative workhorse for the
global optimization problems.

## Acknowledgments

The authors are grateful to Prof. Jonathan Eckstein and Prof. Nick Sahinidis
for their comments and suggestions which greatly improve presentation of this
paper.

## Declarations

Funding: This work was supported in part by Grant 61876199 from National
Natural Science Foundation of China, and Grant YJCB2011003HI from the
Innovation Research Program of Huawei Technologies Co., Ltd..

Conflicts of interest/Competing interests: Not applicable.

Availability of data and material (data transparency): If it is requested, we
will provide the test data.

Code availability (software application or custom code):
https://teacher.bupt.edu.cn/luoxinlong/zh_CN/zzcg/41406/list/index.htm.

## References

  * (1) Abbott, J. P.: _Numerical continuation methods for nonlinear equations and bifurcation problems_ , Ph.D. Thesis, Computer Center, Australian National University, 1977.
  * (2) Andrei, N.: _An unconstrained optimization test functions collection_ , Advanced Modeling and Optimization 10 (2008), 147-161.
  * (3) Adorio, E. P., Diliman, U. P.: _MVF-Multivariate test functions library in C for unconstrained global optimization_ , available at http://www.geocities.ws/eadorio/mvf.pdf, 2005.
  * (4) Allgower, E. L., Georg, K.: _Introduction to Numerical Continuation Methods_ , SIAM, Philadelphia, PA, 2003.
  * (5) Ascher, U. M., Petzold, L. R.: _Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations_ , SIAM, Philadelphia, PA, 1998.
  * (6) Axelsson, O., Sysala, S.: _Continuation Newton methods_ , Comput. Math. Appl. 70 (2015), 2621-2637.
  * (7) Belotti, P., Lee, J., Liberti, L., Margot, F., W$\ddot{a}$chter, A.: _Branching and bounds tightening techniques for non-convex MINLP_ , Optim. Methods Softw. 24 (2009), 597-634.
  * (8) Boender, C. G. E.: _Bayestian stopping rules for multistart global optimization methods_ , Math. Program. 37 (1987), 59-80.
  * (9) Branin, F. H.: _Widely convergent method for finding multiple solutions of simultaneous nonlinear equations_ , IBM J. Res. Dev. 16 (1972), 504-521.
  * (10) Brown, K. M., Gearhart, W. B.: _Deflation techniques for the calculation of further solutions of a nonlinear system_ , Numer. Math. 16 (1971), 334-342.
  * (11) Conn, A. R., Gould, N., Toint, Ph. L.: _Trust-Region Methods_ , SIAM, Philadelphia, PA, 2000.
  * (12) Couenne, _a solver for non-convex MINLP problems_ , available at https://www.coin-or.org/Couenne/, February 2020.
  * (13) CMA-ES, _The covariance matrix adaptation evolution strategy_ , available at http://www.cmap.polytechnique.fr/~nikolaus.hansen/cmaes.m, 2012.
  * (14) Czyzyk, J., Mesnier, M. P., Moré, J. J.: _The NEOS Server_ , IEEE Comput. Sci. Eng. 5 (1998), 68-75.
  * (15) Davidenko, D. F.: _On a new method of numerical solution of systems of nonlinear equations_ (in Russian), Dokl. Akad. Nauk SSSR 88 (1953), 601-602.
  * (16) Deuflhard, P.: _Newton Methods for Nonlinear Problems: Affine Invariance and Adaptive Algorithms_ , Springer-Verlag, Berlin, 2004.
  * (17) Dolan, E.: _The NEOS Server 4.0 Administrative Guide_ , Technical Memorandum ANL/MCS-TM-250, Mathematics and Computer Science Division, Argonne National Laboratory, 2001\.
  * (18) Deuflhard, P., Pesch, H. J., Rentrop, P.: _A modified continuation method for the numerical solution of nonlinear two-point boundary value problems by shooting techniques_ , Numer. Math. 26 (1975), 327-343.
  * (19) Dennis, J. E., Schnabel, R. B.: _Numerical Methods for Unconstrained Optimization and Nonlinear Equations_ , SIAM, Philadelphia, PA, 1996.
  * (20) Dong, H. C., Song, B. W., Dong, Z. M., Wang, P.: _Multi-start space reduction (MSSR) surrogate-based global optimization method_ , Struct. Multidisc. Optim. 54 (2016), 907-926.
  * (21) Elhara, O., Varelas, K., Nguyen, D., Tusar, T., Brockhoff, D., Hansen, N., Auger, A.: _COCO: The large scale black-box optimization benchmarking (bbob-largescale) Test Suite_ , arXiv preprint available at https://arxiv.org/abs/1903.06396, 2019.
  * (22) Gropp, W., Moré, J. J.: _Optimization Environments and the NEOS Server_ , Approximation Theory and Optimization, M. D. Buhmann and A. Iserles, eds., Cambridge University Press, 1997.
  * (23) Gottlieb, S., Shu, C. W.: _Total variation diminishing Runge-Kutta schemes_ , Math. Comput. 67, 73-85 (1998).
  * (24) Gottlieb, S., Shu, C. W., Tadmor, E.: _Strong stability preserving high order time discretization methods_ , SIAM Rev. 43, 89-112 (2001).
  * (25) Golub, G. H., Van Loan, C. F.: _Matrix Computation_ (4th ed.), The John Hopkins University Press, Baltimore, 2013.
  * (26) Hansen, N.: _The CMA evolution strategy: a comparing reivew_ , In: Lozano, J. A., Larranaga, P., Inza, I., Bengoetxea, E. (eds.) Towards a New Evolutionary Computation. Advances on Estimation of Distribution Algorithms, Springer, Berlin, 75-102 (2006).
  * (27) Hansen, N.: _The CMA Evolution Strategy: A Tutorial_ , available at https://arxiv.org/abs/1604.00772, 2010.
  * (28) Hart, W. E.: _Adaptive Global Optimization with Local Search_ , Ph.D. dissertation, University of California, San Diego, CA, USA, 1994.
  * (29) Higham, D. J.: _Trust region algorithms and timestep selection_ , SIAM J. Numer. Anal. 37 (1999), 194-210.
  * (30) Huyer, W., Neumaier, A.: _Global optimization by multilevel coordinate search_ , J. Glob. Optim. 14, 331-355 (1999).
  * (31) Hairer, E., Wanner, G.: _Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems_ , 2nd ed., Springer-Verlag, Berlin, 1996.
  * (32) Kearfott, R. B.: _Rigorous Global Search: Continuous Problems_ , Nonconvex Optimization and Applications, Kluwer Academic, Dordrecht, 1996.
  * (33) Kelley, C. T.: _Solving Nonlinear Equations with Newton’s Method_ , SIAM, Philadelphia, PA, 2003.
  * (34) Kelley, C. T.: _Numerical methods for nonlinear equations_ , Acta Numer. 27 (2018), 207-287.
  * (35) Lavor, C., Maculan, N.: _A function to test methods applied to global minimization of potential energy of molecules_ , Numer. Algorithms 35 (2004), 287-300.
  * (36) Leung, Y.-W., Wang, Y. P.: _An orthogonal genetic algorithm with quantization for global numerical optimization_ , IEEE Trans. Evol. Comput. 5 (2001), 41-53.
  * (37) Liu, S.-T., Luo, X.-L.: _A method based on Rayleigh quotient gradient flow for extreme and interior eigenvalue problems_ , Linear Algebra Appl. 432 (2010), 1851-1863.
  * (38) Luo, X.-L.: _Singly diagonally implicit Runge-Kutta methods combining line search techniques for unconstrained optimization_ , J. Comput. Math. 23, 153-164 (2005)
  * (39) Luo, X.-L., Kelley, C. T., Liao, L.-Z., Tam, H.-W.: _Combining trust region techniques and Rosenbrock methods to compute stationary points_ , J. Optim. Theory Appl. 140 (2009), 265-286.
  * (40) Luo, X.-L.: _A second-order pseudo-transient method for steady-state problems_ , Appl. Math. Comput., 216 (2010), 1752-1762.
  * (41) Luo, X.-L.: _A dynamical method of DAEs for the smallest eigenvalue problem_ , J. Comput. Sci. 3 (2012), 113-119.
  * (42) Luo, X.-L., Lv, J.-H., Sun, G.: _Continuation method with the trusty time-stepping scheme for linearly constrained optimization with noisy data_ , Optim. Eng., published online at published online, http://doi.org/10.1007/s11081-020-09590-z, 1-32, January 2021.
  * (43) Luo, X.-L., Xiao, H., Lv, J.-H.: _Continuation Newton methods with the residual trust-region time-stepping scheme for nonlinear equations_ , Numer. Algorithms, published online at http://doi.org/10.1007/s11075-021-01112-x, 1-25, April 2, 2021.
  * (44) Luo, X.-L., Yao, Y. Y.: _Primal-dual path-following methods and the trust-region strategy for linear programming with noisy data_ , J. Comput. Math., published online at http://doi.org/10.4208/jcm.2101-m2020-0173, or available at http://arxiv.org/abs/2006.07568, 1-21, March 15, 2021.
  * (45) Luo, X.-L., Xiao, H., Lv, J.-H., Zhang, S.: _Explicit pseudo-transient continuation and the trust-region updating strategy for unconstrained optimization_ , Appl. Numer. Math., 165 (2021), 290-302.
  * (46) Luo, X.-L., Xiao, H.: _Generalized continuation Newton methods and the trust-region updating strategy for the underdetermined system_ , J. Sci. comput. 88 (2021), published online at http://doi.org/10.1007/s10915-021-01566-0, 1-22.
  * (47) Luo, X.-L., Xiao, H.: _The regularization continuation method with an adaptive time step control for linearly equality-constrained optimization problems_ , submitted, arXiv preprint available at http://arxiv.org/abs/2106.01122, July 7, 2021.
  * (48) Man, K. F., Tang, K. S., Kwong, S.: _Genetic Algorithms: Concepts and Designs_ , Springer, Berlin, 1999.
  * (49) Macêdo, M. J. F. G., Karas, E. W., Costa, M. F. P., Rocha, A. M. A. C.: _Filter-based stochastic algorithm for global optimization_ , J. Glob. Optim. 77 (2020), 777-805.
  * (50) MATLAB v9.8.0 (R2020a), The MathWorks Inc., http://www.mathworks.com, 2020.
  * (51) MCS, _The multilevel coordinate search_ , available at https://www.mat.univie.ac.at/~neum/software/mcs/, 2000.
  * (52) Mitchell, M.: _An Introduction to Genetic Algorithms_ , MIT press, Cambridge, MA, 1996.
  * (53) Moré, J. J., Garbow, B. S., Hillstrom, K. E.: _Testing unconstrained optimization software_ , ACM Trans. Math. Soft. 7 (1981), 17-41.
  * (54) Moscato, P.: _On evolution, search, optimization, gas and martial arts: Toward memetic algorithms_ , Technical report C3P, Vol. 826, California Institute of Technology, Caltech Concurrent Computing Program, Pasadena, California, 1989.
  * (55) Neumaier, A.: _MCS: Global Optimization by Multilevel Coordinate Search_ , https://www.mat.univie.ac.at/~neum/software/mcs/, 2000.
  * (56) NEOS Server, https://neos-server.org/neos/, 2021.
  * (57) Nocedal, J., Wright, S. J.: _Numerical Optimization_ , Springer-Verlag, Berlin, 1999.
  * (58) Ortega, J. M., Rheinboldt, W. C.: _Iteration Solution of Nonlinear Equations in Several Variables_ , SIAM, Philadelphia, PA, 2000.
  * (59) Regis, R. G., Shoemaker, C. A.: _A quasi-multistart framework for global optimization of expensive functions using response surface models_ , J. Glob. Optim. 56 (2013), 1719-1753.
  * (60) Rios, L. M., Sahinidis, N. V.: _Derivative-free optimization: a review of algorithms and comparison of software implementations_ , J. Glob. Optim. 56 (2013), 1247-1293.
  * (61) Shampine, L. F., Gladwell, I., Thompson, S.: _Solving ODEs with MATLAB_ , Cambridge University Press, Cambridge, 2003.
  * (62) Sahinidis, N. V.: _BARON 21.1.13: Global optimization of mixed-integer nonlinear programs, user’s manual_ , available at https://minlp.com/downloads/docs/baronmanual.pdf, 2021.
  * (63) Surjanovic, S., Bingham, D.: _Virtual library of simulation experiments: test functions and datasets_ , available at http://www.sfu.ca/~ssurjano, January 2020.
  * (64) Sun, J., Garibaldi, J. M., Krasnogor, N., Zhang, Q.: _An intelligent muti-restart memetic algorithm for box constrained global optimisation_ , Evol. Comput. 21 (2014), 107-147.
  * (65) Sergeyev, Y. D., Kvasov, D. E., Mukhametzhanov, M. S.: _On the efficiency of nature-inspired metaheuristics in expensive global optimization with limited budget_ , Sci. Rep. 8 (2018), 1-9.
  * (66) Sun, W.-Y, Yuan, Y.-X.: _Optimization Theory and Methods: Nonlinear Programming_ , Springer, Berlin, 2006.
  * (67) Tanabe, K.: _Continuous Newton-Raphson method for solving an underdetermined system of nonlinear equations_ , Nonlinear Anal. 3 (1979), 495-503.
  * (68) Tawarmalani, M., Sahinidis, N.V.: _A polyhedral branch-and-cut approach to global optimization_ , Math. Program. 103 (2005), 225-249.
  * (69) Teughels, A., Roeck, G. De, Suykens, J. A. K.: _Global optimization by coupled local minimizers and its application to FE model updating_ , Comput. Sturct. 81 (2003), 2337-2351.
  * (70) Ugray, Z., Lasdon, L., Plummer, J., Glover, F., Kelly, J., Marti, R.: _Scatter search and local NLP solvers: A multistart framework for global optimization_ , INFORMS J. Comput. 19 (2007), 328-340.
  * (71) Yuan, Y.-X.: _Recent advances in trust region algorithms_ , Math. Program. 151 (2015), 249-281.


\tfrac{1}{4}\Delta(u_{1}-u_{2})=e^{u_{0}}+3e^{u_{1}}-3e^{u_{2}}+e^{u_{3}},\\\
\tfrac{1}{4}\Delta(u_{2}-u_{3})=-e^{u_{0}}-e^{u_{1}}+3e^{u_{2}}-2e^{u_{3}};\end{cases}$
(A.4) $\displaystyle\text{if $n\geq 4$:
}\begin{cases}\tfrac{1}{4}\Delta(u_{0}-u_{1})=2e^{u_{0}}-2e^{u_{1}},\\\
\tfrac{1}{4}\Delta(u_{1}-u_{2})=e^{u_{0}}+3e^{u_{1}}-3e^{u_{2}}+e^{u_{3}},\\\
\tfrac{1}{4}\Delta(u_{2}-u_{3})=-e^{u_{0}}-e^{u_{1}}+3e^{u_{2}}-3e^{u_{3}}+e^{u_{4}},\\\
\tfrac{1}{4}\Delta(u_{k-1}-u_{k})=-e^{u_{k-2}}+3e^{u_{k-1}}-3e^{u_{k}}+e^{u_{k+1}}\quad(4\leq
k\leq n-1),\\\
\tfrac{1}{4}\Delta(u_{n-1}-u_{n})=-e^{u_{n-2}}+3e^{u_{n-1}}-2e^{u_{n}}.\end{cases}$

We are now ready to carry out the proof of statements (1) and (2).

(1) The assumption $({\boldsymbol{\gamma}}_{1})\prec({\boldsymbol{\beta}})$
implies that
$e^{u_{0}-u_{1}}=\|{\boldsymbol{\beta}}\|^{2}/\|{\boldsymbol{\gamma}}_{1}\|^{2}$
attains its maximum $A_{1}>0$ on a set $P_{1}\subset\Sigma$ not containing the
zeros of ${\boldsymbol{\beta}}$ (note that the condition
$({\boldsymbol{\gamma}}_{1})\leq({\boldsymbol{\beta}})$ is sufficient for
$\|{\boldsymbol{\beta}}\|^{2}/\|{\boldsymbol{\gamma}}_{1}\|^{2}$ to attain
maximum, but in this case a maximal point might be a common zero of
${\boldsymbol{\gamma}}_{1}$ and ${\boldsymbol{\beta}}$ with the same
multiplicity). Since a maximal point of $e^{u_{0}-u_{1}}$ is also a maximal
point of $u_{0}-u_{1}$, by the first equation in (A.3) or (A.4) and Maximum
Principle, we obtain the following inequality, where every function is
evaluated at a given $p\in P_{1}$:

$\text{at $p\in P_{1}$:}\quad
0\geq\tfrac{1}{4}\Delta(u_{0}-u_{1})=2e^{u_{0}}-2e^{u_{1}}=2e^{u_{1}}(e^{u_{0}-u_{1}}-1)=2e^{u_{1}}(A_{1}-1).$

It follows that $A_{1}\leq 1$ because
$e^{u_{1}}=\|{\boldsymbol{\gamma}}_{1}\|^{2}>0$ at $p$.

We now improve this to the strict inequality $A_{1}<1$ by using E. Hopf’s
Strong Maximum Principle (see e.g. [31, Theorem 3.5]), which says that any
function $u$ on a disk in $\mathbb{C}$ satisfying
$\partial_{z\bar{z}}^{2}u\geq cu$ for a constant $c>0$ cannot achieve
nonnegative maximum unless it is constant. Suppose by contradiction that
$A_{1}=1$. Given $p\in P_{1}$, on any disk $U$ centered at $p$ and avoiding
the zeros of ${\boldsymbol{\beta}}$, we have

$\tfrac{1}{h_{0}}\partial_{z\bar{z}}^{2}(u_{0}-u_{1})=\tfrac{1}{4}\Delta(u_{0}-u_{1})=2e^{u_{1}}(e^{u_{0}-u_{1}}-1)\geq
e^{u_{1}}(u_{0}-u_{1}).$

Since the maximum $\log A_{1}=0$ of $u_{0}-u_{1}$ on $U$ is achieved at $p$,
we can use Strong Maximum Principle to infer that $u_{0}-u_{1}$ is identically
$0$ on $U$. Therefore, the nonempty closed set $P_{1}$ is also open, hence is
the whole $\Sigma$. This is impossible because $P_{1}$ does not intersect the
zero locus of ${\boldsymbol{\beta}}$, which is nonempty by assumption. Thus,
we have shown $A_{1}<1$. The required statement follows:

$\tfrac{1}{h_{0}}\partial_{z\bar{z}}^{2}\log
h_{1}=\|{\boldsymbol{\beta}}\|^{2}-\|{\boldsymbol{\gamma}}_{1}\|^{2}=\begin{cases}0&\text{
on $\\{{\boldsymbol{\gamma}}_{1}=0\\}$},\\\
\|{\boldsymbol{\gamma}}_{1}\|^{2}\big{(}\tfrac{\|{\boldsymbol{\beta}}\|^{2}}{\|{\boldsymbol{\gamma}}_{1}\|^{2}}-1\big{)}\leq\|{\boldsymbol{\gamma}}_{1}\|^{2}(A_{1}-1)<0&\text{
on $\\{{\boldsymbol{\gamma}}_{1}\neq 0\\}$}.\end{cases}$

(2) Assuming
$({\boldsymbol{\gamma}}_{n})\prec\cdots\prec({\boldsymbol{\gamma}}_{2})\prec\min\big{\\{}({\boldsymbol{\gamma}}_{1}),({\boldsymbol{\beta}})\big{\\}}$,
we first observe that the divisors $({\boldsymbol{\gamma}}_{1})$ and
$({\boldsymbol{\beta}})$ cannot be both empty: otherwise we have
$({\boldsymbol{\gamma}}_{2})=\cdots=({\boldsymbol{\gamma}}_{n})=\emptyset$ as
well, and hence the holomorphic line bundles
${\mathcal{K}}\\!{\mathcal{L}}_{j}^{-1}{\mathcal{L}}_{j+1}$,
${\mathcal{K}}\\!{\mathcal{L}}_{n}^{-1}$ and
${\mathcal{K}}\\!{\mathcal{L}}_{1}{\mathcal{L}}_{2}$ are all trivial; but the
triviality of the first two implies that
${\mathcal{L}}_{i}={\mathcal{K}}^{n+1-i}$, so
${\mathcal{K}}\\!{\mathcal{L}}_{1}{\mathcal{L}}_{2}={\mathcal{K}}^{2n}$ is not
trivial (as $\Sigma$ has genus $\geq 2$), a contradiction. Therefore, after
possibly replacing $({\mathcal{L}}_{1},{\boldsymbol{h}}_{1})$ by
$({\mathcal{L}}_{1}^{-1},{\boldsymbol{h}}_{1}^{-1})$ and interchanging
${\boldsymbol{\gamma}}_{1}$ with ${\boldsymbol{\beta}}$ (which does not affect
the equations), we may assume that $({\boldsymbol{\beta}})\neq\emptyset$.

We also infer from the assumption that
$e^{u_{0}-u_{2}}+e^{u_{1}-u_{2}}=\tfrac{\|{\boldsymbol{\gamma}}_{1}\|^{2}+\|{\boldsymbol{\beta}}\|^{2}}{\|{\boldsymbol{\gamma}}_{2}\|^{2}}$
attains its maximum $A_{2}>0$ on a set $P_{2}\subset\Sigma$ not containing the
zeros of ${\boldsymbol{\gamma}}_{2}$, while
$e^{u_{k-1}-u_{k}}=\|{\boldsymbol{\gamma}}_{k-1}\|^{2}/\|{\boldsymbol{\gamma}}_{k}\|^{2}$
($3\leq k\leq n$) attains its maximum $A_{k}>0$ on $P_{k}\subset\Sigma$ not
containing the zeros of ${\boldsymbol{\gamma}}_{k}$. We proceed in the
following steps:

Step 1. Use Maximum Principle to produce a set of algebraic inequalities on
$A_{2},\cdots,A_{n}$.

First consider the case $n\geq 4$. Using (A.4) and Maximum Principle, we get

$\displaystyle\text{at $p\in P_{3}$:}\quad 0$
$\displaystyle\geq\tfrac{1}{4}e^{-u_{3}}\Delta(u_{2}-u_{3})$
$\displaystyle=e^{-u_{3}}(-e^{u_{0}}-e^{u_{1}}+3e^{u_{2}}-3e^{u_{3}}+e^{u_{4}})$
$\displaystyle=-e^{u_{2}-u_{3}}(e^{u_{0}-u_{2}}+e^{u_{1}-u_{2}}-1)+2(e^{u_{2}-u_{3}}-1)-(1-e^{u_{4}-u_{3}})$
$\displaystyle\geq-
A_{3}(A_{2}-1)+2(A_{3}-1)-\big{(}1-\tfrac{1}{A_{4}}\big{)},$

where the multiplier $e^{-u_{3}}$ in the first line is allowed because $P_{3}$
does not contain the zeros of ${\boldsymbol{\gamma}}_{3}$. Similarly, we have

$\displaystyle\text{at $p\in P_{k}$ ($4\leq k\leq n-1$):}\quad 0$
$\displaystyle\geq\tfrac{1}{4}e^{-u_{k}}\Delta(u_{k-1}-u_{k})$
$\displaystyle=e^{-u_{k}}(-e^{u_{k-2}}+3e^{u_{k-1}}-3e^{u_{k}}+e^{u_{k+1}})$
$\displaystyle=-e^{u_{k-1}-u_{k}}(e^{u_{k-2}-u_{k-1}}-1)+2(e^{u_{k-1}-u_{k}}-1)-(1-e^{u_{k+1}-u_{k}})$
$\displaystyle\geq-
A_{k}(A_{k-1}-1)+2(A_{k}-1)-\big{(}1-\tfrac{1}{A_{k+1}}\big{)};$
$\displaystyle\text{at $p\in P_{n}$:}\quad 0$
$\displaystyle\geq\tfrac{1}{4}e^{-u_{n}}\Delta(u_{n-1}-u_{n})$
$\displaystyle=e^{-u_{n}}(e^{-u_{n-2}}+3e^{u_{n-1}}-2e^{u_{n}})$
$\displaystyle=-e^{u_{n-1}-u_{n}}(e^{u_{n-2}-u_{n-1}}-1)+2(e^{u_{n-1}-u_{n}}-1)$
$\displaystyle\geq-A_{n}(A_{n-1}-1)+2(A_{n}-1).$

As for $P_{2}$, we use the inequality
$\Delta\log(e^{f}+1)=\tfrac{e^{f}}{e^{f}+1}\Delta
f+\tfrac{e^{f}}{(e^{f}+1)^{2}}\|\nabla f\|^{2}\geq\tfrac{e^{f}}{e^{f}+1}\Delta
f$, which holds for any $f\in C^{\infty}(\Sigma)$, to get

(A.5) $\displaystyle\text{at $p\in P_{2}$:}\quad 0$
$\displaystyle\geq\tfrac{1}{4}\Delta\log\big{(}e^{u_{0}-u_{2}}+e^{u_{1}-u_{2}}\big{)}$
$\displaystyle=\tfrac{1}{4}\Delta\log\big{(}e^{u_{0}-u_{1}}+1\big{)}+\tfrac{1}{4}\Delta(u_{1}-u_{2})$
$\displaystyle\geq\frac{e^{u_{0}-u_{1}}}{e^{u_{0}-u_{1}}+1}\cdot\tfrac{1}{4}\Delta(u_{0}-u_{1})+\tfrac{1}{4}\Delta(u_{1}-u_{2})$
$\displaystyle=\frac{e^{u_{0}-u_{1}}}{e^{u_{0}-u_{1}}+1}\cdot
2(e^{u_{0}}-e^{u_{1}})+(e^{u_{0}}+3e^{u_{1}}-3e^{u_{2}}+e^{u_{3}})$
$\displaystyle=\frac{e^{-u_{1}}(e^{u_{0}}-e^{u_{1}})^{2}}{e^{u_{0}-u_{1}}+1}+(e^{u_{0}}-e^{u_{1}})+(e^{u_{0}}+3e^{u_{1}}-3e^{u_{2}}+e^{u_{3}})$
$\displaystyle\geq(e^{u_{0}}-e^{u_{1}})+(e^{u_{0}}+3e^{u_{1}}-3e^{u_{2}}+e^{u_{3}})$
$\displaystyle=e^{u_{2}}\left[2\big{(}e^{u_{0}-u_{2}}+e^{u_{1}-u_{2}}-1\big{)}-\big{(}1-e^{u_{3}-u_{2}}\big{)}\right]$
$\displaystyle\geq
e^{u_{2}}\left[2(A_{2}-1)-\big{(}1-\tfrac{1}{A_{3}}\big{)}\right].$

Thus, we have obtained

(A.6) $\text{ if $n\geq 4$:
}\begin{cases}2(A_{2}-1)-\big{(}1-\tfrac{1}{A_{3}}\big{)}\leq 0,\\\
-A_{3}(A_{2}-1)+2(A_{3}-1)-\big{(}1-\tfrac{1}{A_{4}}\big{)}\leq 0,\\\
-A_{k}(A_{k-1}-1)+2(A_{k}-1)-\big{(}1-\tfrac{1}{A_{k+1}}\big{)}\leq
0\quad(4\leq k\leq n-1),\\\ -A_{n}(A_{n-1}-1)+2(A_{n}-1)\leq 0.\end{cases}$

When $n=2,3$, the same argument yields

(A.7) $\displaystyle\text{ if $n=2$: }A_{2}-1\leq 0;$ (A.8)
$\displaystyle\text{ if $n=3$:
}\begin{cases}2(A_{2}-1)-\big{(}1-\tfrac{1}{A_{3}}\big{)}\leq 0,\\\
-A_{3}(A_{2}-1)+2(A_{3}-1)\leq 0.\end{cases}$

Step 2. Infer from these inequalities that $A_{2},\cdots,A_{n}$ are either all
equal to $1$ or all less than $1$.

If $n=2$, there is nothing to prove. If $n=3$, we can rewrite (A.8) as
$4\big{(}1-\tfrac{1}{A_{3}}\big{)}\leq 2(A_{2}-1)\leq 1-\tfrac{1}{A_{3}}$,
whence the required conclusion follows easily. Now assume that $n\geq 4$ and
$A_{2},\cdots,A_{n}$ are not all $1$. Put
$B_{k}:=\big{(}1-\tfrac{1}{A_{k}}\big{)}-(A_{k-1}-1)$ for $3\leq k\leq n$.
Then (A.6) can be rewritten as

(A.9) $\begin{cases}A_{2}-1\leq B_{3},\\\ A_{3}B_{3}\leq B_{4},\\\
A_{k}B_{k}\leq B_{k+1}\quad(4\leq k\leq n-1),\\\
A_{n}B_{n}\leq-(A_{n}-1).\end{cases}$

Let us first show $A_{n}\leq 1$. Applying the definition of the $B_{k}$’s
repeatedly, we get

$\displaystyle 1-\tfrac{1}{A_{n}}$
$\displaystyle=B_{n}+A_{n-1}-1=B_{n}+A_{n-1}\big{(}1-\tfrac{1}{A_{n-1}}\big{)}$
$\displaystyle=B_{n}+A_{n-1}(B_{n-1}+A_{n-2}-1)=B_{n}+A_{n-1}B_{n-1}+A_{n-1}A_{n-2}\big{(}1-\tfrac{1}{A_{n-2}}\big{)}$
$\displaystyle=B_{n}+A_{n-1}B_{n-1}+A_{n-1}A_{n-2}(B_{n-2}+A_{n-3}-1)=\cdots$
$\displaystyle=B_{n}+A_{n-1}B_{n-1}+A_{n-1}A_{n-2}B_{n-2}+\ \cdots\
+A_{n-1}\cdots A_{3}B_{3}+A_{n-1}\cdots A_{3}(A_{2}-1)$

The result has $n-1$ terms and each term is no larger than $B_{n}$ by (A.9).
Therefore, we have

$1-\tfrac{1}{A_{n}}\leq(n-1)B_{n}\leq-(n-1)\big{(}1-\tfrac{1}{A_{n}}\big{)},$

which implies $A_{n}\leq 1$, as required.

Next, we show $A_{2}<1$. Suppose by contradiction that $A_{2}\geq 1$. Then the
inequalities in (A.9) imply, successively for $3\leq k\leq n$, that $B_{k}\geq
0$; or equivalently, $1-\tfrac{1}{A_{k}}\geq A_{k-1}-1$. This in turn implies
successively that $A_{k}\geq 1$ for $3\leq k\leq n$. Since we already know
that $A_{n}\leq 1$, it follows that $A_{n}=1$. This allows us to deduce again
from (A.9) that $B_{k}\leq 0$ for $3\leq k\leq n$ (in reversed order).
Therefore, we have $B_{3}=\cdots=B_{n}=0$, which implies
$A_{2}=\cdots=A_{n}=1$, a contradiction.

Finally, we show that $A_{3},\cdots,A_{n}<0$ in a similar way as follows.
Suppose by contradiction that $A_{l-1}<1$ but $A_{l}\geq 1$ for some $3\leq
l\leq n$. Then we have
$B_{l}:=\big{(}1-\tfrac{1}{A_{l}}\big{)}-(A_{l-1}-1)>0$. Similarly as in the
above proof of $A_{2}<1$, we infer that $B_{k}>0$ and hence $A_{k}\geq 1$ for
$l\leq k\leq n$. So again we have $A_{n}=1$, which implies $B_{n}\leq 0$, a
contradiction.

Step 3. Exclude the case $A_{2}=\cdots=A_{n}=1$ by using Stronger Maximum
Principle.

Suppose $A_{2}=\cdots=A_{n}=1$ by contradiction. Observe that all the
inequalities in (A.5) other than the first and the last hold outside of the
zeros of $({\boldsymbol{\beta}})$ or $({\boldsymbol{\gamma}}_{1})$ (namely,
they do not require the functions to be evaluated at $p\in P_{2}$). So we have
the following estimates for
$w:=\log\big{(}e^{u_{0}-u_{2}}+e^{u_{1}-u_{2}}\big{)}=\log\tfrac{\|{\boldsymbol{\gamma}}_{1}\|^{2}+\|{\boldsymbol{\beta}}\|^{2}}{\|{\boldsymbol{\gamma}}_{2}\|^{2}}$
away from the zeros of $({\boldsymbol{\beta}})$ or
$({\boldsymbol{\gamma}}_{1})$:

$\displaystyle\tfrac{1}{4}\Delta w$ $\displaystyle\geq
e^{u_{2}}\left[2\big{(}e^{u_{0}-u_{2}}+e^{u_{1}-u_{2}}-1\big{)}+\big{(}e^{-(u_{2}-u_{3})}-1\big{)}\right]$
$\displaystyle\geq 2e^{u_{2}}\big{(}e^{u_{0}-u_{2}}+e^{u_{1}-u_{2}}-1\big{)}$
$\displaystyle=2e^{u_{2}}(e^{w}-1)\geq 2e^{u_{2}}w,$

where the second inequality is because $e^{u_{2}-u_{3}}\leq A_{3}=1$. Since
$w$ achieves its maximum $\log A_{2}=0$ on $P_{2}$, we may use Strong Maximum
Principle in the same way as in Part (1) to conclude that $P_{2}=\Sigma$, or
in other words,
$\tfrac{\|{\boldsymbol{\gamma}}_{1}\|^{2}+\|{\boldsymbol{\beta}}\|^{2}}{\|{\boldsymbol{\gamma}}_{2}\|^{2}}=1$
throughout $\Sigma$. By the assumption
$({\boldsymbol{\gamma}}_{2})\prec\min\big{\\{}({\boldsymbol{\gamma}}_{1}),({\boldsymbol{\beta}})\\}$,
it follows that ${\boldsymbol{\gamma}}_{1}$ and ${\boldsymbol{\beta}}$ do not
have any common zero and ${\boldsymbol{\gamma}}_{2}$ is nowhere zero.

In order to derive a contradiction, we consider the maximum $A_{2}^{\prime}$
of
$e^{u_{1}-u_{2}}=\|{\boldsymbol{\gamma}}_{1}\|^{2}/\|{\boldsymbol{\gamma}}_{2}\|^{2}$.
By the second equation in (A.3) or (A.4), at any point where $A_{2}^{\prime}$
is achieved (which cannot be a zero of
$e^{u_{2}}=\|{\boldsymbol{\gamma}}_{2}\|^{2}$ because
$({\boldsymbol{\gamma}}_{2})\prec({\boldsymbol{\gamma}}_{1})$), we have

$\displaystyle 0$
$\displaystyle\geq\tfrac{1}{4}e^{-u_{2}}\Delta(u_{1}-u_{2})=e^{-u_{2}}\big{(}e^{u_{0}}+3e^{u_{1}}-3e^{u_{2}}+e^{u_{3}}\big{)}$
$\displaystyle\geq
e^{u_{0}-u_{2}}+e^{u_{1}-u_{2}}+2e^{u_{1}-u_{2}}-3+e^{-(u_{2}-u_{3})}$
$\displaystyle=2e^{u_{1}-u_{2}}-2+e^{-(u_{2}-u_{3})}$ $\displaystyle\geq
2A_{2}^{\prime}-2+\tfrac{1}{A_{3}}=2A_{2}^{\prime}-1.$

It follows that $A_{2}^{\prime}\leq\tfrac{1}{2}$. But this is impossible
because
$\|{\boldsymbol{\gamma}}_{1}\|^{2}/\|{\boldsymbol{\gamma}}_{2}\|^{2}=1$ on the
zero locus of ${\boldsymbol{\beta}}$ (which is nonempty as observed at the
beginning). This contradiction shows that $A_{2}=\cdots=A_{n}=1$ cannot hold.
Therefore, we have $A_{2},\cdots,A_{n}<1$.

Step 4. Deduce the required statement from the inequality
$A_{2},\cdots,A_{n}<1$ that we just obtained: since
$({\boldsymbol{\gamma}}_{n})\leq\cdots\leq({\boldsymbol{\gamma}}_{2})\leq\min\big{\\{}({\boldsymbol{\gamma}}_{1}),({\boldsymbol{\beta}})\big{\\}}$,
we have

$\tfrac{1}{h_{0}}\partial_{z\bar{z}}^{2}\log
h_{2}=\|{\boldsymbol{\gamma}}_{1}\|^{2}-\|{\boldsymbol{\gamma}}_{2}\|^{2}+\|{\boldsymbol{\beta}}\|^{2}\\\
=\begin{cases}0&\text{on $\\{{\boldsymbol{\gamma}}_{2}=0\\}$},\\\
\|{\boldsymbol{\gamma}}_{2}\|^{2}\big{(}\tfrac{\|{\boldsymbol{\gamma}}_{1}\|^{2}+\|{\boldsymbol{\beta}}\|^{2}}{\|{\boldsymbol{\gamma}}_{2}\|^{2}}-1\big{)}\leq\|{\boldsymbol{\gamma}}_{2}\|^{2}(A_{2}-1)<0&\text{on
$\\{{\boldsymbol{\gamma}}_{2}\neq 0\\}$};\end{cases}$
$\tfrac{1}{h_{0}}\partial_{z\bar{z}}^{2}\log
h_{k}=\|{\boldsymbol{\gamma}}_{k-1}\|^{2}-\|{\boldsymbol{\gamma}}_{k}\|^{2}=\begin{cases}0&\text{on
$\\{{\boldsymbol{\gamma}}_{k}=0\\}$},\\\
\|{\boldsymbol{\gamma}}_{k}\|^{2}\big{(}\tfrac{\|{\boldsymbol{\gamma}}_{k-1}\|^{2}}{\|{\boldsymbol{\gamma}}_{k}\|^{2}}-1\big{)}\leq\|{\boldsymbol{\gamma}}_{k}\|^{2}(A_{k}-1)<0&\text{on
$\\{{\boldsymbol{\gamma}}_{k}\neq 0\\}$}.\end{cases}$

Thus, we have proven statement (2) in the case ${\boldsymbol{\beta}}\not\equiv
0$. To adapt this proof to ${\boldsymbol{\beta}}\equiv 0$, we basically only
need to erase all the terms containing $u_{0}$ or
$\|{\boldsymbol{\beta}}\|^{2}$ and take into account some simplications along
the way. For example, when $n\geq 4$, the equations that we started with
become

$\begin{cases}\tfrac{1}{4}\Delta(u_{1}-u_{2})=3e^{u_{1}}-3e^{u_{2}}+e^{u_{3}}\geq
3e^{u_{1}}-3^{u_{2}},\\\
\tfrac{1}{4}\Delta(u_{2}-u_{3})=-e^{u_{1}}+3e^{u_{2}}-3e^{u_{3}}+e^{u_{4}},\\\
\tfrac{1}{4}\Delta(u_{k-1}-u_{k})=-e^{u_{k-2}}+3e^{u_{k-1}}-3e^{u_{k}}+e^{u_{k+1}}\quad(4\leq
k\leq n-1),\\\
\tfrac{1}{4}\Delta(u_{n-1}-u_{n})=-e^{u_{n-2}}+3e^{u_{n-1}}-2e^{u_{n}}.\end{cases}$

Putting $A_{k}:=\max_{\Sigma}e^{u_{k-1}-u_{k}}$ for $k=2,\cdots,n$ (only
$A_{2}$ differs from the above one), we may use Maximum Principle to deduce
almost the same inequalities on $A_{2},\cdots,A_{n}$ as (A.6), except that the
first one is replaced by the simpler inequality $A_{2}-1\leq 0$. Then, in the
same way as above, we may show that $A_{2},\cdots,A_{n}<1$, whence the
required statement follows. ∎

## References

  * [1] D. Alessandrini and B. Collier, The geometry of maximal components of the $\mathrm{PSp}(4,\mathbb{R})$ character variety, Geom. Topol., 23 (2019), pp. 1251–1337.
  * [2] H. Anciaux, Minimal submanifolds in pseudo-Riemannian geometry, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2011.  With a foreword by F. Urbano.
  * [3] M. A. Arroyo, The geometry of $\mathrm{SO}(p,q)$-Higgs bundles, Ph.D thesis, Facultad de Ciencias de la Universidad de Salamanca, (2009).
  * [4] J. C. Baez and J. Huerta, $\mathrm{G}_{2}$ and the rolling ball, Trans. Amer. Math. Soc., 366 (2014), pp. 5257–5293.
  * [5] D. Baraglia, $\mathrm{G}_{2}$ geometry and integrable systems, Ph.D thesis at University of Oxford, (2009).
  * [6] D. Baraglia, Cyclic Higgs bundles and the affine Toda equations, Geom. Dedicata, 174 (2015), pp. 25–42.
  * [7] J. a. L. M. Barbosa, On minimal immersions of $\mathbb{S}^{2}$ into $\mathbb{S}^{2m}$, Trans. Amer. Math. Soc., 210 (1975), pp. 75–106.
  * [8] N. Berline, E. Getzler, and M. Vergne, Heat kernels and Dirac operators, Grundlehren Text Editions, Springer-Verlag, Berlin, 2004.  Corrected reprint of the 1992 original.
  * [9] J. Beyrer and B. Pozzetti, Positive surface group representations in $\mathrm{PO}(p,q)$, arXiv e-prints, (2021), p. arXiv:2106.14725.
  * [10] F. Bonsante and J.-M. Schlenker, Maximal surfaces and the universal Teichmüller space, Invent. Math., 182 (2010), pp. 279–333.
  * [11] F. Bonsante and A. Seppi, Anti-de sitter geometry and Teichmüller theory, in In the tradition of Thurston, Springer, Cham, [2020] ©2020, pp. 545–643.
  * [12] R. L. Bryant, Conformal and minimal immersions of compact surfaces into the $4$-sphere, J. Differential Geometry, 17 (1982), pp. 455–473.
  * [13] , Metrics with exceptional holonomy, Ann. of Math. (2), 126 (1987), pp. 525–576.
  * [14] M. Burger, A. Iozzi, and A. Wienhard, Surface group representations with maximal Toledo invariant, Ann. of Math. (2), 172 (2010), pp. 517–566.
  * [15] E. Calabi, Minimal immersions of surfaces in Euclidean spheres, J. Differential Geometry, 1 (1967), pp. 111–125.
  * [16] S. S. Chern, On the minimal immersions of the two-sphere in a space of constant curvature, in Problems in analysis (Lectures at the Sympos. in honor of Salomon Bochner, Princeton Univ., Princeton, N.J., 1969), 1970, pp. 27–40.
  * [17] T. H. Colding and W. P. Minicozzi, II, A course in minimal surfaces, vol. 121 of Graduate Studies in Mathematics, American Mathematical Society, Providence, RI, 2011.
  * [18] B. Collier, Maximal ${\rm Sp}(4,\mathbb{R})$ surface group representations, minimal immersions and cyclic surfaces, Geom. Dedicata, 180 (2016), pp. 241–285.
  * [19] , $\mathrm{SO}(n,n+1)$-surface group representations and Higgs bundles, Ann. Sci. Éc. Norm. Supér. (4), 53 (2020), pp. 1561–1616.
  * [20] B. Collier and Q. Li, Asymptotics of Higgs bundles in the Hitchin component, Adv. Math., 307 (2017), pp. 488–558.
  * [21] B. Collier, N. Tholozan, and J. Toulisse, The geometry of maximal representations of surface groups into ${\rm SO}_{0}(2,n)$, Duke Math. J., 168 (2019), pp. 2873–2949.
  * [22] K. Corlette, Flat $G$-bundles with canonical metrics, J. Differential Geom., 28 (1988), pp. 361–382.
  * [23] S. Dai and Q. Li, Minimal surfaces for Hitchin representations, J. Differential Geom., 112 (2019), pp. 47–77.
  * [24] , On cyclic Higgs bundles, Math. Ann., 376 (2020), pp. 1225–1260.
  * [25] S. K. Donaldson, Twisted harmonic maps and the self-duality equations, Proc. London Math. Soc. (3), 55 (1987), pp. 127–131.
  * [26] J. Eells, Jr. and J. H. Sampson, Harmonic mappings of Riemannian manifolds, Amer. J. Math., 86 (1964), pp. 109–160.
  * [27] F. Forstneric, Proper superminimal surfaces of given conformal types in the hyperbolic four-space, arXiv e-prints, (2020), p. arXiv:2005.02201.
  * [28] , The Calabi-Yau property of superminimal surfaces in self-dual Einstein four-manifolds, arXiv e-prints, (2020), p. arXiv:2004.03536.
  * [29] T. Friedrich, On superminimal surfaces, Arch. Math. (Brno), 33 (1997), pp. 41–56.
  * [30] O. Garcia-Prada, P. B. Gothen, and I. M. i. Riera, The Hitchin-Kobayashi correspondence, Higgs pairs and surface group representations, arXiv e-prints, (2009), p. arXiv:0909.4487.
  * [31] D. Gilbarg and N. S. Trudinger, Elliptic partial differential equations of second order, Classics in Mathematics, Springer-Verlag, Berlin, 2001\.  Reprint of the 1998 edition.
  * [32] P. B. Gothen, Components of spaces of representations and stable triples, Topology, 40 (2001), pp. 823–850.
  * [33] A. Gray, Minimal varieties and almost Hermitian submanifolds, Michigan Math. J., 12 (1965), pp. 273–287.
  * [34] , Nearly Kähler manifolds, J. Differential Geometry, 4 (1970), pp. 283–309.
  * [35] O. Guichard, An introduction to the differential geometry of flat bundles and of Higgs bundles, in The geometry, topology and physics of moduli spaces of Higgs bundles, vol. 36 of Lect. Notes Ser. Inst. Math. Sci. Natl. Univ. Singap., World Sci. Publ., Hackensack, NJ, 2018, pp. 1–63.
  * [36] O. Guichard, F. Labourie, and A. Wienhard, Positivity and representations of surface groups, arXiv e-prints, (2021), p. arXiv:2106.14584.
  * [37] O. Guichard and A. Wienhard, Positivity and higher Teichmüller theory, in European Congress of Mathematics, Eur. Math. Soc., Zürich, 2018, pp. 289–310.
  * [38] R. D. Gulliver, II, R. Osserman, and H. L. Royden, A theory of branched immersions of surfaces, Amer. J. Math., 95 (1973), pp. 750–812.
  * [39] S. u. Helgason, Differential geometry and symmetric spaces, Pure and Applied Mathematics, Vol. XII, Academic Press, New York-London, 1962.
  * [40] N. J. Hitchin, The self-duality equations on a Riemann surface, Proc. London Math. Soc. (3), 55 (1987), pp. 59–126.
  * [41] , Lie groups and Teichmüller space, Topology, 31 (1992), pp. 449–473.
  * [42] Z. Huang, M. Lucia, and G. Tarantello, Bifurcation for minimal surface equation in hyperbolic 3-manifolds, Ann. Inst. H. Poincaré C Anal. Non Linéaire, 38 (2021), pp. 243–279.
  * [43] S. Karigiannis, Introduction to $G_{2}$ Geometry, in Lectures and surveys on $G_{2}$-manifolds and related topics, vol. 84 of Fields Inst. Commun., Springer, New York, 2020, pp. 3–50.
  * [44] K. Kommerell, Die Krümmung der zweidimensionalen Gebilde im ebenen Raum von vier Dimensionen, Dissertation Tübingen, (1897).
  * [45] F. Labourie, Flat projective structures on surfaces and cubic holomorphic differentials, Pure Appl. Math. Q., 3 (2007), pp. 1057–1099.
  * [46] F. Labourie, Cross ratios, Anosov representations and the energy functional on Teichmüller space, Ann. Sci. Éc. Norm. Supér. (4), 41 (2008), pp. 437–469.
  * [47] , Cross ratios, Anosov representations and the energy functional on Teichmüller space, Ann. Sci. Éc. Norm. Supér. (4), 41 (2008), pp. 437–469.
  * [48] , Cyclic surfaces and Hitchin components in rank 2, Ann. of Math. (2), 185 (2017), pp. 1–58.
  * [49] F. Labourie and J. Toulisse, Quasicircles and quasiperiodic surfaces in pseudo-hyperbolic spaces, arXiv e-prints, (2020), p. arXiv:2010.05704.
  * [50] F. Labourie, J. Toulisse, and M. Wolf, Plateau Problems for Maximal Surfaces in Pseudo-Hyperbolic Spaces, arXiv e-prints, (2020), p. arXiv:2006.12190.
  * [51] J. C. Loftin, Affine spheres and convex $\mathbb{RP}^{n}$-manifolds, Amer. J. Math., 123 (2001), pp. 255–274.
  * [52] V. Marković, Non-uniqueness of minimal surfaces in a product of closed Riemann surfaces, Geom. Funct. Anal., 32 (2022), pp. 31–52.
  * [53] V. Markovic, N. Sagman, and P. Smillie, Unstable minimal surfaces in $\mathbb{R}^{n}$ and in products of hyperbolic surfaces, arXiv e-prints, (2022), p. arXiv:2206.02938.
  * [54] M.-L. Michelsohn, Surfaces minimales dans les sphères, no. 154-155, 1987, pp. 9, 115–130, 350 (1988).  Théorie des variétés minimales et applications (Palaiseau, 1983–1984).
  * [55] K. Sagerschnig, Split octonions and generic rank two distributions in dimension five, Arch. Math. (Brno), 42 (2006), pp. 329–339.
  * [56] C. T. Simpson, Constructing variations of Hodge structure using Yang-Mills theory and applications to uniformization, J. Amer. Math. Soc., 1 (1988), pp. 867–918.
  * [57] M. Spivak, A comprehensive introduction to differential geometry. Vol. IV, Publish or Perish, Inc., Boston, Mass., 1975.
  * [58] C. H. Taubes, Minimal surfaces in germs of hyperbolic 3-manifolds, in Proceedings of the Casson Fest, vol. 7 of Geom. Topol. Monogr., Geom. Topol. Publ., Coventry, 2004, pp. 69–100.
  * [59] K. K. Uhlenbeck, Closed minimal surfaces in hyperbolic $3$-manifolds, in Seminar on minimal submanifolds, vol. 103 of Ann. of Math. Stud., Princeton Univ. Press, Princeton, NJ, 1983, pp. 147–168.
  * [60] A. Wienhard, An invitation to higher Teichmüller theory, in Proceedings of the International Congress of Mathematicians—Rio de Janeiro 2018. Vol. II. Invited lectures, World Sci. Publ., Hackensack, NJ, 2018, pp. 1013–1039.
  * [61] Y. Xin, Minimal submanifolds and related topics, vol. 16 of Nankai Tracts in Mathematics, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2019.  Second edition of [ MR2035469].