Dw-03299 RR247

IRIA
laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B.P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

<figure><img src="image_1.png" /></figure>

OPTIMAL SEARCH
FOR THE ZERO OF
THE (n - 1)st DERIVATIVE

Laurent HYAFIL

Rapport de Recherche N° 247

Août 1977
OPTIMAL SEARCH FOR THE ZERO OF THE (n-1)st DERIVATIVE

Laurent Hyafil

IRIA/LABORIA

Résumé :

La recherche « Fibonacci » introduite par Kiefer pour la recherche du maximum d'une fonction unimodale est une généralisation de la méthode de la dichotomie pour trouver le zéro d'une fonction monotone. Dans cet article nous généralisons le résultat de Kiefer à la recherche du zéro de la dérivée (n-1)ième d'une fonction à partir des seules valeurs de la fonction. On peut utiliser cet algorithme pour trouver tous les zéros d'une fonction dont la dérivée (n-1)ième ne change de si-gne qu'une fois à partir des seules valeurs de la fonction.

Abstract :

The well known method of binary search to find the zero of a monotone function has been extended by Kiefer to the search for the maximum of an unimodal function. He introduced Fibonacci search which be proved to be optimal. In this paper we generalize Kiefer's result to finding a zero of the (n-1)st derivative of a function, using only the values of this function. Applications of this algorithm include root separation and zero-finding algorithms for functions whose (n-1)st derivative changes sign only once.

Keywords : Analytic complexity, zero-finding algorithm.
<figure><img src="image_2.png" /></figure>

# 1. Introduction

The problem of finding zeroes of the $(n-1)^{st}$ derivative of a function without knowing the derivatives has been considered by Brent [3] who gives an algorithm with superlinear convergence assuming the starting values are "sufficiently good". In this paper we give an algorithm with guaranteed convergence assuming the $(n-1)^{st}$ derivative changes sign only once*.

The algorithm is an extension of the Fibonacci search (Kiefer [1]) and uses some generalisations of Fibonacci numbers. The algorithm is proved to be optimal among the algorithms which use only $n$ values of the function to determine the next point of evaluation.

In section 2 we describe the general algorithm and we study its optimality in section 3. Optimal algorithms are provided when $n$ is odd and $n=2$ and 4. For larger odd number the algorithm would require further properties of generalized Fibonacci numbers. The method used to prove optimality is derived from previous work by the author and H.T. Kung (which is best described in [2]) on the zero finding problem for an asynchronous multiprocessor.

# 2. Finding the zero of the $(n-1)^{st}$ derivative of $f$.

We first give a theorem which is used to prove the correctness of the algorithm :

**Theorem 1** : Let $f$ be a function from $[a,b] \subset \mathbb{R}$ to $\mathbb{R}$ such that $f^{(n-1)}$ is continuous, has a unique zero $c$ and $f(a)^{(n-1)} < 0$ and $f(b)^{(n-1)} > 0$.

<footnote>*Note : Throughout this paper when we write "the unique zero of a function" we implicitely assume that this function changes of sign only once.</footnote>
Let $x_1, x_2, \dots, x_n \in [a, b]$ such that $a < x_1 < x_2 < \dots < x_n < b$.

Let $L(x_1, x_2, \dots, x_n, f) = \sum_{i=1}^{n} \left[ \frac{f(x_i)}{\prod_{j=1, n} (x_j - x_i)} \right]^{\frac{1}{n}}$ the divided difference of order n.

if $L(x_1, x_2, \dots, x_n, f) \le 0$ then $c \in [a, x_1]$,

if $L(x_1, x_2, \dots, x_n, f) \ge 0$ then $c \in [x_n, b]$.

**Proof:** In order to prove Theorem 1 we need to prove the following lemma :

**Lemma 1 :** Let $f$ be a function from $[a,b] \subset \mathbb{R}$ to $\mathbb{R}$ such that $\forall x \in [a,b], f^{(n-1)}(x) > 0$, then for every set of distinct $x_i \in [a,b]$, $i=1,n: L(x_1, x_2, \dots, x_n, f) > 0$.

**Proof of lemma 1 :** Let $g(x) = f(x) - \sum_{i=1,n} \left[ f(x_i) \prod_{j \neq i} \frac{(x-x_j)}{(x_i-x_j)} \right]^{\frac{1}{n}}$ then $g(x_i) = 0$ for $i = 1,n$. By repetitive applications of Rolle's theorem we deduce that there exists $d, d \in [a,b]$, $g^{(n-1)}(d) = 0$. The result follows from : $f^{(n-1)}(d) = g^{(n-1)}(d) + (n-1)L(x_1, x_2, \dots, x_n, f)$

Then assume the zero c of $f^{(n-1)}$ is in $[a, x_1]$. Then $f^{(n-1)}$ is strictly positive in $[x_1, b]$. Therefore if $L(x_1, x_2, \dots, x_n, f) \le 0$ then $c \in [a, x_1]$ by lemma 1. A symmetric proof for $[x_n, b]$ completes the proof of theorem 1. []

From Theorem 1 we deduce an iterative algorithm to compute to within $\varepsilon$ the unique zero of the $(n-1)^{st}$ derivative of a function $f$. At each step of the algorithm we select a $n^{th}$ point, compute the value of $f$ at this point and $L(x_1, x_2, \dots, x_n, f)$ then we discard one of the two intervals $[a, x_1]$ or $[x_n, b]$ (or both) and reduce the length of the interval of uncertainty until it becomes smaller than $\varepsilon$.
Algorithm Z_n: zero finding of the (n-1)st derivative of f.

Initialisation :

- select all distinct x₁, x₂, ..., xₙ₋₁ in ]a, b[

- compute f(x₁), ..., f(xₙ₋₁)

- u ← a

- for i=1, ..., n-1 compute $L_i = \frac{f(x_i)}{\prod_{j=1, n-1}^{j \neq i} (x_j - x_i)}$

- L ← $\sum_{i=1, n-1} L_i$

Loop : until |b-a| < ε do

begin - select $x_{n+1}$ in ]a, b[ different from $x_1, x_2, ..., x_n$ and relabel the n $x_i$'s in ]a, b[ in increasing order so that the new index of $x_{n+1}$ is k. Rearrange the corresponding values of the $L_i$.

- for i=1, n i≠k compute $L_i ← L_i \frac{(u - x_i)}{x_k - x_i}$.

- compute $L_k ← \frac{f(x_k)}{\prod_{j=1, n} (x_j - x_k)}$

- compute $L ← \sum_{i=1}^{n} L_i$

- if $L \ge 0$ then $b \leftarrow x_n$; $u \leftarrow x_n$; $L \leftarrow L-L_n$;

else $a \leftarrow x_1$; $u \leftarrow x_1$; $L \leftarrow L-L_1$;

end
In the case where L is null we discard the two extremal intervals and compute the values of f at two new points before evaluating the new value of L.

In order to describe fully algorithm Z_n we must give the algorithm A_n which is used to provide a new point inside the interval of uncertainty.

$$- \underline{Algorithm A}_{2q+1}$$

We maintain the $2q+1$ points inside the interval of uncertainty in positions determined by $\alpha, \beta, T$ and the length of the interval of uncertainty.

- the points are in positions :

$$iT \text{ for } i=1, \dots, \beta-1$$

$$\beta T + iT/2 \quad i=0, \dots, \alpha$$

$$\beta T + \alpha T/2 + iT \quad i=1, \dots, \beta-1$$

- the length of the interval is $2\beta T + \alpha T/2$

When reducing the interval of uncertainty the position of the new point is determined by the positions of the $2q+1$ points given by the above formula where $\alpha, \beta, T$ are computed from the old values of $\alpha, \beta, T$ by :

$$\beta \leftarrow \beta-1$$

$$\alpha \leftarrow \alpha+2$$

$$\underline{if} \quad \alpha = 2q+2 \quad \underline{then} \quad \begin{cases} T \leftarrow T/2 \\ \alpha \leftarrow 0 \\ \beta \leftarrow q+1 \end{cases}$$
Analysis of $Z_{2q+1}$ (with $A_{2q+1}$)

It is clear that using $A_{2q+1}$ the length of the interval of uncertainty is divided by 2 after $q+1$ executions of the loop of $Z_{q+1}$ ($q+1$ evaluations of $f$). So the number of evaluations of $f$ necessary to reduce the length of the interval from $L$ to 1 is bounded from above by $2q+1+(q+1) \lfloor \log_2 L \rfloor$.

Algorithm $A_{2q}$

The algorithm $A_2$ was given by Kiefer [1] and we give $A_4$ since the only known algorithm $A_{2q}$ for $q > 2$ is to use $A_{2q+1}$. To describe $A_4$ we assume the length of the interval of uncertainty is $T_n$ where :

$T_n = \psi^n$, where $\psi$ satisfies $\psi^3 = \psi+1$. Hence $T_n = T_{n-2} + T_{n-3}$.

From this recurrence we derive $T_n - T_{n-1} = T_{n-2} - T_{n-4} = T_{n-5}$ and $T_n - T_{n-4} < T_{n-1}$. Therefore if the distances of the points to the left extremity are $T_{n-1}$, $T_{n-2}$, $T_{n-3}$, and $T_{n-4}$, the distances to the right extremity are $T_{n-5}$, $T_{n-3}$, $T_{n-2}$ and $T'_{n-1} < T_{n-1}$. We conclude that when one of the extremal intervals is discarded it is possible to choose a new point which will insure the same distribution (within a symmetry) and with a value of $n$ decreased by 1; if the right interval is discarded the point is chosen at distance $T_{n-5}$ from the left extremity, else it is chosen at distance $T_{n-4}$ from the right extremity.

Analysis of $Z_4$ (with $A_4$)

Starting from $T_n$ each execution of the loop reduces the size of the interval from less than $T_i$ to less than $T_{i-1}$. Hence the number of evaluations of $f$ necessary to reduce the size of the interval from 1 to 1 is bounded from above by $4 + \lfloor \log_2 L \rfloor$.
### 3. Optimality of the algorithms

In this section we first prove that $A_n$ are optimal for n odd and n=2,4.
Then we establish that $Z_n$ is optimal among the algorithms which memorize
only n values of the function.

**Theorem 2 :** *Algorithm A*<sub>*n*</sub> are optimal for *n* odd and *n*=2 and 4.

*Proof:* We first relate the problem of finding for the unique zero of the (n-1)st derivative to the problem of finding the unique zero of function with an asynchronous n-processor. An asynchronous n-processor is constituted of n processors which obey the following rules :

(a) at a given time at most one processor has completed its allocated job.

(b) once a processor has completed its job it can access common variables (and only at that time) and it is allocated a new job.

(c) at the beginning of the computation a job is allocated to each of the n processors.

Iterative methods using an asynchronous n-processor are studied in [2].

**Lemma 2 :** Let $(Z_n, A_n)$ be an algorithm for finding the unique zero of the $(n-1)$st derivative. There exists an algorithm for finding the unique zero of function with an asynchronous n-processor with a number of function evaluations not exceeding the number of function evaluations of $(Z_n, A_n)$.

*Proof:* The states of an algorithm for finding the unique zero of a function with n processors can be given by n points inside the interval of uncertainty. A point will be said *real* if a processor is busy evaluating the function at this point and *virtual* otherwise. All virtual points are paired with a processor working outside the interval of uncertainty. When a processor outside the interval of uncertainty completes its job the corresponding virtual point becomes *real* and the job of the processor is to evaluate the value of the function at this point.
When a processor completes its job inside the interval of uncertainty, the value of f is known and according to its sign we can eliminate from the interval of uncertainty one of the two intervals I and J this point determines. Call I the eliminated interval and assume that I contains q points (real or virtual). To describe the zero finding algorithm we explicit the choice of q points inside J, then choose one of them arbitrarily to be real, the q-1 others being virtual.

In order to choose these q points we use $A_n$. It is important to notice that the n points inside the interval of uncertainty are always, by this algorithm, configurations of points that have been generated by $A_n$. Assume I is on the right of J, the following algorithm determines the q new positions in J.

**initialisation** : the n initial positions are the n positions of points in the interval of uncertainty (IuJ)

**loop** : **until** no more position is left in I

do take the right-most position in I and replace it by the position given by $A_n$ when the interval of uncertainty IuJ is reduced on its right.

This simulation argument builds the asynchronous algorithm for zero finding using $A_n$. Hence it uses no more evaluations than $Z_n$ with $A_n$. □

To complete the proof of theorem 2 we need lower bounds on the number of evaluations required by asynchronous n processors for the zero finding problem.

**Lemma 3** : Let $F(L)$ be the number of evaluations necessary to reduce the interval of uncertainty where lies the unique zero of a function f from L to a length < 1 with an asynchronous n-processor, then :

if $L > T$, $F(L) > p$ where $T(L)$ satisfies :

$$ \begin{cases} T_i = T_{i-\left\lfloor \frac{(n+1)}{2} \right\rfloor} + T_{i-\left\lceil \frac{(n+1)}{2} \right\rceil} & \text{for } i \ge n \\ T_i = 1 & \text{for } 0 \le i \le n-1 \end{cases} $$
**Proof** : the proof is by induction on p

- if $p \le n$ the proof is obvious (rule (c) of asynchronous n-processors).

- assume $p > n$ and let $L$ be an interval of uncertainty of length $L$.

The n processors are evaluating $f$ at points $P_1, P_2, \dots, P_n$ numbered in increasing order. Assume the processor at point $P_{\lfloor \frac{n}{2} \rfloor}$ gives an answer :

- if the remaining interval is on its right, and if it is larger than $T_{p-\lfloor \frac{n}{2} \rfloor}$, then since $\lfloor \frac{n}{2} \rfloor$ evaluations have already been performed $F_n(L)>p$.

- if it is on its left, and if it is larger than $T_{p-(n+1-\lfloor \frac{n}{2} \rfloor)}$, then $F_n(L)>p$. If $L = T_{p-\lfloor \frac{n}{2} \rfloor} + T_{p-(n+1-\lfloor \frac{n}{2} \rfloor)}$, then $F_n(L)>p$. $\square$

For $n=2q+1$ the recurrence of lemma 3 becomes for $i$ sufficiently large :

$$T_i = 2 T_{i-(q+1)}$$

therefore from lemmas 2 and 3 we deduce that the interval of uncertainty can be at most divided by 2 every $q+1$ evaluations. Hence algorithm $A_{2q+1}$ is optimal. Similarly we deduce that $A_2$ and $A_4$ are optimal. $\square$

**Theorem 3** : $Z_n$ is optimal among the algorithms which search for the unique zero of the $(n-1)^{th}$ derivative of a function $f$ in memorizing only $n$ values of this function.

**Proof** : The proof results from the following facts :

**Fact 1** : The knowledge of the values of the function at less than $n$ points cannot reduce the length of the interval of uncertainty.

**Fact 2** : The knowledge of the values of the function at points $x_1, x_2, \dots, x_n$ reduces this interval by discarding not more than one of the two extremal intervals if $L(x_1, x_2, \dots, x_n, f) \ne 0$ and not more than the two extremal intervals if $L(x_1, x_2, \dots, x_n, f) = 0$.
To prove fact 1 we build a polynomial $Q$ such that $Q(x_i) = y_i$ for $i=1, m$ ($m<n$) and $Q^{(n-1)}(x) < 0$ for $x<\alpha$, $Q^{(n-1)}(\alpha) = 0$ and $Q^{(n-1)}(x) > 0$ for $x>\alpha$. We assume without loss of generality that $\alpha = 0$.

Let $P$ be a polynomial of degree $m-1$ such that $P(x_i) = y_i - x_i^n$ for $i = 1, m$ ($m<n$)

Let $Q(x) = x^n + P(x)$.

$Q(x)$ satisfies :

$$ egin{cases} Q(x_i) = y_i & \text{for } i=1, m. \\ Q^{(n-1)}(x) = n!x. & \end{cases} $$

Since the unique zero of $Q^{(n-1)}$ is $x=0$, the proof of fact 1 is completed.

In order to prove fact 2 we must build a function $f$ such that :

$$ - f(x_i) = y_i \quad \text{for } i=1, n $$

$- f^{(n-1)}(x)$ is continuous and changes sign only at $\alpha$ where :

$$ \alpha < x_n \quad \text{if } L(x_1, x_2, \dots, x_n, f) \ge 0 $$

$$ \alpha > x_1 \quad \text{if } L(x_1, x_2, \dots, x_n, f) \le 0. $$

Without loss of generality we assume $L(x_1, x_2, \dots, x_n, f) \ge 0$ and $\alpha = 0$. Let $i$ such that $x_1, x_2, \dots, x_i < 0$ and $x_{i+1}, \dots, x_n > 0$.

We define $f(x)$ by :

$$ egin{cases} f(x) = \lambda x^n \prod_{k \le i} (x - x_k) + P(x) & \text{for } x \ge 0 \\ f(x) = \mu x^{n'} \prod_{k > i} (x - x_k) + P(x) & \text{for } x < 0 \end{cases} $$

where $n'$ is $n$ or $n+1$ so that $n'-i+1$ is odd and $P$ is a polynomial of degree at most $n-1$. We choose $P$ such that :
$$
\begin{cases}
P(x_j) = y_i - \lambda x_j^n \prod_{k \le i} (x_j - x_k) & \text{if } i < j \le n \\
P(x_j) = y_i - \mu x_j^{n'} \prod_{k > i} (x_j - x_k) & \text{if } 0 < j \le i.
\end{cases}
$$

$\lambda$ and $\mu$ are chosen so that $P^{(n-1)} = 0$ in the following way :

$$
P^{(n-1)} = \sum_{j=1}^{i} \frac{y_j - \mu x_j^{n'} \prod_{k>i} (x_j - x_k)}{\prod_{k \neq j} (x_j - x_k)} + \sum_{j=i+1}^{n} \frac{y_j - \lambda x_j^n \prod_{k \le i} (x_j - x_k)}{\prod_{k \neq i} (x_j - x_k)} \\
= L(x_1, x_2, \dots, x_n, f) - \mu A - \lambda B
$$

$$
\text{where } A = \sum_{j=1}^{i} \frac{x_j^{n}}{\prod_{\substack{k \le i \\ k \ne j}} (x_j - x_k)} = L(x_1, x_2 \dots x_i ; x^n)
$$

$$
\text{and } B = \sum_{j=i+1}^{n} \frac{x_j^{n'}}{\prod_{\substack{k>i \\ k \ne j}} (x_j - x_k)} = L(x_{i+1}, \dots, x_n ; x^{n'}).
$$

Since $[x_i^{n'}]^{(i-1)}$ is strictly negative for $x<0$ and $[x_i^n]^{(n-i-1)}$ is strictly positive for $x>0$, we deduce by lemma 1 that $A<0$ if $i>0$ and $B>0$ if $i<n$.

If $L(x_1 x_2 \dots x_n, f) > 0$ and $i < n$ we can choose $\lambda > 0$ and $\mu > 0$ such that
$P^{(n-1)} = 0$.

If $L(x_1, x_2, \dots, f) = 0$ and $0 < i < n$ we can choose $\lambda > 0$ and $\mu > 0$ such that
$P^{(n-1)} = 0$.

Such choices of $\lambda$ and $\mu$ imply that :

$$
\begin{cases}
f^{(n-1)}(x) < 0 & \text{for } x < 0 \\
f^{(n-1)}(0) = 0 \\
f^{(n-1)}(x) > 0 & \text{for } x > 0. \quad \Box
\end{cases}
$$
## 4. Summary

We have given an algorithm which searches for the unique zero of the $(n-1)^{st}$ derivative of a function from its values. The following questions are still open :

- How may this algorithm be extended to any value of n ?

- Is this algorithm still optimal if we use more than n values of f to compute the new point of iteration ?

## 5. References

[1] Kiefer, J., "Sequential minimax search for a maximum", Proc. Amer. Soc. 4, 1953, 502-506.

[2] Kung, H.T., "Synchronized and asynchronous parallel algorithms for multiprocessors" in Proceedings of Symposium on Algorithms and Complexity : New Directions and Recent Results, J.F. Traub (Editor), Academic Press, 1976.

[3] Brent, R., "Algorithms for minimization without derivatives", Prentice Hall, 1973, pp 19-46.

## 6. Acknowledgments.

Many thanks to Gilles Kahn for his helpful comments and to J.F. Traub for many valuable suggestions on the first draft of this paper.