Algebraic functions, iterated exponentials and polynomial systems: Fixed points of folded polynomial systems, part 1

$$ \newcommand{\bint}{\displaystyle{\int\hspace{-10.4pt}\Large\mathit{8}}} \newcommand{\res}{\displaystyle{\text{Res}}} \newcommand{Log}{\text{Log}} \newcommand{Arg}{\text{Arg}} \newcommand{pLog}{\text{pLog}} $$

A challenge for the interested reader:

After reading parts 1, 2, and 3 of Fixed points of folded polynomial systems, figure out what these two complex numbers represent (note the comma at the end of the second line): $$ \small \begin{align*} &\{0.02690400181196000919147192678327314910722771734900974123698018975335433562331571334721018629546791683229285582+\\ & 0.61402033212285429358434466607324792561123168919695925972343173572457123006580513164177839232470395409586828918i,\\ & 0.88692017760625646551460420340523084603788255282101275892501273453921267501344123646396443605740554814180281570+\\ & 0.31942170123013658964283374753891056582213182638856410603888905608438478484705058992835729232032645503873540847 i\} \end{align*} $$

Background:

Consider a polynomial system $F(x,y)=\{p_1(x,y),p_2(x,y)\}$ where $p_1$ and $p_2$ are polynomials with rational coefficients. A point $(x,y)$ is a fixed-point of $F(x,y)$ if $F(x,y)=(x,y)$. For example, let $$ F(x,y)=\left\{\begin{array}{c} 1/2+3/4x+1/3y^2 \\ 2/5+1/2x^2-4/5xy \end{array}\right\} $$ Then a fixed-point of $F(x,y)$ is a point $(x_f,y_f)$ such that $F(x_f,y_f)={x_f,y_f}$. That is, one in which $1/2+3/4 x_f+1/3 y_f=x_f$ and $2/5+1/2 x_f^2-4/5 x_f y_f=y_f$. This is equivalent to finding the roots of $V_1(x,y)=F(x,y)-\{x,y\}$.

The number of roots of $n$ polynomials in $n$ variables is, in general, given by Bezout's Theorem. The theorem states that in general, the (maximum) number of common zeros equals the product of the degrees of the polynomials. There are exceptions to this maximum number however. In the example above, both polynomials $p_1$ and $p_2$ have degree $2$. So that the maximum number of zeros to $V_1$ would be four.

It's easy to find four solutions in Mathematica via

Mathematica code 1

  NSolve[{1/2 + 3/4 x + 1/3 y^2 - x == 0,2/5 + 1/2 x^2 - x y - y == 0}, {x, y}]
  
    

What happens if we fold $F(x,y)$? That is, the composition $F\circ F=F^{((2))}=F(F(x,y))$? We then have the system $$ \begin{align*} F^{((2))}(x,y)&=\left\{\frac{1}{3} \left(\frac{x^2}{2}-x y+\frac{2}{5}\right)^2+\frac{3}{4} \left(\frac{3 x}{4}+\frac{y^2}{3}+\frac{1}{2}\right)+\frac{1}{2}\right.,\\ &-\left.\left(\frac{x^2}{2}-x y+\frac{2}{5}\right) \left(\frac{3 x}{4}+\frac{y^2}{3}+\frac{1}{2}\right)+\frac{1}{2} \left(\frac{3 x}{4}+\frac{y^2}{3}+\frac{1}{2}\right)^2+\frac{2}{5}\right\} \end{align*} $$ Folding a system of $2$ polynomials in $2$ variables each with degree $d_1$ and $d_2$ a total of $n$ times increases the degree by $(d_1 d_2)^n$. So that the maximum number of zeros of $V_2$ for this example would be $4^2$ or $16$. An easy way to visualize this is to let $u=1/2+3/4 x+1/3 y^2$ and $v=2/5+1/2 x^2-4 x y$ then $F^{((2))}(x,y)=F(u,v)$. Zeros of $V_2$ are easily computed in Mathematica. However, this method has computational limitations when $n$ becomes greater than about $10$ for the quadratic case in two variables and more limited when either the degree of the polynomials or the number of variables is increased.

The central question of this section then is how do we study the fixed-points of $$ F^{((n))}(X)=\underbrace{F(F(F(\cdots F))\cdots)}_{\text{n times}} $$ with $X=\{x_1,x_2,\cdots,x_k\}$ when $n$ is large say $n=100$ and the polynomials are higher than degree $2$ with more than two variables?

A simple system:

We can start by computing the roots of relatively accessible systems. Consider the system: $$ F(z,w)=\left\{\begin{array}{l} -1-1/4 w+z^2 \\ z \end{array} \right\} $$

The following is a timing report to compute the fixed points of $F^{((n))}$ for $n=1$ through $10$ with $25$ digits of precision using Mathematica's $\texttt{NSolve}$ running on a $4.5$ GHz quad-core Intel Zeon machine: $$ \begin{array}{c} \begin{array}{|c|c|} \hline \text{n} & \text{Time (min)} \\ \hline 1 & 0.00110648 \\ 2 & 0.000122265 \\ 3 & 0.00119263 \\ 4 & 0.00112888 \\ 5 & 0.00820462 \\ 6 & 0.0467463 \\ 7 & 0.340628 \\ 8 & 3.32432 \\ 9 & 33.0478 \\ 10 & 394.904 \\ \hline \end{array}\\ \text{Table 1: Timing results for computing roots of $V^{((1))}$ through $V^{((10))}$} \end{array} $$ The first thing to notice about the table is that the fixed points of $F^{((10)}$ took about $6.5$ hours and was about $10$ times the duration of $F^{((9))}$ which was about $10$ times the duration of $F^{((8))}$ and likewise so of $F^{((7))}$. So that a rough prediction of the time needed to compute the fixed points for $F^{((11))}$ would be a little more than $2.5$ days. With this trend, calculation of fixed points of $F^{((12))}$ using $\texttt{NSolve}$ becomes relatively-inaccessible at this speed considering the function uses Groebner elimination algorithms which cannot be parallelized. We need another way to access the roots of higher-folded systems.

Newton Iteration:

One way to find roots of $V_n(z,w)$ is by iteration. The simplest is Newton iteration. Mathematica implements this method in $\texttt{FindRoot}$. But how should we implement the method in complex $\{z,w\}$-space? For example, since this space is four-dimensional, how to decide which seeds to use for the iteration? A simple method is to search the $\{z,w=z\}$-space keeping in mind a seed such as $\{a,a\}$ may be far-removed from any root $\{p,q\}$ and so not lie in any root's basin of attraction,i.e., the method will diverge. Take for example, Newton iteration of $V_4(z,w)$ for the simple example above using seed $\{10+10i,10+10i\}$ and the following (default) root search using $\texttt{FindRoot}$:

Mathematica code 2

 g0[{z_, w_}] = -1 - 1/4 w + z^2;

h0[{z_, w_}] = z;

vectorF[{z_, w_}] = {g0[{z, w}], h0[{z, w}]};

cycleNum = 4;

foldedF[z_, w_] = Nest[vectorF, {z, w}, cycleNum] - {z, w};

theSeed = {10 + 10. I, 10 + 10. I};

FindRoot[foldedF[z, w], {{z, theSeed[[1]]}, {w, theSeed[[2]]}}]

Mathematica returns the following error:

Mathematica output 2

   FindRoot::jsing: Encountered a singular Jacobian at the point {z,w} = {10. +10. I,10. +10. I}. Try perturbing the initial point(s).

$\texttt{FindRoot}$ can encounter a number of other errors. Four will concern us for this study:

FindRoot::lstol: insufficient decrease in merit function
FindRoot::jsing: Jacobian is singular
FindRoot::cvmit: maximum iterations reached without achieving desired accuracy
FindRoot::sszero: step size has become less than tolerance of precision goal

Mathematica can identify and mark these error cases. Suppose we search a $20\times 20$ region in the complex $\{z,z\}$-space centered at $\{0,0\}$ with a resolution of say $\frac{1}{5}$. This lattice of points equates to a $100\times 100$ array with boundaries $10-10i$ to $10+10i$ or a total of $10,000$ data points, and each can serve as a seed for the Newton iteration of $V_4$. We can then code Mathematica to return an error code if any of the above errors are encountered during the iteration. The code below returns a standard $\{\text{NaN},\text{NaN}\}$ (for "not a number") code if one of the four errors are encountered, otherwise returning a root $\{p,q\}$ of $V_4(z,w)$ with a default machine-precision of about $15$ digits.

Coding a Newton iteration of $V_4(z,w)$:

Mathematica code 3

  (* 

  construct 2x2 array of points in region 10-10i, 10+10i centered at origin w/ res=1/5 

  *) 

   testTable =Array[{(#2 + I #1 ), (#2 + I #1)} &, {100,100}, {{10, -10}, {-10, 10}}];

  (* 

  flatten to 1-D array for parallelizing via ParallelTable function 

  *) 

   theSeedList = Flatten[testTable, 1];

We then set up parallel processing using $\texttt{ParallelTable}$ to iterate $V_4$ over the region returning either an error code or a root of $V_4$. $\texttt{ParallelTable}$ basically distributes the $10,000$ seeds equally among available cores on a multi-core machine thus reducing computation time.

Mathematica code 4

cycleNum = 4;
theV[z_, w_] = Nest[vectorF, {z, w}, cycleNum] - {z, w};
newSeeds = SetPrecision[theSeedList, 50];

PrintTemporary["Analyzing seed ", Dynamic@theIndex, "/", totalPoints];

timing1 = AbsoluteTiming[
  iterationRecords = ParallelTable[
     If[Mod[rootNum, 100] == 0,
      theIndex = rootNum;
      ];
     theSeed = newSeeds[[rootNum]];
     Quiet[
      Check[
        {z, w} /. 
        FindRoot[theV[z, w], {{z, theSeed[[1]]}, {w, theSeed[[2]]}}],
          {"NaN", "NaN"},
          {FindRoot::lstol, (* insufficient decrease in merit function *)
           FindRoot::jsing,  (* Jacobian is singular *)
           FindRoot::cvmit,  (* maximum iterations reached without achieving desired accuracy *)
           FindRoot::sszero  (*step size has become less than tolerance of precision goal *)
        }],
        {FindRoot::lstol, FindRoot::jsing, FindRoot::cvmit,FindRoot::sszero}],
     {rootNum, 1, Length@theSeedList}];
  ]

This code (running in parallel) took about $8$ seconds to compute $1084$ convergent points in the array. The plot below is a basin plot of the $z$-plane for the seed array: Black corresponds to seeds $\{a,a\}$ in which the iteration encountered one of the error cases above. The color table in the plot identifies seed regions which converged to a root. For example, the first color in the color table corresponds to regions converging to the first root of $V_4$ in a default sorted list or $\{-1.53 - 0.0642003 I, -0.0104631 + 0.487218 I\}$. The light green color section surrounding the origin, according to the color table, are seeds converging to the the seventh root or $\{-0.5542,-0.5542\}$ or more explicitly: $$ F^{((4))}(-0.5542,-0.5542)\approx \{-0.5542,-0.5542\} $$ The iteration in this case found all $16$ roots of $V_4$.

Iterating $V_{10}$:

And at this point it's instructive to create folded root lists by folding $F(z,w)$ ten times directly using the first and fourth root in the list above. Note in each list below, the first number is the root followed by a list of sequentially folding $F(z,w)$ ten times with the root: $$ \begin{array}{c} \begin{array}{ccc} \begin{array}{|c|c|} \text{z} & \text{w} \\ \hline -1.20643-0.189987 i & -0.213772+0.692465 i \\ 0.47283\, +0.285298 i & -1.20643-0.189987 i \\ -0.556218+0.317292 i & 0.47283\, +0.285298 i \\ -0.909503-0.424291 i & -0.556218+0.317292 i \\ -0.213772+0.692465 i & -0.909503-0.424291 i \\ -1.20643-0.189987 i & -0.213772+0.692465 i \\ 0.47283\, +0.285298 i & -1.20643-0.189987 i \\ -0.556218+0.317292 i & 0.47283\, +0.285298 i \\ -0.909503-0.424291 i & -0.556218+0.317292 i \\ -0.213772+0.692465 i & -0.909503-0.424291 i \\ -1.20643-0.189987 i & -0.213772+0.692465 i \\ \end{array} & \text{ } & \begin{array}{|c|c|} \text{z} & \text{w} \\ \hline -0.102735-0.872071 i & -1.03771+0.345885 i \\ -1.49052+0.0927125 i & -0.102735-0.872071 i \\ 1.23875\, -0.058363 i & -1.49052+0.0927125 i \\ 0.903731\, -0.167773 i & 1.23875\, -0.058363 i \\ -0.521106-0.288652 i & 0.903731\, -0.167773 i \\ -1.0377+0.34278 i & -0.521106-0.288652 i \\ 0.0896016\, -0.639242 i & -1.0377+0.34278 i \\ -1.14118-0.200249 i & 0.0896016\, -0.639242 i \\ 0.239785\, +0.61685 i & -1.14118-0.200249 i \\ -1.03771+0.345885 i & 0.239785\, +0.61685 i \\ -0.102735-0.872071 i & -1.03771+0.345885 i \\ \end{array} \end{array} \\ \text{Table 3: Folding roots one and four} \end{array} $$

The first thing to note about both lists above is the value of $F(z,w)$ cycles after $10$ foldings. However there is a more important feature of both lists. In the first case the root cycles after folding five times (note the sixth entry in the first list), and then cycles again after the tenth folding. And in the second, the root cycles only after folding ten times. The first case is a cycle-$5$ root and the other is a cycle-$10$ root. And it's important to note that the list of nine roots actually represents $71$ distinct roots of $V_{10}$ since each distinct value in a folded root list is also a root of $V_{10}$.

Cycles of iterated polynomial systems:

In the folded root lists, we encountered cycle-$5$ and cycle-$10$ roots of $V_{10}$. This makes sense since if $\{p,q\}$ is a fixed point of $F^{((r))}$ then it is likewise a fixed point of $F^{((nr))}$ where $n$ is a positive integer. This along with Bezout's Theorem allows us to compute the number of fixed points for each cycle of $F^{((n))}$. If we let $C_i$ equal to the list of cycle-$i$ points and $N(C_i)$ the number of cycle-$i$ points then we may write for a polynomial system of total degree $2$ as in the above example: $$ \begin{align*} N(C_1)&=2 \\ N(C_2)&=2^2-N(C_1) \\ N(C_3)&=2^3-N(C_1) \\ N(C_4)&=2^4-\left(N(C_1)+N(C_2)\right)\\ N(C_5)&=2^5-N(C_1)\\ N(C_6)&=2^6-\left(N(C_1)+N(C_2)+N(C_3)\right)\\ \vdots \\ N(C_n)&=2^n-\sum_{i\vert n}N(C_i). \end{align*} $$ where the notation $i\vert n$ means the positive divisors of $n$ with $i\lt n$. For example: $$ N(C_{10})=2^{10}-\left(N(C_1)+N(C_2)+N(C_5)\right)=990, $$ and these are the cycles experimentally found by analyzing the roots of $V_{10}$ computed for the timing study above: $$ \begin{array}{c} \begin{array}{|c|c|} \hline \text{Cycle Type} & \text{Size} \\ \hline 1 & 2 \\ 2 & 2 \\ 5 & 30 \\ 10 & 990 \\ \hline \end{array} \\ \text{Table 4: Cycles of $F^{((10))}$} \\ \end{array} $$

Algebraic functions, iterated exponentials and polynomial systems

Menu

Fixed points of folded polynomial systems, part 1

A challenge for the interested reader:

Background:

Mathematica code 1

A simple system:

Newton Iteration:

Mathematica code 2

Mathematica output 2

Coding a Newton iteration of $V_4(z,w)$:

Mathematica code 3

Mathematica code 4

Iterating $V_{10}$:

Cycles of iterated polynomial systems:

No comments:

Post a Comment

Blog Archive