GPU Accelerated Newton for Taylor Series Solutions
of Polynomial Homotopies in Multiple Double Precision^††thanks: Supported by the National Science Foundation under grant DMS 1854513.

The two main concerns are performance and convergence. For performance, the input must be sufficiently large, but still well conditioned enough to allow for Newton’s method to converge. In examining the scalability we address the first key question: how much of the overhead can be compensated by GPU acceleration? Prior work showed that teraflop performance was achieved in the acceleration of the convolutions to evaluate and differentiate polynomials in several variables at power series [42] and in the acceleration of the blocked Householder QR [43] to solve linear systems in the least squares sense in multiple double precision. The differentiation adapts the reverse mode of algorithmic differentiation [14] with power series arithmetic. The acceleration of the blocked Householder QR [8] is explained in [22], and also addressed in [6] and [45]. Additional related work on multicore and accelerated QR can be found in [5], [2], [25], [35], and [36]. The acceleration of the back substitution algorithm applies the formulas of [15], developed further using the ideas of [31].

The second question concerns the combination of various kernels in the linearization of the power series. Of the various different types of kernels that are launched, which types require the most amount of time?

A multiple double number is an unevaluated sum of nonoverlapping doubles. The renormalization and arithmetical operations are explained in [30].

MPLAPACK [28] supports quad double arithmetic and provides implementations of arbitrary multiprecision linear algebra operations. A recent application to matrix-matrix multiplication is in [39]. In [24], the 3M and 4M methods for complex multiplication are applied in linear algebra operations on complex matrices.

In this paper, for multiple double precision, the software libraries QDlib [17] and CAMPARY [20] are applied, customized as follows. Instead of working with an array of double double numbers, two arrays of doubles are used: the first array for the most significant doubles and the second one for the least significant doubles, as this memory layout benefits memory coalescing, especially for complex quad double and octo double numbers. The GPU version of the QDlib library [27] uses the double2 and double4 types of the CUDA SDK, which promote good memory access for double doubles and quad doubles, but are not longer suited for complex quad doubles or octo doubles.

Multiple double precision is not true multiprecision in the sense that one cannot select any number of bits for the fraction. The other drawback is the limited size of the exponents (limited to the 11 bits of the 64-bit hardware doubles), which will prohibit the computation with infinitesimal values. In the context of GPU acceleration, recent work of [19] makes an interesting comparison with double double arithmetic: “The double double arithmetic of CAMPARY performs best for the problem of matrix-vector multiplication.” Concerning quad double precision, the authors of [19] write: “the CAMPARY library is faster than our implementation; however as the precision increases the execution time of CAMPARY also increases significantly.” The advantage of multiple double arithmetic is that simple counts of the number of floating-point operations quantify the cost overhead precisely and the flops metrics for performance are directly applicable.

The goal of applying GPU acceleration to offset the cost of multiple double precision arithmetic is related to the recent trend of mixed precision [1], [18], [21].

Multiprecision is often deemed meaningless on problems with noisy input data, but the problems of interest are solution paths highly sensitive to small changes in the continuation parameter, corresponding to Taylor series with a small convergence radius. Reconditioning of the Taylor series is proposed in [44].

The convergence concern is closely related to the numerical conditioning of the problem of power series solutions. Consider the following classical result, applied in [33] to detect nearby singularities.

In a numerical interpretation of this theorem, observe that the smaller convergence radius, the larger the growth in the coefficients, because then $|c_{d+1}|>|c_{d}|$. For example if $|z|$ is 1/2, then for sufficiently large $d$, $|c_{d+1}|\approx 2|c_{d}|$. Thus, for series of order 64, we could observe coefficients of magnitude $2^{64}\approx 1.8\times 10^{19}$. Therefore, for convergence, it is best to generate examples which have a unique power series solution and use decaying coefficients as in the series developments of exponential functions.

In this paper Newton’s method is applied to compute Taylor series. An alternative is to apply Fourier series methods as done in analytic continuation, see e.g. [13]. These methods are very sensitive to a good choice of the step size for taking samples of the function to be differentiated. As explained in [13] a smaller step size benefits the lower order coefficients while a larger step size may be needed to compute the higher order coefficients accurately.

The main experimental result of this research is that on systems of 1024 equations in 1024 variables, when doubling the precision from quad double to octo double, the increase in wall clock time is much less than what can be predicted from the cost overhead factors of the multiple double arithmetical operations as the increase in the time occupied by the kernels is significant. As expected, with multiple double arithmetic, the arithmetic intensity of the computations increases and the GPU becomes fully occupied.

The next section contains a high level description of Newton’s method on power series, using linearizations. The setup of the test problems with the levels of precision is justified in the third section. As defined in the fourth section, Newton’s method is executed using a staggered progression of the order of the power series. To predict the performance of the accelerated code, in section five the arithmetic intensities of different kernels are computed. Section six contains the results of the computational experiments.

In the computation of the arithmetic intensity of convolutions, or equivalently, the number of floating-point computations per double, it suffices to consider one monomial. For example, take $n=4$ and let $f=x_{1}x_{2}x_{3}x_{4}$ be the monomial we evaluate and differentiate. Using the reverse mode of algorithmic differentiation, the computations are organized in three columns as follows:

where each $\star$ indicates a new multiplication. In the three columns we count respectively $n-1=3$, $n-2=2$, $n-2=2$, for a total of $(n-1)+(n-2)+(n-2)=3n-5$ multiplications, for $n$ inputs. If the inputs were doubles, then the arithmetic intensity would be

Each input is a power series of order $d$. To avoid thread divergence, the coefficients of the second series in each product are padded with zeros. For example, for $d=3$:

where coefficients with negative indices are zero. Ignoring the additions, we count $d^{2}$ multiplications, so we now have $(3n-5)d^{2}$ multiplications for $nd$ inputs. If the coefficients of the power series were doubles, then the arithmetic intensity would be

The coefficients of the power series are multiple doubles. For double double, quad double, and octo double, the number of doubles in the inputs are respectively $2nd$, $4nd$, and $8nd$. Working with complex coefficients doubles the size of the input. Doubling the precision doubles the size of the input, but increase the arithmetical cost significantly, as illustrated by Table 2 (with data from [41]).

Then the number of floating-point operations per double for the evaluation and differentiation of a product of $n$ power series of order $d$ are respectively for double doubles, quad doubles, and octo doubles:

where the corresponding multiplication factors $23/2$, $336/4$, and $1742/8$ evaluate respectively to 11.5, 84, and 217.75.

With each doubling of the precision, the arithmetic intensity increases by a factor of 11.5 (from double to double double), by a factor of $7.30\approx 84/11.5$ (from double double to quad double), and by a factor of $2.59\approx 217.75/84$ (from quad double to octo double). As the evaluation and differentiation is needed at every stage of Newton’s method, the high numbers of floating-point operations per double are promising indicators for the success of GPU acceleration.

The analysis in the previous section concerned real numbers. On complex data, the computational intensity increases, and moreover, offers more parallelism.

Writing the real and imaginary parts of ${\bf x}$ as ${\cal R}({\bf x})$ and ${\cal I}({\bf x})$, respectively, and applying the same notation to ${\bf y}$, the product of ${\bf x}$ with ${\bf y}$ is then executed as

All $\star$’s in (17) can be executed simultaneously first, with the $-$ and $+$ executed next. The formulas in (17) correspond to what is called the 4M method in [24].

The vectorization of the complex arithmetic introduces more parallelism. Alternatively, the computation can be arranged such that one block of threads simultaneously collaborate on the real and imaginary parts of one product. This alternative is beneficial for vector shorter than the typical warp size of 32 threads.

With the QR decomposition of $A$, solving $A{\bf x}={\bf b}$ in the least squares sense is reduced to $R{\bf x}=Q^{H}{\bf b}$, to the multiplication of $Q^{H}$ with ${\bf b}$, followed by a back substitution.

Assuming the leading coefficient vector ${\bf x}_{0}$ of the series has an accuracy of at least half the working precision, the QR decomposition happens only once in the first stage of Newton’s method, at a cost of $O(n^{3})$. While every stage involves the solution of $R{\bf x}=Q^{H}{\bf b}$, the cost of computing $Q^{H}{\bf b}$ and the back substitution is both $O(n^{2})$. Even if the acceleration of the QR decomposition works better than the acceleration of $R{\bf x}=Q^{H}{\bf b}$, the factor $n$ in the cost overhead of QR over $R{\bf x}=Q^{H}{\bf b}$ makes that the proportion of the QR decomposition will still dominate all $Q^{H}{\bf b}$ computations and all back substitutions, as $n$ equals 1024 and the number of stages is capped to 24.

The updates of the right hand side vectors involve many different matrices. For example, the updates to ${\bf b}_{3}$ happen as ${\bf b}_{3}:={\bf b}_{3}-A_{3}{\bf x}_{0}$, ${\bf b}_{3}:={\bf b}_{3}-A_{2}{\bf x}_{1}$, ${\bf b}_{3}:={\bf b}_{3}-A_{1}{\bf x}_{2}$, each time with different matrices $A_{3}$, $A_{2}$, and $A_{1}$ which cannot remain all in the main memory of the device. Even as the cost of these computations is $O(n^{2})$ we may expect the updates to occupy a significant portion of the total execution times.

The same arguments apply to the computation of the residuals, when measuring the accuracy of the computed updates $\Delta{\bf x}$. However, one could significantly reduce the cost by selecting only one or a few random equations instead of computing the residuals for all equations.

For comparison with the arithmetic intensities of the convolutions, consider the matrix-vector product, for an $n$-by-$n$ matrix. Executed on doubles, $n^{2}$ multiplications are performed on $(n+1)n$ doubles. The ratio $n^{2}/(n^{2}+n)$ improves on power series of order $d$, as convolutions with padding take $d^{2}$ multiplications, while the size of the data is multiplied by $d$. Restricting to multiplications, the arithmetic intensity then is

Compared to the arithmetic intensity of convolutions, the leading terms of the numerator of (15) divided by the denominator evaluates to $3d$. So the arithmetic intensity of convolutions is three times more than that of matrix-vector products.

The obvious storage of a multiple double is as an array of many doubles, stored at consecutive locations. But this storage is unfavorable for coalesced memory accesses when working with vectors of multiple doubles. Therefore, all arithmetical operations on multiple doubles are defined by functions which take as input arguments not arrays of doubles, but sequences of doubles. For example, the function to add two double doubles has prototype

     void ddf_add
      ( double a_hi, double a_lo,
        double b_hi, double b_lo,
        double *c_hi, double *c_lo );

where high doubles of the numbers $a$, $b$, and $c$ are respectively given as a_hi, b_hi, c_hi and their low doubles as a_lo, b_lo, c_lo, for the addition $c:=a+b$. The GPU counterpart to the ddf_add above has the same arguments but starts with

  __device__ __forceinline__  void ddg_add

forcing the inlining of the kernels. All kernels (in the .cu files) are compiled with nvcc -O3.

Every supported precision level has dedicated functions, once defined as regular C functions (for correctness comparisons) and once defined as kernel functions to be executed on the device.

A vector of double doubles is stored in main memory on the device as two vectors, one vector with the high doubles, and another vector with the low doubles. All threads in one block load the high and low doubles in vectors in memory shared by all threads in the block. As adjacent threads read into adjacent memory locations, this staging the multiple doubles benefits memory coalescing.

The inlined arithmetical kernels are applied to numbers in registers. Before executing the arithmetical operations, the data is copied from shared memory into local variables, variables that are local to each thread in a block. The additional instruction cycles allow to hide the latencies of the data transfers.

The code was developed for the “Volta” V100 NVIDIA GPU, and tested on the “Pascal” P100 and RTX 2080, 4080 NVIDIA GPUs. The “Ampere” A100 was added most recently. Table 3 lists the main characteristics of the GPUs.

The double precision peak performance of the P100 is 4.7 TFLOPS. At 7.8 TFLOPS, the V100 is 1.66 times faster than the P100. The theoretical peak performance of the A100 is 8.7 TFLOPS. To evaluate the algorithms, compare the ratios of the wall clock times on the P100 over V100 with the factor 1.66. For every kernel, the number of arithmetical operations is accumulated. The total number of double precision operations is computed using the cost overhead multipliers.

The double precision peak performance of the RTX 2080 Max Q-design is 201.5 GFLOPS, which corresponds to 761.5 GFLOPS for the more recent RTX 4080. Both the RTX 2080 and 4080 are housed in a laptop computer.

The operating system of the host of the P100 and V100 is CentOS 7. For the A100, the host runs Rocky Linux. The code on the host is compiled with g++ 4.8.5 (Red Hat), with optimization flag ’-O3’, using version 8.0 of the CUDA Software Development Kit. The RTX 2080 and RTX 4080 are housed in laptops running Windows 11. The code on the host of RTX 2080 is compiled with the 2019 Community Visual Studio compiler, using CUDA 10.2. On the host of RTX 4080, the compiler is the 2022 Community Visual Studio version, using CUDA 12.1. On Windows, the optimization flag for the compiler is ’-02’.

Table 4 lists six different types of kernels.

Visualizing the data in Table 5, Figure 1 shows the percentages of the kernels on one column of monomials defined by a triangular exponent matrix of dimension 1024 to compute series of order 64 in octo double precision, done on the V100.

The largest portion of the time goes to the residual computations, for all equations in the system. The residuals are important to measure the convergence and must be computed in multiple double precision. One optimization could be to select at random one or a couple of equations and compute the residuals for those selected equations instead of for all equations. Figure 2 shows that for a 2-column monomial system, the time spent on convolutions dominates.

Figure 1 shows that the convolutions occupy a substantial part. For what orders of the series do we observe teraflop performance? Consider Table 6.

On one column of monomials, triangular exponent matrix of ones, $n=1024$, performance of the evaluation and differentiation, in octo double precision, for increasing orders of the series, visualizing the data in Table 6, Figure 3 shows that teraflop performance is observed after order 40 on the V100. Teraflop is observed after order 40 on the V100, but on the P100 only half a teraflop is reached at order 64. On the A100, teraflop performance is observed at order 27.

In the implementation, in the convolution of two power series each thread is responsible for one coefficient of the result. Threads are launched in blocks of the size that match the number of coefficients. An implementation better suited for series of lower order would employ a finer granularity and have several threads collaborate to compute one coefficient of a convolution of two series.

As the V100 is theoretically 1.66 faster than the P100, multiplying 554.654 (on the last row of Table 6) by 1.66 yields 920.726, which is less than the observed 1,658.382 on the V100. Thus, for higher levels of precision (octo double), newer GPUs scale better than older ones. The observed 1.7 teraflop performance on the V100 is only at 22% of the theoretical peak performance. The performance improves on series of higher order.

To investigate how much of the cost overhead can be compensated by the acceleration consider the wall clock times and the elapsed times spend by all kernels when the precision is doubled.

Visualizing Table 7, Figure 4 shows the 2-logarithms of the times of 24 steps with Newton’s method on one column of monomials defined by a triangular exponent matrix of ones of dimension 1024, on the V100. Doubling the precision less than doubles the wall clock time and increases the time spent by all kernels. The reason for this is that problems become more computationally intensive: the amount of data doubles but the amount of arithmetical computations on average is multiplied by factors larger than two, as explained in section 5.1, see also Table 2.

On the problem of applying the block Householder QR on real matrices of octo double numbers, the RTX 4080 reaches a performance of 413.2 GFLOPS on dimension 1,024 (8 blocks of size 128), and 719.5 GFLOPS (very close to the theoretical peak performance) is reached on dimension 2,048 ($8~{}\times~{}256$).

The last experiments concern the RTX 2080, on 16 steps with Newton’s method on one column of monomials defined by a triangular exponent matrix of ones of dimension 512. The results are summarized in Figure 5, visualizing the data in Table 8.

Table 9 compares the performance of the newer RTX 4080 versus the older RTX 2080, illustrating the scalability, on one column of 512 monomials, on series of order 64. Although the RTX 4080 is at least twice as powerful as the RTX 2080, the order 64 is unfavorable for the RTX 4080 which has 128 cores per multiprocessor. Relative to the P100 and V100, Figure 5 looks better than the corresponding Figure 4 because the RTX 2080 is already better occupied, already at a column of 512 monomials.