Tight and efficient enclosure of matrix multiplication by using optimized BLAS

Abstract

This paper is concerned with the tight enclosure of matrix multiplication AB for two floating-point matrices A and B. The aim of this paper is to compute component-wise upper and lower bounds of the exact result C of the matrix multiplication AB by floating-point arithmetic. Namely, an interval matrix enclosing C is obtained. In this paper, new algorithms for enclosing C are proposed. The proposed algorithms are designed to mainly exploit the level 3 operations in BLAS. Although the proposed algorithms take around twice as much costs as a standard algorithm promoted by Oishi and Rump, the accuracy of the result by the proposed algorithms is better than that of the standard algorithm. At the end of this paper, we present numerical examples showing the efficiency of the proposed algorithms.We define the radius of an interval [C ij , C ij ] by (C ij −C ij )/2. For a dominant approach for this purpose, a method promoted by Oishi and Rump [7,9,12] is fast and useful. Throughout this paper, we call it Oishi-Rump's method. It exploits switches of rounding mode defined in IEEE 754, in particular, rounding upward and rounding downward are used. It is particularly worth noting that only matrix operations are required. Concretely speaking, two matrix products and only two switches of rounding mode are required to enclose the matrix multiplication (1).There is a so-called optimized BLAS (basic linear algebra subprograms) whose performance for the matrix multiplication is nearly peak. For example, GotoBLAS [13], Intel Math Kernel Library and ATLAS [14] are well known as the optimized BLAS. Moreover, matrix multiplication routines in these optimized BLAS can automatically be parallelized on several symmetric multi-processing environments. As an advantage of Oishi-Rump's method, such convenient routines for the matrix multiplication can be exploited for the enclosure. The Oishi-Rump's method is known to be wellbalanced between the tightness of the resultant interval matrix and the computational performance. This method has already been implemented on INTLAB developed by Rump [15]. Here, INTLAB is a fast and useful interval computation toolbox for MATLAB.To obtain an interval result whose radius is tighter than that calculated by Oishi-Rump's method, the use of some multi-precision libraries are seen as a means of achieving this purpose. For instance, MPFR [16] is a fast tool for multi-precision floating-point arithmetic. As it supports directed rounding in its arithmetic as in IEEE 754, it is possible to obtain enclosure of the matrix multiplication by MPFR. It is expected that the radius of the matrix given by MPFR becomes tight in proportion to computational precision.Recently, we have investigated an accurate algorithm for the matrix multiplication [17,18]. By specializing the algorithm, we develop efficient algorithms for obtaining tight enclosure of the matrix multiplication. As the proposed algorithms mainly exploit the matrix multiplication routines in BLAS, it receives much benefit from the optimized BLAS in terms of th...

Improvement of error-free splitting for accurate matrix multiplication

“…If the extra matrix multiplication by the binary32 arithmetic is accepted, then the accuracy of the result becomes more better. We confirm that the new error-free splitting algorithm can efficiently be applied into the interval enclosure of the matrix multiplication as in [14]. Algorithm 7.…”

“…In the following subsection, we compare accuracy of numerical results and computing times by each algorithm. In Section 4.2, we introduce the application of our algorithm to interval computations for matrix multiplication as in [14]. (rand(n) − 0.5).…”

“…It is shown in [14] that the error-free splitting is useful for interval arithmetic [24]. If (3) is satisfied, then we can develop efficient algorithms for interval enclosure of AB.…”

“…Moreover, even if we evaluate A (1) B (1) by the floating-point arithmetic, no rounding error occurs. Therefore, it is shown in [14,15] that the accuracy of the computed result by A (1) B (1) + A (1) B (2) + A (2) B is better than that by the pure floating-point evaluation of AB in many cases. In this paper, we improve an algorithm which produces A (1) , A (2) , B (1) and B (2) from A and B.…”

“…1 Recently, we have developed accurate algorithms for floating-point matrix multiplication [14,15]. Our algorithms involve several floating-point matrix products which can be computed by functions supported in BLAS (Basic Linear Algebra Subprograms).…”