The Approximate Irreducible Factorization - ISSAC 2009 - KIAS
The
Approximate Irreducible Factorization of a Univariate Polynomial. Revisited Zhonggang Zeng Northeastern Illinois University
ISSAC 2009, KIAS, Seoul, July 31, 2009 (supported in part by NSF under Grant DMS-0715137)
This talk revisits the work
An algorithm for accurate computation of multiple roots and multiplicities using floating point arithmetic even if the coefficients are perturbed. 1
A two-staged algorithm proposed in ISSAC 2003 For the polynomial
2.5
2
1.5
15
10
5
with (inexact ) coefficients in machine precision
1
imaginary part
( x − 1) ( x − 2) ( x − 3) ( x − 4) 20
0.5
0
−0.5
−1
−1.5
−2
−2.5
Stage I results: The backward error: Computed roots 1.000000000000353 2.000000000030904 3.000000000176196 4.000000000109542
0
1
2
3 real part
4
5
6
Stage II results: 6.05 x 10-10
The backward error:
multiplicities 20 15 10 5
Software package MultRoot:
Computed roots 1.000000000000000 1.999999999999997 3.000000000000011 3.999999999999985
Z. Zeng, ACM TOMS 2004
6.16 x 10-16 multiplicities 20 15 10 5
2
Example of a new method: For polynomial ( x − 1)80 ( x − 2)60 ( x − 3) 40 ( x − 4) 20 with (inexact ) coefficients in hardware precision
3
Example: For polynomial ( x − 1)80 ( x − 2)60 ( x − 3) 40 ( x − 4)20 with (inexact ) coefficients in hardware precision Exact factorization
Approximate factorization: >> [F,res,fcnd]
=
uvFactor(f,1e-10,1);
THE CONDITION NUMBER: THE BACKWARD ERROR: THE ESTIMATED FORWARD ROOT ERROR:
914.329 5.71e-015 1.04e-011
FACTORS ( ( ( (
x x x x
-
3.999999999999990 3.000000000000008 1.999999999999998 1.000000000000000
)^20 )^40 )^60 )^80
A significant advancement in robustness but not really the point for a revisit
Question: What problem are we really solving? 4
The Approximate Irreducible Factorization (also known as Root-finding) Problem:
Given
p(x) = x + .6667 x – 2.333 x - 1.333 x 6
S
tu r e ll p a m
(
4
3
+ 1.667 x2 + 0.6667 x – 0.3333
io n t a rb
pˆ (x) = (x −1)2 (x +1)3 x − 1 3
5
Conventional Factorization
)
1
(x+1.0189+0.0034i)(x+1.019-0.00336i) (x+0.9621)(x-0.3332) (x-1+0.0144i) (x-1-0.0144i)
= aˆ0 (aˆ1x + bˆ1)2 (aˆ2 x + bˆ2 )3 (aˆ3 x + bˆ3 )1
≈ ≈ ≈ ≈
1. Match multiplicities
AIF
≈ ≈
~ 2. aˆi − a~i + bˆi − bi = O( p − ~ p)
well-posed?
≈
≈
~ 2 ~ ~ 3 ~ ~ 1 ~ ~ ~ p ( x) = a0 (a1 x + b1 ) (a2 x + b2 ) (a3 x + b3 ) = 0.9999(1.0 x − 1.00001) 2 (1.0 x + 1.000009) 3 (1.0 x − 0.3334)1
5
A well-posed problem: (Hadamard, 1923) the solution satisfies
•
• •
existence uniqueness continuity w.r.t data
An ill-posed problem is infinitely sensitive to perturbation tiny perturbation
Î
huge error
A frontier in scientific computing Though frequently needed in application, the adequate handling of such ill-posed … problems is hardly ever touched upon in numerical analysis textbooks. --- Arnold Neumaier, SIAM Review
6
Challenge in solving ill-posed problems: Can we recover the lost solution when the problem is inexact?
P
Solution
P
Data
P : Data Æ Solution
P 7
Geometry of AIF (simplified view) (x – t)3
Π(3)
= –t3 + ( 3 t2) x + (–3 t) x2 + x3
⎡− t 3 ⎤ ⎢ ⎥ F (t ) = ⎢ 3t 2 ⎥ ⎢− 3t ⎥ ⎣ ⎦
F
F-1 diffeomorphism
t
Π (1,2) (x – u)1 (x – v)2
= –uv2 + (v2+ 2uv) x + (–2v–u) x2 + x3
⎡ − uv 2 ⎤ ⎢ 2 ⎥ G (u, v ) = ⎢v + 2uv ⎥ ⎢ − 2v − u ⎥ ⎣ ⎦
v
G
u
G-1 diffeomorphism
Polynomials form (factorization) manifolds
8
Proposition 1:
Polynomials
c0 + c1 x + c2 x 2 + L + cm x m = a0 (a1 x + bˆ1 ) k1 (a2 x + b2 ) k2 L (an x + bn ) kn
form a differentiable manifold Π[k
1 ,…, kn]
of codimension
codim (Π[k1 ,…, kn] ) = m - n
Number of factors dimension of the polynomial vector space (≥ polynomial degree) 9
Are ill-posed problems really sensitive? Kahan: It is a misconception. W. Kahan’s observation (1972) • Problems form a “pejorative manifolds”
Plot of pejorative manifolds of degree 3 polynomials with multiple roots
• Ill-posedness: a tiny perturbation pushes the problem out of the manifold
• A problem is not sensitive at all if it stays on the manifold.
10
Stratification of factorization manifolds of degree 3 monic polynomials
Π (1,1,1) = {p( x ) = ( x − α ) ( x − β ) ( x − γ ) α ≠ β ≠ γ } 1
1
1
Π (1,2) = {p( x ) = ( x − α )1 ( x − β )2 α ≠ β }
Π (3) = {p( x ) = ( x − α )3 α ∈ C }
Π (3) ⊂ Π (1,2) ⊂ Π (1,1,1) = C 3 Codimensions: 2
1
0
Factorization manifold stratification of degree 4 polynomials:
Π ( 2,2) Π (1,1,2)
Π ( 4)
Π (1,1,1,1)
= C4
Π (1,3) Codimensions: 3
2
1
0
11
Factorization manifolds and their stratification
{
Π[ k1k 2 Lk n ] = a0 ( a1 x + b1 ) k1 ( a2 x + b2 )k 2 L( an x + bn ) k n
ai , bi ∈ C , ai b j ≠ a j bi , ∀i ≠ j
}
⊂ Cm [ x ] = { c0 + c1 x + L + cm x m ci ∈ C }
p ∈ ∏[2,2]
Proposition 3:
⇄
dist(p, ∏[2,2] ) = dist(p, ∏[2,1,1] ) = dist(p, ∏[1,1,1,1] ) = 0
p ∈ ∏[k1… kn]
if and only if
codim(∏[k1… kn] ) = max { codim(∏) | dist(p,∏ ) = 0 } 12
Formulation of the approximate irreducible factorization
Π
p The approximate factorization of - the exact factorization of -
p
pˆ
is
~ p
~ p
~ p
lies in the nearby manifold ∏ of the highest codimension
~ p from
is the nearest polynomial on ∏
p
31
A “three-strikes” principle for formulating an approximate irreducible factorization: • Backward nearness: The AIF is the exact factorization of a nearby polynomial • Maximum codimension: The AIF is the exact factorization of a polynomial in the nearby factorization manifold of the highest codimension.
• Minimum distance: The AIF is the exact factorization of the nearest polynomial in the nearby factorization manifold of the highest codimension.
In comparison: Symbolic computation: Backward nearness with distance = 0
Numerical computation: (straightforward) Backward nearness with minimal distance
Finding the AIF is (apparently) a well-posed problem The AIF is a generalization of exact factorization.
13
pˆ ( x) = aˆ0 (aˆ1 x + bˆ1 ) k1 L (aˆ n x + bˆn ) k n ∈ Π k
Theorem 1 Assume Then,
p − pˆ
is small
∃ an interval I and ∀ε∈ I
∃ a unique AIF within ε ~ ~ ~ p ( x) = a~0 (a~1 x + b1 ) k1 L (a~n x + bn ) k n ∈ Π k
Πκ
p pˆ
such that
~ (a~i x + bi ) − (aˆi x + bˆi ) = O ( p − pˆ
)
~ p
Moreover, the AIF is continuous w.r.t. p
17
The two-staged algorithm after formulating the AIF of p within ε
Stage I: Find the factorization manifold Π of the highest dimension s.t.
dist ( P, Π k ) < ε
Πκ
p pˆ ~ p
Stage II: Find/solve problem Q such that
p− ~ p = min p − q q∈Π k
18
Stage I: Example:
Identify the AIF manifold by a squarefree factorization
p15 p23 p33 p4
= ( p1 p2 p3 p4 )( p1 p2 p3 )( p1 p2 p3 )( p1 )( p1 )
= ( p 4 ) (1 ) ( p 2 p 3 ) (1 ) ( p1 ) 1
2
3
4
5
--- flat SFF
--- staircase SFF
A new staircase SFF algorithm:
p ≈ u11 u22 L pkk
a0 (a1 x + b1 ) k1 (a2 x + b2 ) k 2 L (an x + bn ) k n
19
Stage II:
Minimize the distance to the AIF manifold
a0 (a1 x + b1 ) k1 (a2 x + b2 ) k 2 L (an x + bn ) k n = p α1a1 + β1b1 = γ1
α 2 a2 + β 2b2
= γ2 O
M M
α n an + β nbn = γ n
G(a0, a1 ,…,an, b1 ,…,bn) = d Nonlinear least squares problem solved by the Gauss-Newton iteration
20
Example: For polynomial ( x − 1)80 ( x − 2)60 ( x − 3) 40 ( x − 4)20 with (inexact ) coefficients in hardware precision
Exact factorization
Approximate factorization: >> [F,res,fcnd]
=
uvFactor(f,1e-10,1);
THE CONDITION NUMBER: THE BACKWARD ERROR: THE ESTIMATED FORWARD ROOT ERROR:
914.329 5.71e-015 1.04e-011
FACTORS ( ( ( (
x x x x
-
3.999999999999990 3.000000000000008 1.999999999999998 1.000000000000000
)^20 )^40 )^60 )^80 21
Summary: • Factorizations are sensitive because polynomials form manifolds of positive codimensions in strata. • An AIF can be formulated as an exact factorization of the nearest polynomial on a nearby manifold of the highest codimension. • The AIF approximates the factorization of the (hidden) underlying polynomial from the perturbed data. • The AIF can be computed by an (improved) algorithm in two stages.
Software is available in the package ApaTools (google apatools)
Approximate Irreducible Factorization of a Univariate Polynomial. Revisited Zhonggang Zeng Northeastern Illinois University
ISSAC 2009, KIAS, Seoul, July 31, 2009 (supported in part by NSF under Grant DMS-0715137)
This talk revisits the work
An algorithm for accurate computation of multiple roots and multiplicities using floating point arithmetic even if the coefficients are perturbed. 1
A two-staged algorithm proposed in ISSAC 2003 For the polynomial
2.5
2
1.5
15
10
5
with (inexact ) coefficients in machine precision
1
imaginary part
( x − 1) ( x − 2) ( x − 3) ( x − 4) 20
0.5
0
−0.5
−1
−1.5
−2
−2.5
Stage I results: The backward error: Computed roots 1.000000000000353 2.000000000030904 3.000000000176196 4.000000000109542
0
1
2
3 real part
4
5
6
Stage II results: 6.05 x 10-10
The backward error:
multiplicities 20 15 10 5
Software package MultRoot:
Computed roots 1.000000000000000 1.999999999999997 3.000000000000011 3.999999999999985
Z. Zeng, ACM TOMS 2004
6.16 x 10-16 multiplicities 20 15 10 5
2
Example of a new method: For polynomial ( x − 1)80 ( x − 2)60 ( x − 3) 40 ( x − 4) 20 with (inexact ) coefficients in hardware precision
3
Example: For polynomial ( x − 1)80 ( x − 2)60 ( x − 3) 40 ( x − 4)20 with (inexact ) coefficients in hardware precision Exact factorization
Approximate factorization: >> [F,res,fcnd]
=
uvFactor(f,1e-10,1);
THE CONDITION NUMBER: THE BACKWARD ERROR: THE ESTIMATED FORWARD ROOT ERROR:
914.329 5.71e-015 1.04e-011
FACTORS ( ( ( (
x x x x
-
3.999999999999990 3.000000000000008 1.999999999999998 1.000000000000000
)^20 )^40 )^60 )^80
A significant advancement in robustness but not really the point for a revisit
Question: What problem are we really solving? 4
The Approximate Irreducible Factorization (also known as Root-finding) Problem:
Given
p(x) = x + .6667 x – 2.333 x - 1.333 x 6
S
tu r e ll p a m
(
4
3
+ 1.667 x2 + 0.6667 x – 0.3333
io n t a rb
pˆ (x) = (x −1)2 (x +1)3 x − 1 3
5
Conventional Factorization
)
1
(x+1.0189+0.0034i)(x+1.019-0.00336i) (x+0.9621)(x-0.3332) (x-1+0.0144i) (x-1-0.0144i)
= aˆ0 (aˆ1x + bˆ1)2 (aˆ2 x + bˆ2 )3 (aˆ3 x + bˆ3 )1
≈ ≈ ≈ ≈
1. Match multiplicities
AIF
≈ ≈
~ 2. aˆi − a~i + bˆi − bi = O( p − ~ p)
well-posed?
≈
≈
~ 2 ~ ~ 3 ~ ~ 1 ~ ~ ~ p ( x) = a0 (a1 x + b1 ) (a2 x + b2 ) (a3 x + b3 ) = 0.9999(1.0 x − 1.00001) 2 (1.0 x + 1.000009) 3 (1.0 x − 0.3334)1
5
A well-posed problem: (Hadamard, 1923) the solution satisfies
•
• •
existence uniqueness continuity w.r.t data
An ill-posed problem is infinitely sensitive to perturbation tiny perturbation
Î
huge error
A frontier in scientific computing Though frequently needed in application, the adequate handling of such ill-posed … problems is hardly ever touched upon in numerical analysis textbooks. --- Arnold Neumaier, SIAM Review
6
Challenge in solving ill-posed problems: Can we recover the lost solution when the problem is inexact?
P
Solution
P
Data
P : Data Æ Solution
P 7
Geometry of AIF (simplified view) (x – t)3
Π(3)
= –t3 + ( 3 t2) x + (–3 t) x2 + x3
⎡− t 3 ⎤ ⎢ ⎥ F (t ) = ⎢ 3t 2 ⎥ ⎢− 3t ⎥ ⎣ ⎦
F
F-1 diffeomorphism
t
Π (1,2) (x – u)1 (x – v)2
= –uv2 + (v2+ 2uv) x + (–2v–u) x2 + x3
⎡ − uv 2 ⎤ ⎢ 2 ⎥ G (u, v ) = ⎢v + 2uv ⎥ ⎢ − 2v − u ⎥ ⎣ ⎦
v
G
u
G-1 diffeomorphism
Polynomials form (factorization) manifolds
8
Proposition 1:
Polynomials
c0 + c1 x + c2 x 2 + L + cm x m = a0 (a1 x + bˆ1 ) k1 (a2 x + b2 ) k2 L (an x + bn ) kn
form a differentiable manifold Π[k
1 ,…, kn]
of codimension
codim (Π[k1 ,…, kn] ) = m - n
Number of factors dimension of the polynomial vector space (≥ polynomial degree) 9
Are ill-posed problems really sensitive? Kahan: It is a misconception. W. Kahan’s observation (1972) • Problems form a “pejorative manifolds”
Plot of pejorative manifolds of degree 3 polynomials with multiple roots
• Ill-posedness: a tiny perturbation pushes the problem out of the manifold
• A problem is not sensitive at all if it stays on the manifold.
10
Stratification of factorization manifolds of degree 3 monic polynomials
Π (1,1,1) = {p( x ) = ( x − α ) ( x − β ) ( x − γ ) α ≠ β ≠ γ } 1
1
1
Π (1,2) = {p( x ) = ( x − α )1 ( x − β )2 α ≠ β }
Π (3) = {p( x ) = ( x − α )3 α ∈ C }
Π (3) ⊂ Π (1,2) ⊂ Π (1,1,1) = C 3 Codimensions: 2
1
0
Factorization manifold stratification of degree 4 polynomials:
Π ( 2,2) Π (1,1,2)
Π ( 4)
Π (1,1,1,1)
= C4
Π (1,3) Codimensions: 3
2
1
0
11
Factorization manifolds and their stratification
{
Π[ k1k 2 Lk n ] = a0 ( a1 x + b1 ) k1 ( a2 x + b2 )k 2 L( an x + bn ) k n
ai , bi ∈ C , ai b j ≠ a j bi , ∀i ≠ j
}
⊂ Cm [ x ] = { c0 + c1 x + L + cm x m ci ∈ C }
p ∈ ∏[2,2]
Proposition 3:
⇄
dist(p, ∏[2,2] ) = dist(p, ∏[2,1,1] ) = dist(p, ∏[1,1,1,1] ) = 0
p ∈ ∏[k1… kn]
if and only if
codim(∏[k1… kn] ) = max { codim(∏) | dist(p,∏ ) = 0 } 12
Formulation of the approximate irreducible factorization
Π
p The approximate factorization of - the exact factorization of -
p
pˆ
is
~ p
~ p
~ p
lies in the nearby manifold ∏ of the highest codimension
~ p from
is the nearest polynomial on ∏
p
31
A “three-strikes” principle for formulating an approximate irreducible factorization: • Backward nearness: The AIF is the exact factorization of a nearby polynomial • Maximum codimension: The AIF is the exact factorization of a polynomial in the nearby factorization manifold of the highest codimension.
• Minimum distance: The AIF is the exact factorization of the nearest polynomial in the nearby factorization manifold of the highest codimension.
In comparison: Symbolic computation: Backward nearness with distance = 0
Numerical computation: (straightforward) Backward nearness with minimal distance
Finding the AIF is (apparently) a well-posed problem The AIF is a generalization of exact factorization.
13
pˆ ( x) = aˆ0 (aˆ1 x + bˆ1 ) k1 L (aˆ n x + bˆn ) k n ∈ Π k
Theorem 1 Assume Then,
p − pˆ
is small
∃ an interval I and ∀ε∈ I
∃ a unique AIF within ε ~ ~ ~ p ( x) = a~0 (a~1 x + b1 ) k1 L (a~n x + bn ) k n ∈ Π k
Πκ
p pˆ
such that
~ (a~i x + bi ) − (aˆi x + bˆi ) = O ( p − pˆ
)
~ p
Moreover, the AIF is continuous w.r.t. p
17
The two-staged algorithm after formulating the AIF of p within ε
Stage I: Find the factorization manifold Π of the highest dimension s.t.
dist ( P, Π k ) < ε
Πκ
p pˆ ~ p
Stage II: Find/solve problem Q such that
p− ~ p = min p − q q∈Π k
18
Stage I: Example:
Identify the AIF manifold by a squarefree factorization
p15 p23 p33 p4
= ( p1 p2 p3 p4 )( p1 p2 p3 )( p1 p2 p3 )( p1 )( p1 )
= ( p 4 ) (1 ) ( p 2 p 3 ) (1 ) ( p1 ) 1
2
3
4
5
--- flat SFF
--- staircase SFF
A new staircase SFF algorithm:
p ≈ u11 u22 L pkk
a0 (a1 x + b1 ) k1 (a2 x + b2 ) k 2 L (an x + bn ) k n
19
Stage II:
Minimize the distance to the AIF manifold
a0 (a1 x + b1 ) k1 (a2 x + b2 ) k 2 L (an x + bn ) k n = p α1a1 + β1b1 = γ1
α 2 a2 + β 2b2
= γ2 O
M M
α n an + β nbn = γ n
G(a0, a1 ,…,an, b1 ,…,bn) = d Nonlinear least squares problem solved by the Gauss-Newton iteration
20
Example: For polynomial ( x − 1)80 ( x − 2)60 ( x − 3) 40 ( x − 4)20 with (inexact ) coefficients in hardware precision
Exact factorization
Approximate factorization: >> [F,res,fcnd]
=
uvFactor(f,1e-10,1);
THE CONDITION NUMBER: THE BACKWARD ERROR: THE ESTIMATED FORWARD ROOT ERROR:
914.329 5.71e-015 1.04e-011
FACTORS ( ( ( (
x x x x
-
3.999999999999990 3.000000000000008 1.999999999999998 1.000000000000000
)^20 )^40 )^60 )^80 21
Summary: • Factorizations are sensitive because polynomials form manifolds of positive codimensions in strata. • An AIF can be formulated as an exact factorization of the nearest polynomial on a nearby manifold of the highest codimension. • The AIF approximates the factorization of the (hidden) underlying polynomial from the perturbed data. • The AIF can be computed by an (improved) algorithm in two stages.
Software is available in the package ApaTools (google apatools)
