Convergence of matrix solution for base e

[jaydfox]

I've taken the time to analyze the rate of convergence of various solutions of truncations of the matrix method (Andrew's solution, essentially the Abel matrix solution).

I started with the accelerated version, mainly because I'm confident that it converges on the same solution as the natural solution, but fast enough to be useful for numerical analysis.

For all the graphs that follow, I found the accelerated solutions of the following matrix sizes: 16, 24, 32, 48, 64, 88, 128, 184, 256, 360, 512. Note that these are half-powers of 2, from 2^4 to 2^9, and the non-integer powers (e.g., 2^6.5=90.51) have been rounded to a multiple of 8.

So, first things first, here's a graph of the accumulated absolute error in the coefficients. Note that this will greatly exceed the maximum error for the slog of any point on the unit disc, so this isn't a measure of error per se, but it's a good indicator, because the maximum error should be within one or two orders of magnitude of the accumulated absolute error.

   

A few things to mention. The scale on the y-axis is in bits. The red lines are integer powers of 2, and the blue lines are the rounded half-powers. The red line at the top is for the 16x16 solution, and the red line at the bottom is for the 512x512 solution.

Essentially, each time you double the system size (double the rows and the columns, so yes, four times as many matrix elements), you get an additional 6.8 bits of precision, give or take.

I find it interesting, because 6.8 is very close to the periodicity of the terms of the singularities I use for the accelerated solution. (Actually, the periodicity is half that, about 3.36 or so, but if you look at doubling the number of matrix elements, then the relationship holds).

On top of that, the precision lost when solving the accelerated system is about 1.375 bits per row or column, and 1.375 is very close to 1.374557, the distance of the primary singularities from the origin. Coincidence? I'd like to know!

Anyway, to make it easier to see that the relationship is almost linear (1 bit of rows for 6.8 bits of precision), here's a chart with logarithmic scales on both axes. So, 4 bits of rows means 16 rows, and 9 bits means 512 rows.

   

[jaydfox]

I suppose I should explain how I got these. I compared each solution to a 1200x1200 accelerated solution. Assuming that we can expect 6.5-7 bits of precision to be gained with each doubling, there should be more than enough precision in the 1200x1200 solution for these graphs to be very precise. On the other hand, if I had gone all the way up to a 1024x1024 solution, there would only be about 1.5 bits difference, which would most likely cause a significant deviation from the linear relationship seen in the log-log comparision graph.

Later this week I hope to repeat this analysis with various natural solutions, which I should be able to take up 1024x1024 (perhaps even 1448x1448, depending on how much precision I can get out of the matrix solver). I expect to see an additional 2.5-4.0 bits for each doubling, based on initial observations I made a week or two ago.

Edit: changed "2.5-3.0 bits" to "2.5-4.0 bits", to be a bit more conservative on my estimate for the natural solution.

[andydude]

So let me get this straight. On the first graph you are plotting arg-vs-log(err) and on the second graph you are plotting log(size)-vs-log(err), is that right? I think that would mean the direct quantities would have a power relationship, if \( \log(E) = -6.8\log(n) \) then \( E = n^{-6.8} < n^{-6} \) which means you may have found an upper bound for the error.

Wow!

Andrew Robbins

[jaydfox]

Hmm, I'm not quite sure if we mean the same thing about the first graph. The first graph has, for the x-axis, the coefficient index "n" of the power series, and for the y-axis, the log of the cumulative absolute error for all coefficients from 1 to n. In other words, assuming the x-axis is linear in n, then the y-axis is:

\( \log_2\left(\sum_{k=1}^{n}\left|a^{'}_k-a_k\right|\right) \)

Here, \( a^{'}_k \) is the k-th coefficient of a particular solution of a truncated matrix, and \( a_k \) is the k-th coefficient of the "theoretical" infinite solution. (In fact, in these graphs, it's the 1200-term truncated solution.)

As you can see, for any particular power series of a finite truncation, the cumulative error rises rapidly for the first dozen or so coefficients, but tapers off quickly, approaching some value asymptotically. This is because the terms decrease exponentially (per the radius of convergence), so the absolute error in the later terms is negligible. This is assuming we evaluated at some point on the unit circle, and took the absolute value of each term in the power series, so this is only a useful bound within the unit disc.


BTW, if I were to do the same calculation, but assuming a radius of 4/3 (1.333...), then the total error would rise higher and continue to rise much further into the series. It would be interesting to see, as it would dictate a very conservative upper bound on error within the disc with radius 4/3. I could extend this logic up to, but not including, the radius of convergence, at which point, the total error for any finite truncation would in theory be infinite, making comparison...difficult .

As for the error relationship, I think yes, we can say we have found a very conservative upper bound for error, but only within the unit disc. I need to try this for various radii and see if we can find an upper bound as a function of n and radius. That would be quite nice!

[jaydfox]

I could extend this logic up to, but not including, the radius of convergence, at which point, the total error for any finite truncation would in theory be infinite, making comparison...difficult

After further reflection, I think the infinity issue is only a problem for the finite natural solutions. The accelerated solutions will not suffer from this problem.

I have two reasons for coming to this conclusion.

The first is conceptual: the infinity arises from the fact that there is a singularity in the natural solution. However, after deriving the power series for the "residue" as I call it, by subtracting the logarithms from the accelerated solution, we are left with a removable singularity. I.e., the value of the "residue" at a primary fixed point is the difference of two infinities, so it is technically a singularity, but it has a distinct finite value in the limit, much like \( \frac{sin(z)}{z} \) evaluated at 0.

As such, all values of the residue along the radius of convergence are finite.

The second reason for me to come to this conclusion is numerical: testing with my 1200-term accelerated solution indicates that, even at the radius of convergence, the magnitude of the terms of the power series decreases as we go further and further into the series. The decrease is not quite exponential, but it at least appears that the error is bounded asymptotically on the radius of convergence. Outside the radius, I would of course expect divergence, but it would take at least tens of thousands of terms to notice it appreciably for points very close to the radius, and I only have 1200.

Anyway, with all this said, I think I can calculate the cumulative absolute error at the radius of convergence, and assuming I see a similar linear plot on the log-log graph, we should be able to come up with a very conservative upper bound for the error on the entire domain of the power series. But only for accelerated solutions. So all this presupposes that the accelerated and natural solutions converge on the same infinite solution.

[jaydfox]

I threw together some graphs of the cumulative absolute error of the accelerated solution, at the radius of convergence:

The first graph is the accumulated error by coefficient index:

   

The bolder section of each line represents the coefficients actually calculated by a matrix solution, while the pale section represents the coefficients of the logarithms used to accelerate the solution. Notice that the error rises quite rapidly, then seems to get basically flat once we reach the end of the matrix solution. The flatness is due to the fact that the error has grown so large by the final term of the truncated matrix, that the remaining terms (from the power series for the logarithms) add negligible error, especially on this logarithmic scale. Thus, with larger systems, the initial error is very small, allowing for a slow and steady climb, which flattens off when we've reached peak error. I'll post more graphs later showing the term by term error.

The second graph is the plot of log(error) vs. log(size).

   

Once again, the plot appears linear, though we only get 3.6 additional bits of accuracy for each additional doubling of matrix size. With 36 bits for a 512x512 system, we're looking at about 40 bits for my 1200x1200 solution, perhaps as much as 42-45 bits in practice, because the errors in the coefficients will not all be oriented in the same direction (in the complex plane). So I don't have enough for double precision over the entire domain yet. But excluding values very close to the singularities, the accuracy should be enough for double precision math. (I'll need to investigate what regions can expect double precision, based on solution size.)

[jaydfox]

I have finally gotten around to analyzing the natural solutions. The 1024x1024 matrix solution took about 4 hours, with base_ring=RealField(2560) for the curious.

First I show the cumulative absolute error in the coefficients, which sets a conservative upper bound on error in the unit disk. The lines are the same as before, with red lines being solutions of systems with size being an integral power of 2. The blue lines are the half-powers, rounded to the nearest multiple of 8. Notice that there is an extra pair of lines, compared to the accelerated solutions I analyzed. They represent system sizes of 728 and 1024.

   

Next I show a plot of log(error) vs. log(size). The logarithms are binary logarithms, giving error in bits versus system size in bits.

   

And there you have it. The natural systems get about 3.5 bits of accuracy per doubling, roughly half that of the accelerated solutions. Therefore, I now feel confident in my original conjecture that my 1200x1200 accelerated solution is more accurate than a 1,000,000x1,000,000 natural solution.