ITP method

This article has multiple issues. Please help improve it or discuss these issues on the talk page. (Learn how and when to remove these messages) This article relies excessively on references to primary sources. Please improve this article by adding secondary or tertiary sources. Find sources: "ITP method" – news · newspapers · books · scholar · JSTOR (November 2024) (Learn how and when to remove this message) This article needs additional citations for verification. Please help improve this article by adding citations to reliable sources. Unsourced material may be challenged and removed. Find sources: "ITP method" – news · newspapers · books · scholar · JSTOR (January 2021) (Learn how and when to remove this message) The topic of this article may not meet Wikipedia's general notability guideline. Please help to demonstrate the notability of the topic by citing reliable secondary sources that are independent of the topic and provide significant coverage of it beyond a mere trivial mention. If notability cannot be shown, the article is likely to be merged, redirected, or deleted. Find sources: "ITP method" – news · newspapers · books · scholar · JSTOR (May 2024) (Learn how and when to remove this message) (Learn how and when to remove this message)

In numerical analysis, the ITP method (Interpolate Truncate and Project method) 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 perturbs/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 κ 1 ∈ ( 0 , ∞ ) , κ 2 ∈ [ 1 , 1 + ϕ ) {\displaystyle \kappa _{1}\in (0,\infty ),\kappa _{2}\in \left[1,1+\phi \right)} and n 0 ∈ [ 0 , ∞ ) {\displaystyle n_{0}\in [0,\infty )} where ϕ {\displaystyle \phi } is the golden ratio 1 2 ( 1 + 5 ) {\displaystyle {\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.^[a]

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

Problem Definition: Find x ^ {\displaystyle {\hat {x}}} such that | x ^ − x ∗ | ≤ ϵ {\displaystyle |{\hat {x}}-x^{*}|\leq \epsilon } , where x ∗ {\displaystyle x^{*}} satisfies f ( x ∗ ) = 0 {\displaystyle 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 ≡ ⌈ log 2 ⁡ ( ( b 0 − a 0 ) / 2 ϵ ) ⌉ {\displaystyle n_{1/2}\equiv \lceil \log _{2}((b_{0}-a_{0})/2\epsilon )\rceil } iterations when initiated on an interval [ a 0 , b 0 ] {\displaystyle [a_{0},b_{0}]} .

Given κ 1 ∈ ( 0 , ∞ ) , κ 2 ∈ [ 1 , 1 + ϕ ) {\displaystyle \kappa _{1}\in (0,\infty ),\kappa _{2}\in \left[1,1+\phi \right)} , n 1 / 2 ≡ ⌈ log 2 ⁡ ( ( b 0 − a 0 ) / 2 ϵ ) ⌉ {\displaystyle n_{1/2}\equiv \lceil \log _{2}((b_{0}-a_{0})/2\epsilon )\rceil } and n 0 ∈ [ 0 , ∞ ) {\displaystyle n_{0}\in [0,\infty )} where ϕ {\displaystyle \phi } is the golden ratio 1 2 ( 1 + 5 ) {\displaystyle {\tfrac {1}{2}}(1+{\sqrt {5}})} , in each iteration j = 0 , 1 , 2 … {\displaystyle j=0,1,2\dots } the ITP method calculates the point x ITP {\displaystyle x_{\text{ITP}}} following three steps:

1. [Interpolation Step] Calculate the bisection and the regula falsi points: x 1 / 2 ≡ a + b 2 {\displaystyle x_{1/2}\equiv {\frac {a+b}{2}}} and x f ≡ b f ( a ) − a f ( b ) f ( a ) − f ( b ) {\displaystyle x_{f}\equiv {\frac {bf(a)-af(b)}{f(a)-f(b)}}}  ;
2. [Truncation Step] Perturb the estimator towards the center: x t ≡ x f + σ δ {\displaystyle x_{t}\equiv x_{f}+\sigma \delta } where σ ≡ sign ( x 1 / 2 − x f ) {\displaystyle \sigma \equiv {\text{sign}}(x_{1/2}-x_{f})} and δ ≡ min { κ 1 | b − a | κ 2 , | x 1 / 2 − x f | } {\displaystyle \delta \equiv \min\{\kappa _{1}|b-a|^{\kappa _{2}},|x_{1/2}-x_{f}|\}}  ;
3. [Projection Step] Project the estimator to minmax interval: x ITP ≡ x 1 / 2 − σ ρ k {\displaystyle x_{\text{ITP}}\equiv x_{1/2}-\sigma \rho _{k}} where ρ k ≡ min { ϵ 2 n 1 / 2 + n 0 − j − b − a 2 , | x t − x 1 / 2 | } {\displaystyle \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 ITP ) {\displaystyle 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 following algorithm (written in pseudocode) assumes the initial values of y a {\displaystyle y_{a}} and y b {\displaystyle y_{b}} are given and satisfy y a < 0 < y b {\displaystyle y_{a}<0<y_{b}} where y a ≡ f ( a ) {\displaystyle y_{a}\equiv f(a)} and y b ≡ f ( b ) {\displaystyle y_{b}\equiv f(b)} ; and, it returns an estimate x ^ {\displaystyle {\hat {x}}} that satisfies | x ^ − x ∗ | ≤ ϵ {\displaystyle |{\hat {x}}-x^{*}|\leq \epsilon } in at most n 1 / 2 + n 0 {\displaystyle n_{1/2}+n_{0}} function evaluations.

Input: 
  
    
      
        a
        ,
        b
        ,
        ϵ
        ,
        
          κ
          
            1
          
        
        ,
        
          κ
          
            2
          
        
        ,
        
          n
          
            0
          
        
        ,
        f
      
    
    {\displaystyle a,b,\epsilon ,\kappa _{1},\kappa _{2},n_{0},f}
  
 
    Preprocessing: 
  
    
      
        
          n
          
            1
            
              /
            
            2
          
        
        =
        ⌈
        
          log
          
            2
          
        
        ⁡
        
          
            
              
                b
                −
                a
              
              
                2
                ϵ
              
            
          
        
        ⌉
      
    
    {\displaystyle n_{1/2}=\lceil \log _{2}{\tfrac {b-a}{2\epsilon }}\rceil }
  
, 
  
    
      
        
          n
          
            max
          
        
        =
        
          n
          
            1
            
              /
            
            2
          
        
        +
        
          n
          
            0
          
        
      
    
    {\displaystyle n_{\max }=n_{1/2}+n_{0}}
  
, and  
  
    
      
        j
        =
        0
      
    
    {\displaystyle j=0}
  
;
    While ( 
  
    
      
        b
        −
        a
        >
        2
        ϵ
      
    
    {\displaystyle b-a>2\epsilon }
  
 )
  
        Calculating Parameters:
            
  
    
      
        
          x
          
            1
            
              /
            
            2
          
        
        =
        
          
            
              
                a
                +
                b
              
              2
            
          
        
      
    
    {\displaystyle x_{1/2}={\tfrac {a+b}{2}}}
  
, 
  
    
      
        r
        =
        ϵ
        
          2
          
            
              n
              
                max
              
            
            −
            j
          
        
        −
        (
        b
        −
        a
        )
        
          /
        
        2
      
    
    {\displaystyle r=\epsilon 2^{n_{\max }-j}-(b-a)/2}
  
, 
  
    
      
        δ
        =
        
          κ
          
            1
          
        
        (
        b
        −
        a
        
          )
          
            
              κ
              
                2
              
            
          
        
      
    
    {\displaystyle \delta =\kappa _{1}(b-a)^{\kappa _{2}}}
  
;
        Interpolation:
            
  
    
      
        
          x
          
            f
          
        
        =
        
          
            
              
                
                  y
                  
                    b
                  
                
                a
                −
                
                  y
                  
                    a
                  
                
                b
              
              
                
                  y
                  
                    b
                  
                
                −
                
                  y
                  
                    a
                  
                
              
            
          
        
      
    
    {\displaystyle x_{f}={\tfrac {y_{b}a-y_{a}b}{y_{b}-y_{a}}}}
  
;
        Truncation:
            
  
    
      
        σ
        =
        
          sign
        
        (
        
          x
          
            1
            
              /
            
            2
          
        
        −
        
          x
          
            f
          
        
        )
      
    
    {\displaystyle \sigma ={\text{sign}}(x_{1/2}-x_{f})}
  
;
            If 
  
    
      
        δ
        ≤
        
          |
        
        
          x
          
            1
            
              /
            
            2
          
        
        −
        
          x
          
            f
          
        
        
          |
        
      
    
    {\displaystyle \delta \leq |x_{1/2}-x_{f}|}
  
 then 
  
    
      
        
          x
          
            t
          
        
        =
        
          x
          
            f
          
        
        +
        σ
        δ
      
    
    {\displaystyle x_{t}=x_{f}+\sigma \delta }
  
,
            Else 
  
    
      
        
          x
          
            t
          
        
        =
        
          x
          
            1
            
              /
            
            2
          
        
      
    
    {\displaystyle x_{t}=x_{1/2}}
  
;
        Projection:
            If 
  
    
      
        
          |
        
        
          x
          
            t
          
        
        −
        
          x
          
            1
            
              /
            
            2
          
        
        
          |
        
        ≤
        r
      
    
    {\displaystyle |x_{t}-x_{1/2}|\leq r}
  
 then 
  
    
      
        
          x
          
            ITP
          
        
        =
        
          x
          
            t
          
        
      
    
    {\displaystyle x_{\text{ITP}}=x_{t}}
  
,
            Else 
  
    
      
        
          x
          
            ITP
          
        
        =
        
          x
          
            1
            
              /
            
            2
          
        
        −
        σ
        r
      
    
    {\displaystyle x_{\text{ITP}}=x_{1/2}-\sigma r}
  
;
        Updating Interval:
            
  
    
      
        
          y
          
            ITP
          
        
        =
        f
        (
        
          x
          
            ITP
          
        
        )
      
    
    {\displaystyle y_{\text{ITP}}=f(x_{\text{ITP}})}
  
;
            If 
  
    
      
        
          y
          
            ITP
          
        
        >
        0
      
    
    {\displaystyle y_{\text{ITP}}>0}
  
 then 
  
    
      
        b
        =
        
          x
          
            I
            T
            P
          
        
      
    
    {\displaystyle b=x_{ITP}}
  
 and 
  
    
      
        
          y
          
            b
          
        
        =
        
          y
          
            ITP
          
        
      
    
    {\displaystyle y_{b}=y_{\text{ITP}}}
  
,
            Elseif 
  
    
      
        
          y
          
            ITP
          
        
        <
        0
      
    
    {\displaystyle y_{\text{ITP}}<0}
  
 then 
  
    
      
        a
        =
        
          x
          
            ITP
          
        
      
    
    {\displaystyle a=x_{\text{ITP}}}
  
 and 
  
    
      
        
          y
          
            a
          
        
        =
        
          y
          
            ITP
          
        
      
    
    {\displaystyle y_{a}=y_{\text{ITP}}}
  
,
            Else 
  
    
      
        a
        =
        
          x
          
            ITP
          
        
      
    
    {\displaystyle a=x_{\text{ITP}}}
  
 and 
  
    
      
        b
        =
        
          x
          
            ITP
          
        
      
    
    {\displaystyle b=x_{\text{ITP}}}
  
;
            
  
    
      
        j
        =
        j
        +
        1
      
    
    {\displaystyle j=j+1}
  
;
Output: 
  
    
      
        
          
            
              x
              ^
            
          
        
        =
        
          
            
              
                a
                +
                b
              
              2
            
          
        
      
    
    {\displaystyle {\hat {x}}={\tfrac {a+b}{2}}}
  

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

Iteration	a n {\displaystyle a_{n}}	b n {\displaystyle b_{n}}	c n {\displaystyle c_{n}}	f ( c n ) {\displaystyle 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.

The main advantage of the ITP method is that it is guaranteed to require no more iterations than the bisection method when n 0 = 0 {\displaystyle n_{0}=0} . And so its average performance is guaranteed to be better than the bisection method even when interpolation fails. Furthermore, if interpolations do not fail (smooth functions), then it is guaranteed to enjoy the high order of convergence as interpolation based methods.

Because the ITP method projects the estimator onto the minmax interval with a n 0 {\displaystyle n_{0}} slack, it will require at most n 1 / 2 + n 0 {\displaystyle n_{1/2}+n_{0}} iterations (Theorem 2.1 of ^[3]). This is minmax optimal like the bisection method when n 0 {\displaystyle n_{0}} is chosen to be n 0 = 0 {\displaystyle n_{0}=0} .

Because it does not take more than n 1 / 2 + n 0 {\displaystyle n_{1/2}+n_{0}} iterations, the average number of iterations will always be less than that of the bisection method for any distribution considered when n 0 = 0 {\displaystyle n_{0}=0} (Corollary 2.2 of ^[3]).

If the function f ( x ) {\displaystyle f(x)} is twice differentiable and the root x ∗ {\displaystyle x^{*}} is simple, then the intervals produced by the ITP method converges to 0 with an order of convergence of κ 2 {\displaystyle {\sqrt {\kappa _{2}}}} if n 0 ≠ 0 {\displaystyle n_{0}\neq 0} or if n 0 = 0 {\displaystyle n_{0}=0} and ( b − a ) / ϵ {\displaystyle (b-a)/\epsilon } is not a power of 2 with the term ϵ 2 n 1 / 2 b − a {\displaystyle {\tfrac {\epsilon 2^{n_{1/2}}}{b-a}}} not too close to zero (Theorem 2.3 of ^[3]).

• The itp^[4] contributed package in R.

• Bisection method
• Ridders' method
• Regula falsi
• Brent's method

• An Improved Bisection Method, by Kudos