Standard Gaussian Elimination Research Articles

On the complexity of the [formula omitted] Gröbner basis algorithm

We study the complexity of Gröbner bases computation, in particular in the generic situation where the variables are in simultaneous Noether position with respect to the system.We give a bound on the number of polynomials of degree d in a Gröbner basis computed by Faugère's F5 algorithm (2002) in this generic case for the grevlex ordering (which is also a bound on the number of polynomials for a reduced Gröbner basis, independently of the algorithm used). Next, we analyse more precisely the structure of the polynomials in the Gröbner bases with signatures that F5 computes and use it to bound the complexity of the algorithm.Our estimates show that the version of F5 we analyse, which uses only standard Gaussian elimination techniques, outperforms row reduction of the Macaulay matrix with the best known algorithms for moderate degrees, and even for degrees up to the thousands if Strassen's multiplication is used. The degree being fixed, the factor of improvement grows exponentially with the number of variables.

Iterative convergence of boundary-volume integral equation method

The boundary-volume integral equation numerical technique can be a powerful tool for piecewise heterogeneous media, but it is limited to small problems or low frequencies because of great computational cost. Therefore, a restarted GMRES method is applied to solve large-scale boundary-volume scattering problems in this paper to overcome the computational barrier. The iterative method is firstly applied to responses of dimensionless frequency to a semicircular alluvial valley filled with sediments, compared with the standard Gaussian elimination method. Then the method is tested by a heterogeneous multilayered model to show its applicability. Numerical experiments indicate that the preconditioned GMRES method can significantly improve computational efficiency especially for large Earth models and high frequencies, but with a faster convergence for the left diagonal preconditioning.

The Cache-Oblivious Gaussian Elimination Paradigm: Theoretical Framework, Parallelization and Experimental Evaluation

We consider triply-nested loops of the type that occur in the standard Gaussian elimination algorithm, which we denote by GEP (or the Gaussian Elimination Paradigm). We present two related cache-oblivious methods I-GEP and C-GEP, both of which reduce the number of cache misses incurred (or I/Os performed) by the computation over that performed by standard GEP by a factor of $\sqrt{M}$, where M is the size of the cache. Cache-oblivious I-GEP computes in-place and solves most of the known applications of GEP including Gaussian elimination and LU-decomposition without pivoting and Floyd-Warshall all-pairs shortest paths. Cache-oblivious C-GEP uses a modest amount of additional space, but is completely general and applies to any code in GEP form. Both I-GEP and C-GEP produce system-independent cache-efficient code, and are potentially applicable to being used by optimizing compilers for loop transformation. We present parallel I-GEP and C-GEP that achieve good speed-up and match the sequential caching performance cache-obliviously for both shared and distributed caches for sufficiently large inputs. We present extensive experimental results for both in-core and out-of-core performance of our algorithms. We consider both sequential and parallel implementations, and compare them with finely-tuned cache-aware BLAS code for matrix multiplication and Gaussian elimination without pivoting. Our results indicate that cache-oblivious GEP offers an attractive trade-off between efficiency and portability.

Efficient method for analyzing leaky cavities in two-dimensional photonic crystals

An efficient numerical method is developed for analyzing leaky cavities in two-dimensional photonic crystals, assuming that the cavities are located in the vicinity of a photonic crystal waveguide or close to a homogeneous medium. A leaky cavity mode is characterized by its complex eigenfrequency and by the outgoing wave behavior of the field in the nearby waveguide or homogeneous medium. In our method, the eigenfrequency of a leaky cavity mode is determined from a condition formulated on one single edge of a defect unit cell. Such a condition is obtained by using rigorous boundary conditions to terminate photonic crystal waveguides and homogeneous media, Dirichlet-to-Neumann maps of the unit cells to reduce the problem to the cell edges, and standard Gaussian elimination to further reduce the problem to one single edge. The reduction process can be performed efficiently. Accurate solutions can be obtained when a small number of discrete points are used on each edge of the unit cells.

Efficient Algorithms for Reconstructing Zero-Recombinant Haplotypes on a Pedigree Based on Fast Elimination of Redundant Linear Equations

Computational inference of haplotypes from genotypes has attracted a great deal of attention in the computational biology community recently, partially driven by the international HapMap project. In this paper, we study the question of how to efficiently infer haplotypes from genotypes of individuals related by a pedigree, assuming that the hereditary process was free of mutations (i.e., the Mendelian law of inheritance) and recombinants. The problem has recently been formulated as a system of linear equations over the finite field of $F(2)$ and solved in $O(m^3n^3)$ time by using standard Gaussian elimination, where $m$ is the number of loci (or markers) in a genotype and $n$ the number of individuals in the pedigree. We give a much faster algorithm with running time $O(mn^2+n^3\log^2n\log\log n)$. The key ingredients of our construction are (i) a new system of linear equations based on some spanning tree of the pedigree graph and (ii) an efficient method for eliminating redundant equations in a system of $O(mn)$ linear equations over $O(n)$ variables. Although such a fast elimination method is not known for general systems of linear equations, we take advantage of the underlying pedigree graph structure and recent progress on low-stretch spanning trees.

Fully implicit parallel simulation of single neurons

When a multi-compartment neuron is divided into subtrees such that no subtree has more than two connection points to other subtrees, the subtrees can be on different processors and the entire system remains amenable to direct Gaussian elimination with only a modest increase in complexity. Accuracy is the same as with standard Gaussian elimination on a single processor. It is often feasible to divide a 3-D reconstructed neuron model onto a dozen or so processors and experience almost linear speedup. We have also used the method for purposes of load balance in network simulations when some cells are so large that their individual computation time is much longer than the average processor computation time or when there are many more processors than cells. The method is available in the standard distribution of the NEURON simulation program.

Fully implicit parallel simulation of single neurons

When tree topology matrices are divided into subtrees where each subtree is on a different cpu and with the constraint that other subtrees are not connected to a given subtree at more than two distinct points (defining a backbone path on that subtree), the entire system remains amenable to direct gaussian elimination. The complexity increase is twice the number of divisions and four times the number of multiplications normally required along the backbones due to the necessity, during the triangularization phase, of transforming the tridiagonal backbone into an N topology matrix. In addition, each subtree is required to send its root diagonal and right hand side element, or, in the case of a subtree with a backbone, the 2 × 2 matrix and right hand sides of the backbone end points, to one of the cpus where that information is added together to form a reduced tree matrix of rank equal to the number of split points on the cell. The reduced tree matrix equation is solved, giving the voltages at the split points, and this information is sent back to the appropriate subtrees on the other cpus. Those subtrees with backbones can then use the N topology to quickly compute the voltages along the backbone and everyone can complete the back substitution phase of their gaussian elimination. Accuracy is the same as with standard gaussian elimination on a single cpu and any quantitative differences are attributed to accumulated round off error due to different ordering of subtrees containing backbones. With this method, it is often feasible to divide a 3-d reconstructed neuron model into a dozen or so pieces and experience almost linear speedup. We have used the method for purposes of load balance in network simulations when some cells are very much larger than the average cell and there are more cpus than cells. The method is available in the current standard distribution of the NEURON simulation environment.

Parallel Complexity of Numerically Accurate Linear System Solvers

We prove a number of negative results about practical (i.e., work efficient and numerically accurate) algorithms for computing the main matrix factorizations. In particular, we prove that the popular Householder and Givens methods for computing the QR decomposition are P-complete, and hence presumably inherently sequential, under both real and floating point number models. We also prove that Gaussian elimination (GE) with a weak form of pivoting, which aims only at making the resulting algorithm nondegenerate, is likely to be inherently sequential as well. Finally, we prove that GE with partial pivoting is P-complete over GF(2) or when restricted to symmetric positive definite matrices, for which it is known that even standard GE (no pivoting) does not fail. Altogether, the results of this paper give further formal support to the widespread belief that there is a tradeoff between parallelism and accuracy in numerical algorithms.

Improved T-matrix computations for large, nonabsorbing and weakly absorbing nonspherical particles and comparison with geometrical-optics approximation

We show that the use of a matrix inversion scheme based on a special lower triangular-upper triangular factorization rather than on the standard Gaussian elimination significantly improves the numerical stability of T-matrix computations for nonabsorbing and weakly absorbing nonspherical particles. As a result, the maximum convergent size parameter for particles with small or zero absorption can increase by a factor of several and can exceed 100. We describe an improved scheme for evaluating Clebsch-Gordon coefficients with large quantum numbers, which allowed us to extend the analytical orientational averaging method developed by Mishchenko [J. Opt. Soc. Am. A 8, 871 (1991)] to larger size parameters. Comparisons of T-matrix and geometrical optics computations for large, randomly oriented spheroids and finite circular cylinders show that the applicability range of the ray-tracing approximation depends on the imaginary part of the refractive index and is different for different elements of the scattering matrix.

On the ordering of sparse linear systems

In this paper we consider the algorithms for transforming an n × n sparse matrix A into another matrix B such that Gaussian elimination applied to B takes time asymptotically less than n 3. These algorithms take the sparse matrix A as input, and return a pair of permutation matrices P, Q such that B = PAQ has a small bandwidth, or some other desirable form. We study the average effectiveness of these algorithms by using random matrices with Θ( n) nonzero elements. We prove that with high probability these algorithms cannot produce a reduction of the asymptotic cost of the standard Gaussian elimination algorithm. We also study the effectiveness of these algorithms for ordering very sparse matrices. We show that there exist matrices with 3 n nonzeros for which reordering rows and columns does not reduce the asymptotic cost of Gaussian elimination. We also prove that each matrix with at most two nonzeros in each row and in each column, can be transformed into a banded matrix with bandwidth five.

Perturbation results for nearly uncoupled Markov chains with applications to iterative methods

The standard perturbation theory for linear equations states that nearly uncoupled Markov chains (NUMCs) are very sensitive to small changes in the elements. Indeed, some algorithms, such as standard Gaussian elimination, will obtain poor results for such problems. A structured perturbation theory is given that shows that NUMCs usually lead to well conditioned problems. It is shown that with appropriate stopping, criteria, iterative aggregation/disaggregation algorithms will achieve these structured error bounds. A variant of Gaussian elimination due to Grassman, Taksar and Heyman was recently shown by O'Cinneide to achieve such bounds.

A parallel elimination algorithm for the solution of dense linear systems

For the solution of an N × N linear system Ax = d, based on the folding method of Evans and Hatzopoulos [4] and by partitioning the system into r 2 blocks each of size n × n (N = r n), in [3] we had described a parallel elimination method suitable for an r-processor machine. The serial arithmetical operations count for this algorithm is 0(1/3(3 − l/r)N 3) thereby the algorithm could perform with efficiency E r = 2/(3 − 1/r); note that 2/3 ≤ Er ≤ 4/5. In the present paper we consider an alternative elimination procedure together with the method of partitioning to obtain a parallel algorithm for the solution of dense linear systems. The serial arithmetical operations count for the present algorithm is 0(2/3N 3), which is of the same order as that for the standard Gaussian elimination method; thus the present algorithm could perform with efficiency close to 1 on a multiprocessor machine.

Efficiency of different equation solvers in cokriging

When the system of equations for cokriging is written in matrix form the sample-sample covariance matrix may be considered either as an mn × mn matrix of scalar entries, where n is the number of sample locations and m is the number of variables, or as an n × n matrix whose entries are m × m matrices. Similarly, the point-sample covariance matrix may be considered as m column vectors or as a single column whose entries are m × m matrices. The formulation in the original program assumed that the submatrix structure should be preserved, but this is not necessary. The scalar matrix formulation allows for the use of a standard Gaussian elimination to reduce the matrix to diagonal form or for reduction to upper triangular form together with back substitution. Both methods result in significant reductions in computing time.

An improved finite element iterative method for three-dimensional fluid flow problems

A new iterative method is presented for the solution of three-dimensional fluid flow problems. The method is based on a preconditioned Uzawa algorithm and uses an incomplete Choleski factorization of the diffusion matrix as preconditioning. Convergence properties and computational efficiency have been investigated. It is shown that the CPU requirements are O(NEQ 1.65) compared to O(NEQ 2.33) for a standard Gaussian elimination, where NEQ is the number of equations. Numerical results are obtained for the excess entrance pressure loss in a sudden contraction for Newtonian and power-law fluids, and compared to available data for code validation.

A Method to Calculate Steady-State Distributions of Large Markov Chains by Aggregating States

This paper develops an efficient iterative algorithm to calculate the steady-state distribution of nearly all irreducible discrete-time Markov chains. Computational experiences suggest that, for large Markovian systems (more than 130 states), the proposed algorithm can be ten times faster than standard Gaussian elimination in finding solutions to an accuracy of 0.1%. The proposed algorithm is developed in three stages. First, we develop a very efficient algorithm for determining steady-state distributions of a restricted class of Markovian systems. A second result establishes a relationship between a general irreducible Markovian system and a system in the restricted class of Markovian systems. Finally, we combine the two results to produce an efficient, iterative algorithm to solve Markov systems. The paper concludes with a discussion of the observed performance of the algorithm.

A boundary element solution to the problem of interacting a.c. fields in parallel conductors

The a.c. field in electrically insulated conductors interact through the surrounding electromagnetic fields. The pertinent field equations reduce to the Helmholtz equation inside each conductor (interior problem), and to the Laplace equation outside the conductors (exterior problem). These equations are transformed to integral equations, with the magnetic vector potential and its normal derivative on the boundaries as unknowns. The integral equations are then approximated by sets of algebraic equations. The interior problem involves only unknowns on the boundary of each conductor, while the exterior problem couples unknows from several conductors. The interior and the exterior problem are coupled through the field continuity conditions. The full set of equations is solved by standard Gaussian elimination. We also show how the total current and the dissipated power within each conductor can be expressed as boundary integrals. Finally, computational results for a sample problem are compared with a finite difference solution.

Inverse matrix representation with one triangular array

We describe a technique that permits the representation of the inverse of a matrix A with only one additional triangular array. Let L∗A = U, with L lower and U upper triangular arrays of order N. Algorithms are presented that use A and L to compute the matrix-vector products A -1∗b and b T∗A -1 with N 2 multiplications and additions. The array L can be computed, with N 3/3 multiplications, with a technique that avoids the computation of U. Standard Gaussian elimination simultaneously computes L and U as follows: start with I∗A = A, where I is the identity matrix; perform identical linear combinations of rows on I and on the right hand side array A; gradually transform I into L and A into U. At an intermediate stage, where A has not yet been fully triangularized, we have L′∗A = U′.L′ and Ú represent one of the pairs of arrays present before each linear combination of rows is performed. The key observation is that we only need two elements of Ú to compute each linear combination of rows of Ĺ. Compute them with a scalar product of the appropriate rows of Ĺ and columns of A. Instead of storing the arrays Ú, recompute their few needed elements whenever necessary.

Application of Vector Processors To Solve Finite Difference Equations

Abstract As computer technology approaches limitations imposed by the speed of light, increased emphasis is placed on vector processors. These have the ability to increase greatly the speed of arithmetic even without improvements in such basic computer characteristics as memory cycle time. This paper deals with solving systems of finite difference equations on the STAR 100 and the CYBER 203, two Control Data Corp. computers with built-in vector processors. Systems of three-dimensional finite difference equations having from 2,000 to 8,000 unknowns were solved by means of Gaussian elimination and line successive overrelaxation (LSOR). On these machines, the D4 Gaussian elimination technique reduced computer time by factors as large as 4.6 relative to standard Gaussian elimination. Vectorization of the D4 code on the STAR 100 reduced computer times relative to scalar results by factors as large as 26, despite nonoptimal coding. LSOR was vectorized successfully with computer time reduction factors of 35 to 43 on the STAR 100. On the CYBER 203, run times were reduced by factors of 45 to 54, relative to the scalar performance of the STAR 100. On an 8,000-block problem, average processing speed for a complete LSOR solution was approximately 25 million floating operations per second (megaflops).