SLIDE 1
Numerical Rootfinding in a Compact Region Suzanna Stephenson - - PowerPoint PPT Presentation
Numerical Rootfinding in a Compact Region Suzanna Stephenson - - PowerPoint PPT Presentation
Numerical Rootfinding in a Compact Region Suzanna Stephenson Brigham Young University January 26, 2019 Introduction Optimization problems are extremely important Many problems have lots of variables, especially in machine learning Most
SLIDE 2
SLIDE 3
The Problem
Find the common roots of f1, f2, . . . , fn in a given compact set C ⊂ Rn. (i.e. x ∈ C such that f1(x) = f2(x) = · · · = fn(x) = 0.) f (x, y) = sin (xy) + x log (y + 3) − x2 + 1 y − 4 g(x, y) = cos (3xy) + exp ( 3y x − 2) − x − 6
SLIDE 4
Introduction
This problem is surprisingly difficult!
◮ Bad Complexity - Scales Poorly ◮ Numerical Instability
Roots of 20
i=1(x − i) = 0 in black.
Roots given when perturbing coefficients by 10−10 in red.
SLIDE 5
Current Methods
◮ Newton’s Method ◮ Homotopy continuation (Bertini) ◮ Spectral Methods ◮ Subdivision
SLIDE 6
Our Algorithm
Subdivision Method System of general functions
- 1. Subdivide
- 2. Chebyshev approximate
- 3. Eliminate sub-domains with no roots
- 4. For each remaining sub-domain,
◮ Solve Chebyshev approximations with modified TVB ◮ Eliminate extra roots
SLIDE 7
Companion Matrix
For a single univariate polynomial p(x), construct a matrix A with p(x) as its characteristic polynomial. Then use numerical methods to find the eigenvalues of A: the roots of p(x).
SLIDE 8
Companion Matrix
If p(x) = xn + an−1xn−1 + ... + a1x + a0, then A = 1 · · · 1 · · · . . . . . . . . . ... . . . · · · 1 −a0 −a1 −a2 · · · −an−1 has p(x) as its characteristic polynomial.
SLIDE 9
Generalized Companion Matrix
Can the companion matrix be generalized to systems of multivariate polynomials?
SLIDE 10
Generalized Companion Matrix
”multiplication-by-x” operator on quotient algebra C[x]/p(x) 1 → x x → x2 . . . xn−2 → xn−1 xn−1 → xn ≡ −an−1xn−1 − · · · − a1x − a0 mod p(x) Companion matrix: representation of this operator (standard basis) 1 · · · 1 · · · . . . . . . . . . ... . . . · · · 1 −a0 −a1 −a2 · · · −an−1
SLIDE 11
M¨
- ller-Stetter Matrix
◮ A := C[x1, . . . , xn]/f1, . . . , fn ◮ dim(A) = # of common roots of f1, . . . , fn.
Mh : g → hg for some h ∈ C[x1, ..., xn]/f1, ..., fn
◮ Eigenvalues : h(x1), . . . , h(xN) ◮ Left-eigenvectors : [b1(xi), . . . , bN(xi)]T ◮ Common roots can be extracted from the
eigenvalues/eigenvectors
SLIDE 12
M¨
- ller-Stetter Matrix
◮ A := C[x1, . . . , xn]/f1, . . . , fn ◮ dim(A) = # of common roots of f1, . . . , fn.
Mh : g → hg for some h ∈ C[x1, ..., xn]/f1, ..., fn
◮ Eigenvalues : h(x1), . . . , h(xN) ◮ Left-eigenvectors : [b1(xi), . . . , bN(xi)]T ◮ Common roots can be extracted from the
eigenvalues/eigenvectors
SLIDE 13
M¨
- ller-Stetter Matrix
◮ A := C[x1, . . . , xn]/f1, . . . , fn ◮ dim(A) = # of common roots of f1, . . . , fn.
Mh : g → hg for some h ∈ C[x1, ..., xn]/f1, ..., fn
◮ Eigenvalues : h(x1), . . . , h(xN) ◮ Left-eigenvectors : [b1(xi), . . . , bN(xi)]T ◮ Common roots can be extracted from the
eigenvalues/eigenvectors
SLIDE 14
M¨
- ller-Stetter Matrix
Example: A = C[x, y, z]/f1, f2, f3 Eigenvalues Method:
◮ Eigenvalues of Mx are x-coordinates of roots ◮ Make Mx, My, Mz
EIgenvectors Method:
◮ Include 1, x, y, z in basis for A ◮ ∀h, eigenvectors of Mh will look like
[c, cx, cy, cz, ...]T
◮ Compute the roots as
(x, y, z) = (cx c , cy c , cz c )
SLIDE 15
Challenge
Compute a basis for A and a matrix representation of Mh in higher dimensions
SLIDE 16
Macaualy Matrix
f1 : x1x2 − 2x2 = 0 f2 : x2 − 3 = 0
x3
1
x2
1x2
x1x2
2
x3
2
x2
1
x1x2 x2
2
x1 x2 1 1 −3 x2
1f2
1 −3 x1x2f2 1 −3 x2
2f2
1 −3 x1f2 1 −3 x2f2 1 −3 f2 1 −2 x1f1 1 −2 x2f1 1 −2 f1
SLIDE 17
Telen-Van Barel Method
◮ Reduce the Macaulay Matrix to Echelon form ◮ Monomials of free columns form a basis for A ◮ Every row is in f1, ..., fn ◮ Relations to find Mh: g → hg for some
h ∈ C[x1, ..., xn]/f1, ..., fn
◮ Reducing with pivoting make the reduction more stable
SLIDE 18
Subdivision Method
Alex Townsend
◮ Subdivide the region ◮ Chebyshev Approximate ◮ Bezout Resultant ◮ Unstable in 3+ dimensions
SLIDE 19
Chebyshev
◮ Orthogonal Basis of Polynomials ◮ Rootfinding Well-Conditioned ◮ Approximation can be found quickly with FFT ◮ Coefficients drop to 0 exponentially as degree increases
SLIDE 20
Putting it together
Subdivision Method System of general functions
- 1. Subdivide
- 2. Chebyshev approximate
- 3. Eliminate sub-domains with no roots
- 4. For each remaining sub-domain,
◮ Solve Chebyshev approximations with modified TVB ◮ Eliminate extra roots
SLIDE 21
Numerical Instability
◮ Coefficient of highest degree monomials very small for
Chebyshev approximations
◮ Highest degree monomials can’t be in vector basis for building
Mh
◮ Must be a pivot column of the matrix. ◮ Thus computing Mh is numerically impossible!
SLIDE 22
Division Matrix
Solution - M¨
- ller Stetter Matrix using division (Mx−1
i
) instead of multiplication (Mxi)
◮ Increases Numerical Stability for Chebyshev approximations ◮ Numerically, x−1 i
exists
◮ Now terms with no xi must be pivot ◮ Much more stable numerically
SLIDE 23
Results
f (x, y) = sin (xy) + x log (y + 3) − x2 + 1 y − 4 g(x, y) = cos (3xy) + exp ( 3y x − 2) − x − 6
SLIDE 24
Results
f (x, y) = sin(30x − y 30) + y g(x, y) = cos( x 30 − 30y) − x
SLIDE 25
Results
f (x, y) = sin (20x + y) g(x, y) = cos (x2 + xy) − 1 4
SLIDE 26
Results
Nick Trefethen’s Hundred-dollar, Hundred-digit Challenge
f (x, y) = esin(50x)+sin(60ey)+sin(70 sin(x))+sin(sin(80y))−sin(10(x+y))+1/4(x2+y2)
SLIDE 27
Results
SLIDE 28
Results
f (x, y, z) = sin(5x + y + z) g(x, y, z) = sin(xyz) h(x, y, z) = x2 + y2 − z2 − 1
SLIDE 29
Timing Results
SLIDE 30