A multiple double is an unevaluated sum of doubles. An NVIDIA tensor core is a specialized high performance compute core for matrix multiplication. The Ampere A100, released in 2020, introduced tensor cores capable of 64-bit floating-point arithmetic. Every multiple double arithmetical operation requires renormalization, which involves branching, for which tensor cores are unsuited. To solve this problem caused by renormalization, we apply a solution similar to the Ozaki scheme [Ozaki et al, Numerical Algorithms, 2012]. Our software is available under the GPU GPL license on github.

In many applications, the precision by the available hardware arithmetic is insufficient to guarantee accurate results. Multiword arithmetic is a special type of multiprecision arithmetic where a multiple double is an unevaluated sum of 64-bit doubles, or where a multiple integer is an unevaluated sum of 64-bit integers. Parallel computing is applied to compensate for the cost overhead of multiword arithmetic. This type of arithmetic exploits naturally the optimized hardware, allows for efficient type conversions, memory layouts, all favorable for parallel computing. For example, storing a multiword in registers rather than arrays is beneficial to parallel computing by tasking and acceleration by graphics processing units. Code for multiword arithmetic is available in the software PHCpack, written mainly in Ada, publicly available at github, and as an Alire crate, released under the GNU GPL v3.0 license.

Neural networks are important tools in machine learning. Representing piecewise linear activation functions with tropical arithmetic enables the application of tropical geometry. Algorithms are presented to compute regions where the neural networks are linear maps. Through computational experiments, we provide insights on the difficulty to train neural networks, in particular on the problems of overfitting and on the benefits of skip connections.

A polynomial homotopy is a family of polynomial systems in one parameter, which defines solution paths starting from known solutions and ending at solutions of a system that has to be solved. We consider paths leading to isolated singular solutions, to which the Taylor series converges logarithmically. Whether or not extrapolation algorithms manage to accelerate the slowly converging series depends on the proximity of poles close to the disk of convergence of the Taylor series.

A polynomial homotopy is a family of polynomial systems, typically in one parameter t. Our problem is to compute power series expansions of the coordinates of the solutions in the parameter t, accurately, using multiple double arithmetic. One application of this problem is the location of the nearest singular solution in a polynomial homotopy, via the theorem of Fabry. Power series serve as input to construct Padé approximations. Exploiting the massive parallelism of Graphics Processing Units capable of performing several trillions floating-point operations per second, the objective is to compensate for the cost overhead caused by arithmetic with power series in multiple double precision. The application of Newton's method for this problem requires the evaluation and differentiation of polynomials, followed by solving a blocked lower triangular linear system. Experimental results are obtained on NVIDIA GPUs, in particular the RTX 2080, RTX 4080, P100, V100, and A100. Code generated by the CAMPARY software is used to obtain results in double double, quad double, and octo double precision. The programs in this study are self contained, available in a public github repository under the GPL-v3.0 License.

Convolutional neural networks (CNNs) are a popular choice of model for tasks in computer vision. When CNNs are made with many layers, resulting in a deep neural network, skip connections may be added to create an easier gradient optimization problem while retaining model expressiveness. In this paper, we show that arbitrarily complex, trained, linear CNNs with skip connections can be simplified into a single-layer model, resulting in greatly reduced computational requirements during prediction time. We also present a method for training nonlinear models with skip connections that are gradually removed throughout training, giving the benefits of skip connections without requiring computational overhead during prediction time. These results are demonstrated with computational examples on Residual Networks (ResNet) architecture.

PHCpack is a software package for polynomial homotopy continuation, which provides a robust path tracker [Telen, Van Barel, Verschelde, SISC 2020]. This tracker computes the radius of convergence of Newton's method, estimates the distance to the nearest path, and then applies Padé approximants to predict the next point on the path. A priori step size control is less sensitive to finely tuned tolerances than a posteriori step size control, and is therefore robust. The Python interface phcpy is extended with a new step-by-step tracker and is applied to experiment with extrapolation methods to accurately locate the singular points at the end of solution paths.

A polynomial homotopy is a family of polynomial systems, where the systems in the family depend on one parameter. If for one value of the parameter we know a regular solution, then what is the nearest value of the parameter for which the solution in the polynomial homotopy is singular? For this problem we apply the ratio theorem of Fabry. Richardson extrapolation is effective to accelerate the convergence of the ratios of the coefficients of the series expansions of the solution paths defined by the homotopy. For numerical stability, we recondition the homotopy. To compute the coefficients of the series we propose the quaternion Fourier transform. We locate the closest singularity computing at a regular solution, avoiding numerical difficulties near a singularity.

This paper describes the application of the code generated by the CAMPARY software to accelerate the solving of linear systems in the least squares sense on Graphics Processing Units (GPUs), in double double, quad double, and octo double precision. The goal is to use accelerators to offset the cost overhead caused by multiple double precision arithmetic. For the blocked Householder QR and the back substitution, of interest are those dimensions at which teraflop performance is attained. The other interesting question is the cost overhead factor that appears each time the precision is doubled. Experimental results are reported on five different NVIDIA GPUs, with a particular focus on the P100 and the V100, both capable of teraflop performance. Thanks to the high Compute to Global Memory Access (CGMA) ratios of multiple double arithmetic, teraflop performance is already attained running the double double QR on 1,024-by-1,024 matrices, both on the P100 and the V100. For the back substitution, the dimension of the upper triangular system must be as high as 17,920 to reach one teraflops on the V100, in quad double precision, and then taking only the times spent by the kernels into account. The lower performance of the back substitution in small dimensions does not prevent teraflop performance of the solver at dimension 1,024, as the time for the QR decomposition dominates. In doubling the precision from double double to quad double and from quad double to octo double, the observed cost overhead factors are lower than the factors predicted by the arithmetical operation counts. This observation correlates with the increased performance for increased precision, which can again be explained by the high CGMA ratios.

The objective is to demonstrate the making of Ada software available to Python and Julia programmers using GPRbuild. GPRbuild is the project manager of the GNAT toolchain. With GPRbuild the making of shared object files is fully automated and the software can be readily used in Python and Julia. The application is the build process of PHCpack, a free and open source software package to solve polynomial systems by homotopy continuation methods, written mainly in Ada, with components in C++, available at github at https://github.com/janverschelde/PHCpack.

The Macaulay2 package NumericalSchubertCalculus provides methods for the numerical computation of Schubert problems on Grassmannians. It implements both the Pieri homotopy algorithm and the Littlewood-Richardson homotopy algorithm. Each algorithm has two independent implementations in this package. One is in the scripting language of Macaulay2 using the package NumericalAlgebraicGeometry, and the other is in the compiled code of PHCpack.

The problem is to evaluate a polynomial in several variables and its gradient at a power series truncated to some finite degree with multiple double precision arithmetic. To compensate for the cost overhead of multiple double precision and power series arithmetic, data parallel algorithms for general purpose graphics processing units are presented. The reverse mode of algorithmic differentiation is organized into a massively parallel computation of many convolutions and additions of truncated power series. Experimental results demonstrate that teraflop performance is obtained in deca double precision with power series truncated at degree 152. The algorithms scale well for increasing precision and increasing degrees.

We develop the Littlewood-Richardson homotopy algorithm, which uses numerical continuation to compute solutions to Schubert problems on Grassmannians and is based on the geometric Littlewood-Richardson rule. One key ingredient of this algorithm is our new optimal formulation of Schubert problems in local Stiefel coordinates as systems of equations. Our implementation can solve problem instances with tens of thousands of solutions.

We consider the problem of tracking one solution path defined by a polynomial homotopy on a parallel shared memory computer. Our robust path tracker applies Newton's method on power series to locate the closest singular parameter value. On top of that, it computes singular values of the Hessians of the polynomials in the homotopy to estimate the distance to the nearest different path. Together, these estimates are used to compute an appropriate adaptive step size. For n-dimensional problems, the cost overhead of our robust path tracker is O(n), compared to the commonly used predictor-corrector methods. This cost overhead can be reduced by a multithreaded program on a parallel shared memory computer.

Hardware double precision is often insufficient to solve large scientific problems accurately. Computing in higher precision defined by software causes significant computational overhead. The application of parallel algorithms compensates for this overhead. Newton's method to develop power series expansions of algebraic space curves is the use case for this application.

We propose a new algorithm for numerical path tracking in polynomial homotopycontinuation. The algorithm is “robust” in the sense that it is designed toprevent path jumping, and in many cases it can be used in (only) doubleprecision arithmetic. It is based on an adaptive stepsize predictor that usesPadé techniques to detect local difficulties for function approximation anddanger for path jumping. We show the potential of the new path trackingalgorithm through several numerical examples and compare it with existingimplementations.

Polynomial homotopies define solution paths of polynomial systems. By random complex constants in the homotopy, the paths are free of singular solutions and do not diverge, except perhaps at the end when the target system has singular solutions or solutions at infinity. Numerical continuation methods apply adaptive steplength control algorithms [7] to track the solution paths. Path jumping (identified in [9] as a reliability issue) occurs when the numerical approximations of one path jump onto another path. Although the paths themselves are free from singularities, path jumping occurs most likely when paths come near singularities. We developed a new adaptive steplength control algorithm to track the solution paths using estimates of the distance to the nearest value of the continuation parameter where singularities or solutions at infinity occur. Our point of departure is Fabry’s ratio theorem [6]. With this theorem we can detect singular points based on the coefficients of the Taylor series [5]. In particular, if for the series

The solutions of a system of polynomials in several variables are often needed, e.g.: in the design of mechanical systems, and in phase-space analyses of nonlinear biological dynamics. Reliable, accurate, and comprehensive numerical solutions are available through PHCpack, a FOSS package for solving polynomial systems with homotopy continuation. This paper explores new developments in phcpy, a scripting interface for PHCpack, over the past five years. For instance, phcpy is now available online through a JupyterHub server featuring Python2, Python3, and SageMath kernels. As small systems are solved in real-time by phcpy, they are suitable for interactive exploration through the notebook interface. Meanwhile, phcpy supports GPU parallelization, improving the speed and quality of solutions to much larger polynomial systems. From various model design and analysis problems in STEM, certain classes of polynomial system frequently arise, to which phcpy is well-suited.

A numerical irreducible decomposition for a polynomial system provides representations for the irreducible factors of all positive dimensional solution sets of the system, separated from its isolated solutions. Homotopy continuation methods are applied to compute a numerical irreducible decomposition. Load balancing and pipelining are techniques in a parallel implementation on a computer with multicore processors. The application of the parallel algorithms is illustrated on solving the cyclic $n$-roots problems, in particular for $n = 8, 9$, and~12.

We consider the extension of the method of Gauss-Newton from complex floating-point arithmetic to the field of truncated power series with complex floating-point coefficients. With linearization we formulate a linear system where the coefficient matrix is a series with matrix coefficients, and provide a characterization for when the matrix series is regular based on the algebraic variety of an augmented system. The structure of the linear system leads to a block triangular system. In the regular case, solving the linear system is equivalent to solving a Hermite interpolation problem. We show that this solution has cost cubic in the problem size. In general, at singular points, we rely on methods of tropical algebraic geometry to compute Puiseux series. With a few illustrative examples, we demonstrate the application to polynomial homotopy continuation.