BACKTRACKING NEW Q-NEWTON’S METHOD, NEWTON’S FLOW, VORONOI’S DIAGRAM AND STOCHASTIC ROOT FINDING

. A new variant of Newton’s method - named Back-tracking New Q-Newton’s method (BNQN) - which has strong theoretical guarantee, is easy to implement, and has good experimental performance, was recently introduced by the third author. Experiments performed previously showed some remarkable properties of the basins of attractions for finding roots of polynomials and meromorphic functions, with BNQN. In general, they look more smooth than that of Newton’s method. In this paper, we continue to experimentally explore in depth this remarkable phenomenon, and connect BNQN to Newton’s flow and Voronoi’s diagram. This link poses a couple of challenging puzzles to be explained. Experiments also indicate that BNQN is more robust against random perturbations than Newton’s method and Random Relaxed Newton’s method.


Introduction
This paper mainly concerns the problem of finding roots of a polynomial or meromorphic function in 1 complex variable f (z) = 0.
Newton's method is a well known method to solve equations and optimization problems.In the case of a function f (z) concerns here, the method amounts to the update rule: One can also apply the higher dimensional version of Newton's method for systems of equations, to the system G(x, y) = (Re(f (x+iy)), Im(f (x+ iy))) = 0 of real and imaginary parts of f , where x, y ∈ R, i.e. to consider the iterative method z n+1 = z n − JG(z n ) −1 .G(z n ), Date: January 9, 2024.where JG is the Jacobian matrix of G.However, in the case where f is a meromorphic function, this iterative method is the same as that of the usual Newton's method, thanks to Cauchy-Riemann's equations.
A new variant of Newton's method -named Backtracking New Q-Newton's method (BNQN) -which has strong theoretical guarantee, is easy to implement, and has good experimental performance, was recently introduced by the third author [19].It is a variant of Newton's method for optimization.For details about the method, also comparisons between it and some well known variants of Newton's method (including Newton's method itself), the readers can see [19] [20].In Section 2 we will recall essential information on BNQN, enough for the purpose of this paper.
Experiments performed previously in [19], indicate a remarkable phenomenon: The basins of attraction, when BNQN is applied to find roots of meromorphic functions, seem to be more smooth than that of Newton's method.In particular, these basins seem not to be fractal, as is often the case for Newton's method [16].Moreover, when finding roots of a polynomial of degree 2 in 1 complex variable, the picture one obtains for basins of attraction of BNQN looks exactly the same as that in the classical Schröder's theorem for Newton's method [17] [1].
Which leads to a natural and interesting question: Can we explain this phenomenon rigorously, or is it just a lucky occurrence of numerical calculations?
After presentations of BNQN at ICIAM 2023 (Waseda University) and a satellite conference, a couple of candidates, which have a lot of connections in mathematics and physics, for explaining the pictures obtained in [19] have emerged: -These seem relevant to Newton's flow.
-These seem relevant to Voronoi's diagrams.More details about Newton's flow and Voronoi's diagrams will be recalled in Section 2. Note that Voronoi's diagram for a pair of distinct points in the plane is the same as what produced in Schröder's theorem for Newton's method.Also, Newton's flow for a polynomial of degree 2 produces exactly the same picture.
In a previous paper [8], we proved that indeed BNQN applied to finding roots of a polynomial f (z) of degree 2 will produce the same picture as in Schröder's theorem for Newton's method.Even so, BNQN is indeed more smooth than Newton's method in this case: On the boundary line, the dynamics of Newton's method is chaotic, while the dynamics of BNQN has a global attraction.Newton's method for optimization, on the other hand, produces a picture with many black blobs of initial points z 0 for which the constructed sequence {z n } does not converge to the roots of f (z) but instead to the critical point of f , and the basins of attraction even seem to be fractal!In this paper, we will explore in depth the possible links mentioned above between BNQN, Newton's flow and Voronoi's diagrams, when finding roots of meromorphic functions.We verify with many different geometric configurations of the roots, that indeed the link is very apparent and well above a random coincidence.We work with three different versions of Newton's flow, with two versions of Newton's method, as well as Random Relaxed Newton's method.For the readers's convenience, we will recall in enough detail on Newton's method and Random Relaxed Newton's method in the next Section.It seems that the picture one obtains from BNQN is a better version of what one obtains with Newton's method for Optimization, and while BNQN is not a flow method it somehow attains many good features of flow methods.From these experiments, another interesting new phenomenon seems to occur when one applies Newton's flow to f /f ′ , even if one only wants to find roots of f and even if f has only roots of multiplicity 1.
An interesting possible application of the study in this paper is as follows.In the literature, there is no unique construction of Voronoi's diagrams when two or more points meet together (limits of Voronoi's diagrams) [10].In such a case, we expect that a curved version of a "canonical" Voronoi's diagram can be obtained by using Newton's flow or BNQN.For example, working with the polynomial P (z) = (z − a) 2  Given that when using computers one cannot avoid dealing with random errors/perturbations, and that solving equations with random perturbations is interesting itself, we also present some experiments concerning this issue.The experiments indicate that BNQN is more robust compared to Newton's method and Random Relaxed Newton's method.
Plan of the paper: In Section 2 we recall relevant information concerning the methods, as well as Newton's flow and Voronoi's diagrams, as well as stochastic root finding.In Section 3, we introduce the setting of our experiments and then present the experimental results.In the final section, we draw some conclusions and some challenging puzzles and directions for future's research.
Acknowledgments: Mi Hu and Tuyen Trung Truong are partially supported by Young Research Talents grant 300814 from Research Council of Norway.Takayuki Watanabe is partially supported by JSPS Grant-in-Aid for Early-Career Scientists Grant Number JP 23K13000.Some ideas of the work were discussed/carried out in the inspiring environments of the conferences ICIAM 2023 (Waseda University) and "Birational geometry and algebraic dynamics" 2023 (University of Tokyo), for which we would like to thank the organisers and participants, in particularly Mark Comerford for suggesting the connection to Voronoi's diagrams.

Preliminaries
In this section, we briefly review some previous results and algorithms, with the emphasis on properties which will be used later.
2.1.Newton's method.Besides the main version in the introduction, here we recall Newton's method in the optimization setting (Newton's method for optimization) also, to easily compare later to BNQN.Assume that one wants to find (local) minima of an objective function F : R m → R. Let ∇F be the gradient of F , and ∇ 2 F be the Hessian of F .Then Newton's method for the function F is the following iterative algorithm: Choose z 0 ∈ R m an initial point, and define: In this paper, Newton's method for optimization will be applied to the function F (x, y) = |f (x + iy)| 2 /2, for a meromorphic function f in 1 complex variable z = x + iy.
Both versions of Newton's method are invariant under a linear change of coordinates.However, their behaviour can be very different, see the experimental results for detail.
Concerning finding roots of polynomials, we note that Newton's method gives rise to an algebraic dynamical system on the Riemann sphere.As such, [11] shows that Newton's method does not guarantee to always find roots of polynomials of degree at least 4, see the cited paper for the precise statement.For Newton's method for optimization, it has the tendency to converge to the nearest (non-degenerate) critical point of F , and hence also has no guarantee for finding roots.See the experimental results for detail.

2.2.
Random relaxed Newton's method.A direct variant of Newton's method is Relaxed Newton's method.It works as follows.We choose a nonzero complex number α, and consider the iterative procedure: Again, this gives rise to an algebraic dynamical system on the Riemann sphere, and hence by [11] does not have guarantee to finding roots of polynomials.However, see the discussion in the subsection on Newton's flow about how the basins of attraction for Relaxed Newton's method behave when α goes to 0. A surprising fact is that if one adds randomness into the design, then the method can find roots of polynomials.The iterative procedure, named Random Relaxed Newton's method, in this case is the following: where α n is randomly chosen in an appropriate manner.
We recall here the relevant result in [18] on convergence of Random Relaxed Newton's method when applied to finding roots of polynomials.
Let α n be randomly chosen from the uniform distribution on the set {α ∈ C : Then the Random Relaxed Newton's method, applied to finding roots of the polynomial P (z), will converge to a root of P (z) for all -except a finite number of exceptions -initial points z 0 .

2.3.
Backtracking New Q-Newton's method.In this section, we recall enough details on both theoretical and practical aspects of BNQN.

2.3.1.
Heuristics and some delicate issues.A tendency of Newton's method for optimization, applied to find minima of a function F : R m → R, is that if the initial point z 0 is close to a (non-degenerate) critical point z * of F (z), then the sequence constructed by Newton's method will converge to z * .Hence, if z * is a saddle point or local maximum, then Newton's method is not desirable for finding (local) minima.In the special case where F (z) is a quadratic function whose Hessian is invertible, then for every initial point z 0 , the point z 1 constructed by Newton's method will be the unique critical point z * = 0.
To overcome this undesirable behaviour of Newton's method, a new variant called New Q-Newton's method (NQN) was recently proposed in [20].This consists of the following two main ideas: -Add a perturbation δ||∇F (z)|| τ Id to the Hessian ∇ 2 F (z), where τ > 0 is a constant, and δ is appropriately chosen within a previously chosen set {δ 0 , . . ., δ m } ⊂ R. Thus, we consider a matrix A(z) = ∇ 2 F (z) + δ||∇F (z)|| τ Id, instead of ∇ 2 F (z) as in the vanilla Newton's method.This has two benefits.First, it avoids the (minor) difficulty one encounters in Newton's method if ∇ 2 F (z) is a singular matrix.Second, more importantly, it turns out that if the δ 0 , . . ., δ m are randomly chosen, then NQN can avoid saddle points.Since the perturbation δ||∇F (z)|| τ Id is negligble near non-degenerate local minima, the rate of convergence of NQN is the same as that of Newton's method there.
-If one mimics Newton's method, in defining z n+1 = z n −A(z n ) −1 .∇F(z), then the constructed sequence still has the same tendency of convergence to the nearest critical point of F (z).This is remedied in NQN by the following idea: We let B(z) be the matrix with the same eigenvectors as A(z), but whose eigenvalues are all absolute values of the corresponding eigenvalues of A(z).The update rule of NQN is The precise definition of the algorithm NQN is given in the next subsection.If F (z) is a C 3 function, then NQN applied to F (z) can avoid saddle points, and has the same rate of convergence near nondegenerate local minima.
However, NQN has no global convergence guarantee.In [19], a variation of NQN, that is BNQN, was introduced and shown to keep the same good theoretical guarantees as NQN, with the additional bonus that global convergence can be proven for BNQN for very general classes of functions F (z): functions which have at most countably many critical points, or functions which satisfy a Lojasiewicz gradient inequality type.The main idea is to incorporate Armijo's Backtracking line search [2] into NQN.For the readers' convenience, we recall here the idea of Armijo's Backtracking line search.Let F : R m → R be a C 1 function.Let z, w ∈ R such that < ∇F (z), w > is strictly positive.Then, there exists a positive real number γ for which F (x − γw) − F (x) ≤ − < ∇F (x), γw > /3.If we choose γ by a backtracking manner (that is, start γ from a positive number, and then reduce it exponentially until the above mentioned inequality is satisfied), then the procedure is called Armijo's Backtracking line search.
We notice that there are some delicate issues when doing this incorporation between NQN and Armijo's Backtracking line search, which rely crucially on the fact that we are using a perturbation of the Hessian matrix here.First, in BNQN, the choice of δ from among {δ 0 , . . ., δ m } is more complicated than NQN.Second, the analog Backtracking line search for Gradient descent is not yet known to be able to avoid saddle points, even though there is a slight variant which can avoid saddle points, and experiments support that the method should be able to avoid saddle point.Third, a priori the learning rate one finds by Armijo's Backtracking line search can be smaller than 1, and hence the rate of convergence can a priori slower than being quadratic.
The precise definition of BNQN is given in the next subsection.
2.3.2.Algorithms: NQN and BNQN.Here we present the basic versions of NQN and BNQN.Many more variations can be found in [19].
Let A : R m → R m be an invertible symmetric square matrix.In particular, it is diagonalisable.Let V + be the vector space generated by eigenvectors of positive eigenvalues of A, and V − the vector space generated by eigenvectors of negative eigenvalues of A. Then pr A,+ is the orthogonal projection from R m to V + , and pr A,− is the orthogonal projection from R m to V − .As usual, Id means the m × m identity matrix.
BNQN includes a more sophisticated choice of δ in NQN, together with a combination of Armijo's Backtracking line search.For a symmetric, square real matrix A, we define: sp(A) = the maximum among |λ|'s, where λ runs in the set of eigenvalues of A, this is usually called the spectral radius in the Linear Algebra literature; and minsp(A) = the minimum among |λ|'s, where λ runs in the set of eigenvalues of A, this number is non-zero precisely when A is invertible.
One can easily check the following more familiar formulas: sp(A) = max ∥e∥=1 ∥Ae∥ and minsp(A) = min ∥e∥=1 ∥Ae∥, using for example the fact that A is diagonalisable.
We recall that a function F has compact sublevels if for all C ∈ R the set {z : Algorithm 2 has two different versions depending whether the objective function F has compact sublevels or not.In [8], we introduced a new variant, named BNQN New Variant, with a new parameter, which includes both these versions as special cases.Indeed, θ = 0 in BNQN New Variant is the version of BNQN for functions with compact sublevels, and θ = 1 in BNQN New Variant is the version of BNQN for general functions.
A main theoretical result for finding roots of meromorphic functions by BNQN.As mentioned previously, BNQN has strong theoretical guarantees for some big classes of functions in any dimension.However, to keep the presentation concise, here we present only one main result relevant to the question pursued in this paper, that of finding roots of meromorphic functions in 1 complex variable.For more general results, the readers can consult [19], [20].
Given an initial point z 0 ∈ C, which is not a pole of g, we let {z n } be the sequence constructed by BNQN New Variant applied to the function F with initial point z 0 .If the objective function has compact sublevels (i.e. for all C ∈ R the set {(x, y) ∈ R 2 : F (x, y) ≤ C} is compact), we choose θ ≥ 0, while in general we choose θ > 0.
1) Any critical point of F is a root of g(z)g ′ (z) = 0.
2) If z * is a cluster point of {z n } (that is, if it is the limit of a subsequence of {z n }), then z * is a critical point of F .Moreover, in this case the whole sequence {z n } converges to z * .
3) If F has compact sublevels, then {z n } converges.4) Assume that the parameters δ 0 , δ 1 , δ 2 in BNQN New Variant are randomly chosen.Assume also that g(z) is generic, in the sense that {z ∈ C : g(z)g"(z) = g ′ (z) = 0} = ∅.There exists an exceptional set E ⊂ C of zero Lebesgue measure so that if z 0 ∈ C\E, then {z n } must satisfy one of the following two options: Option 1: {z n } converges to a root z * of g(z), and if γ 0 = 1 in the algorithm then the rate of convergence is quadratic.
Option 2: For BNQN, Theorem 2.2 is proven in [20].The stated version for BNQN New Variant is from [8].An example for which part 3 of Theorem 2.2 can apply is when g is a polynomial or g = P/Q where P, Q are polynomials and P has bigger degree than Q (a special case is when Q = P ′ , in which case the zeros of g are exactly that of P , with the advantage that they all have multiplicity 1).[19].The implementation is flexible in that one does not need precise values of the gradient and Hessian of the function F , but approximate values are good enough.Experiments also show that the performance of BNQN is very stable, with respect to its parameters and to the values of the objective function and its gradient and Hessian matrix.Some experiments later in this paper will illustrate this.

Newton's flow.
A general strategy when studying an iterative method is that one also studies its flow counterpart.Usually, the flow counterpart is a smooth version of the (discrete) iterative method.In this subsection, we briefly review several variants of Newton's flow.Newton's flow means solutions of certain ordinary differential equations which we will define later.
The authors observed that the basins of attractions for BNQN method look smoother than that for Newton's method or Backtracking line search for Gradient Descent.This observation inspired us to compare the basin structure for BNQN with that for Newton's flow.In the following, we consider three types of continuous-time dynamical systems.
2.4.1.Newton's flow for f .The first Newton-type flow is the ordinary differential equation of complex-valued function z(t) = x(t) + iy(t) of one real variable t: with initial value z(0) = z 0 .Here, f (z) is a holomorphic function.Note that the right-hand side of (1) has singularity at {f ′ = 0} \ {f = 0}.
If a global solution z(t) for the ODE above exists, then it satisfies Thus, f (z(t)) = f (z 0 )e −t and f (z(t)) → 0 as t → ∞.By this calculation, we expect that the solution of the ODE (1) approximates a zero of the target function f.Neuberger et al. showed that we can interpret the ODE (1) as the Sobolev gradient of F (z) = ∥f (z)∥ 2 /2 relative to the natural Riemannian metric.See Appendix A of [15].
From our point of view, it is worth noting that the Euler discretization of the ODE (1) gives the Relaxed Newton's method.Let α > 0 be a constant time step.The Euler discretization gives us with initial condition z(0) = z 0 .Thus, we have the recursive formula z(t + α) = z(t) − αf (z(t))/f ′ (z(t)), which is the relaxed Newton's method.In particular, if the time step α is 1, then the Euler discretization is precisely the same as the usual Newton's method.It is shown in [9][12] that if f is a rational function, then the union of basins of attraction for Newton's flow of all the roots of f has full measure.It is also shown that in case f is a polynomial, there is a small value of α for which Relaxed Newton's method has convergence guarantee to find all roots of f .Note that the choice of α heavily depends on the polynomial f .See [4] for a good review on some results on this direction, as well as dynamics of meromorphic functions in general.
Compared to the original Newton's method, the local convergence rate of Newton's flow is slower.Namely, if z = a is a simple zero of f , then we have |z(t) − a| = O(e −t ) as t → ∞.Thus, Newton's flow converges to a simple root with order 1.On the other hand, the original Newton's method converges to a simple root with order 2 as is well-known.
From now on, this method will be called Newton's flow.
2.4.2.Newton's flow for f /f ′ .The second Newton-type flow is the ordinary differential equation of complex-valued function z(t) of one real variable t: with initial value z(0) = z 0 .Here, g(z) = f (z)/f ′ (z) is a meromorphic function and f (z) is a holomorphic function.The roots of g are exactly the roots of f , with the advantage that their multiplicities are all 1.
The ODE ( 2) is obtained by substituting g = f /f ′ for f in the ODE (1).By a similar calculation, we can show that g(z(t)) = g(z 0 )e −t → 0 as t → ∞ if a global solution z(t) exists.Since the zeros of g is same as the zeros of f, we expect the flow z(t) will converge to a root of f.
The conceptual difference between two ODEs (1) and ( 2) is the singular sets.The singular set for (1) From now on, this version will be named Newton's flow vFraction.
2.4.3.Newton's flow for the optimization version.The last Newton's flow is applicable to Newton's method v2, the optimization version of Newton's method.This Newton-type flow is the ordinary differential equations related to the function z(t) = x(t) + iy(t) of one real variable t: d dt with initial value z(0) = z 0 .Here, F (z) = ∥f (z)∥ 2 /2 and f (z) is a holomorphic function.Also, by identifying z = x + iy, we denote by where F x = ∂F ∂x , F y = ∂F ∂y and so on.The right-hand side of ( 3) is defined outside the singular set {det ∇ 2 F = 0}.
If this ODE has a global solution z(t), then we have Therefore, the global solution satisfies ∇F (z(t)) = e −t ∇F (z 0 ), and hence we have ∇F (z(t)) → 0 as t → ∞.This shows that we can expect to obtain an approximation of points where ∇F vanishes by solving the ODE.
We are now interested in the local behavior near simple zeros of f .Suppose f (a) = 0. We can assume that a = 0 by changing the coordinate.Since f is holomorphic, we can write f (z) = cz + g(z) where c ̸ = 0 is a complex number and g(z) is a holomorphic function defined near z = 0 which satisfies g(0) = g ′ (0) = 0 and g ′′ (0) ̸ = 0. Let us denote f (z) = u(z) + iv(z) using real functions.Then u(0) = 0, v(0) = 0, and Therefore, the ODE (3) near z = 0 is reduced to be d dt Thus, the flow that starts near a simple root has a time-global solution and will converge to the simple root.
From now on, this Newton's flow is named Newton's flow vOptimization.
2.4.4.Implementation details.In this subsubsection, we present our implementation of three Newton-type flows.There are various known methods for solving differential equations numerically, but in this paper, we use the Explicit Runge-Kutta method of order 5(4).This method is the default method used in scipy.integrate.solveivp function.
The Explicit Runge-Kutta method of order 5(4) is discussed in [7].This embedded method automatically adjusts the time step size to enable fast and accurate numerical computation, but for simplicity, we review the classical Runge-Kutta method here.
We want to numerically solve the autonomous ODE ẏ = g(y), y(0) = y 0 with a given initial condition.Here, y and g are R n -valued, or they are C-valued functions.Fix a small time step h > 0, and define Y (0) = y 0 .
If we already define the value Y (t), then define where It is known that the truncation error is proportional to h 5 for every time step (if the solution y(t) is sufficiently smooth).Thus, if we calculate until t = T, then the error |y(T ) − Y (T )| is proportional to h 5 • T /h = T h 4 .See Section 5.4 and following of [5] for more details.We are interested in the global structure of basins for Newton's flow.Therefore, we performed the following experiment for every three ODEs above.For several initial values, we numerically integrated the trajectories using scipy.integrate.solveivp.The time step was h = 0.01 and the endpoint of integration was T = 100.The same color was applied to the initial values whose trajectory reached within a certain distance from the same root of f .If the trajectories are away from all the roots, then we colored their initial values black.An illustration is shown in Fig. 1.

2.5.
Voronoi's diagrams.Voronoi's diagrams are a classical object [21] [22], with an ever increasing domain of applications: geophysics, surface metrology, hydrology, archaeology, dialectometry, political sciences, biology, ecology, ethology, computational chemistry, astrophysics, computational fluid dynamics, computational physics, medical diagnosis, epidemiology, polymer physics, materials science, aviation, architecture, urban planing, mining, robotics, wireless network, computer graphics, machine learning, and a lot more.For example, the skin pattern of a giraffe has Voronoi-type diagrams.A good place to start with Voronoi's diagrams and applications is [23] and its references.
There are many versions of Voronoi's diagrams.In this paper, we consider the most basis version, that of a finite set of points z 1 , . . ., z m in R 2 , with the Euclidean distance d E (., .).If z ∈ R 2 is such that there is a unique i ∈ {1, . . ., m} for which d E (z, z i ) < min j̸ =i d E (z, z j ), then we say that z belongs to the Voronoi's cell of the point z i .The boundary consists of line intervals where the maximum of min j d E (z, z j ) is attained at two or more points z i 's.
2.6.Stochastic root finding.Because of various reasons, it is necessary and interesting to consider root finding problem for not only the determnistic case f (z) = 0, but also the stochastic version g(z, ξ), where ξ is a random variable.
We assume that the expectation of g, with respect to the variable ξ, is f , that is E(g(z, .))= f (z).
In stochastic root finding, we aim to test the robustness of a root finding method IM in a stochastic environment.When given a function h in the variable z, the root finding method IM will construct a sequence defined by the formula z n+1 = IM (h(z), z n ).More precisely, we will follow the following procedure: -Generate a sequence {ξ n } from the distribution for the random variable ξ.
-Choose an initial point z 0 and construct a sequence: We hope that {z n } will converge to a root of f (z).To ensure the stability, we need the variation of g is small: V ar(g(z, .)) is small.

Experimental results
Here we discuss the settings for our experiments, and then present the experiments.
3.1.Settings.We present the settings of our experiments.
1) The functions: We consider many functions the geometric configurations of whose roots vary (for example, some functions have roots which are vertices of convex polygons, while some functions have some roots inside the convex hull of other roots, some functions have three roots on the same line, and some functions have simple roots while other have non-simple roots).The list of the functions is as follows: The function f 23 (z) has infinitely many roots.In the region −10 ≤ x, y ≤ 10, it has 8 (approximate) roots: 0.01453348 + 0.24577632i, −1.79690338 − 0.16311646i, 2.65293461 − 2.52795741i, 2.70778504 + 2.4386467i, −7.27782023−4.1230358i,−7.26685729+4.13462414i,9.62682067− 4.62305718i, and 9.63392763 + 4.61683271i.
In some experiments, we also consider functions of the form g = f /f ′ , where f is one of the above functions.
In stochastic root finding experiments, we consider functions of the form g(z, ξ) = f (z) + ϵξ(z 3 + 2z − 5).Here f (z) is one of the above functions, ξ is a random variable with normal distribution (E(ξ) = 0, V ar(ξ) = 1), ϵ > 0 is a small number.We can check that E(g(z, .))= f (z).However, in experiments involving Newton's method vOptimization and BNQN we need to work with the functions . Therefore, we should choose ϵ small so that E(G(z, .)) is close enough to F (z).
2) The experiments: We will compute basins of attraction when using different methods (iterative or flow) to certain complex-valued functions in 1 complex variable (listed in the previously).We also draw pictures of (reduced) Voronoi's diagrams, where the multiplicity of a point is disregarded.
3) Colours and interpretations of the results: If the factorisation of a function f is written as (note that the order z 1 , z 2 , . . . is important) . ., then we will colour the basins of attraction as follows: that for the first root has the green colour, that for the second root has the yellow colour, that for the third root has the blue colour, that for the fourth root has the red colour, that for the fifth root has the pink colour, that for the sixth root has the cyan colour, that for the seventh root has the orange colour, that for the eighth root has the purple colour.Other points (e.g. points whose corresponding sequence goes to infinity or does not converge to one of the roots) will have the black colour.
For example, for the function f 10 (z) = z 2 (z − i)(z − 1 − i), by our convenience, the first root is 0, the second root is i and the third root is 1 + i.Therefore, the basin of attraction for 0 will have the green colour, the basin of attraction for i will have the yellow colour, and the basin of attraction for 1 + i has the blue colour.
In Voronoi's diagrams, we also follow the same convenience.4) Implementation of the codes: Implementation of BNQN and Newton's flows are as discussed in Section 2. Implementation of Newton's method and Random Relaxed Newton's method is straightforward from their descriptions.
For BNQN, most of the case we will run BNQN New Variant with θ = 0, and this will be reported in the experimental results as BNQN.In some cases, we also run with BNQN New Variant with θ = 1, and report it as BNQN v2.
To avoid overflow errors, we use the mpmath library [14].

5) Experimental procedures:
We will choose a grid of 240 × 240 points in [−10, 10] × [−10, 10], whose center is randomly chosen in a small neighbourhood of (0, 0).Each of the points z 0 in the grid will be the initial point to run the concerned iterative method.We will run a maximum of 10000 iterations, and can stop earlier if either the function value is smaller than a threshold ϵ or the gradient of the objective cost function is smaller than the threshold (for BNQN and Newton's method vOptimization).If z n is the last constructed point, we declare that z n belongs to the basin of attraction for a root z * if the Euclidean distance d(z n , z * ) is smaller than a threshold.We choose the threshold to be 10 −6 or smaller.
This manner of determining the basins of attraction is reasonable.First, we know that if z n is close enough to a root z * , then the considered iterative methods will converge to the root z * .Second, in the stochastic root finding setting, as explained in Section 3, E(G(z, .)) is not the same as F (z), and hence we do not expect that z n converges precisely to a root of f (z) (or more generally, a minimum of F (z)).However, if E(G(z, .))− F (z) is small, then we expect that the minizers of E(G(z, .))are within a small neighbourhood of that of F (z). Since we expect that z n will converge to minimizers of E(G(z, .)), it follows that z n will converge to within a small neighbourhood of a minimizer of F (z). Explicitly, in g(z, ξ) = f (z) + ϵξ(z 3 + 2z − 5), we choose ϵ to be 10 −4 .In this case, we reduce the threshold to declare that the point z n belongs to the basin of attraction of a root z * to be 10ϵ, which is 10 −3 .

Conclusions and further directions
From experiments, we draw the following conclusions, when finding the roots of f : -Newton's method does not reflect well the geometric configuration of the roots of f (i.e.Voronoi's diagrams).Newton's flow method, while more smooth than Newton's method, also does not reflect well the geometric configuration of the roots.The basins of attraction for Random Relaxed Newton's method are distorted versions of Newton's method.
Both Newton's method and Random Relaxed Newton's method seem not work well (take too long time to run) with functions of the form f /f ′ , where f is transcendental.A reason could be that in these cases, it can take many iterates for these methods to converge, or they are more likely to diverge, in both cases it takes more time to verify the STOP conditions.
-Newton's method for Optimization performs poorly.Newton's flow for Optimization, except for some small open sets which are basins of attraction of the critical points of f , produces pictures which match well with Voronoi's diagrams.This is in big contrast to its discrete version.However, the method is quite heavy, and takes too long time for the transcendental examples.
-Newton's flow for f /f ′ , on the other hand, reflects well Voronoi's diagrams, in contrast to Newton's flow for f .The reason could be as follows: The RHS in Newton's flow for f /f ′ is which involves the second derivative and hence reflects better the curvature of the function landscape.
-BNQN reflects well Voronoi's diagram.In particular, BNQN v2 (where we choose θ = 1) is more similar to Newton's flow for f /f ′ .-If some of the roots are contained inside the convex hull of other roots, then the similarity between these iterative methods and Voronoi's diagrams becomes less apparent.A reason could be that indeed one needs to use a more appropriate metric than the usual Euclidean metric when constructing Voronoi's diagrams.-BNQN is more robust against stochastic root finding than Newton's method, Random Relaxed Newton's method and Newton's method for Optimization.
The above conclusions suggest a couple of natural follow ups, which are beyond the current paper and are left for future explorations: -While Newton's flow for f /f ′ and BNQN both involve the second derivatives of f , how they depend on the second derivatives is different.In BNQN, we need to change the signs of negative eigenvalues, and we also need to run an Armijo's Backtracking line search.What is the reason that they produce similar pictures?
-Is it true that in BNQN New Variant, if we choose bigger values for the parameter θ, then we will obtain more smooth pictures?-What should be the metric to be chosen in constructing Voronoi's diagrams, so that we will obtain more similar pictures to basins of attraction?
-For the case of functions with multiple roots, are there corresponding Voronoi's diagrams to the basins of attraction?

1
arXiv:2401.01393v2 [math.OC] 8 Jan 2024 (z − b)(z − c) can provide one some insights into how Voronoi's diagram for three points a (with multiplicity 2), b (with multiplicity 1) and c (with multiplicity 1) should approximately look like.In the experiments, we have several such examples.

Figure 1 .
Figure 1.Basins of attraction for the ODEs (1), (2), and (3) from left to right, applied to a polynomial of degree 3 (the function f 1 listed in Section 3).

Figure 2 .
Figure 2. Streamline plots for the ODE (3), from Figure 1.The figure on the right is an enlarged image of the figure on the left, whose range is 0 ≤ Rez, Imz ≤ 5. We can see a stable equilibrium point near z = 1.966755 + 1.516588i which does not correspond to the roots of f 1 , but a critical point of f 1 .

Figure 3 .
Figure 3. Basins of attraction for finding roots of the function f 1 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 1 .

Figure 4 .
Figure 4. Basins of attraction for finding roots of the function f 2 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 2 .

Figure 5 .
Figure 5. Basins of attraction for finding roots of the function f 3 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 3 .

Figure 6 .
Figure 6.Basins of attraction for finding roots of the function f 4 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 4 .

Figure 7 .
Figure 7. Basins of attraction for finding roots of the function f 5 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 5 .

Figure 8 .
Figure 8. Basins of attraction for finding roots of the function f 6 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points for Newton's flow vOptimization are those in the basin of attraction of critical points of f 6 which are not roots, while the black points for Newton's method for Optimization seem to due to the fact that the cost function is highly degenerate.

Figure 9 .
Figure 9. Basins of attraction for finding roots of the function f 7 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 7 .

Figure 10 .
Figure 10.Basins of attraction for finding roots of the function f 8 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 8 .

Figure 11 .
Figure 11.Basins of attraction for finding roots of the function f 9 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 9 .

Figure 12 .
Figure 12.Basins of attraction for finding roots of the function f 10 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 10 .

Figure 13 .
Figure 13.Basins of attraction for finding roots of the function f 11 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 11 .

Figure 14 .
Figure 14.Basins of attraction for finding roots of the function f 12 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 12 .

Figure 15 .
Figure 15.Basins of attraction for finding roots of the function f 13 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 13 .

Figure 16 .
Figure 16.Basins of attraction for finding roots of the function f 14 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 14 .

Figure 17 .
Figure 17.Basins of attraction for finding roots of the function f 15 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 15 .

Figure 18 .
Figure 18.Basins of attraction for finding roots of the function f 16 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 16 .

Figure 19 .
Figure 19.Basins of attraction for finding roots of the function f 17 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 17 .

Figure 20 .
Figure 20.Basins of attraction for finding roots of the function f 18 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 18 .Department of Mathematics, NTNU, Norway Email address: fornaess@gmail.comDepartment of Mathematics, University of Oslo, Norway Email address: humihqu@gmail.comDepartment of Mathematics, University of Oslo, Norway Email address: tuyentt@math.uio.noChubu University Academy of Emerging Sciences/Center for Mathematical Science and Artificial Intelligence, Chubu University, JapanEmail address: takawatanabe@isc.chubu.ac.jp

Figure 21 .
Figure 21.Basins of attraction for finding roots of the function f 19 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 19 .

Figure 22 .
Figure 22.Basins of attraction for finding roots of the function f 20 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 20 .

Figure 23 .
Figure 23.Basins of attraction for finding roots of the function f 21 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 21 .

Figure 24 .
Figure 24.Basins of attraction for finding roots of the function f 22 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is Voronoi's diagram, central picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, right picture is for BNQN.Row 3: left picture is for Newton's flow, central picture is for Newton's flow vFraction, right picture is for Newton's flow vOptimization.The black points in some of these pictures are those in the basin of attraction of critical points of f 22 .

Figure 25 .
Figure 25.Basins of attraction for finding roots of the function f 23 by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: Voronoi's diagram (using only 8 roots inside the domain (−10, 10) × (−10, 10)).Row 2: left picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 3: left picture is for BNQN, right picture is for BNQN v2.Row 4: left picture is for Newton's flow, right picture is for Newton's flow vFraction.The method Newton's method vOptimization encounters errors, while Newton's flow vOptimization takes too long time to finish.

Figure 26 .
Figure 26.Basins of attraction for finding roots of the function f 7 e z by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for BNQN, right picture is for BNQN v2.Row 3: left picture is for Newton's flow, right picture is for Newton's flow vFraction.Row 4: Voronoi's diagram.The method Newton's method vOptimization encounters errors, while Newton's flow vOptimization takes too long time to finish.

Figure 27 .
Figure 27.Basins of attraction for finding roots of the function f 17 e z by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for BNQN, right picture is for BNQN v2.Row 3: left picture is for Newton's flow, right picture is for Newton's flow vFraction.The method Newton's method vOptimization encounters errors, while Newton's flow vOptimization takes too long time to finish.

Figure 28 .
Figure 28.Basins of attraction for finding roots of the function f 17 /f ′ 17 by different methods.Left picture is for Newton's method, central picture is for Random Relaxed Newton's method, right picture is for BNQN.

Figure 29 .
Figure 29.Basins of attraction for finding roots of the function f 18 /f ′ 18 by different methods.Left picture is for Newton's method, central picture is for Random Relaxed Newton's method, right picture is for BNQN.

Figure 30 .
Figure 30.Basins of attraction for finding roots of the function f 19 /f ′ 19 by different methods.Left picture is for Newton's method, central picture is for Random Relaxed Newton's method, right picture is for BNQN.

Figure 31 .
Figure 31.Basins of attraction for finding roots of the function f 20 /f ′ 20 by different methods.Left picture is for Newton's method, central picture is for Random Relaxed Newton's method, right picture is for BNQN.

Figure 32 .
Figure 32.Basins of attraction for finding roots of the function f 21 /f ′ 21 by different methods.Left picture is for Newton's method, central picture is for Random Relaxed Newton's method, right picture is for BNQN.

Figure 33 .
Figure 33.Basins of attraction for finding roots of the function f 22 /f ′ 22 by different methods.Left picture is for Newton's method, central picture is for Random Relaxed Newton's method, right picture is for BNQN.

Figure 34 .
Figure 34.Basins of attraction for finding roots of the function f 23 /f ′ 23 , f 24 /f ′ 24 and f 25 /f ′ 25 by BNQN.For this case, Newton's method and Random Relaxed Newton's method either encounter errors or take very long time to finish.

Figure 35 .
Figure 35.Basins of attraction for finding roots of the function f 24 /f ′ 24 and f 25 /f ′ 25 by BNQN v2, which are more smooth than that by BNQN in Figure 34.

Figure 36 .
Figure 36.Basins of attraction for finding roots of the stochastic function f 1 + ϵξ(z 3 + 2z − 5) by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, central picture is for BNQN, right picture is for BNQN v2.

Figure 37 .
Figure 37. Basins of attraction for finding roots of the stochastic function f 5 + ϵξ(z 3 + 2z − 5) by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, central picture is for BNQN, right picture is for BNQN v2.

Figure 38 .
Figure 38.Basins of attraction for finding roots of the stochastic function f 14 +ϵξ(z 3 +2z −5) by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, central picture is for BNQN, right picture is for BNQN v2.

Figure 39 .
Figure 39.Basins of attraction for finding roots of the stochastic function f 21 +ϵξ(z 3 +2z −5) by different methods.Pictures are referenced to from top to bottom, from left to right.Row 1: left picture is for Newton's method, right picture is for Random Relaxed Newton's method.Row 2: left picture is for Newton's method vOptimization, central picture is for BNQN, right picture is for BNQN v2.