Fast algorithms for floating-point interval matrix multiplication

We discuss several methods for real interval matrix multiplication. First, earlier studies of fast algorithms for interval matrix multiplication are introduced: naive interval arithmetic, interval arithmetic by midpoint-radius form by Oishi-Rump and its fast variant by Ogita-Oishi. Next, three new and fast algorithms are developed. The proposed algorithms require one, two or three matrix products, respectively. The point is that our algorithms quickly predict which terms become dominant radii in interval computations. We propose a hybrid method to predict which algorithm is suitable for optimizing performance and width of the result. Numerical examples are presented to show the efficiency of the proposed algorithms.     

Software to support numerical analysis teaching

This paper describes a Fortran90 library designed to support the teaching of numerical analysis and its applications. As well as covering traditional material it introduces recent and important ideas in numerical computation such as interval arithmetic and automatic differentiation. The library rests on a module realpac which provides real arithmetic in a range of precisions with a choice of rounding strategies. This, in turn, supports the implementation of an interval arithmetic module intpac. Derived data types and overloaded operations help inexperienced users to interface with unfamiliar data types such as intervals. The library also includes more conventional modules such as lepac for solving linear systems and minpac for nonlinear optimization. These, however, can be enhanced by being linked to more sophisticated tools for sparse matrix handling and automatic differentiation. As well as showing the main structure and scope of the software, the paper mentions some exercises that have successfully been performed by students.     

An interval component for continuous constraints

Local consistency techniques for numerical constraints over interval domains combine interval arithmetic, constraint inversion and bisection to reduce variable domains. In this paper, we study the problem of integrating any specific interval arithmetic library in constraint solvers. For this purpose, we design an interface between consistency algorithms and arithmetic. The interface has a two-level architecture: functional interval arithmetic at low-level, which is only a specification to be implemented by specific libraries, and a set of operations required by solvers, such as relational interval arithmetic or bisection primitives. This work leads to the implementation of an interval component by means of C++ generic programming methods. The overhead resulting from generic programming is discussed.     

Intrinsic Compiler Support for Interval Arithmetic

Interval arithmetic provides an efficient method for monitoring errors in numerical computations and for solving problems that cannot be efficiently solved with floating-point arithmetic. To support interval arithmetic, several software tools have been developed including interval arithmetic libraries, extended scientific programming languages, and interval-enhanced compilers. The main disadvantage of these software tools is their speed, since interval operations are implemented using function calls. In this paper, compiler support for interval arithmetic is investigated. In particular, the performance benefits of having the compiler inline interval operations to eliminate function call overhead is researched. Interval operations are inlined with the GNU gcc compiler and the performance of interval arithmetic is evaluated on a superscalar architecture. To implement interval operations with compiler support, the compiler produces sequences of instructions that use existing floating point hardware. Simulation results show that the compiler implementation of interval arithmetic is approximately 4 to 5 times faster than a functionally equivalent interval arithmetic software implementation with function call overhead and approximately 1.2 to 1.5 times slower than a dedicated interval arithmetic hardware implementation.     

One-modulus residue arithmetic algorithm to solve linear equations exactly

The paper presents an error-free algorithm to solve linear equations using the residue arithmetic. Simultaneously with solving linear equation system, the exact value of determinant of the system matrix is also calculated. The algorithm removes roundoff errors and according to this kind of errors ensures stability of the solution. It is suitable for implementation for computers with possibility of vector operations.     

On the Non‐Hierarchical Preconditioner Constructions

In this paper a new cyclic matrix representation of the H^–1/2 norm is presented. Its application as Schur complement preconditioning matrix requires only matrix‐vector multiplication. The computational cost of this matrix‐vector multiplication is O (N · log(N)) arithmetic operations, where N is the number of unknowns. The efficiency of the construction to elliptic problems has been verified by numerical tests. (© 2004 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)     

Applying fuzzy GERT with approximate fuzzy arithmetic based on the weakest t-norm operations to evaluate repairable reliability

In general, the fuzzy Graphical Evaluation and Review Technique (GERT) usually evaluates/analyzes variables with interval arithmetic (α-cut arithmetic) operations, especially those with complicated fuzzy systems. Thus the interval arithmetic operations may occur accumulating phenomenon of fuzziness in complicated systems, and the accumulating phenomenon of fuzziness may make decision-maker that cannot effectively evaluate problems/systems under vague environment. In order to overcome the accumulating phenomenon of fuzziness or credibly reduce fuzzy spreads, this study adopts approximate fuzzy arithmetic operations under the weakest t-norm arithmetic operations (T_ω) to evaluate fuzzy reliability models based on fuzzy GERT simulation technology. The approximate fuzzy arithmetic operations employ principle of interval arithmetic under the weakest t-norm arithmetic operations. Therefore, the novel fuzzy arithmetic operations may obtain fitter decision values, which have smaller fuzziness accumulating, under vague environment. In numerical examples the approximate fuzzy arithmetic operations has evidenced that it can successfully calculate results of fuzzy operations as interval arithmetic, and can more effectively reduce fuzzy spreads. In the real fuzzy repairable reliability model the performance also shows that the approximate fuzzy arithmetic operations successfully analyze the reliability problem and obtain more confident fuzzy results.     

Using successive approximations for improving the convergence of GMRES method

In this paper, our attention is concentrated on the GMRES method for the solution of the system (I–T)x=b of linear algebraic equations with a nonsymmetric matrix. We perform m pre-iterations y _l+1 =T _yl +b before starting GMRES and put y _m for the initial approximation in GMRES. We derive an upper estimate for the norm of the error vector in dependence on the mth powers of eigenvalues of the matrix T Further we study under what eigenvalues lay-out this upper estimate is the best one. The estimate shows and numerical experiments verify that it is advisable to perform pre-iterations before starting GMRES as they require fewer arithmetic operations than GMRES. Towards the end of the paper we present a numerical experiment for a system obtained by the finite difference approximation of convection-diffusion equations.     

一致性区间数判断矩阵的性质及排序

讨论了具有序传递性的区间数判断矩阵的性质,给出了一致性区间数判断矩阵具有序传递性的一个充要条件.在此基础上,证明了一致性区间效判断矩阵的权值可行域中所有向量的分量在不劣于的次序关系下具有完全相同的排列次序.因此,用权值可行域的顶点的算术平均或几何平均作为一致性区间数判断矩阵的排序向量是合理的.     

Using Chebyshev polynomials and approximate inverse triangular factorizations for preconditioning the conjugate gradient method

In order to precondition a sparse symmetric positive definite matrix, its approximate inverse is examined, which is represented as the product of two sparse mutually adjoint triangular matrices. In this way, the solution of the corresponding system of linear algebraic equations (SLAE) by applying the preconditioned conjugate gradient method (CGM) is reduced to performing only elementary vector operations and calculating sparse matrix-vector products. A method for constructing the above preconditioner is described and analyzed. The triangular factor has a fixed sparsity pattern and is optimal in the sense that the preconditioned matrix has a minimum K-condition number. The use of polynomial preconditioning based on Chebyshev polynomials makes it possible to considerably reduce the amount of scalar product operations (at the cost of an insignificant increase in the total number of arithmetic operations). The possibility of an efficient massively parallel implementation of the resulting method for solving SLAEs is discussed. For a sequential version of this method, the results obtained by solving 56 test problems from the Florida sparse matrix collection (which are large-scale and ill-conditioned) are presented. These results show that the method is highly reliable and has low computational costs.     

Optimal accumulation of Jacobian matrices by elimination methods on the dual computational graph

The accumulation of the Jacobian matrix F of a vector function can be regarded as a transformation of its linearized computational graph into a subgraph of the directed complete bipartite graph K_n,m. This transformation can be performed by applying different elimination techniques that may lead to varying costs for computing F. This paper introduces face elimination as the basic technique for accumulating Jacobian matrices by using a minimal number of arithmetic operations. Its superiority over both edge and vertex elimination methods is shown. The intention is to establish the conceptual basis for the ongoing development of algorithms for optimizing the computation of Jacobian matrices.     

Eigenvectors of interval matrices over max–plus algebra

The behaviour of a discrete-event dynamic system is often conveniently described using a matrix algebra with operations max and plus. Such a system moves forward in regular steps of length equal to the eigenvalue of the system matrix, if it is set to operate at time instants corresponding to one of its eigenvectors. However, due to imprecise measurements, it is often unappropriate to use exact matrices. One possibility to model imprecision is to use interval matrices. We show that the problem to decide whether a given vector is an eigenvector of one of the matrices in the given matrix interval is polynomial, while the complexity of the existence problem of a universal eigenvector remains open. As an aside, we propose a combinatorial method for solving two-sided systems of linear equations over the max–plus algebra.     

An error-free algorithm to solve a linear system of polynomial equations

The paper presents an error-free algorithm to solve a system of linear equations with polynomial coefficients. Modular arithmetic in residual polynomial class and in residual numeric class is employed. The algorithm is iterative and well suited for implementation for computers with vector operations and fast and error-free convolutors.     

A hybrid iterative method for symmetric indefinite linear systems

Hybrid iterative methods that combine a conjugate direction method with a simpler iteration scheme, such as Chebyshev or Richardson iteration, were first proposed in the 1950s. The ease with which Chebyshev and Richardson iteration can be implemented efficiently on a large variety of computer architectures has in recent years lead to renewed interest in iterative methods that use Chebyshev or Richardson iteration. This paper presents a new hybrid iterative method for the solution of linear systems of equations with a symmetric indefinite matrix. Our method combines the conjugate residual method with Richardson iteration. Special attention is paid to the determination of two real intervals, one on each side of the origin, that contain most of the eigenvalues of the matrix. These intervals are used to compute suitable iteration parameters for Richardson iteration. We also discuss when to switch between the methods. The hybrid scheme typically uses the Richardson method for most iterations, and this reduces the number of arithmetic vector operations significantly compared with the number of arithmetic vector operations required when only the conjugate residual method is used. Computed examples illustrate the competitiveness of the hybrid scheme.     

Numerical enclosure for multiple eigenvalues of an Hermitian matrix whose graph is a tree

In this paper we consider a numerical enclosure method for multiple eigenvalues of an Hermitian matrix whose graph is a tree. If an Hermitian matrix A whose graph is a tree has multiple eigenvalues, it has the property that matrices which are associated with some branches in the undirected graph of A have the same eigenvalues. By using this property and interlacing inequalities for Hermitian matrices, we show an enclosure method for multiple eigenvalues of an Hermitian matrix whose graph is a tree. Since we do not generally know whether a given matrix has exactly a multiple eigenvalue from approximate computations, we use the property of interlacing inequalities to enclose some eigenvalues including multiplicities.In this process, we only use the enclosure of simple eigenvalues to enclose a multiple eigenvalue by using a computer and interval arithmetic.     

Tropical Arithmetic and Matrix Algebra

This article introduces a new structure of commutative semiring, generalizing the tropical semiring, and having an arithmetic that modifies the standard tropical operations, i.e., summation and maximum. Although our framework is combinatorial, notions of regularity and invertibility arise naturally for matrices over this semiring; we show that a tropical matrix is invertible if and only if it is regular.     

A look-ahead algorithm for the solution of general Hankel systems

Summary The solution of systems of linear equations with Hankel coefficient matrices can be computed with onlyO(n ²) arithmetic operations, as compared toO(n ³) operations for the general cases. However, the classical Hankel solvers require the nonsingularity of all leading principal submatrices of the Hankel matrix. The known extensions of these algorithms to general Hankel systems can handle only exactly singular submatrices, but not ill-conditioned ones, and hence they are numerically unstable. In this paper, a stable procedure for solving general nonsingular Hankel systems is presented, using a look-ahead technique to skip over singular or ill-conditioned submatrices. The proposed approach is based on a look-ahead variant of the nonsymmetric Lanczos process that was recently developed by Freund, Gutknecht, and Nachtigal. We first derive a somewhat more general formulation of this look-ahead Lanczos algorithm in terms of formally orthogonal polynomials, which then yields the look-ahead Hankel solver as a special case. We prove some general properties of the resulting look-ahead algorithm for formally orthogonal polynomials. These results are then utilized in the implementation of the Hankel solver. We report some numerical experiments for Hankel systems with ill-conditioned submatrices.The research of the first author was supported by DARPA via Cooperative Agreement NCC 2-387 between NASA and the Universities Space Research Association (USRA).The research of the second author was supported in part by NSF grant DRC-8412314 and Cooperative Agreement NCC 2-387 between NASA and the Universities Space Research Association (USRA).     

Chebyshev model arithmetic for factorable functions

This article presents an arithmetic for the computation of Chebyshev models for factorable functions and an analysis of their convergence properties. Similar to Taylor models, Chebyshev models consist of a pair of a multivariate polynomial approximating the factorable function and an interval remainder term bounding the actual gap with this polynomial approximant. Propagation rules and local convergence bounds are established for the addition, multiplication and composition operations with Chebyshev models. The global convergence of this arithmetic as the polynomial expansion order increases is also discussed. A generic implementation of Chebyshev model arithmetic is available in the library MC++. It is shown through several numerical case studies that Chebyshev models provide tighter bounds than their Taylor model counterparts, but this comes at the price of extra computational burden.     

Factorization of analytic functions by means of Koenig's theorem and Toeplitz computations

Summary. By providing a matrix version of Koenig's theorem we reduce the problem of evaluating the coefficients of a monic factor r(z) of degree h of a power series f(z) to that of approximating the first h entries in the first column of the inverse of an Toeplitz matrix in block Hessenberg form for sufficiently large values of n. This matrix is reduced to a band matrix if f(z) is a polynomial. We prove that the factorization problem can be also reduced to solving a matrix equation for an matrix X, where is a matrix power series whose coefficients are Toeplitz matrices. The function is reduced to a matrix polynomial of degree 2 if f(z) is a polynomial of degreeN and . These reductions allow us to devise a suitable algorithm, based on cyclic reduction and on the concept of displacement rank, for generating a sequence of vectors that quadratically converges to the vector having as components the coefficients of the factor r(z). In the case of a polynomial f(z) of degree N, the cost of computing the entries of given is arithmetic operations, where is the cost of solving an Toeplitz-like system. In the case of analytic functions the cost depends on the numerical degree of the power series involved in the computation. From the numerical experiments performed with several test polynomials and power series, the algorithm has shown good numerical properties and promises to be a good candidate for implementing polynomial root-finders based on recursive splitting strategies. Applications to solving spectral factorization problems and Markov chains are also shown. Received September 9, 1998 / Revised version received November 14, 1999 / Published online February 5, 2001