Currently FORTRAN versions are available which compare execution with LINPACK's dgefa/dgesl. If you would like to use this source code, please contact Stott Parker.
D. Stott Parker and Dinh Le
Gaussian elimination is probably the best known and most widely used method for solving linear systems, computing determinants, and finding matrix decompositions. While the basic elimination procedure is simple to state and implement, it becomes more complicated with the addition of a pivoting procedure, which handles degenerate matrices having zero elements on the diagonal. Pivoting can significantly complicate the algorithm, increase data movement, and reduce speed, particularly on high-performance computers.In this paper we propose an alternative scheme for performing Gaussian elimination that first preconditions the input matrix by multiplying it with random matrices, whose inverses can be applied subsequently. At the expense of these multiplications, and making the linear system dense if it was not already, this approach makes the system `nondegenerate' --- subsystems have full rank --- with probability 1. This preconditioning has the effect of (almost certainly) eliminating the need for pivoting.
The randomization approach also suggests many interesting developments beyond the elimination of pivoting. Generally, randomization is a new technique for developing matrix computations that are suitable for high-performance computers.
D. Stott Parker
Theory and practice of computational linear algebra differ over the issue of degeneracy. Block matrix decompositions are used heavily in theory, but less in practice, since even when a matrix is nondegenerate (has full rank) its block submatrices can be degenerate. The potential degeneracy of block submatrices can completely prevent practical use of block matrix algorithms.Gaussian elimination is an important example of an algorithm affected by the possibility of degeneracy. While the basic elimination procedure is simple to state and implement, it becomes more complicated with the addition of a pivoting procedure, which handles degenerate matrices having zeros on the diagonal. Pivoting can significantly complicate the algorithm, increase data movement, and reduce speed, particularly on high-performance computers.
We propose a randomization scheme that preconditions an input matrix by multiplying it with random matrices, where this multiplication can be performed efficiently. At the expense of these multiplications, and making the linear system dense if it were not already, this approach makes the system nondegenerate with probability 1. We show that this preconditioning (almost certainly) eliminates the need for pivoting.
In this paper we study the properties of Random Butterfly Transformations (RBT), an interesting class of random matrices. RBTs resemble the Fast Fourier Transform (FFT). RBTs have three important benefits. First, they have provable randomization properties. Second, RBTs can be performed efficiently, and in particular should be useful on parallel and vector machines with architectures that support FFT-like computations. Third, RBTs deliver important numerical properties when restricted to unitary matrices, and even when restricted to well-behaved, non-unitary matrices.
The RBT approach also suggests many interesting developments beyond the elimination of pivoting. Randomization offers new possibilities for linear systems solution that are well-suited to high-performance computers.
D. Stott Parker
We present a new randomization scheme that preconditions an input linear system by multiplying it with nonsingular random matrices. At the expense of these multiplications, and making the linear system dense if it was not already, this approach makes blocks of the resulting system nonsingular with probability 1.Specifically, in this paper we consider random matrices that are `random butterflies', and call their application a Random Butterfly Transformation (RBT). RBTs can be performed efficiently, and in particular should be useful on parallel and vector machines with architectures that support FFT-like computations.
Block Gaussian elimination is an important example of an algorithm with which this RBT can be used. Gaussian elimination is complicated by pivoting, which handles degenerate matrices having zero elements on the diagonal. Pivoting can significantly complicate the algorithm, increase data movement, and reduce speed, particularly on high-performance computers. We show that this preconditioning has the effect of eliminating the need for pivoting (with probability 1).
The RBT approach also suggests many interesting developments beyond the elimination of pivoting. Randomization appears to offer new possibilities for linear systems solution that are well-suited to high-performance computers.
D. Stott Parker
Gaussian Elimination is probably the best known and most widely used method for solving linear systems, computing determinants, and decomposing matrices. It is an important tool in the kernel of most computational linear algebra packages today. While the basic elimination procedure is simple to state and implement, it is complex enough that it is nontrivial to reason about.In this paper we resurrect a characterization of the behavior of Gaussian Elimination, namely that every one of its intermediate results is a ratio of determinants of submatrices of the input matrix, and show its equivalence to another characterization in terms of Schur complements. These characterizations seem to be unfamiliar to numerical analysts, but deserve wider attention since they lead to a deep appreciation of Gaussian Elimination.
The characterizations can be used in obtaining bounds for roundoff error with Gaussian Elimination. The characterizations also help in better understanding the effects of roundoff error and pivoting.
D. Stott Parker
In this paper we use Sylvester's determinantal identity to derive two reformulations of the Gaussian Elimination algorithm. One reformulation obviates the need for division operations, and another produces intermediate results that are always specific determinants of submatrices of the input matrix. Both algorithms find L D^{-1} U decompositions of the input matrix A, and L D^{-1} L^T decompositions when A is symmetric. Although they perform quite different computations than the traditional algorithm, they differ formally only by an underlying change of variables.
D. Stott Parker
For two decades Schur complements have seen increasing applications in linear algebra, often as abstractions of Gaussian elimination. It is known that they obey certain nontrivial identities, such as Crabtree and Haynsworth's quotient property. We began this work asking if there were a theory for deciding their properties in general.Lambek's Categorial Grammar is a deductive system introduced in 1958 by Lambek as a mathematical foundation for a syntactic calculus of English sentences. We show that Categorial Grammar gives a deductive system for deciding properties LU- and UL-decompositions, and about Gaussian elimination. It also can be used to derive identities obeyed by Schur complements.
At first impression this seems to be a strange result, connecting two unrelated topics. In retrospect, though, it is a consequence of the way both use quotients. It may have applications in developing both grammatical formalisms and numerical algorithms.
D. Stott Parker and Dinh Le
With the quadtree representation of matrices popularized by Wise, many matrix algorithms are naturally implemented in a recursive functional style. This style is useful in exposing parallelism, identifying mappings for matrix computations onto particular computer architectures, and more generally in exposing new ways to reason about matrix algorithms.The world has not beaten a path to quadtree-based algorithms, unfortunately. It seems there are a number of reasons that quadtree algorithms and architectures have not reached their potential in practice. Important among these are that quadtree representions are not unique, and that straightforward implementations of block-matrix algorithms fail for degenerate matrices (e.g., with singular quadrants). Proposed extensions of quadtree-based algorithms to deal with multiple representations, or degenerate matrices (such as resorting to pivoting), have not really removed the objections: they lose either numerical stability or the benefits in elegance, parallelism, and functional style that encouraged the use of quadtrees in the first place.
We survey the situation in this paper, and present new approaches for handling degenerate quadrants in quadtree algorithms: First, we rely on canonical quadtree representations. Second, we randomize the input matrix in a way that guarantees nondegeneracy of the quadrants, and also permits derandomization after the algorithm. This approach should make the quadtree representation more successful in real matrix computations.
D. Stott Parker, Brad Pierce
Recently, randomizing input transformations have been developed for streamlining Gaussian elimination. Multiplying by special sparse random matrices, these transformations produce an equivalent linear system which, with probability 1, Gaussian elimination can solve without any augmentation by pivoting.In this paper we study properties of the unitary Randomizing Discrete Fourier Transform (RDFT), which scales rows by random values from the complex unit circle and then applies the conventional DFT. When the problem size is a power of 2, the Randomizing Fast Fourier Transform (RFFT) algorithm can efficiently implement the RDFT, using existing FFT code, on the many parallel and vector machines whose architectures support FFT-like computations.
Dinh Le, D. Stott Parker, Brad Pierce
Many standard matrix algorithms are difficult to implement as systolic arrays, because they involve data movement that cannot be determined a priori, resulting in high timing complexity. In the case of dense matrix inversion, Gaussian elimination (GE) with pivoting has largely been supplanted by alternative schemes, such as Givens rotations and the Gram-Schmidt method, which are costly and complex, and the simpler GE with pairwise pivoting and GE without pivoting, which can break down on ``degenerate'' inputs, i.e., invertible matrices with noninvertible submatrices.A new input randomization technique efficiently transforms the linear problem Ax=b into the randomized problem VAx = Vb, where the matrix V is chosen from a special class of random matrices. If A is nonsingular, then, with probability 1, GE without pivoting can be successfully applied to VA. As is demonstrated by an extended tutorial example, this simple algorithm is amenable to current-generation matrix algorithm compilers, such as UCLA's MAMACG system.
Dinh Le, D. Stott Parker
Recursive block decomposition algorithms (also known as quadtree algorithms when the blocks are all square) have been proposed to solve well-known problems such as matrix addition, multiplication, inversion, determinant computation, block LDU decomposition, discrete Fourier transform, and Cholesky and QR factorization.Until now, such algorithms have often been viewed as impractical for building general tools. Specifically, recursive block decomposition methods are impractical for matrix inversion and Gaussian elimination, since both require leading submatrices of the input matrix to be nonsingular (which is not always guaranteed).
We show how to use randomization to guarantee that submatrices meet these requirements, and to make recursive block decomposition methods practical on well-conditioned input matrices. The resulting algorithms are elegant, and we show the recursive programs can perform well for both dense and sparse matrices, although with randomization dense computations seem most practical.
Brad Pierce, D. Stott Parker
A block matrix algorithm for linear systems solution is derived using Haynsworth's quotient formula for Schur complements. The algorithm is essentially a generalization of Gauss-Jordan elimination.
Brad Pierce, D. Stott Parker
Recently we have reported on new LU factorization and block matrix methods that use randomized linear preconditioning in place of dynamic strategies such as partial or pairwise pivoting. We elaborate on the theory underlying these a priori methods, and apply it to several examples, including a new method for the calculation of Moore-Penrose pseudoinverses.