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 methods only find local min/max ◮ Finding the common zeros of the partial derivatives gives global min/max

The Problem Find the common roots of f 1 , f 2 , . . . , f n in a given compact set C ⊂ R n . (i.e. x ∈ C such that f 1 ( x ) = f 2 ( x ) = · · · = f n ( x ) = 0.) 1 f ( x , y ) = sin ( xy ) + x log ( y + 3) − x 2 + y − 4 g ( x , y ) = cos (3 xy ) + exp ( 3 y x − 2) − x − 6

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.

Current Methods ◮ Newton’s Method ◮ Homotopy continuation (Bertini) ◮ Spectral Methods ◮ Subdivision

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

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 ).

Companion Matrix If p ( x ) = x n + a n − 1 x n − 1 + ... + a 1 x + a 0 , then   0 1 0 · · · 0 0 0 1 · · · 0   . . . .  ...  . . . . A =  . . . .      0 0 0 · · · 1   − a 0 − a 1 − a 2 · · · − a n − 1 has p ( x ) as its characteristic polynomial.

Generalized Companion Matrix Can the companion matrix be generalized to systems of multivariate polynomials?

Generalized Companion Matrix ”multiplication-by- x ” operator on quotient algebra C [ x ] / � p ( x ) � 1 �→ x x �→ x 2 . . . x n − 2 �→ x n − 1 x n − 1 �→ x n ≡ − a n − 1 x n − 1 − · · · − a 1 x − a 0 mod p ( x ) Companion matrix: representation of this operator (standard basis)   0 1 0 · · · 0 0 0 1 · · · 0    . . . .  ... . . . .   . . . .     0 0 0 · · · 1   − a 0 − a 1 − a 2 · · · − a n − 1

M¨ oller-Stetter Matrix ◮ A := C [ x 1 , . . . , x n ] / � f 1 , . . . , f n � ◮ dim ( A ) = # of common roots of f 1 , . . . , f n . M h : g �→ hg for some h ∈ C [ x 1 , ..., x n ] / � f 1 , ..., f n � ◮ Eigenvalues : h ( x 1 ) , . . . , h ( x N ) ◮ Left-eigenvectors : [ b 1 ( x i ) , . . . , b N ( x i )] T ◮ Common roots can be extracted from the eigenvalues/eigenvectors

M¨ oller-Stetter Matrix ◮ A := C [ x 1 , . . . , x n ] / � f 1 , . . . , f n � ◮ dim ( A ) = # of common roots of f 1 , . . . , f n . M h : g �→ hg for some h ∈ C [ x 1 , ..., x n ] / � f 1 , ..., f n � ◮ Eigenvalues : h ( x 1 ) , . . . , h ( x N ) ◮ Left-eigenvectors : [ b 1 ( x i ) , . . . , b N ( x i )] T ◮ Common roots can be extracted from the eigenvalues/eigenvectors

M¨ oller-Stetter Matrix ◮ A := C [ x 1 , . . . , x n ] / � f 1 , . . . , f n � ◮ dim ( A ) = # of common roots of f 1 , . . . , f n . M h : g �→ hg for some h ∈ C [ x 1 , ..., x n ] / � f 1 , ..., f n � ◮ Eigenvalues : h ( x 1 ) , . . . , h ( x N ) ◮ Left-eigenvectors : [ b 1 ( x i ) , . . . , b N ( x i )] T ◮ Common roots can be extracted from the eigenvalues/eigenvectors

M¨ oller-Stetter Matrix Example: A = C [ x , y , z ] / � f 1 , f 2 , f 3 � Eigenvalues Method: ◮ Eigenvalues of M x are x-coordinates of roots ◮ Make M x , M y , M z EIgenvectors Method: ◮ Include 1 , x , y , z in basis for A ◮ ∀ h , eigenvectors of M h will look like [ c , cx , cy , cz , ... ] T ◮ Compute the roots as ( x , y , z ) = ( cx c , cy c , cz c )

Challenge Compute a basis for A and a matrix representation of M h in higher dimensions

Macaualy Matrix � f 1 : x 1 x 2 − 2 x 2 = 0 f 2 : x 2 − 3 = 0 x 3 x 2 x 1 x 2 x 3 x 2 x 2 1 x 2 x 1 x 2 x 1 x 2 1 1 2 2 1 2 x 2 1 f 2 0 1 0 0 − 3 0 0 0 0 0   x 1 x 2 f 2 0 0 1 0 0 − 3 0 0 0 0   x 2 2 f 2 0 0 0 1 0 0 − 3 0 0 0     x 1 f 2 0 0 0 0 0 1 0 − 3 0 0     x 2 f 2  0 0 0 0 0 0 1 0 − 3 0    f 2  0 0 0 0 0 0 0 0 1 − 3    x 1 f 1  0 1 0 0 0 − 2 0 0 0 0    0 0 1 0 0 0 − 2 0 0 0 x 2 f 1   0 0 0 0 0 1 0 0 − 2 0 f 1

Telen-Van Barel Method ◮ Reduce the Macaulay Matrix to Echelon form ◮ Monomials of free columns form a basis for A ◮ Every row is in � f 1 , ..., f n � ◮ Relations to find M h : g �→ hg for some h ∈ C [ x 1 , ..., x n ] / � f 1 , ..., f n � ◮ Reducing with pivoting make the reduction more stable

Subdivision Method Alex Townsend ◮ Subdivide the region ◮ Chebyshev Approximate ◮ Bezout Resultant ◮ Unstable in 3+ dimensions

Chebyshev ◮ Orthogonal Basis of Polynomials ◮ Rootfinding Well-Conditioned ◮ Approximation can be found quickly with FFT ◮ Coefficients drop to 0 exponentially as degree increases

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

Numerical Instability ◮ Coefficient of highest degree monomials very small for Chebyshev approximations ◮ Highest degree monomials can’t be in vector basis for building M h ◮ Must be a pivot column of the matrix. ◮ Thus computing M h is numerically impossible!

Division Matrix Solution - M¨ oller Stetter Matrix using division ( M x − 1 ) instead of i multiplication ( M x i ) ◮ Increases Numerical Stability for Chebyshev approximations ◮ Numerically, x − 1 exists i ◮ Now terms with no x i must be pivot ◮ Much more stable numerically

Results 1 f ( x , y ) = sin ( xy ) + x log ( y + 3) − x 2 + y − 4 g ( x , y ) = cos (3 xy ) + exp ( 3 y x − 2) − x − 6

Results f ( x , y ) = sin(30 x − y 30) + y g ( x , y ) = cos( x 30 − 30 y ) − x

Results f ( x , y ) = sin (20 x + y ) g ( x , y ) = cos ( x 2 + xy ) − 1 4

Results Nick Trefethen’s Hundred-dollar, Hundred-digit Challenge f ( x , y ) = e sin(50 x ) +sin(60 e y )+sin(70 sin( x ))+sin(sin(80 y )) − sin(10( x + y ))+1 / 4( x 2 + y 2 )

Results

Results f ( x , y , z ) = sin (5 x + y + z ) g ( x , y , z ) = sin ( xyz ) h ( x , y , z ) = x 2 + y 2 − z 2 − 1

Timing Results

Future Work We have cool ideas to make it even better! ◮ Adaptive Subdivision ◮ Dimension Reduction ◮ Better Interval Checks Look up our code on GitHub! ◮ https://github.com/tylerjarvis/RootFinding