Finding real roots of polynomials of degree 1 000 000 Guillaume - - PowerPoint PPT Presentation

▶
finding real roots of polynomials of degree 1 000 000
SMART_READER_LITE
LIVE PREVIEW

Finding real roots of polynomials of degree 1 000 000 Guillaume - - PowerPoint PPT Presentation

May 26, 2023 155 likes •392 views

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

slide-1
SLIDE 1

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

slide-2
SLIDE 2

New software: Clenshaw

f(x) := a0 + a1x + · · · + an−1xn−1 + anxn = 0 ai ∼ 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

slide-3
SLIDE 3

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 ?

slide-4
SLIDE 4

State of the art samples

Numerical factoring [Pan 02] Worst case complexity O(n2(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

slide-5
SLIDE 5

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

slide-6
SLIDE 6

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

slide-7
SLIDE 7

Subdivision with evaluation of f and f ′

Monomial approximation over [-1, 1] xn ≃ 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

slide-8
SLIDE 8

Conversion to Chebyshev basis

f(x) = a0 + a1x + · · · + an−1xn−1 + anxn g(cos(θ)) = c0 + c1 cos(θ)

T1(x)

+ · · · + cn−1 cos((n − 1)θ)

  • Tn−1(x)

+ cn cos(nθ)

  • Tn(x)

Number of arithmetic operations Complexity O(n log n) [Pan 1998] Main idea f(1 2(z + 1 z )) where z = eiθ 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

slide-9
SLIDE 9

Conversion to Chebyshev basis

p0 1 p1

1 2

p2

1 4

2

1

  • 1

4

p3

1 8

3

1

  • 1

8

p4

1 16

4

2

  • 1

16

4

1

  • 1

16

p5 . . .

1 32

5

2

  • .

. .

1 32

5

1

  • 1

32

. . . . . . . . . pn

1 2n

k

n/2

  • 1

2n

  • n

n/2−1

  • 1

2n

  • n

n/2−2

  • 1

2n

  • n

n/2−k

  • · · ·

1 2n

pn+1

1 2n+1

n+1

n/2

  • 1

2n+1

n+1

n/2−1

  • 1

2n+1

n+1

n/2−2

  • 1

2n+1

n+1

n/2−k

  • · · ·

1 2n+1

c0 c1 c2 c3 c4 c5 c6 · · · ck ck+1

Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 9 / 23

slide-10
SLIDE 10

Conversion to Chebyshev basis

Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 10 / 23

slide-11
SLIDE 11

Conversion to Chebyshev basis

p0 1 p1

1 2

p2

1 4

2

1

  • 1

4

p3

1 8

3

1

  • 1

8

p4

1 16

4

2

  • 1

16

4

1

  • 1

16

p5 . . .

1 32

5

2

  • .

. .

1 32

5

1

  • 1

32

. . . . . . . . . pn

1 2n

k

n/2

  • 1

2n

  • n

n/2−1

  • 1

2n

  • n

n/2−2

  • 1

2n

  • n

n/2−k

  • · · ·

1 2n

pn+1

1 2n+1

n+1

n/2

  • 1

2n+1

n+1

n/2−1

  • 1

2n+1

n+1

n/2−2

  • 1

2n+1

n+1

n/2−k

  • · · ·

1 2n+1

c0 c1 c2 c3 c4 c5 c6 · · · ck ck+1

Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 11 / 23

slide-12
SLIDE 12

Conversion to Chebyshev basis

p0 1 p1

1 2

p2

1 4

2

1

  • 1

4

p3

1 8

3

1

  • 1

8

p4

1 16

4

2

  • 1

16

4

1

  • 1

16

p5 . . .

1 32

5

2

  • .

. .

1 32

5

1

  • 1

32

. . . . . . . . . pn

1 2n

k

n/2

  • 1

2n

  • n

n/2−1

  • 1

2n

  • n

n/2−2

  • 1

2n

  • n

n/2−k

  • · · ·

1 2n

pn+1

1 2n+1

n+1

n/2

  • 1

2n+1

n+1

n/2−1

  • 1

2n+1

n+1

n/2−2

  • 1

2n+1

n+1

n/2−k

  • · · ·

1 2n+1

c0 c1 c2 c3 c4 c5 c6 · · · ck ck+1 Tails of the binomial distribution

≤ e− k2

n+1

≤ e−2 k2

n

[Hoeffding, 1963]

Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 12 / 23

slide-13
SLIDE 13

Evaluation in Chebyshev basis

g(cos(θ)) = c0 + c1 cos(θ) + · · · + cn−1 cos((n − 1)θ) + cn cos(nθ)

Clenshaw algorithm

Input: x = cos(θ) Return: g(x) un+1(x) ← 0 un(x) ← an For k from n − 1 to 1:

uk(x) ← 2xuk+1(x) − uk+2(x) + ck

u0(x) ← xu1(x) − u2(x) + a0 Return u0(x)

Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 13 / 23

slide-14
SLIDE 14

Evaluation with interval arithmetic

uk(x) = 2xuk+1(x) − uk+2(x) + ck Using interval arithmetic on [ 1

2 − ε, 1 2 + ε]

  • uk+1 = mk+1 ± ek+1

uk+2 = mk+2 ± ek+2 = ⇒ uk ± ek = 2(1

2 ± ε)(mk+1 ± ek+1) − (mk+2 ± ek+2) + ck

= (mk+1 − mk+2) ± (εmk+1 + ek+1 + ek+2) = ⇒ ek ≥ ek+1 + ek+2 Error is exponential in n

Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 14 / 23

slide-15
SLIDE 15

Evaluation with backward error analysis [Elliott 1968]

uk(x) = 2xuk+1(x) − uk+2(x) + ck Recurrence step with x = m + ε uk(x) = 2(m + ε)uk+1(x) − uk+2(x) + ck = 2muk+1 − uk+2 + ck −

ek

2εuk+2

  • ck

Evaluation g(x) = c0 + c1T1(x) + · · · + cn−1Tn−1(x) + cnTn(x) =

  • g(m) =

c0 + c1T1(m) + · · · + cn−1Tn−1(m) + cnTn(m) = ⇒ g(x) = g(m) ⊂ g(m) ±

  • |ek|

Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 15 / 23

slide-16
SLIDE 16

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

slide-17
SLIDE 17

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

slide-18
SLIDE 18

The curious case of the resultant

f(x) = Resy(P(x, y), Q(x, y)) Degree Numpy MPSolve Sage ANewDsc Clenshaw 52 3.1e-04 3.4e-03 4.5e-03 1.0e-02 5.2e-03 62 5.2e-04 5.1e-03 6.3e-03 8.5e-03 1.0e-02 72 9.8e-04 7.8e-03 7.8e-03 1.0e-02 1.0e-02 82 1.7e-03 1.1e-02 2.7e-01 1.3e-02 7.7e-02 92 2.5e-03 1.7e-02 1.4e-02 2.5e-02 2.5e+00 102 9.5e-03 2.5e-02 8.7e-03 4.2e-02 7.4e-02 112 1.4e-02 3.4e-02 2.1e-01 4.0e-02 > 5 122 2.3e-02 4.7e-02 3.7e-01 4.6e-02 132 4.1e-02 6.2e-02 2.4e-01 8.6e-02 142 5.5e-02 8.3e-02 2.7e-01 6.4e-02 152 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

slide-19
SLIDE 19

The curious case of the resultant

f(x) = Resy(P(x, y), Q(x, y)) % Time Line Contents ===================== 0.7 chebyshev(Vf, f, err_f) [...] 88.3 solve(f, df, err_f, err_df) [...] 0.1 chebyshev(Vf, f, err_f) [...] 6.5 solve(f, df, err_f, err_df)

Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 19 / 23

slide-20
SLIDE 20

Condition number

Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 20 / 23

slide-21
SLIDE 21

Conclusion

Code Prototype with machine precision accessible with GPL licence http://gitlab.inria.fr/gmoro/clenshaw Open questions Complexity for random case in O(n√n) ? Approximated conversion to Chebyshev basis in O(n) ? Generalization to multivariate ?

Guillaume Moroz Finding real roots of polynomials of degree 1 000 000 March 5, 2020 21 / 23

slide-22
SLIDE 22

Thank you

slide-23
SLIDE 23