Computational Solutions of Matrix Problems Over an Integral Domain

Recent methods for handling matrix problems over an integraldomain are investigated from a unifying point of view. Emphasizedare symbolic matrix inversion and numerically exact methodsfor solving Ax = b. New proofs are given for the theory of themultistep method. A proof for the existence and an algorithmfor the exact solution of Tx = b, where T is a finite Toeplitzmatrix, is given. This algorithm reduces the number of requiredsingle precision multiplications by a factor of order n overthe corresponding Gaussian elimination method. The use of residuearithmetic is enhanced by a new termination process. The matrixinversion problem with elements in the ring of polynomials isreduced to operations over a Galois field. It is shown thatinterpolation methods are equivalent to congruence methods withlinear modulus and that the Chinese remainder theorem over GF(x-p_k)is the Lagrange interpolation formula. With regard to the numerical problem of exact matrix inversion,the One- and Two-step Elimination methods are critically comparedwith the methods using modular or residue arithmetic. Formulasfor estimating maximum requirements for storage and timing ofthe salient parts of the algorithms are developed. The resultsof a series of recent tests, using existing codes, standardmatrices and matrices with random elements are reported andsummarized in tabular form. The paper concludes that the two-stepelimination method be used for the inversion problem of numericmatrices, and in particular when a black-box approach to thematrix inversion problem is attempted such as in commercialtime sharing systems. It is recommended that the inversion problemof matrices with elements over the polynomial ring be reducedto the numeric inversion problem with subsequent interpolation.An extensive Reference list is added.     

On the Complexity of Computing Error Bounds

We consider the cost of estimating an error bound for the computed solution of a system of linear equations, i.e., estimating the norm of a matrix inverse. Under some technical assumptions we show that computing even a coarse error bound for the solution of a triangular system of equations costs at least as much as testing whether the product of two matrices is zero. The complexity of the latter problem is in turn conjectured to be the same as matrix multiplication, matrix inversion, etc. Since most error bounds in practical use have much lower complexity, this means they should sometimes exhibit large errors. In particular, it is shown that condition estimators that: (1) perform at least one operation on each matrix entry; and (2) are asymptotically faster than any zero tester, must sometimes over or underestimate the inverse norm by a factor of at least , where n is the dimension of the input matrix, k is the bitsize, and where either or grows faster than any polynomial in n . Our results hold for the RAM model with bit complexity, as well as computations over rational and algebraic numbers, but not real or complex numbers. Our results also extend to estimating error bounds or condition numbers for other linear algebra problems such as computing eigenvectors. September 10, 1999. Final version received: August 23, 2000.     

Error analysis of algorithms for matrix multiplication and triangular decomposition using Winograd's identity

Summary The number of multiplications required for matrix multiplication, for the triangular decomposition of a matrix with partial pivoting, and for the Cholesky decomposition of a positive definite symmetric matrix, can be roughly halved if Winograd's identity is used to compute the inner products involved. Floating-point error bounds for these algorithms are shown to be comparable to those for the normal methods provided that care is taken with scaling.     

Uniform Convergence of Variational Methods for Elliptic Problems in N Dimensions

Earlier work in variational methods in the plane is extendedto N dimensions. Smoothness theorems are constructed and usedot convert means convergence to uniform convergence, establishingpointwise error bounds. A constrained variational method isused so that a priori bounds for the derivatives of the solutionmay be imposed on the approximating functions. In the non-linearproblem of compressible flow such a bound is provided by anassumption that the solution is subsonic.     

Fast Structured Matrix Computations: Tensor Rank and Cohn–Umans Method

We discuss a generalization of the Cohn–Umans method, a potent technique developed for studying the bilinear complexity of matrix multiplication by embedding matrices into an appropriate group algebra. We investigate how the Cohn–Umans method may be used for bilinear operations other than matrix multiplication, with algebras other than group algebras, and we relate it to Strassen’s tensor rank approach, the traditional framework for investigating bilinear complexity. To demonstrate the utility of the generalized method, we apply it to find the fastest algorithms for forming structured matrix–vector product, the basic operation underlying iterative algorithms for structured matrices. The structures we study include Toeplitz, Hankel, circulant, symmetric, skew-symmetric, f-circulant, block Toeplitz–Toeplitz block, triangular Toeplitz matrices, Toeplitz-plus-Hankel, sparse/banded/triangular. Except for the case of skew-symmetric matrices, for which we have only upper bounds, the algorithms derived using the generalized Cohn–Umans method in all other instances are the fastest possible in the sense of having minimum bilinear complexity. We also apply this framework to a few other bilinear operations including matrix–matrix, commutator, simultaneous matrix products, and briefly discuss the relation between tensor nuclear norm and numerical stability.     

Bounds for the Singular Values of a Matrix

A lower-bound theorem is developed for the singular values ofa matrix A (and therefore for the eigenvalues of B = A^HA). Itis found that there is always a unique column scaling of A whichproduces the optimum bound. However, sharper bounds still maysometimes be obtained by taking advantage of matrix partitioning.It is shown that the resulting bounds may often (but not always)be better than those obtained by applying Gerschgorin's theoremto B. The equivalent upper-bound theorem is found to be weak.     

Error Bounds for Polynomial Evaluation and Complex Arithmetic

Explicit and computable strict bounds of a posteriori type areconstructed for the evaluation of a polynomial and its derivativesby nested multiplication, using floating-point arithmetic. Thepolynomial may be real or complex, and may contain errors inits coefficients and argument. Some new error bounds for arithmetic operations with complexnumbers are included.     

A new non-conforming Petrov-Galerkin finite-element method with triangular elements for a singularly perturbed advection-diffusion problem

A new non-conforming exponentially fitted Petrov-Galerkin finite-elementmethod based on Delaunay triangulation and its Dirichlet tessellationis constructed for the numerical solution of singularly perturbedstationary advectiondiffusion problems with a singular perturbationparameter . The method is analyzed mathematically and its stabilityis shown to be independent of . The error estimate in an -independentdiscrete energy norm for the approximate solution is shown todepend on first-order seminorms of the flux and the zero-orderterm of the equation, the sup norm of the exact solution, thefirst-order seminorm of the coefficient of the advection term,and the approximation error of the inhomogeneous term. Sincethe first two seminorms are not bounded uniformly in , the -uniformconvergence of the method still remains an open question. Noassumption is required that the angles in the triangulationare all acute. Since the system matrix for this method is identicalto that for the exponentially fitted box method, the theoreticalresults also provide an analysis of that box method. The newmethod also contains the central-difference and upwind methodsas two limiting cases. It can be regarded as a weighted finite-differencemethod on a triangular mesh. Numerical results are presentedto show the superior performance of the method in comparisonwith the upwind and central-difference methods for a small increasein the computation cost. Present address: School of Mathematics, The University of NewSouth Wales, Kensington, NSW 2033, Australia.     

On the Accuracy of Quadrature Laplace Transform Inversion and the Solution of State Space Equations

A quadrature formula is shown to be an approximation of thepower-series method of inverting Laplace transforms. This togetherwith the properties of the constants derived from a power-seriesexpansion of the Pad? approximation to exp (s) yield an importantupper limit on t which is quite sharp in determining the breakdownpoint up to and after which the approximation is accurate andinaccurate respectively. The solution of state space equationsusing the quadrature inversion formula is also discussed.     

A fully discrete stabilized finite-element method for the time-dependent Navier-Stokes problem

A fully discrete stabilized finite-element method is presentedfor the two-dimensional time-dependent Navier–Stokes problem.The spatial discretization is based on a finite-element spacepair (X_h, M_h) for the approximation of the velocity and thepressure, constructed by using the Q₁ – P₀ quadrilateralelement or the P₁ – P₀ triangular element; the time discretizationis based on the Euler semi-implicit scheme. It is shown thatthe proposed fully discrete stabilized finite-element methodresults in the optimal order error bounds for the velocity andthe pressure.     

Partition and composition matrices

This paper introduces two matrix analogues for set partitions. A composition matrix on a finite set X is an upper triangular matrix whose entries partition X, and for which there are no rows or columns containing only empty sets. A partition matrix is a composition matrix in which an order is placed on where entries may appear relative to one-another.We show that partition matrices are in one-to-one correspondence with inversion tables. Non-decreasing inversion tables are shown to correspond to partition matrices with a row ordering relation. Partition matrices which are s-diagonal are classified in terms of inversion tables. Bidiagonal partition matrices are enumerated using the transfer-matrix method and are equinumerous with permutations which are sortable by two pop-stacks in parallel.We show that composition matrices on X are in one-to-one correspondence with (2+2)-free posets on X. Also, composition matrices whose rows satisfy a column-ordering relation are shown to be in one-to-one correspondence with parking functions. Finally, we show that pairs of ascent sequences and permutations are in one-to-one correspondence with (2+2)-free posets whose elements are the cycles of a permutation, and use this relation to give an expression for the number of (2+2)-free posets on {1,…,n}.     

A Posteriori Error Bounds for Linear Elliptic Equations of Potential Type

Second-order Linear elliptic partial differential equations of potential type with Dirichlet (type 1) or Neumann (type II)boundary conditions on a simply-connected two-dimensional domainare considered. Conjugate problems, that is, a pair of one type1 and one type II problem, are introduced along with an auxiliaryeppiptic system of two equations in such a way that the energiesof the given problem, its conjugate problem, and the auxiliarysystem add to a known constant. There result two-sided boundsfor the energy of the given problem and, as a consequence, aposteriori error bounds for the norm of the difference of anapproximate solution and the exact solution of the problem.A method by which the amount of computation required to obtainthe a posteriori error bounds can be almost halved in many casesof practical interest is given. A posteriori error bounds forapproximate solutions of the auxiliary system are also given.     

Accuracy of the Kogbetliantz method for scaled diagonally dominant triangular matrices

The paper proves that Kogbetliantz method computes all singular values of a scaled diagonally dominant triangular matrix, which can be well scaled from both sides symmetrically, to high relative accuracy. Special attention is paid to deriving sharp accuracy bounds for one step, one batch and one sweep of the method. By a simple numerical test it is shown that the methods based on bidiagonalization are generally not accurate on that class of well-behaved matrices.     

On the Application of a Generalization of Toeplitz Matrices to the Numerical Solution of Integral Equations with Weakly Singular Convolution Kernels

This paper discusses the numerical solution by product integrationof weakly singular Fredholm integral equations of the secondkind with symmetric difference kernels. The product integrationmethod uses a piecewise polynomial, in general, at most, continuousat the knots. A main result of the paper is to show that, owingto the difference kernel, a highly patterned linear system ofequations arises if the knots are equally spaced. Specificallythe order-N coefficient matrix is block-Toeplitz or a generalization,and centrosymmetric. An algorithm to solve this linear systemin O(N²) operations is presented. Unfortunately the asymptotic rate of convergence of the productintegration solution is limited if a uniform mesh is used. Amethod to improve the rate of convergence, while retaining apatterned coefficient matrix, is described. This method involvessubtraction of the singularities in the solution induced bythe weakly singular kernel. A higher rate of convergence forthe modified product integration method is shown. These resultsare illustrated by application to an integral equation arisingin outdoor sound propagation.     

On the Structure of the Moving Finite-element Equations

Address from 1st April 1985, School of Mathematics, Universityof Bristol, University Walk, Bristol BS8 1TW. The morning finite-element method for evolutionary partial differentialequations leads to a coupled non-linear system of ordinary differentialequations in time, with a coefficien matrix A, say, for thetime derivaties, We show for linear elements in any number ofdimensions, A can be written in the form M^TCM, where the matrixC depends solely on the mesh geometry and the matrix M on thegradient of the section, As a simple consequence we show thatA is singular only in the cases (i) element degeneracy () and (ii) collinearity of nodes (M not out of fullrank). We give constructions for the inversion of A in all cases. In one dimension, if A is non-singular, it has a simple explicitinverse. If A is singular we replace it by reduced matrix A^*.It can be shown that every case the spectral radius of the Jacobiiteration matrix ia ?and that A or A^* can be efficiently invertedby conjugate gradient methods. Finally, we discuss the applicability of these arguments tosystem of equations in any number of dimensions.     

A Mixed Stochastic-Deterministic Method for Systems of Linear Equations

A Monte Carlo procedure is presented with a view to solvingthe system of linear equations x = Hx+b, where H is a substochasticmatrix, by inversion of some principal submatrices of H. Themethod is described on the basis of the directed graph representationof matrix H and it is shown how it is possible to constructa sampling model which, compared with the usual Monte Carlomethod, avoids transitions between the rows (or columns) ofthese submatrices. Systems with more general matrices are alsoconsidered.     

Some Condition-Number bounds for the Gaussian Elimination Process

It is well-known that Gaussian Elimination is equivalent toLU decomposition. This paper shows that if Gaussian Eliminationis applied to a well-conditioned matrix, and if an increasein element size is observed during the elimination process,then both the triangular factors will be badly conditioned.It is further shown that the effect of partial pivoting is toplace an upper bound upon the condition number of the lowertriangular factor.     

Toeplitz型矩阵的逆矩阵的快速三角分解算法

针对有关“型”矩阵的三角分解问题 ,提出了一种 Toeplitz型矩阵的逆矩阵的快速三角分解算法 .首先假设给定 n阶非奇异矩阵 A,利用一组线性方程组的解 ,得到 A- 1的一个递推关系式 ,进而利用该关系式得到 A- 1的一种三角分解表达式 ,然后从 Toeplitz型矩阵的特殊结构出发 ,利用上述定理的结论 ,给出了Toeplitz型矩阵的逆矩阵的一种快速三角分解算法 ,算法所需运算量为 O( mn2 ) .最后 ,数值计算表明该算法的可靠性 .     

Non-autonomous Riccati-type matrix differential equations: existence interval, construction of continuous numerical solutions and error bounds

In this paper, non-autonomous Riccati-type matrix differentialequations with twice continuously differentiable coefficientsare considered. First of all an existence interval in termsof the data for an initial-value problem is determined. Secondly,a discrete numerical solution is constructed at a net of pointsof the predetermined existence interval using linear one-stepmatrix methods. Finally, using linear B-spline matrix functions,a continuous numerical solution with error bounds is obtained.Given an admissible error >0 we construct a continuous numericalsolution whose error is uniformly upper bounded by in the existencedomain.     

A Fast Direct Method for the Least Squares Solution of Slightly Overdetermined Sets of Linear Equations

Noble (1969) has described a method for the solution of N+Mlinear equations in N unknowns, which is based on an initialpartitioning of the matrix A, and which requires only the solutionof square sets of equations. He assumed rank (A) = N. We describehere an efficient implementation of Noble's method, and showthat it generalizes in a simple way to cover also rank deficientproblems. In the common case that the equation is only slightlyoverdetermined (M << N) the resulting algorithm is muchfaster than the standard methods based on M.G.S. or Householderreduction of A, or on the normal equations, and has a very similaroperation count to the algorithm of Cline (1973). Slightly overdetermined systems arise from Galerkin methodsfor non-Hermitian partial differential equations. In these systems,rank (A) = N and advantage can be taken of the structure ofthe matrix A to yield a least squares solution in (N²) operations.