ITP method

In numerical analysis, the ITP method, short for Interpolate Truncate and Project, is the first root-finding algorithm that achieves the superlinear convergence of the secant method^[1] while retaining the optimal^[2] worst-case performance of the bisection method.^[3] It is also the first method with guaranteed average performance strictly better than the bisection method under any continuous distribution.^[3] In practice it performs better than traditional interpolation and hybrid based strategies (Brent's Method, Ridders, Illinois), since it not only converges super-linearly over well behaved functions but also guarantees fast performance under ill-behaved functions where interpolations fail.^[3]

The ITP method follows the same structure of standard bracketing strategies that keeps track of upper and lower bounds for the location of the root; but it also keeps track of the region where worst-case performance is kept upper-bounded. As a bracketing strategy, in each iteration the ITP queries the value of the function on one point and discards the part of the interval between two points where the function value shares the same sign. The queried point is calculated with three steps: it interpolates finding the regula falsi estimate, then it perturbes/truncates the estimate (similar to Regula falsi § Improvements in regula falsi) and then projects the perturbed estimate onto an interval in the neighbourhood of the bisection midpoint. The neighbourhood around the bisection point is calculated in each iteration in order to guarantee minmax optimality (Theorem 2.1 of ^[3]). The method depends on three hyper-parameters $\kappa _{1}\in (0,\infty ),\kappa _{2}\in \left[1,1+\phi \right)$ and $n_{0}\in [0,\infty )$ where $\phi$ is the golden ratio ${\tfrac {1}{2}}(1+{\sqrt {5}})$ , the first two control the size of the truncation and the third is a slack variable that controls the size of the interval for the projection step.

Root finding problem

Given a continuous function $f$ defined from $[a,b]$ to $\mathbb {R}$ such that $f(a)f(b)\leq 0$ , where at the cost of one query one can access the values of $f(x)$ on any given $x$ . And, given a pre-specified target precision $\epsilon >0$ , a root-finding algorithm is design to solve the following problem with the least amount of queries as possible:

Problem Definition: Find ${\hat {x}}$ such that $|{\hat {x}}-x^{*}|\leq \epsilon$ , where $x^{*}$ satisfies $f(x^{*})=0$ .

This problem is very common in numerical analysis, computer science and engineering; and, root-finding algorithms are the standard approach to solve it. Often, the root-finding procedure is called by more complex parent algorithms within a larger context, and, for this reason solving root problems efficiently is of extreme importance since an inefficient approach might come at a high computational cost when the larger context is taken into account. This is what the ITP method attempts to do by simultaneously exploiting interpolation guarantees as well as minmax optimal guarantees of the bisection method that terminates in at most $n_{1/2}\equiv \lceil (b_{0}-a_{0})/2\epsilon \rceil$ iterations when initiated on an interval $[a_{0},b_{0}]$ .

The method

Given $\kappa _{1}\in (0,\infty ),\kappa _{2}\in \left[1,1+\phi \right)$ , $n_{1/2}\equiv \lceil \log _{2}((b_{0}-a_{0})/2\epsilon )\rceil$ and $n_{0}\in [0,\infty )$ where $\phi$ is the golden ratio ${\tfrac {1}{2}}(1+{\sqrt {5}})$ , in each iteration $j=0,1,2...$ the ITP method calculates the point $x_{\text{ITP}}$ following three steps:

[Interpolation Step] Calculate the bisection and the regula falsi points: $x_{1/2}\equiv {\frac {a+b}{2}}$ and $x_{f}\equiv {\frac {bf(a)-af(b)}{f(a)-f(b)}}$ ;
[Truncation Step] Perturb the estimator towards the center: $x_{t}\equiv x_{f}+\sigma \delta$ where $\sigma \equiv {\text{sign}}(x_{1/2}-x_{f})$ and $\delta \equiv \min\{\kappa _{1}|b-a|^{\kappa _{2}},|x_{1/2}-x_{f}|\}$ ;
[Projection Step] Project the estimator to minmax interval: $x_{\text{ITP}}\equiv x_{1/2}-\sigma \rho _{k}$ where $\rho _{k}\equiv \min \left\{\epsilon 2^{n_{1/2}+n_{0}-j}-{\frac {b-a}{2}},|x_{t}-x_{1/2}|\right\}$ .

The value of the function $f(x_{\text{ITP}})$ on this point is queried, and the interval is then reduced to bracket the root by keeping the sub-interval with function values of opposite sign on each end.

The algorithm

The following algorithm (written in pseudocode) assumes the initial values of $y_{a}$ and $y_{b}$ are given and satisfy $y_{a}<0<y_{b}$ where $y_{a}\equiv f(a)$ and $y_{b}\equiv f(b)$ ; and, it returns an estimate ${\hat {x}}$ that satisfies $|{\hat {x}}-x^{*}|\leq \epsilon$ in at most $n_{1/2}+n_{0}$ function evaluations.

Input:  $a,b,\epsilon ,\kappa _{1},\kappa _{2},n_{0},f$  
    Preprocessing:  $n_{1/2}=\lceil \log _{2}{\tfrac {b-a}{2\epsilon }}\rceil$ ,  $n_{\max }=n_{1/2}+n_{0}$ , and   $j=0$ ;
    While (  $b-a>2\epsilon$  )
  
        Calculating Parameters:
             $x_{1/2}={\tfrac {a+b}{2}}$ ,  $r=\epsilon 2^{n_{\max }-j}-(b-a)/2$ ,  $\delta =\kappa _{1}(b-a)^{\kappa _{2}}$ ;
        Interpolation:
             $x_{f}={\tfrac {y_{b}a-y_{a}b}{yb-ya}}$ ;
        Truncation:
             $\sigma ={\text{sign}}(x_{1/2}-x_{f})$ ;
            If  $\delta \leq |x_{1/2}-x_{f}|$  then  $x_{t}=x_{f}+\sigma \delta$ ,
            Else  $x_{t}=x_{1/2}$ ;
        Projection:
            If  $|x_{t}-x_{1/2}|\leq r$  then  $x_{\text{ITP}}=x_{t}$ ,
            Else  $x_{\text{ITP}}=x_{1/2}-\sigma r$ ;
        Updating Interval:
             $y_{\text{ITP}}=f(x_{\text{ITP}})$ ;
            If  $y_{\text{ITP}}>0$  then  $b=x_{ITP}$  and  $y_{b}=y_{\text{ITP}}$ ,
            Elseif  $y_{\text{ITP}}<0$  then  $a=x_{\text{ITP}}$  and  $y_{a}=y_{\text{ITP}}$ ,
            Else  $a=x_{\text{ITP}}$  and  $b=x_{\text{ITP}}$ ;
             $j=j+1$ ;
Output:  ${\hat {x}}={\tfrac {a+b}{2}}$

Example: Finding the root of a polynomial

Suppose that the ITP method is used to find a root of the polynomial $f(x)=x^{3}-x-2\,.$ Using $\epsilon =0.0005,\kappa _{1}=0.1,\kappa _{2}=2$ and $n_{0}=1$ we find that:

Iteration	$a_{n}$	$b_{n}$	$c_{n}$	$f(c_{n})$
1	1	2	1.43333333333333	-0.488629629629630
2	1.43333333333333	2	1.52713145056966	0.0343383329048983
3	1.43333333333333	1.52713145056966	1.52009281150978	-0.00764147709265051
4	1.52009281150978	1.52713145056966	1.52137899116052	-4.25363464540141e-06
5	1.52137899116052	1.52713145056966	1.52138301273268	1.96497878177659e-05
6	1.52137899116052	1.52138301273268	<-- Stopping Criteria Satisfied

This example can be compared to Bisection method § Example:finding the root of a polynomial. The ITP method required less than half the number of iterations than the bisection to obtain a more precise estimate of the root with no cost on the minmax guarantees. Other methods might also attain a similar speed of convergence (such as Ridders, Brent etc.) but without the minmax guarantees given by the ITP method.

Análise

A principal vantagem do método ITP é que é garantido que não exija mais interações do que o método da bissecção quando $n_{0}=0$ . E assim, seu desempenho médio é garantido ser melhor do que o método de bissecção, mesmo quando a interpolação falha. Além disso, se as interpolações não falham (funções suaves), então é garantido desfrutar da alta ordem de convergência como métodos baseados em interpolação.

Pior caso de desempenho

O método ITP projeta o estimador no intervalo minmax com um $n_{0}$ folga, vai exigir no máximo $n_{1/2}+n_{0}$ interações (Teorema 2.1 de ^[3]). Este é minmax ideal como o método de bissecção quando $n_{0}$ é escolhido para ser $n_{0}=0$ .

Desempenho médio

Porque não é preciso mais do que $n_{1/2}+n_{0}$ interações, o número médio de interações será sempre menor do que o método da bissecção para qualquer distribuição considerada quando $n_{0}=0$ (Corolário 2.2 de ^[3]).

Desempenho assintótico

Se a função $f(x)$ é duas vezes diferenciável e a raiz $x^{*}$ é simples, então os intervalos produzidos pelo método ITP convergem para 0 com uma ordem de convergência de ${\sqrt {\kappa _{2}}}$ e se $n_{0}\neq 0$ ou se $n_{0}=0$ e $(b-a)/\epsilon$ não é uma potência de 2 com o termo ${\tfrac {\epsilon 2^{n_{1/2}}}{b-a}}$ não muito perto de zero (Teorema 2.3 de ^[3]).

References

^ Argyros, I. K.; Hernández-Verón, M. A.; Rubio, M. J. (2019). "On the Convergence of Secant-Like Methods". Current Trends in Mathematical Analysis and Its Interdisciplinary Applications: 141–183. doi:10.1007/978-3-030-15242-0_5. ISBN 978-3-030-15241-3.
^ Sikorski, K. (1982-02-01). "Bisection is optimal". Numerische Mathematik. 40 (1): 111–117. doi:10.1007/BF01459080. ISSN 0945-3245. S2CID 119952605.
^ ^a ^b ^c ^d ^e ^f ^g Oliveira, I. F. D.; Takahashi, R. H. C. (2020-12-06). "An Enhancement of the Bisection Method Average Performance Preserving Minmax Optimality". ACM Transactions on Mathematical Software. 47 (1): 5:1–5:24. doi:10.1145/3423597. ISSN 0098-3500.

External links

An Improved Bisection Method, by Kudos

[1] Argyros, I. K.; Hernández-Verón, M. A.; Rubio, M. J. (2019). "On the Convergence of Secant-Like Methods". Current Trends in Mathematical Analysis and Its Interdisciplinary Applications: 141–183. doi:10.1007/978-3-030-15242-0_5. ISBN 978-3-030-15241-3.

[2] Sikorski, K. (1982-02-01). "Bisection is optimal". Numerische Mathematik. 40 (1): 111–117. doi:10.1007/BF01459080. ISSN 0945-3245. S2CID 119952605.

[:0-3] ^ ^a ^b ^c ^d ^e ^f ^g Oliveira, I. F. D.; Takahashi, R. H. C. (2020-12-06). "An Enhancement of the Bisection Method Average Performance Preserving Minmax Optimality". ACM Transactions on Mathematical Software. 47 (1): 5:1–5:24. doi:10.1145/3423597. ISSN 0098-3500.

[1]

[2]

[3]