Finding real roots of polynomials of degree 1 000 000 Guillaume Moroz March 5, 2020 Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 1 / 23
New software: Clenshaw f ( x ) := a 0 + a 1 x + · · · + a n − 1 x n − 1 + a n x n = 0 a i ∼ Uniform ( − 100 , 100) Degree Numpy MPSolve Sage ANewDsc Clenshaw 50 4.2 e-04 3.3 e-03 3.2 e-03 5.3 e-03 1.2e-03 100 4.3 e-03 2.8 e-02 5.5 e-03 5.9 e-03 1.3e-03 500 2.7 e-01 5.9 e-01 3.3 e-01 6.4 e-02 2.3e-03 1000 1.1 2.5 9.3 e-01 3.8 e-01 3.0e-03 2000 >5 3.0 >5 3.3 3.7e-03 3000 >5 3.4 5.1e-03 4000 >5 5.2e-03 5000 6.0e-03 10000 1.1e-02 100000 1.2e-01 1000000 2.9 Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 2 / 23
Finding real roots of polynomials of degree 1 000 000 Rise and fall of a numerical approach 1 Why is it fast ? 2 Wait ... is it actually fast ?
State of the art samples Numerical factoring [Pan 02] Worst case complexity � O ( n 2 ( n + ℓ )) Not implemented Companion matrix [Moler 1991] Computes the eigen values of the companion matrix Implemented in Matlab and Numpy Aberth-Ehrlich method [Aberth 1967] Newton on the vector of roots Implemented in MPSolve [Bini, Fiorentino 2000, Bini, Robol 2014] Hybrid subdivision/Newton [Sagraloff 2012] Achieves state-of-the-art complexity Implemented in ANewDsc [Kobel, Rouillier, Sagraloff 2016] Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 4 / 23
Subdivision methods = ⇒ = ⇒ Inclusion/exclusion criteria Descartes rule [1637], available in Sage, Maple, ANewDsc Budan theorem [1807], available in Mathematica Sturm’s theorem [1829] Evaluation of f and f ′ Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 5 / 23
Subdivision methods = ⇒ = ⇒ Main algorithm Input: polynomial f and a list L of one interval I Returns: partition of I While L is not empty: I ← pop L If 0 / ∈ f ( I ) : add I to partition and continue If 0 / ∈ f ′ ( I ) : add I to partition and continue Otherwise subdivide I and add the two halves to L Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 6 / 23
Subdivision with evaluation of f and f ′ Monomial approximation over [-1, 1] � x n ≃ polynomial of degree O ( n log(1 /ε )) [Newman, Rivlin 1976] Best approximation with truncated Chebyshev series [Powell 1967] Polynomial with random coefficients ∼ Normal(0, 1) Number of real solutions : O (log n ) Main algorithm Input: polynomial f and a list L of one interval I Returns: partition of I g, h ← convert f ( x ) and f ′ ( x ) to Chebyshev basis and truncate While L is not empty: I ← pop L If 0 / ∈ g ( I ) : add I to partition and continue If 0 / ∈ h ( I ) : add I to partition and continue Otherwise subdivide I and add the two halves to L Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 7 / 23
Conversion to Chebyshev basis a n − 1 x n − 1 a n x n f ( x ) = a 0 + a 1 x + · · · + + g (cos( θ )) = c 0 + c 1 cos( θ ) + · · · + c n − 1 cos(( n − 1) θ ) + c n cos( nθ ) � �� � � �� � � �� � T 1 ( x ) T n − 1 ( x ) T n ( x ) Number of arithmetic operations Complexity O ( n log n ) [Pan 1998] Main idea f (1 2( z + 1 where z = e iθ z )) Taylor shift: f ( x + 1) O ( n log n ) arithmetic operations [Aho, Steiglitz, Ullman 1975] � O ( n ( n + ℓ )) bit operations [Gathen, Gerhard 1997] Size of the output: O ( n ( n + ℓ )) Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 8 / 23
Conversion to Chebyshev basis 1 p 0 1 p 1 2 � 2 1 � 1 p 2 4 1 4 � 3 1 � 1 p 3 8 1 8 � 4 � 4 1 � 1 � 1 p 4 16 2 16 1 16 . . . . 1 � 5 1 � 5 1 � � . . p 5 32 2 32 1 32 . . . . . . . . . � k 1 1 n 1 n 1 n 1 � � � � � � � p n · · · 2 n n/ 2 2 n n/ 2 − 1 2 n n/ 2 − 2 2 n n/ 2 − k 2 n � n +1 � n +1 � n +1 � n +1 1 � 1 � 1 � 1 � 1 p n +1 · · · 2 n +1 n/ 2 2 n +1 n/ 2 − 1 2 n +1 n/ 2 − 2 2 n +1 n/ 2 − k 2 n +1 c 0 c 1 c 2 c 3 c 4 c 5 c 6 c k c k +1 · · · Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 9 / 23
Conversion to Chebyshev basis Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 10 / 23
Conversion to Chebyshev basis 1 p 0 1 p 1 2 � 2 1 � 1 p 2 4 1 4 � 3 1 � 1 p 3 8 1 8 � 4 � 4 1 � 1 � 1 p 4 16 2 16 1 16 . . . . 1 � 5 1 � 5 1 � � . . p 5 32 2 32 1 32 . . . . . . . . . � k 1 1 n 1 n 1 n 1 � � � � � � � p n · · · 2 n n/ 2 2 n n/ 2 − 1 2 n n/ 2 − 2 2 n n/ 2 − k 2 n � n +1 � n +1 � n +1 � n +1 1 � 1 � 1 � 1 � 1 p n +1 · · · 2 n +1 n/ 2 2 n +1 n/ 2 − 1 2 n +1 n/ 2 − 2 2 n +1 n/ 2 − k 2 n +1 c 0 c 1 c 2 c 3 c 4 c 5 c 6 c k c k +1 · · · Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 11 / 23
Conversion to Chebyshev basis 1 p 0 1 p 1 2 � 2 1 � 1 p 2 4 1 4 � 3 1 � 1 p 3 8 1 8 � 4 � 4 1 � 1 � 1 p 4 16 2 16 1 16 . . . . 1 � 5 1 � 5 1 Tails of the binomial distribution . � . � p 5 32 2 32 1 32 [Hoeffding, 1963] . . . . . . ≤ e − 2 k 2 . . . n � k 1 1 n 1 n 1 n 1 � � � � � � � · · · p n 2 n n/ 2 2 n n/ 2 − 1 2 n n/ 2 − 2 2 n n/ 2 − k 2 n � n +1 � n +1 � n +1 � n +1 1 � 1 � 1 � 1 � 1 · · · p n +1 2 n +1 n/ 2 2 n +1 n/ 2 − 1 2 n +1 n/ 2 − 2 2 n +1 n/ 2 − k 2 n +1 ≤ e − k 2 n +1 · · · c 0 c 1 c 2 c 3 c 4 c 5 c 6 c k c k +1 Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 12 / 23
Evaluation in Chebyshev basis g (cos( θ )) = c 0 + c 1 cos( θ ) + · · · + c n − 1 cos(( n − 1) θ ) + c n cos( nθ ) Clenshaw algorithm Input: x = cos( θ ) Return: g ( x ) u n +1 ( x ) ← 0 u n ( x ) ← a n For k from n − 1 to 1 : u k ( x ) ← 2 xu k +1 ( x ) − u k +2 ( x ) + c k u 0 ( x ) ← xu 1 ( x ) − u 2 ( x ) + a 0 Return u 0 ( x ) Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 13 / 23
Evaluation with interval arithmetic u k ( x ) = 2 xu k +1 ( x ) − u k +2 ( x ) + c k Using interval arithmetic on [ 1 2 − ε, 1 2 + ε ] � u k +1 = m k +1 ± e k +1 u k +2 = m k +2 ± e k +2 = 2( 1 = ⇒ u k ± e k 2 ± ε )( m k +1 ± e k +1 ) − ( m k +2 ± e k +2 ) + c k = ( m k +1 − m k +2 ) ± ( εm k +1 + e k +1 + e k +2 ) = ⇒ e k ≥ e k +1 + e k +2 Error is exponential in n Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 14 / 23
Evaluation with backward error analysis [Elliott 1968] u k ( x ) = 2 xu k +1 ( x ) − u k +2 ( x ) + c k Recurrence step with x = m + ε u k ( x ) = 2( m + ε ) u k +1 ( x ) − u k +2 ( x ) + c k e k � �� � = 2 mu k +1 − u k +2 + c k − 2 εu k +2 � �� � � c k Evaluation g ( x ) = c 0 + c 1 T 1 ( x ) + · · · + c n − 1 T n − 1 ( x ) + c n T n ( x ) = � g ( m ) = � c 0 + � c 1 T 1 ( m ) + · · · + � c n − 1 T n − 1 ( m ) + � c n T n ( m ) � = ⇒ g ( x ) = � g ( m ) ⊂ g ( m ) ± | e k | Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 15 / 23
Rounding up in 2020 Triple-check the CPU Arm does not always support rounding up for floating points operations Triple-check your compiler Gcc, Clang don’t guarantee correctness for change of rounding mode https://gcc.gnu.org/wiki/FloatingPointMath, http://lists.llvm.org/pipermail/llvm-dev/2018-May/123529.html Triple-check your parallel computation library OpenMP doesn’t copy the FPU control word from the main thread Triple check your code Misplaced rounding modes are hard to detect Triple-check Wikipedia Wrong bound on tail of binomial distribution. Correction rejected. Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 16 / 23
Implementation Optimizations: C code with Python interface Order operations for cache performance Amortized latency of arithmetic operations Vectorized instructions Multi-threaded Profiling case n = 1 000 000 % Time Line Contents [...] 49.1 chebyshev(Vf, f, err_f) [...] 2.5 solve(f, df, err_f, err_df) [...] 46.1 chebyshev(Vf, f, err_f) [...] 1.6 solve(f, df, err_f, err_df) Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 17 / 23
The curious case of the resultant f ( x ) = Res y ( P ( x, y ) , Q ( x, y )) Degree Numpy MPSolve Sage ANewDsc Clenshaw 5 2 3.1e-04 3.4e-03 4.5e-03 1.0e-02 5.2e-03 6 2 5.2e-04 5.1e-03 6.3e-03 8.5e-03 1.0e-02 7 2 9.8e-04 7.8e-03 7.8e-03 1.0e-02 1.0e-02 8 2 1.7e-03 1.1e-02 2.7e-01 1.3e-02 7.7e-02 9 2 2.5e-03 1.7e-02 1.4e-02 2.5e-02 2.5e+00 10 2 9.5e-03 2.5e-02 8.7e-03 4.2e-02 7.4e-02 11 2 1.4e-02 3.4e-02 2.1e-01 4.0e-02 > 5 12 2 2.3e-02 4.7e-02 3.7e-01 4.6e-02 13 2 4.1e-02 6.2e-02 2.4e-01 8.6e-02 14 2 5.5e-02 8.3e-02 2.7e-01 6.4e-02 15 2 7.8e-02 1.0e-01 5.0e-01 7.7e-02 Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 18 / 23
Recommend
More recommend
