A radix-independent error analysis of the Cornea-Harrison-Tang method

International audienceAssuming floating-point arithmetic with a fused multiply-add operation and rounding to nearest, the Cornea-Harrison-Tang method aims to evaluate expressions of the form $ab+cd$with high relative accuracy. In this paper we provide a rounding error analysis of this method,which unlike previous studiesis not restricted to binary floating-point arithmetic but holds for any radix $\beta$.We show first that an asymptotically optimal bound on the relative error of this method is$2u + O(u^2)$, where $u= \frac{1}{2}\beta^{1-p}$ is the unit roundoff in radix $\beta$ and precision $p$.Then we show that the possibility of removing the $O(u^2)$ term from this bound is governed bythe radix parity andthe tie-breaking strategy used for rounding: if $\beta$ is odd or rounding is \emph{to nearest even}, then the simpler bound $2u$ is obtained,while if $\beta$ is even and rounding is \emph{to nearest away}, then there exist floating-point inputs $a,b,c,d$ that lead to a relative error larger than $2u + \frac{2}{\beta} u^2 - 4u^3$.All these results hold provided underflows and overflows do not occurand under some mild assumptions on $p$ satisfied by IEEE 754-2008 formats

A (hopefully) friendly introduction to the complexity of polynomial matrix computations

This paper aims at a friendly introduction to the field of fast algorithms for polynomial matrices, and surveys the results of the ISSAC 2003 paper 'On the Complexity of Polynomial Matrix Computations' by Pascal Giorgi, Claude-Pierre Jeannerod, and Gilles Villard

Exploiting structure in floating-point arithmetic

Invited paper - MACIS 2015 (Sixth International Conference on Mathematical Aspects of Computer and Information Sciences)International audienceThe analysis of algorithms in IEEE floating-point arithmetic is most often carried out via repeated applications of the so-called standard model, which bounds the relative error of each basic operation by a common epsilon depending only on the format. While this approach has been eminently useful for establishing many accuracy and stability results, it fails to capture most of the low-level features that make floating-point arithmetic so highly structured. In this paper, we survey some of those properties and how to exploit them in rounding error analysis. In particular, we review some recent improvements of several classical, Wilkinson-style error bounds from linear algebra and complex arithmetic that all rely on such structure properties

Computing specified generators of structured matrix inverses

International audienceThe asymptotically fastest known divide-and-conquer methods for inverting dense structured matrices are essentially variations or extensions of the Morf/Bitmead-Anderson algorithm. Most of them must deal with the growth in length of intermediate generators, and this is done by incorporating various generator compression techniques into the algorithms. One exception is an algorithm by Cardinal, which in the particular case of Cauchy-like matrices avoids such growth by focusing on well-specied, already compressed generators of the inverse. In this paper, we extend Cardinal's method to a broader class of structured matrices including those of Vandermonde, Hankel, and Toeplitz types. Besides, some rst experimental results illustrate the practical interest of the approach

Fast Computation of Minimal Interpolation Bases in Popov Form for Arbitrary Shifts

We compute minimal bases of solutions for a general interpolation problem, which encompasses Hermite-Pad\'e approximation and constrained multivariate interpolation, and has applications in coding theory and security. This problem asks to find univariate polynomial relations between $m$ vectors of size $\sigma$; these relations should have small degree with respect to an input degree shift. For an arbitrary shift, we propose an algorithm for the computation of an interpolation basis in shifted Popov normal form with a cost of $\mathcal{O}\tilde{~}(m^{\omega-1} \sigma)$ field operations, where $\omega$ is the exponent of matrix multiplication and the notation $\mathcal{O}\tilde{~}(\cdot)$ indicates that logarithmic terms are omitted. Earlier works, in the case of Hermite-Pad\'e approximation and in the general interpolation case, compute non-normalized bases. Since for arbitrary shifts such bases may have size $\Theta(m^2 \sigma)$, the cost bound $\mathcal{O}\tilde{~}(m^{\omega-1} \sigma)$ was feasible only with restrictive assumptions on the shift that ensure small output sizes. The question of handling arbitrary shifts with the same complexity bound was left open. To obtain the target cost for any shift, we strengthen the properties of the output bases, and of those obtained during the course of the algorithm: all the bases are computed in shifted Popov form, whose size is always $\mathcal{O}(m \sigma)$. Then, we design a divide-and-conquer scheme. We recursively reduce the initial interpolation problem to sub-problems with more convenient shifts by first computing information on the degrees of the intermediate bases.Comment: 8 pages, sig-alternate class, 4 figures (problems and algorithms

On relative errors of floating-point operations: optimal bounds and applications

International audienceRounding error analyses of numerical algorithms are most often carried out via repeated applications of the so-called standard models of floating-point arithmetic. Given a round-to-nearest function fl and barring underflow and overflow, such models bound the relative errors E 1 (t) = |t − fl(t)|/|t| and E 2 (t) = |t − fl(t)|/|fl(t)| by the unit roundoff u. This paper investigates the possibility and the usefulness of refining these bounds, both in the case of an arbitrary real t and in the case where t is the exact result of an arithmetic operation on some floating-point numbers. We show that E 1 (t) and E 2 (t) are optimally bounded by u/(1 + u) and u, respectively, when t is real or, under mild assumptions on the base and the precision, when t = x ± y or t = xy with x, y two floating-point numbers. We prove that while this remains true for division in base β > 2, smaller, attainable bounds can be derived for both division in base β = 2 and square root. This set of optimal bounds is then applied to the rounding error analysis of various numerical algorithms: in all cases, we obtain significantly shorter proofs of the best-known error bounds for such algorithms, and/or improvements on these bounds themselves

Simultaneous floating-point sine and cosine for VLIW integer processors

Accepted for publication in the proceedings of the 23rd IEEE International Conference on Application-specific Systems, Architectures and Processors (ASAP 2012).International audienceGraphics and signal processing applications often require that sines and cosines be evaluated at a same floating-point argument, and in such cases a very fast computation of the pair of values is desirable. This paper studies how 32-bit VLIW integer architectures can be exploited in order to perform this task accurately for IEEE single precision. We describe software implementations for sinf, cosf, and sincosf over [-pi/4,pi/4] that have a proven 1-ulp accuracy and whose latency on STMicroelectronics' ST231 VLIW integer processor is 19, 18, and 19 cycles, respectively. Such performances are obtained by introducing a novel algorithm for simultaneous sine and cosine that combines univariate and bivariate polynomial evaluation schemes

Improved error bounds for inner products in floating-point arithmetic

International audienceGiven two floating-point vectors $x,y$ of dimension $n$ and assuming rounding to nearest, we show that if no underflow or overflow occurs, any evaluation order for inner product returns a floating-point number $\hat r$ such that $|{\hat r}-x^Ty| \le nu|x|^T|y|$ with $u$ the unit roundoff. This result, which holds for any radix and with no restriction on $n$, can be seen as a generalization of a similar bound given in~\cite{Rump12} for recursive summation in radix $2$, namely $|{\hat r}- x^Te| \le (n-1)u|x|^Te$ with $e=[1,1,\ldots,1]^T$. As a direct consequence, the error bound for the floating-point approximation $\hat C$ of classical matrix multiplication with inner dimension $n$ simplifies to $|\hat{C}-AB|\le nu|A||B|$

Further remarks on Kahan summation with decreasing ordering

We consider Kahan's compensated summation of $n$ floating-point numbers ordered as $|x_1| \ge \cdots \ge |x_n|$,and show that in IEEE arithmetic a large relative error can occur for a dimension as small as $n=4$. This answers a question raised in particular by Priest and Higham