Root finding via subdivision/IA

  • 6 Replies
  • 168 Views

0 Members and 1 Guest are viewing this topic.

Offline marcm200

  • *
  • Fractal Freak
  • **
  • Posts: 759
« on: July 11, 2020, 09:12:57 AM »
A subdivision IA-based algorithm to find roots of a simple polynomial based on the Krawczyk operator (follows Subdivision Algorithms for Complex Root Isolation: Empirical Comparisons, N. Kamath, master thesis, 2010.).

\[
K(A) = m - Y\cdot f(m) - (I - Y\cdot J(A))\cdot(A-m)
 \]

where

  • f is the polynomial function to find zeros for
  • A is the input complex square to judge if it contains a root
  • m is an arbitrarily chosen point inside A (I use the center)
  • Y is an abritrary invertible matrix - real valued or real interval-valued (can be the inverse of the Jacobian, I use {{2,-0.25},{-0.5,1.5}}
  • J(A) is the Jacobian of f in its interval extension at square A

Algebraic operations involving different types (A-m e.g.) are enhanced appropriately (m is treated as a degenerate complex interval with zero width in its components).

The algorithm starts with a list of complex squares A (here the entire 8-square) and goes over it repeatedly, computing the operator K(A), then using one test:

  • exclusion predicate If \( K(A) \cap A = \emptyset \) => A definitely does not contain a root and is discarded.
  • inclusion If \( K(A) \subseteq A \) => A definitely contains a root and is subdivided into 2x2 squares put into the list.
  • else nothing can be deduced and A is subdivided as above.

Subdivision occurs for a specified number of refinement steps and the union of all nondiscarded squares after the last refinement cover all roots (provided one starts with a sufficiently large box). Usually there are more boxes left than there are roots.

The image below is a first example for \( z^3+iz+(-2+2.5i) \) for level 50 refinement (for the output list squares, not the image) - and it resulted in only a dozen possible root regions. White regions do not contain a root, the circumferencing red rectangle indicates the size at which those regions could be discarded. Their area decrease leads visually the way to the possible root positions - single gray pixels in a white desert. The right image is a point-sampled classic Newton control, using a modified universal set of starting points forum link.

It was nice to see that the subdivision algorithm also judges the two bottom roots as slightly not on an axis parallel line. I like those subtle features as a plausibility test.

All calculations were done error-free using my custInt library and outward rounding with double precision at interval end points. The image itself is however only a visualization and not error-controlled (pixels can therefore be off by +-1 in position).

The exclusion predicate works quite well, but so far I haven't had a positive inclusion result. I am currently exploring different Y matrices and try to find out whether the outward rounding (computing K(A) involves a lot of operations) always enlarges K(A) to such an extent that it hardly can get smaller than A.

Has anyone experience with subdivision/IA and especially inclusion tests? For my purposes, speed is not a concern, neither is finding all roots, but my subsequent analysis relies on the information that a region in the complex plane definitely contains a root.


Linkback: https://fractalforums.org/fractal-mathematics-and-new-theories/28/root-finding-via-subdivisionia/3658/

Offline marcm200

  • *
  • Fractal Freak
  • **
  • Posts: 759
« Reply #1 on: July 31, 2020, 09:51:49 AM »
"Newton rules!"

Using the Newton interval operator and the mean-value-form for the partial derivatives to compute the inverse of the interval Jacobian led to the detection of boxes that definitely contain a root.

\[
N(B) = M - J^{-1}(B)\cdot f(M)
 \]

where M is the midpoint of the box B.

The left image below shows the subdivision of a random degree-6 polynomial with interval coefficients of width 2^-18. Each black box B contains definitely a root (N(B) is a subset of B). The roots were found at different refinements indicated by their differing size. The unit circle is drawn in blue.

The right image is showing my intended goal - a very simple result image for a rather complicated polynomial. It is also degree-6 and describes the period-5 hyperbolic components of the quadratic Mandelbrot set forum link, where the coefficients were calculated for the complex interval c=[a..a+w] + i*[b..b+w] with a=-131/256, b=-283/512 and w=4/2048 which is around -0.5117,-0.5527 and a numerical estimate of a region deep within the south west period-5 bulb.

The polynomial itself is (numbers are rounded outwards for display purposes for at most 2 decimals):


\[
f(z)= z^6 + ([-177.03..-177.02]+i[35.78..35.79])\cdot z^5 + ([10210.00..10210.18]+i[-1892.33..-1892.15])\cdot z^4 + ([-481051..-481028]+i[61052..61078])\cdot z^3\\ + ([14590943..14593060]+i[-3515628..-3513511])\cdot z^2 + ([-193307003..-193251224]+i[44146251..44202030])\cdot z + ([30388361..36183842]+i[50842801..56638282])
 \]


The intervals of the coefficients are quite large at lower degree terms but still result in a positive outcome.

As this interval polynomial has definitely a root within the unit circle, all individual complex numbers in the initial c interval have, so c resides in the interior of a period-5 component (I do not know where the other roots lie). As a negative control I used a small c-interval around the center of the period-2 bulb and this returned a white blank image, so definitely no root inside the unit circle.

The next step is to evaluate how far I can come with double precision at the end points of the interval number type.

Technical details
  • The midpoint M can be a small interval in case its coordinates need outward rounding during calculation (0.5*x_left+0.5*x_right).
  • Subdividing a box stops as soon as it definitely contains a root.
  • All calculations were done using the custInt library with double precision endpoints.
  • The C++ code for the mean-value form and the calculation path for the coefficient computation were generated in an automated manner through a symbolic parser. Numbers used therein were validated by maxima to fit into double precision.
  • Inverting the Jacobian was done numerically, first evaluating the partial derivatives to arrive at a 2x2 matrix whose entries were real-valued intervals.



Offline Adam Majewski

  • *
  • Fractal Flamingo
  • ****
  • Posts: 346
« Reply #2 on: July 31, 2020, 04:28:56 PM »
see also :
Empirical Study of an Evaluation-Based Subdivision Algorithm for Complex Root Isolation
June 2012 DOI: 10.1145/2331684.2331710 by Narayan KamathIrina VoiculescuIrina VoiculescuChee Keng YapChee Keng Yap

Do you use core library ?

Offline marcm200

  • *
  • Fractal Freak
  • **
  • Posts: 759
« Reply #3 on: July 31, 2020, 04:45:51 PM »
No, I don't use core. As I enjoy the implementation and optimization process as much as using the final software to arrive at my initial goal (here a lower estimate on the area of the Mandelbrot set), I'm coding everything myself (except for some things like interval transcendental functions which I only need a few times).

The article seems to be a follow-up of the master thesis of the first author (post #1) as they use the same interval predicates. As for me the information "root there" is sufficient, the clever approach of CEVAL is not necessary (arc crossings, quite interesting though), but I have to see how often the pure interval predicates fail.

Offline marcm200

  • *
  • Fractal Freak
  • **
  • Posts: 759
« Reply #4 on: August 02, 2020, 11:50:00 AM »
I uploaded the code for the subdivision agorithm to:

https://github.com/marcm200/subdiv-core

A general degree-10 polynomial in all 11 (complex interval) coefficents is implemented.

The main output is a list of definitely root-containing regions of adjustable width (please see README for command-line parameters). In addition, an overview subdivision process image is created.

Example polynomials _z2.txt and _p5.txt are provided. A local version of the custInt library is uploaded.

Quick start:
Code: [Select]
subdiv-core file=_p5.txt range=-2.5,2.25,-2.5,2.25 depth=16 stop=width,0.0625

Exemplary output for z2-2i=0:
Code: [Select]
    #0 DEFINITE:
    exact [9007199201053893*2^-53..4503599635759099*2^-52]+i*[9007199201053893*2^-53..4503599635759099*2^-52])
    approx [0.99999999403953466..1.000000001862644] x [0.99999999403953466..1.000000001862644]
    ||squared=[9007199147366794*2^-52....4503599644147703*2^-51] (approx [2..2]


Offline Adam Majewski

  • *
  • Fractal Flamingo
  • ****
  • Posts: 346
« Reply #5 on: August 02, 2020, 07:00:42 PM »

Offline marcm200

  • *
  • Fractal Freak
  • **
  • Posts: 759
« Reply #6 on: August 02, 2020, 08:15:49 PM »
For reliable computing there's (afaik) box interval, ball arithmetics and affine arithmetics as an underlying method. I started with boxes as the TrueShape Julia set algorithm by Figueiredo uses a square tiling of the plane, so that's where I have the most experience. I haven't worked much with ball arithmetics as I currently do not see how the complex plane can be tiled without overlapping balls, which then makes drawing reliable Julia set images needing a different definition of what the image shows.

I like the simple 2x2 subdivision rule of boxes: just take the vertical and horizontal middle line (or anything therein) and use those 4 adjacent and touching, but not overlapping boxes for further root location.

But I may resort to ball arithmetics in the parabolic case, esp. z^2-0.75, as there box tiling the plane does not lead to interior detection. I'm (slowly) working on a triangular tiling and might try ball tiling there as well.

As specifically for root finding, I haven't used balls so far. Is there a particular advantage over boxes, maybe in terms of convergence or overestimation?