Readers are encouraged to read Fixed points of folded polynomial systems, part 1 first.
Background:
Better seeds
Finding roots of $V_n(z,w)$ with $V_{n-1}(z,w)$ seeds, Case 1:
Roots of $V_n(z,w)$, Case 2:
Next consider iterating $V_{10}$ with greater working precision and larger iteration limit. In the code below we set the precision of the seeds to a working precision of $100$ digits and increase the maximum number of iterations to $500$:
Improving Newton's method, Case 3: Step size control
In the first part of this section, we studied fixed points of a folded system $$ F^{((n))}(z,w)=\underbrace{F(F(F(\cdots F))\cdots)}_{\text{n times}} $$ where $$ \begin{equation} F(z,w)=\left\{\begin{array}{l} 1-1/4 w+z^2 \\ z \end{array} \right\} \end{equation} $$ or equivalently the roots of $V_n(z,w)=F^{((n))}-\{z,w\}$. And we found the time to compute the fixed points of $F^{((10))}$ using $\texttt{NSolve}$ was about $6.5$ hours on a $4.5$ GHz quad-core Intel Zeon machine. And based on the trend of computing $F^{((7))}$ through $F^{((10))}$, estimated the time to compute $F^{((11))}$ to be about $2.5$ days and about $25$ days for $F^{((12))}$. We also developed a method to blanket-search a region in the $\{z,w=z\}$ plane using Netwon's method and Mathematica's $\texttt{FindRoot}$ function to compute cycle-$10$ roots but this approach produced a small number of roots for the test cases studied and would quickly become computationally expensive if the region were more finely partitioned. Clearly we would like a more efficient means of computing the fixed points to higher folded polynomial systems.
One reason for the low number of cycle-10 roots found by the blanket-search Newton iteration was the choice of seeds. Recall, a root is given by two numbers $\{p,q\}\in \mathbb{C}$. And the seeds were of the form $\{a,a\}$. However in general a root $\{p,q\}$ may be far-removed from any $\{a,a\}$ seed. For example, take the cycle-10 root $\{1.54409 -0.00374401 I,-1.59651-0.0411066 I\}$. A seed $\{a,a\}$ closest to this root would lie in a neighborhood of $\{0,0\}$, a distance of about $1.5$ from the root. And with the method used for the search, too far-removed from the root and so missed by the algorithm. Now, it is known a basin plot for a Newton iteration is often fractal which means there are many basins of attraction for reach root, but it is well-known the closer the seed is to the root, the greater chance the seed will lie in a basin of attraction for the root and so lead to a convergent sequence under Newton iteration. What we would like therefore, is to find a set of seeds that are as close to the roots as possible and thus increase the probability and therefore the number of seeds converging to a root.
Using system (1) above, suppose we plot, for the set of cycle-$n$ roots $\big\{\{z_1,w_1\},\{z_2,w_2\},\cdots,\{z_k,w_k\}\big\}$, the points $\big\{\{Re(z_1),Im(z_1)\},\{Re(z_2),Im(z_2)\},\cdots,\{Re(z_k),Im(z_k)\}\big\}$ for cycles $5$ through $10$. We obtain the plots shown in Figure 1. And in this particular example, similar plots are obtained for the $w$-terms.
There is an obvious trend in the figures: the roots appear to increasingly populate a well-defined constrained region. In particular, consider the plots for cycle-$9$ and cycle-$10$. Suppose for each cycle-$10$ root, we compute the nearest cycle-$9$ root and create a bar-chart of the minimum distance for each of the cycle-$10$ roots. The results are shown in Figure 2.
So that if we used the cycle-$9$ roots as seeds for the cycle-$10$ iteration, the maximum difference between a seed and the nearest root would be about $0.12$ and in most other cases, much less than this. How might that affect a Newton iteration of $V_{10}$ using as seeds the roots of $V_{9}$? And more importantly, could we use this to improve the estimated time of $2.5$ days to compute $V_{11}$ roots using $V_{10}$ roots as seeds? And of higher-folded systems?
Recall in Part 1 we computed the roots of $V_9$ to $25$ digits of accuracy using $\texttt{NSolve}$. We would like to use these roots as seeds to compute the roots of $V_{10}$. The following code implements $\texttt{FindRoot}$ with default settings to compute the roots of $V_{10}$ using $V_9$ roots as seeds.
Mathematica code 1
cycleNum = 10;
theV[z_, w_] = Nest[vectorF, {z, w}, cycleNum] - {z, w};
newSeeds=cycle9Roots;
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}];
]
Mathematica output 1
Total divergent seeds: 22 Total roots found: 490 Root precision: MachinePrecision Root accuracy: 15.6983 Total time: 0.915275
The code took about a second to find 490 roots with a default precision of 16 digits. Note only $22$ seeds did not converge to a root. The divergent seeds generated $\texttt{FindRoot::lstol}$ and $\texttt{FindRoot::cvmit}$ errors indicating insufficient working precision to achieve the default accuracy of $16$ digits using the default maximum number of iterations at $100$. If we now attempt to create folded root lists, the process of folding such low precision roots $10$ times through $F(z,w)$ reduces the precision down to $10$ which is below machine precision. An attempt to generate a cycle report of the resulting lists produces extremely poor results: $$ \begin{array}{c} \begin{array}{|c|c|c|c|} \hline \text{Cycle Type} & \text{Size} &\text{Roots Found} & \text{Difference} \\ \hline 1 & 2 & 261 & -259\\ 2 & 2 & 8 & -6\\ 5 & 30 & 154 & -124\\ 10 & 990 & 4499 & -3509\\ \hline \end{array} \\ \text{Table 1: Computed Cycles of $F^{((10))}$, Case 1} \\ \end{array} $$
The problem with this case of course is the lack of sufficient working precision and illustrates the need to fine-tune the $\texttt{FindRoot}$ algorithm.
Mathematica code 2
workingPrecision=100;
maxIterations=500;
theSeedList = SetPrecision[theCurrentRoots, workingPrecision];
timing1 = AbsoluteTiming[
iterationRecords = ParallelTable[
If[Mod[rootNum, 100] == 0,
theIndex = rootNum;
];
theSeed = theSeedList[[rootNum]];
Quiet[
Check[ {z, w} /.
FindRoot[theV[z, w], {{z, theSeed[[1]]}, {w, theSeed[[2]]}},
WorkingPrecision -> workingPrecision,
MaxIterations -> maxIterations],
{"NaN", "NaN"}, {FindRoot::lstol, FindRoot::jsing,
FindRoot::cvmit}], {FindRoot::lstol, FindRoot::jsing,
FindRoot::cvmit}],
{rootNum, 1, Length@theSeedList}];
];
Mathematica output 2
Total divergent seeds: 2 Total roots found: 510 Root precision: 100. Root accuracy: 99.7437 Total time: 3.23674
Note with a higher precision and increasing the maximum number of iterations reduced the number of divergent seeds to two. And more importantly, the roots were computed to an accuracy of of $99$ digits. Now when we fold the roots and compute the cycle report we obtain: $$ \begin{array}{c} \begin{array}{|c|c|c|c|} \hline \text{Cycle Type} & \text{Size} &\text{Roots Found} & \text{Difference} \\ \hline 1 & 2 & 2 & 0\\ 2 & 2 & 2 & 0\\ 5 & 30 & 30 & 0\\ 10 & 990 & 970 & 20\\ \hline \end{array} \\ \text{Table 2: Computed Cycles of $F^{((10))}$, Case 2} \\ \end{array} $$
So by increasing the precision and maximum iterations, and including the post-processing time to delete duplicate roots, fold the roots, and compute the cycle report, we were able to find all but $20$ roots of $V_{10}$ in about $4$ seconds using the cycle-$9$ roots as seeds.
A good step-size control of Newton's method will prevent oscillations about a root or escapes from the root. One method which can avoid these two possibilities is "Line search" method. In this method, the direction of a step is given by Newton's iteration but the size of the step is altered by internally solving a particular minimization problem which improves convergence. Background on this method can be found at the following Mathematica site: Unconstrained optimization: Step Control
The following code employs an explicit step-control method to search for the cycle-$10$ roots. Values for $\texttt{DecreaseFactor}$ and $\texttt{MaxRelativeStepSize}$ were determined heuristically.
Mathematica code 3
workingPrecision=100;
maxIterations=500;
theSeedList = SetPrecision[theCurrentRoots, workingPrecision];
timing1 = AbsoluteTiming[
iterationRecords = ParallelTable[
If[Mod[rootNum, 100] == 0,
theIndex = rootNum;
];
theSeed = theSeedList[[rootNum]];
Quiet[
Check[ {z, w} /.
FindRoot[theV[z, w], {{z, theSeed[[1]]}, {w, theSeed[[2]]}},
Method -> {"Newton",
"StepControl" -> {"LineSearch",
"CurvatureFactor" -> Automatic,
"DecreaseFactor" -> 1/10000,
"MaxRelativeStepSize" -> 1/2, Method -> Automatic}},
WorkingPrecision -> workingPrecision,
MaxIterations -> maxIterations], {"NaN",
"NaN"}, {FindRoot::lstol, FindRoot::jsing,
FindRoot::cvmit}], {FindRoot::lstol, FindRoot::jsing,
FindRoot::cvmit}],
{rootNum, 1, Length@theSeedList}];
];
Mathematica output 3
Total divergent seeds: 2 Total roots found: 510 Root precision: 100. Root accuracy: 99.7437 Total time: 3.04331
And in this case all roots were computed in about $4$ seconds. Note in particular the number of convergent seeds was the same as Case 2: $510$. However, in this case, there were seeds which converged to different roots than that of Case 2 and by folding these different roots, contributed the missing $20$ roots identified in Case 3. $$ \begin{array}{c} \begin{array}{|c|c|c|c|} \hline \text{Cycle Type} & \text{Size} &\text{Roots Found} & \text{Difference} \\ \hline 1 & 2 & 2 & 0\\ 2 & 2 & 2 & 0\\ 5 & 30 & 30 & 0\\ 10 & 990 & 990 & 0\\ \hline \end{array} \\ \text{Table 3: Final cycle report for roots of $V_{10}$} \\ \end{array} $$
So given the $33$ minutes to compute the cycle-$9$ roots reported in Part 1, only $4$ seconds of computation was needed to compute the cycle-$10$ roots. In Part 3, we use the techniques above and other methods to compute roots to higher-folded systems.
No comments:
Post a Comment