Finite Precision Arithmetic Research Articles

On recurrences converging to the wrong limit in finite precision and some new examples

In 1989, Jean-Michel Muller gave a famous example of a recurrence where, for particular initial values, the iteration over real numbers converges to a repellent fixed point, whereas finite precision arithmetic produces a different result, the attracting fixed point. We analyze recurrences in that spirit and remove a gap in previous arguments in the literature, that is, the recursion must be well defined. The latter is known as the Skolem problem. We identify initial values producing a limit equal to the repellent fixed point, show that in every e-neighborhood of such initial values the recurrence is not well defined, and characterize initial values for which the recurrence is well defined. We give some new examples in that spirit. For example, the correct real result, i.e., the repellent fixed point, may be correctly computed in bfloat, half, single, double, formerly extended precision (80 bit format), binary128 as well as many formats of much higher precision. Rounding errors may be beneficial by introducing some regularizing effect.

Square-Root Approach for Chandrasekhar-Based Maximum Correntropy Kalman Filtering

The first two Chandrasekhar recursions for the maximum correntropy criterion (MCC) Kalman filter (KF) have been recently derived for constant discrete-time linear systems. Their key feature is a mathematical re-formulation of the underlying MCC-based Riccati-type difference equation in terms of the involved error covariance matrix increment. Thus, the Chandrasekhar recursion-based solution is proved to yield a significant reduction of the computational complexity. This letter discusses the existence of a stable square-root solution for Chandrasekhar-type MCC-KF estimators, i.e. their computational reliability issue in a finite precision arithmetic. Two square-root solutions are proposed in terms of covariance quantities, namely within the Cholesky factorization and singular value decomposition.

Architect: Arbitrary-Precision Hardware With Digit Elision for Efficient Iterative Compute

Many algorithms feature an iterative loop that converges to the result of interest. The numerical operations in such algorithms are generally implemented using finite-precision arithmetic, either fixed- or floating-point, most of which operate least-significant digit first. This results in a fundamental problem: if, after some time, the result has not converged, is this because we have not run the algorithm for enough iterations or because the arithmetic in some iterations was insufficiently precise? There is no easy way to answer this question, so users will often over-budget precision in the hope that the answer will always be to run for a few more iterations. We propose a fundamentally new approach: with the appropriate arithmetic able to generate results from most-significant digit first, we show that fixed compute-area hardware can be used to calculate an arbitrary number of algorithmic iterations to arbitrary precision, with both precision and approximant index increasing in lockstep. Consequently, datapaths constructed following our principles demonstrate efficiency over their traditional arithmetic equivalents where the latter’s precisions are either under- or over-budgeted for the computation of a result to a particular accuracy. Use of most-significant digit-first arithmetic additionally allows us to declare certain digits to be stable at runtime, avoiding their recalculation in subsequent iterations and thereby increasing performance and decreasing memory footprints. Versus arbitrary-precision iterative solvers without the optimizations we detail herein, we achieve up-to $16\times $ performance speedups and $1.9\times $ memory savings for the evaluated benchmarks.

Sound and robust solid modeling via exact real arithmetic and continuity

Algorithms for solid modeling, i.e., Computer-Aided Design (CAD) and computer graphics, are often specified on real numbers and then implemented with finite-precision arithmetic, such as floating-point. The result is that these implementations do not soundly compute the results that are expected from their specifications. We present a new library, StoneWorks, that provides sound and robust solid modeling primitives. We implement StoneWorks in MarshallB, a pure functional programming language for exact real arithmetic in which types denote topological spaces and functions denote continuous maps, ensuring that all programs are sound and robust. We developed MarshallB as an extension of the Marshall language. We also define a new shape representation, compact representation ( K-rep ), that enables constructions such as Minkowski sum and analyses such as Hausdorff distance that are not possible with traditional representations. K-rep is a nondeterminism monad for describing all the points in a shape. With our library, language, and representation together, we show that short StoneWorks programs can specify and execute sound and robust solid modeling algorithms and tasks.

Novel consistent approach in controllability evaluations of point reactor kinetics models

The controllability of linearized equations of the equivalent normalized version of the reactor point kinetics equations with six delayed neutron groups and without neutron source at an equilibrium power is analyzed. There exists inconsistency between the results from two kinds of approaches. The equations analyzed are controllable according to the controllability evaluation by analytical derivation while they are uncontrollable by numerical calculation. The reason of this inconsistency is that the large difference between the time effect of prompt neutrons and that of delayed neutrons in reactors reduces the controllability matrix to become a rank-deficient matrix under the conditions of using finite-precision arithmetic. In order to improve the analysis approach and address the divergence in outcomes, a novel definition of degree of controllability is proposed by means of applying matrix norm and singular value decomposition (SVD). This definition resolves the inconsistent conclusions, i.e., the equations analyzed are controllable but their degree of controllability is very low. The calculation of the degree of controllability shows that in the multiple representations of point reactor kinetics, the representations of few-group models are more propitious to analyze and design the control system of nuclear reactor than the representation of six-group model.

Factored-form Kalman-like implementations under maximum correntropy criterion

Abstract The maximum correntropy criterion (MCC) methodology is recognized to be a robust filtering strategy with respect to outliers and shown to outperform the classical Kalman filter (KF) for estimation accuracy in the presence of non-Gaussian noise. However, the numerical stability of the newly proposed MCC-KF estimators in finite precision arithmetic is seldom addressed. In this paper, a family of factored-form (square-root) algorithms is derived for the MCC-KF and its improved variant, respectively. The family traditionally consists of three factored-form implementations: (i) Cholesky factorization-based algorithms, (ii) modified Cholesky, i.e. UD-based methods, and (iii) the recently established SVD-based filtering. All these strategies are commonly recognized to enhance the numerical robustness of conventional filtering with respect to roundoff errors and, hence, they are the preferred implementations when solving applications with high reliability requirements. Previously, only Cholesky-based IMCC-KF algorithms have been designed. This paper enriches a factored-form family by introducing the UD- and SVD-based methods as well. A special attention is paid to array algorithms that are proved to be the most numerically stable and, additionally, suitable for parallel implementations. The theoretical properties are discussed and numerical comparison is presented for determining the most reliable implementations.

Computing primary solutions of equations involving primary matrix functions

The matrix equation f(X)=A, where f is an analytic function and A is a square matrix, is considered. Some results on the classification of solutions are provided. When f is rational, a numerical algorithm is proposed to compute all solutions that can be written as a polynomial of A. For real data, the algorithm yields the real solutions using only real arithmetic. Numerical experiments show that the algorithm performs in a stable fashion when run in finite precision arithmetic.

A stabilized normal form algorithm for generic systems of polynomial equations

We propose a numerical linear algebra based method to find the multiplication operators of the quotient ring C [ x ] ∕ I associated to a zero-dimensional ideal I generated by n C -polynomials in n variables. We assume that the polynomials are generic in the sense that the number of solutions in C n equals the Bézout number. The main contribution of this paper is an automated choice of basis for C [ x ] ∕ I , which is crucial for the feasibility of normal form methods in finite precision arithmetic. This choice is based on numerical linear algebra techniques and it depends on the given generators of I .

Risk of Hearing Loss Caused by Multiple Acoustic Impulses in the Framework of Biovariability

We consider the hearing loss injury among subjects in a crowd with a wide spectrum of individual intrinsic injury probabilities due to biovariability. For multiple acoustic impulses, the observed injury risk of a crowd vs the effective combined dose follows the logistic dose-response relation. The injury risk of a crowd is the average fraction of injured. The injury risk was measured in experiments as follows: each subject is individually exposed to a sequence of acoustic impulses of a given intensity and the injury is recorded; results of multiple individual subjects were assembled into data sets to mimic the response of a crowd. The effective combined dose was adjusted by varying the number of impulses in the sequence. The most prominent feature observed in experiments is that the injury risk of the crowd caused by multiple impulses is significantly less than the value predicted based on assumption that all impulses act independently in causing injury and all subjects in the crowd are statistically identical. Previously, in the case where all subjects are statistically identical (i.e., no biovariability), we interpreted the observed injury risk caused by multiple impulses in terms of the immunity effects of preceding impulses on subsequent impulses. In this study, we focus on the case where all sound exposure events act independently in causing injury regardless of whether one is preceded by another (i.e., no immunity effect). Instead, we explore the possibility of interpreting the observed logistic dose-response relation in the framework of biovariability of the crowd. Here biovariability means that subjects in the crowd have their own individual injury probabilities. That is, some subjects are biologically less or more susceptible to hearing loss injury than others. We derive analytically the distribution of individual injury probability that produces the observed logistic dose-response relation. For several parameter values, we prove that the derived distribution is mathematically a proper density function. We further study the asymptotic approximations for the density function and discuss their significance in practical numerical computation with finite precision arithmetic. Our mathematical analysis implies that the observed logistic dose-response relation can be theoretically explained in the framework of biovariability in the absence of immunity effect.

Exact Lookup Tables for the Evaluation of Trigonometric and Hyperbolic Functions

Elementary mathematical functions are pervasively used in many applications such as electronic calculators, computer simulations, or critical embedded systems. Their evaluation is always an approximation, which usually makes use of mathematical properties, precomputed tabulated values, and polynomial approximations. Each step generally combines error of approximation and error of evaluation on finite-precision arithmetic. When they are used, tabulated values generally embed rounding error inherent to the transcendence of elementary functions. In this article, we propose a general method to use error-free values that is worthy when two or more terms have to be tabulated in each table row. For the trigonometric and hyperbolic functions, we show that Pythagorean triples can lead to such tables in little time and memory usage. When targeting correct rounding in double precision for the same functions, we also show that this method saves memory and floating-point operations by up to 29 and 42 percent, respectively.

Stability and sensitivity of a cross-flow-dominated Falkner–Skan–Cooke boundary layer with discrete surface roughness

With the motivation of determining the critical roughness size, a global stability and sensitivity analysis of a three-dimensional Falkner–Skan–Cooke (FSC) boundary layer with a cylindrical surface roughness is performed. The roughness size is chosen such that breakdown to turbulence is initiated by a global version of traditional secondary instabilities of the cross-flow (CF) vortices instead of an immediate flow tripping at the roughness. The resulting global eigenvalue spectra of the systems are found to be very sensitive to numerical parameters and domain size. This sensitivity to numerical parameters is quantified using the $\unicode[STIX]{x1D700}$-pseudospectrum, and the dependency on the domain is analysed through an impulse response, structural sensitivity analysis and an energy budget. It is shown that while the frequencies remain relatively unchanged, the growth rates increase with domain size, which originates from the inclusion of stronger CF vortices in the baseflow. This is reflected in a change in the rate of advective energy transport by the baseflow. It is concluded that the onset of global instability in a FSC boundary layer as the roughness height is increased does not correspond to an immediate flow tripping behind the roughness, but occurs for lower roughness heights if sufficiently long domains are considered. However, the great sensitivity results in an inability to accurately pinpoint the exact parameter values for the bifurcation, and the large spatial growth of the disturbances in the long domains eventually becomes larger than can be resolved using finite-precision arithmetic.

Regularization of Ill-Posed Point Neuron Models

Point neuron models with a Heaviside firing rate function can be ill-posed. That is, the initial-condition-to-solution map might become discontinuous in finite time. If a Lipschitz continuous but steep firing rate function is employed, then standard ODE theory implies that such models are well-posed and can thus, approximately, be solved with finite precision arithmetic. We investigate whether the solution of this well-posed model converges to a solution of the ill-posed limit problem as the steepness parameter of the firing rate function tends to infinity. Our argument employs the Arzelà–Ascoli theorem and also yields the existence of a solution of the limit problem. However, we only obtain convergence of a subsequence of the regularized solutions. This is consistent with the fact that models with a Heaviside firing rate function can have several solutions, as we show. Our analysis assumes that the vector-valued limit function v, provided by the Arzelà–Ascoli theorem, is threshold simple: That is, the set containing the times when one or more of the component functions of v equal the threshold value for firing, has zero Lebesgue measure. If this assumption does not hold, we argue that the regularized solutions may not converge to a solution of the limit problem with a Heaviside firing function.

Towards a Compiler for Reals

Numerical software, common in scientific computing or embedded systems, inevitably uses a finite-precision approximation of the real arithmetic in which most algorithms are designed. In many applications, the roundoff errors introduced by finite-precision arithmetic are not the only source of inaccuracy, and measurement and other input errors further increase the uncertainty of the computed results. Adequate tools are needed to help users select suitable data types and evaluate the provided accuracy, especially for safety-critical applications. We present a source-to-source compiler called Rosa that takes as input a real-valued program with error specifications and synthesizes code over an appropriate floating-point or fixed-point data type. The main challenge of such a compiler is a fully automated, sound, and yet accurate-enough numerical error estimation. We introduce a unified technique for bounding roundoff errors from floating-point and fixed-point arithmetic of various precisions. The technique can handle nonlinear arithmetic, determine closed-form symbolic invariants for unbounded loops, and quantify the effects of discontinuities on numerical errors. We evaluate Rosa on a number of benchmarks from scientific computing and embedded systems and, comparing it to the state of the art in automated error estimation, show that it presents an interesting tradeoff between accuracy and performance.

Automatic source‐to‐source error compensation of floating‐point programs: code synthesis to optimize accuracy and time

SummaryNumerical programs with IEEE 754 floating‐point computations may suffer from inaccuracies because finite precision arithmetic is an approximation of real arithmetic. Solutions that reduce the loss of accuracy are available, such as compensated algorithms or double‐double precision floating‐point arithmetic. Our goal is to automatically improve the numerical quality of a numerical program with the smallest impact on its performance. We define and implement source code transformations in order to derive automatically compensated programs. We present several experimental results to compare the transformed programs and existing solutions. The transformed programs are as accurate and efficient as the implementations of compensated algorithms when the latter exist. Furthermore, we propose some transformation strategies allowing us to partially improve the accuracy of programs and to tune the impact on execution time. Trade‐offs between accuracy and performance are assured by code synthesis. Experimental results show that user‐defined trade‐offs are achievable in a reasonable amount of time, with the help of the tools we present in the paper. Copyright © 2016 John Wiley & Sons, Ltd.

Design Criteria for FIR-Based Echo Cancellers

The most relevant impairment experienced by on-channel repeaters (OCRs) in single frequency networks scenarios is the presence of a coupling-channel between the transmitting and receiving antennas, that generates unwanted echoes. This phenomenon causes a degradation of the repeated signal and, above all, could lead to the instability of OCRs, owing to the positive feedback that could result. For this reason, OCRs are usually equipped with echo canceller units, aimed at removing the coupling contributions. In this paper, we consider a realistic OCR setup and we analytically derive proper design criteria for echo-cancellers, showing the role of system parameters and implementation aspects on their performance. Here, in particular, we investigate the joint effect of the estimation noise and the finite precision arithmetic of digital systems, the system sensitivity to different design parameters and the relation between the echo-cancelling performance and the probability of instability.

Stable evaluation of Green's functions in cylindrically stratified regions with uniaxial anisotropic layers

We present a robust algorithm for the computation of electromagnetic fields radiated by point sources (Hertzian dipoles) in cylindrically stratified media where each layer may exhibit material properties (permittivity, permeability, and conductivity) with uniaxial anisotropy. Analytical expressions are obtained based on the spectral representation of the tensor Green's function based on cylindrical Bessel and Hankel eigenfunctions, and extended for layered uniaxial media. Due to the poor scaling of these eigenfunctions for extreme arguments and/or orders, direct numerical evaluation of such expressions can produce numerical instability, i.e., underflow, overflow, and/or round-off errors under finite precision arithmetic. To circumvent these problems, we develop a numerically stable formulation through suitable rescaling of various expressions involved in the computational chain, to yield a robust algorithm for all parameter ranges. Numerical results are presented to illustrate the robustness of the formulation including cases of practical interest.

Dispersion Compensation FIR Filter With Improved Robustness to Coefficient Quantization Errors

In this paper, we propose a new finite impulse response (FIR) filter for chromatic dispersion compensation which is given in closed form. We identify a relation between the out-of-band gain and the in-band error when the filter is implemented with finite-precision arithmetic. In particular, a large out-of-band gain makes the filter more sensitive to coefficient quantization errors due to finite precision digital signal processing. To improve robustness to coefficient quantization errors, our proposed filter is designed based on confining the out-of-band gain. By means of simulations, we show that our filter outperforms other existing FIR filters. The performance gain improves with increasing modulation order and decreasing number of bits used to represent the filter taps.

Variants of the groupwise update strategy for short-recurrence Krylov subspace methods

Krylov subspace methods often use short-recurrences for updating approximations and the corresponding residuals. In the bi-conjugate gradient (Bi-CG) type methods, rounding errors arising from the matrix---vector multiplications used in the recursion formulas influence the convergence speed and the maximum attainable accuracy of the approximate solutions. The strategy of a groupwise update has been proposed for improving the convergence of the Bi-CG type methods in finite-precision arithmetic. In the present paper, we analyze the influence of rounding errors on the convergence properties when using alternative recursion formulas, such as those used in the bi-conjugate residual (Bi-CR) method, which are different from those used in the Bi-CG type methods. We also propose variants of a groupwise update strategy for improving the convergence speed and the accuracy of the approximate solutions. Numerical experiments demonstrate the effectiveness of the proposed method.

Ill-Posed Point Neuron Models.

We show that point-neuron models with a Heaviside firing rate function can be ill posed. More specifically, the initial-condition-to-solution map might become discontinuous in finite time. Consequently, if finite precision arithmetic is used, then it is virtually impossible to guarantee the accurate numerical solution of such models. If a smooth firing rate function is employed, then standard ODE theory implies that point-neuron models are well posed. Nevertheless, in the steep firing rate regime, the problem may become close to ill posed, and the error amplification, in finite time, can be very large. This observation is illuminated by numerical experiments. We conclude that, if a steep firing rate function is employed, then minor round-off errors can have a devastating effect on simulations, unless proper error-control schemes are used.

Solving Generalized CDT Problems via Two-Parameter Eigenvalues

We consider solving a nonconvex quadratic minimization problem with two quadratic constraints, one of which being convex. This problem is a generalization of the Celis--Denis--Tapia (CDT) problem and thus we refer to it as GCDT (Generalized CDT). The CDT problem has been widely studied, but no polynomial-time algorithm was known until Bienstock's recent work. His algorithm solves the CDT problem in polynomial time with respect to the number of bits in data and $\log\epsilon^{-1}$ by admitting an $\epsilon$ error in the constraints. The algorithm, however, appears to be difficult to implement. In this paper, we present another algorithm for GCDT, which is guaranteed to find a global solution for almost all GCDT instances (and slightly perturbed ones in some exceptionally rare cases), in exact arithmetic (including eigenvalue computation). Our algorithm is based on the approach proposed by Iwata, Nakatsukasa, and Takeda (2015) for computing the signed distance between overlapping ellipsoids. Our algorithm computes all the Lagrange multipliers of GCDT by solving a two-parameter linear eigenvalue problem, obtains the corresponding KKT points, and finds a global solution as the KKT point with the smallest objective value. In practice, in finite precision arithmetic, our algorithm requires $O(n^6\log\log u^{-1})$ computational time, where $n$ is the number of variables and $u$ is the unit roundoff. Although we derive our algorithm under the unrealistic assumption that exact eigenvalues can be computed, numerical experiments illustrate that our algorithm performs well in finite precision arithmetic.