Transcript A numerical Approach toward Approximate Algebraic Computatition Zhonggang Zeng Northeastern Illinois University, USA Oct.

A numerical Approach toward
Approximate Algebraic Computatition
Zhonggang Zeng
Northeastern Illinois University, USA
Oct. 18, 2006, Institute of Mathematics and its Applications
What would happen
when we try
on
numerical computation
algebraic problems?
A numerical analyst got a surprise 50 years ago
on a deceptively simple problem.
1
James H. Wilkinson (1919-1986)
Start of project:
Completed:
Add time:
Input/output:
Memory size:
Memory type:
Technology:
Floor space:
Project leader:
1948
1950
1.8 microseconds
cards
352 32-digit words
delay lines
800 vacuum tubes
12 square feet
J. H. Wilkinson
Britain’s Pilot Ace
2
The Wilkinson polynomial
p(x) = (x-1)(x-2)...(x-20)
= x20 - 210 x19 + 20615 x18 + ...
 ( x  0.99651)1 ( x  2.3145)2 ( x  5.7222)5 ( x 11.98)7 ( x 18.379)5
Wilkinson wrote in 1984:
Speaking for myself I regard it as the most traumatic
experience in my career as a numerical analyst.
3
Matrix rank problem
5
Factoring a multivariate polynomial:
A factorable polynomial
approximation
irreducible
6
Solving polynomial systems:
Example: A distorted cyclic four system:
Translation: There are two 1-dimensional solution set:
 z   6 , z  t, z   3
 1
2
3
t
t

z4   t


3
7
Distorted Cyclic Four system in floating point form:
approximation
1-dimensional solution set
Isolated solutions
8
tiny perturbation
in data
(< 0.00000001)
huge error
In solution
( >105 )
9
What could happen in approximate algebraic computation?
• “traumatic” error
• dramatic deformation of solution structure
• complete loss of solutions
• miserable failure of classical algorithms
•
•
•
•
•
Polynomial division
Euclidean Algorithm
Gaussian elimination
determinants
……
10
So, why bother with approximation in algebra?
1. You may have no choice
(e.g. Abel’s Impossibility Theorem)
Either
or
All subsequent computations become approximate
11
So, why bother with approximation solutions?
1. You may have no choice
2. Approximate solutions are better!
true image
Application: Image restoration (Pillai & Liang)
blurred image
blurred image
p( x, y)
p( x, y)( x, y)   ( x, y)
 f ( x, y )
p( x, y )( x, y)   ( x, y )
 g ( x, y)
GCD( f ( x, y ), g ( x, y ))  1
Application: Image restoration (Pillai & Liang)
blurred image
p( x, y)( x, y)   ( x, y)
 f ( x, y )
true image
blurred image
p( x, y )( x, y)   ( x, y )
p( x, y)
 g ( x, y)
AGCD( f ( x, y ), g ( x, y ))  ~
p( x, y )
restored image
Approximate solution is better than exact solution!
13
Perturbed Cyclic 4
Exact solutions by Maple:
16 isolated (codim 4) solutions
Approximate solutions
by Bertini
(Courtesy of
Bates, Hauenstein,
Sommese, Wampler)
Perturbed Cyclic 4
Exact solutions by Maple:
16 isolated (codim 4) solutions
Or, by an experimental approximate elimination combined with approximate GCD
Approximate solutions are better than exact ones , arguably
15
So, why bother with approximation solutions?
1. You may have no choice
2. Approximate solutions are better
3. Approximate solutions (usually) cost less
Example: JCF computation
Special case:
r



 




1
r
1
r
s
1
s








t
r  2,
s  3,
t 5
Maple takes 2 hours
On a similar 8x8 matrix, Maple and Mathematica run out of memory
16
Pioneer works in numerical algebraic computation (incomplete list)
• Homotopy method for solving polynomial systems
(Li, Sommese, Wampler, Verschelde, …)
• Numerical Algebraic Geometry
(Sommese, Wampler, Verschelde, …)
• Numerical Polynomial Algerba
(Stetter)
17
What is an “approximate solution”?
To solve
x2  2x 1  0
exact computation
108  0.00000001
( x  0.9999)(x  1.0001)  0

( x  1) 2  (10 4 ) 2  0
x 1
104  0.0001
axact solution
forward error
backward error
x2  2x 1  0
with 8 digits precision:
x  0.9999, 1.0001,
backward error:
0.00000001
-- method is good
forward error:
0.0001
-- problem is bad
18
The condition number
[Forward error] < [Condition number] [Backward error]
From problem
From numerical method
A large condition number
<=> The problem is sensitive or, ill-conditioned
An infinite condition number <=> The problem is ill-posed
19
Wilkinson’s Turing Award contribution:
Backward error analysis
• A numerical algorithm solves a “nearby” problem
• A “good” algorithm may still get a “bad” answer,
if the problem is ill-conditioned (bad)
20
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
Ill-posed problems are common in applications
- image restoration
- IVP for stiction damped oscillator
- some optimal control problems
- air-sea heat fluxes estimation
……
-
deconvolution
inverse heat conduction
electromagnetic inverse scatering
the Cauchy prob. for Laplace eq.
21
Ill-posed problems are common
in algebraic computing
- Multiple roots
- Polynomial GCD
- Factorization of multivariate polynomials
- The Jordan Canonical Form
- Multiplicity structure/zeros of polynomial systems
- Matrix rank
22
If the answer is highly sensitive to perturbations, you
have probably asked the wrong question.
Maxims about numerical mathematics, computers,
science and life, L. N. Trefethen. SIAM News
Does that mean:
(Most) algebraic problems are wrong problems?
A numerical algorithm
seeks the exact solution
of a nearby problem
Ill-posed problems
are infinitely sensitive
to data perturbation
Conclusion: Numerical computation is incompatible
with ill-posed problems.
Solution: Formulate the right problem.
23
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
24
Are ill-posed problems really sensitive to perturbations?
William Kahan:
This is a misconception
Kahan’s discovery in 1972:
Ill-posed problems are sensitive to arbitrary perturbation,
but insensitive to structure preserving perturbation.
25
Why are ill-posed problems infinitely sensitive?
W. Kahan’s observation (1972)
• Problems with certain solution structure form a “pejorative manifold”
Plot of pejorative manifolds of degree 3 polynomials with multiple roots
• The solution structure is lost when the problem leaves the manifold due to an
arbitrary perturbation
• The problem may not be sensitive at all if the problem stays on the manifold,
unless it is near another pejorative manifold
26
Geometry of ill-posed algebraic problems
* Rank r matrices M rmn  A  C mn | rank( A)  r
- - codimM rmn   (m  r)(n  r)
M 0mn  M1mn    M nmn
* Polynomi
al pairs Prm,n  ( p, q) | deg(GCD( p, q))  r
- - codimPrmn   r
Pnmn  Pnm1n    P0mn
Similar manifold stratification exists for problems like
factorization, JCF, multiple roots …
27
Manifolds of 4x4 matrices defined by Jordan structures
(Edelman, Elmroth and Kagstrom 1997)
e.g. {2,1} {1}
is the structure of 2 eigenvalues in 3 Jordan blocks of sizes 2, 1 and 1
28
Illustration of pejorative manifolds
codimsion  0
codimsion  3
?
B
?
A
codimensio n  2
codimsion  1
Problem A
perturbation
Problem B
The “nearest” manifold may not be the answer
The right manifold is of highest codimension within a certain distance
29
A “three-strikes” principle for formulating
an “approximate solution” to an ill-posed problem:
• Backward nearness: The approximate solution is the exact solution of a nearby problem
• Maximum codimension: The approximate solution is the exact solution of a problem
on the nearby pejorative manifold of the highest codimension.
• Minimum distance: The approximate solution is the exact solution of the nearest problem
on the nearby pejorative manifold of the highest codimension.
Finding approximate solution is (likely) a well-posed problem
Approximate solution is a generalization of exact solution.
30
Continuity of the approximate solution:
Formulation of the approximate rank /kernel:
 A  C mn and    0
Backward nearness: app-rank of A is
The approximate rank of A within  :
rank ( A)  min rank( B)
B  A 
the exact rank of certain matrix B within .
Maximum codimension: That matrix B
is on the pejorative manifold P possessing the
highest co-dimension and intersecting the
neighborhood of A.
The approximate kernel of A within  :
Ker ( A)  Ker( B)
B A2 
min
with
rank ( C )  rank  ( A)
CA2
Minimum distance: That B is the nearest
matrix on the pejorative manifold P.
• An exact rank is the app-rank within sufficiently small .
• App-rank is continuous (or well-posed)
31
=4
kernel
nullity = 2
Rank
basis
+ E = 6
nullity = 0
After reformulating the rank:
0      4.98
=4
nullity = 2
Rank
+ E = 4
dist ( Ker ( A)  Ker ( A  E )) 
61.26

1 
nullity = 2
Ill-posedness is removed successfully.
App-rank/kernel can be computed by SVD and other rank-revealing algorithms
(e.g. Li-Zeng, SIMAX, 2005)
32
 ( f , g )  C[ x1,, xl ],    0, deg( f )  m, deg( g )  n
Formulation of the approximate GCD
Pjm,n  ( p, q)  C[ x1 ,, xl ]
deg( p)  m, deg( q)  n, degGCD( p, q)   j

 j ( f , g )  inf ( f , g )  ( p, q ) ( p, q )  Pjm ,n
Pkm1,n
( f , g)

( p, q)
Pkm,n
Pkm1,n
The AGCD within  :
AGCD ( f , g )  EGCD( p, q)
( f , g )  ( p, q)  minm ,n ( f , g )  (u, v)   k ( f , g )

( u ,v )Pk

k  max codim(Pjm,n )
 j ( f , g )
• Finding AGCD is well-posed if k(f,g) is sufficiently small
• EGCD is an special case of AGCD for sufficiently small 
(Z. Zeng, Approximate GCD of inexact polynomials, part I&II)
33
Similar formulation strikes out ill-posedness in problems such as
• Approximate rank/kernel
(Li,Zeng 2005, Lee, Li, Zeng 2006)
• Approximate multiple roots/factorization (Zeng 2005)
• Approximate GCD
(Zeng-Dayton 2004, Gao-Kaltofen-May-Yang-Zhi 2004)
• Approximate Jordan Canonical Form (Zeng-Li 2006)
• Approximate irreducible factorization (Sommesse-Wampler-Verschelde 2003,
Gao et al 2003, 2004, in progress)
• Approximate dual basis and multiplicity structure
(Dayton-Zeng 05, Bates-Peterson-Sommese ’06)
• Approximate elimination ideal (in progress)
34
The two-staged algorithm
after formulating the approximate solution to problem P within 
Stage I: Find the pejorative manifold
P of the highest dimension s.t.
dist(P, P)  
P
P
S
Q
Stage II: Find/solve problem Q such that
P Q
 min P  R
RP
Exact solution of Q is the approximate solution of P within 
which approximates the solution of S where P is perturbed from
35
Case study: Univariate approximate GCD:
 ( f , g )  C[ x],    0, deg( f )  m, deg( g )  n
Stage I: Find the pejorative manifold


dist ( f , g ), Pkm,n     min Sk ( f , g )   m  n
where Sk ( f , g ) is thematrix for (v,w)  f  w  g  v
with deg(v)  m  k , deg( w)  n  k
( f  u  v and g  u  w  f  w  g  v  0)
Stage II: solve the (overdetermined) quadratic system
 uv

u  w
 (u )


f

g

1
F (u, v, w)  b( f , g )
for a least squares solution (u,v,w) by Gauss-Newton iteration
(key theorem: The Jacobian of F(u,v,w) is injective.)
36
Start: k = n
Max-codimension
k := k-1
k := k-1
no
Is AGCD of degree k
possible?
no
probably
Refine with
G-N Iteration
Successful?
yes
Min-distance
nearness
Output GCD
Univariate AGCD algorithm
37
Case study: Multivariate approximate GCD:


 ( f , g )  C[ x1,, xl ],    0, deg( f )  m, deg( g )  n
Stage I: Find the max-codimension pejorative manifold by
applying univariate AGCD algorithm on each variable xj
 f  u  v and g  u  w
 f ( , x j , )  u ( , x j , ) v( , x j , )
g ( , x j , )  u ( , x j , ) w( , x j , )
Stage II: solve the (overdetermined) quadratic system
 uv

u  w
 (u )


f

g

1
F (u, v, w)  b( f , g )
for a least squares solution (u,v,w) by Gauss-Newton iteration
(key theorem: The Jacobian of F(u,v,w) is injective.)
38
Case study: univariate factorization:
 f  C[ x],    0, deg( f )  n
Stage I: Find the max-codimension pejorative manifold by
applying univariate AGCD algorithm on (f, f’ )
 f ( x)  ( x  z1 ) m1  ( x  zk ) m k
 f ' ( x)  ( x  z1 ) m1 1  ( x  zk ) m k 1 q( x)

AGCD( f , f ' )  ( x  z1 ) m1 1  ( x  zk ) m k 1
Stage II: solve the (overdetermined) polynomial system
F(z1 ,…,zk )=f
(   z1 )m1 (   zk )mk  f ()
(in the form of coefficient vectors)
for a least squares solution (z1 ,…,zk ) by Gauss-Newton iteration
(key theorem: The Jacobian is injective.)
39
Case study: Finding the nearest matrix with a Jordan structure
l
   
   
Ax1 , x2 , x3 , x4   x1 , x2 , x3 , x4  J
J=
1
l
1
l
l
Segre characteristic = [3,1]
B
   
   
Au1, u2 , u3 , u4   u1, u2 , u3 , u4  (lI  S )
A
codim = -1 + 3 + 3(1) = 5
l
lI+S=
l
s13 s14
s23 s24
l s34
P
 
ui  u j   ij
l
Wyre characteristic = [2,1,1]
3
1
2
   
   
 Au1 , u2 , u3 , u4   u1 , u2 , u3 , u4  (lI  S )  0
    T    


u
 0

1 , u2 , u3 , u4  u1 , u2 , u3 , u4   I



b T u1
 0

   
F ( A, l, u1, u2 , u3, u4 , s13 , s23 , s14 , s24 , s34 )  0
1
1
Equations determining the manifold
Ferrer’s diagram
F ( A, l ,U , S )  0
40
Case study: Finding the nearest matrix with a Jordan structure
Equations determining the manifold
   
   
 Au1 , u2 , u3 , u4   u1 , u2 , u3 , u4  (lI  S )  0

u1, u2 , u3 , u4 Tu1, u2 , u3 , u4   I
 0

T

b
u1
 0

B
A
P
For B not on the manifold, we can still solve
F ( B, l ,U , S )  0
2
2
F ( B, l,U , S ) 2
for a least squares solution : F ( B, lˆ,Uˆ , Sˆ ) 2  min
l ,U ,S
When
BU  U (lI  S ) 2
is minimized, so is
B  U (lI  S )U T
The crucial requirement: The Jacobian J ( A,,,) of
2
F ( A,,,)
 B A2
is injective.
(Zeng & Li, 2006)
41
Solving G(z) = a
a
Solve
G(z0)+J(z
G( z0)(
) =z -az0 ) = a
for nonlinear
for linearleast
leastsquares
squaressolution
solutionz=z
z =*z1
G(z
- z0 0)))=-- aa]
a]
J(z
-[J(z
z0 ) 00=))(+-]z[G(z
0)(
z1 =
z0z0-)+J(z
[G(z
0
42
Stage I: Find the nearby max-codim manifold
Stage II: Find/solve the nearest problem on the manifold
G(z)=a
. ||G(z*)-a||=minz ||G(z)-a||
via solving an overdetermined system
for a least squares solution z* s.t
by the Gauss-Newton iteration
zk 1  zk  J ( zk ) G( zk )  a, k  0,1,2,
Key requirement: Jacobian J(z*) of G(z) at z* is injective
(i.e. the pseudo-inverse exists)
z  zˆ 
1
G( z )  G( zˆ)  h.o.t

J ( zˆ)
G (zˆ )
G (z )
condition number
(sensitivity measure)
43
Summary:
• An (ill-posed) algebraic problem can be formulated using the three-strikes
principle (backward nearness, maximum-codimension, and minimum distance)
to remove the ill-posedness
• The re-formulated problem can be solved by numerical computation in two stages
(finding the manifold, solving least squares)
• The combined numerical approach leads to Matlab/Maple toolbox ApaTools for
approximate polynomial algebra. The toolbox consists of
univariate/multivariate GCD
matrix rank/kernel
dual basis for a polynomial ideal
univariate factorization
irreducible factorization
elimination ideal
……
(to be continued in the workshop
next week)
44