Backward error analysis for linear systems associated with inverses ofH-matrices

In this paper, bounds on the growth factors resulting from Gaussian elimination applied to inverses ofH-matrices are developed and investigated. These bounds are then used in the error analysis for solving linear systemsAx =b whose coefficient matricesA are of this type. For each such system our results show that the Gaussian elimination without pivoting can proceed safely provided that the elements of the inverse of a certainM-matrix (associated with the coefficient matrixA) are not excessively large. We exhibit a particularly satisfactory situation for the special case whenA itself is an inverse of anM-matrix. Part of the first section of this paper is devoted to a discussion on some results of de Boor and Pinkus for the stability of triangular factorizations of systemsAx =b, whereA is a nonsingular totally nonnegative matrix, and to the explanation of why the analysis of de Boor and Pinkus is not applicable to the case when the coefficient matrixA is an inverse of anM-matrix.Research supported in part by NSF Grant MCS-8102114.Research supported in part by the U.S. Army Research Office under contract No. DAAG-29-81-K-0132 and in part by NSF Grant MCS-8219500.     

Optimal preconditioners of a given sparsity pattern

We consider the problem of finding the symmetric positive definite preconditionerM of a given form- e.g., having nonzero elements only in specified positions — which minimizes the ratio of the largest to smallest eigenvalue ofM ^–1 A, for a given symmetric positive definitive matrixA. We show how this problem can be expressed as one of minimizing a convex function and how an optimization code can be used to solve the problem numerically. Results are presented showing optimal preconditioners of various sparsity patterns and comparing these to preconditioners that have been proposed in the literature. Several conjectures are made, based on results from the optimization code.This work was supported by the Advanced Research Projects Agency of the Department of Defense under contract F49620-87-C-0065 and by the Applied Mathematical Sciences Program of the U.S. Department of Energy under contract DE-AC02-76ER03077.     

The Euler Maclaurin expansion for the Cauchy Principal Value integral

Summary The classical Euler Maclaurin Summation Formula expresses the difference between a definite integral over [0, 1] and its approximation using the trapezoidal rule with step lengthh=1/m as an asymptotic expansion in powers ofh together with a remainder term. Many variants of this exist some of which form the basis of extrapolation methods such as Romberg Integration. in this paper a variant in which the integral is a Cauchy Principal Value integral is derived. The corresponding variant of the Fourier Coefficient Asymptotic Expansion is also derived. The possible role of the former in numerical quadrature is discussed.This work was supported by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy, under contract W-31-109-Eng-38     

Adaptive polynomial preconditioning for hermitian indefinite linear systems

This paper explores the use of polynomial preconditioned CG methods for hermitian indefinite linear systems,Ax=b. Polynomial preconditioning is attractive for several reasons. First, it is well-suited to vector and/or parallel architectures. It is also easy to employ, requiring only matrix-vector multiplication and vector addition. To obtain an optimum polynomial preconditioner we solve a minimax approximation problem. The preconditioning polynomial,C(), is optimum in that it minimizes a bound on the condition number of the preconditioned matrix,C(A)A. We also characterize the behavior of this minimax polynomial, which makes possible a thorough understanding of the associated CG methods. This characterization is also essential to the development of an adaptive procedure for dynamically determining the optimum polynomial preconditioner. Finally, we demonstrate the effectiveness of polynomial preconditioning in a variety of numerical experiments on a Cray X-MP/48. Our results suggest that high degree (20–50) polynomials are usually best.This research was supported in part by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Dept. of Energy, by Lawrence Livermore National Laboratory under contract W-7405-ENG-48.This research was supported in part by the Dept. of Energy and the National Science Foundation under grant DMS 8704169.This research was supported in part by U.S. Dept. of Energy grant DEFG02-87ER25026 and by National Science Foundation grant DMS 8703226.     

A note on the computation of an orthonormal basis for the null space of a matrix

A highly regarded method to obtain an orthonormal basis,Z, for the null space of a matrixA ^T is theQR decomposition ofA, whereQ is the product of Householder matrices. In several optimization contextsA(x) varies continuously withx and it is desirable thatZ(x) vary continuously also. In this note we demonstrate that thestandard implementation of theQR decomposition doesnot yield an orthonormal basisZ(x) whose elements vary continuously withx. We suggest three possible remedies. Work supported in part by the Applied Mathematical Sciences Research Program (KC-04-02) of the Office of Energy Research of the U.S. Department of Energy under Contract W-31-109-Eng-38.     

A reinforcement-depletion urn model: A contiguity case

Summary An urn contains balls ofs different colors. The problem of the reinforcement of a specified color and random depletion of balls has been considered by Bernard (1977,Bull. Math. Biol.,39, 463–470) and Shenton (1981,Bull. Math. Biol.,43, 327–340), (1983,Bull. Math. Biol.,45, 1–9). Here we consider a special relation between a reinforcement and depletion, leading to a hypergeometric distribution. Research sponsored in part by the Applied Mathematical Sciences Research Program, Office of Energy Research, U.S. Department of Energy under contract DE-AC05-840R21400 with the Martin Marietta Energy Systems, Inc.     

Properties of a representation of a basis for the null space

Given a rectangular matrixA(x) that depends on the independent variablesx, many constrained optimization methods involve computations withZ(x), a matrix whose columns form a basis for the null space ofA ^T(x). WhenA is evaluated at a given point, it is well known that a suitableZ (satisfyingA ^T Z = 0) can be obtained from standard matrix factorizations. However, Coleman and Sorensen have recently shown that standard orthogonal factorization methods may produce orthogonal bases that do not vary continuously withx; they also suggest several techniques for adapting these schemes so as to ensure continuity ofZ in the neighborhood of a given point.This paper is an extension of an earlier note that defines the procedure for computingZ. Here, we first describe howZ can be obtained byupdating an explicit QR factorization with Householder transformations. The properties of this representation ofZ with respect to perturbations inA are discussed, including explicit bounds on the change inZ. We then introduceregularized Householder transformations, and show that their use implies continuity of the full matrixQ. The convergence ofZ andQ under appropriate assumptions is then proved. Finally, we indicate why the chosen form ofZ is convenient in certain methods for nonlinearly constrained optimization.The research of the Stanford authors was supported by the U.S. Department of Energy Contract DE-AM03-76SF00326, PA No. DE-AT03-76ER72018; the National Science Foundation Grants MCS-7926009 and ECS-8312142; the Office of Naval Research Contract N00014-75-C-0267; and the U.S. Army Research Office Contract DAAG29-84-K-0156.The research of G.W. Stewart was supported by the Air Force Office of Scientific Research Contract AFOSR-82-0078.     

Ergodic automorphisms ofT n are Bernoulli shifts

An automorphism of then-dimensional torusT ⁿ, none of whose eigenvalues is a root of unity includes on the canonical measure space ofT ⁿ a measure preserving transformation which is isomorphic to a Bernoulli shift. This research was supported in part by the National Science Foundation grant GP-18884 and by the European Research Office of the U.S. Army contract DAJA-37-70-C-0701.     

Parallel sparse Gaussian elimination with partial pivoting

In this article, we consider the factorization of a sparse nonsymmetric matrix using Gaussian elimination with partial pivoting on a multiprocessor having a globally-shared memory. The parallel algorithm makes use of a static data structure developed by George, Liu and Ng in [17]. Some numerical experiments on a Sequent Balance 8000 are presented to demonstrate the efficiency of the parallel implementation.Research supported in part by the Applied Mathematical Sciences Research Program, Office of Energy Research, U.S. Department of Energy under contract DE-AC05-84OR21400 with Martin Marietta Energy Systems, Inc.     

Multiple discrete solutions of the incompressible steady-state Navier-Stokes equations

Summary This paper discusses the computation of multiple solutions of various discretizations of the steady state incompressible Navier-Stokes equations. Solution paths (,R) satisfying the discrete system of equationsH(,R)=0, where represents the discrete flow field andR is the Reynolds number, are computed using a pseudo arc-length continuation procedure. For flows over the back end of an axially symmetric body with a cusped tail in a coaxial circular cylinder, the solution paths often exhibit hairpin turning points. Dependence of the paths on the mesh spacing and various selections for the discretizations of the convective and diffusive terms are presented.Dedicated to Professor Ivo Babuka on the occasion of his 60th birthdayThis research was supported in part by the Naval Surface Warfare Center Independent Research Fund, the Naval Sea Systems Command, and the Office for Naval Research under Contracts N001484WR24012 and N0001486R24209, and was performed in part under the auspices of the U.S. Department of Energy by the Lawrence Livermore National Laboratory under contract number W-7405-ENG-48. Partial support under contract Number W-7405-ENG-48 was provided by the Applied Mathematical Sciences Program of the Office of Energy Research     

Approximation algorithms for indefinite quadratic programming

We consider-approximation schemes for indefinite quadratic programming. We argue that such an approximation can be found in polynomial time for fixed andt, wheret denotes the number of negative eigenvalues of the quadratic term. Our algorithm is polynomial in 1/ for fixedt, and exponential int for fixed. We next look at the special case of knapsack problems, showing that a more efficient (polynomial int) approximation algorithm exists.Part of this work was done while the author was visiting Sandia National Laboratories, Albuquerque, New Mexico, supported by the U.S. Department of Energy under contract DE-AC04-76DP00789. Part of this work was also supported by the Applied Mathematical Sciences Program (KC-04-02) of the Office of Energy Research of the U.S. Department of Energy under grant DE-FG02-86ER25013.A000 and in part by the National Science Foundation, the Air Force Office of Scientific Research, and the Office of Naval Research, through NSF grant DMS 8920550.     

On computingINV block preconditionings for the conjugate gradient method

The INV(k) and MINV(k) block preconditionings for the conjugate gradient method require generation of selected elements of the inverses of symmetric matrices of bandwidth 2k+1. Generalizing the previously describedk=1 (tridiagonal) case tok=2, explicit expressions for the inverse elements of a symmetric pentadiagonal matrix in terms of Green's matrix of rank two are given. These expressions are found to be seriously ill-conditioned; hence alternative computational algorithms for the inverse elements must be used. Behavior of thek=1 andk=2 preconditionings are compared for some discretized elliptic partial differential equation test problems in two dimensions.Presented by the first author at the Joint U.S.-Scandinavian Symposium on Scientific Computing and Mathematical Modeling, January 1985, in honor of Germund Dahlquist on the occasion of his 60th birthday. This work was supported in part by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy under contract DE AC03-76SF00098.     

On optimal triangular meshes for minimizing the gradient error

Summary Construction of optimal triangular meshes for controlling the errors in gradient estimation for piecewise linear interpolation of data functions in the plane is discussed. Using an appropriate linear coordinate transformation, rigorously optimal meshes for controlling the error in quadratic data functions are constructed. It is shown that the transformation can be generated as a curvilinear coordinate transformation for anyC data function with nonsingular Hessian matrix. Using this transformation, a construction of nearly optimal meshes for general data functions is described and the error equilibration properties of these meshes discussed. In particular, it is shown that equilibration of errors is not a sufficient condition for optimality. A comparison of meshes generated under several different criteria is made, and their equilibrating properties illustrated.This work was supported by the Natural Sciences and Engineering Research Council of Canada, by the Information Technology Research Centre, which is funded by the Province of Ontario, by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy under contract DE-AC05-84OR21400 with Martin Marietta Energy Systems, Inc., and through an appointment to the U.S. Department of Energy Postgraduate Research Program administered by Oak Ridge Associated Universities     

An analysis of reduced Hessian methods for constrained optimization

We study the convergence properties of reduced Hessian successive quadratic programming for equality constrained optimization. The method uses a backtracking line search, and updates an approximation to the reduced Hessian of the Lagrangian by means of the BFGS formula. Two merit functions are considered for the line search: the ₁ function and the Fletcher exact penalty function. We give conditions under which local and superlinear convergence is obtained, and also prove a global convergence result. The analysis allows the initial reduced Hessian approximation to be any positive definite matrix, and does not assume that the iterates converge, or that the matrices are bounded. The effects of a second order correction step, a watchdog procedure and of the choice of null space basis are considered. This work can be seen as an extension to reduced Hessian methods of the well known results of Powell (1976) for unconstrained optimization.This author was supported, in part, by National Science Foundation grant CCR-8702403, Air Force Office of Scientific Research grant AFOSR-85-0251, and Army Research Office contract DAAL03-88-K-0086.This author was supported by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy, under contracts W-31-109-Eng-38 and DE FG02-87ER25047, and by National Science Foundation Grant No. DCR-86-02071.     

On the zeros and poles of Padé approximants toe z

Summary In this paper, we study the location of the zeros and poles of general Padé approximats toe ^z. The location of these zeros and poles is useful in the analysis of stability for related numerical methods for solving systems of ordinary differential equations.Research supported in part by the Air Force Office of Scientific Research under Grant AFOSR-74-2688, and by the University of South Fla. Research Council.Research supported in part by the Air Force Office of Scientific Research under Grant AFOSR-74-2729, and by the Atomic Energy Commission under Grant AT(11-1)-2075.     

An algorithm to compute a sparse basis of the null space

Summary LetA be a realm×n matrix with full row rankm. In many algorithms in engineering and science, such as the force method in structural analysis, the dual variable method for the Navier-Stokes equations or more generally null space methods in quadratic programming, it is necessary to compute a basis matrixB for the null space ofA. HereB isn×r, r=n–m, of rankr, withAB=0. In many instancesA is large and sparse and often banded. The purpose of this paper is to describe and test a variation of a method originally suggested by Topcu and called the turnback algorithm for computing a banded basis matrixB. Two implementations of the algorithm are given, one using Gaussian elimination and the other using orthogonal factorization by Givens rotations. The FORTRAN software was executed on an IBM 3081 computer with an FPS-164 attached array processor at the Triangle Universities Computing Center and on a CYBER 205 vector computer. Test results on a variety of structural analysis problems including two- and three-dimensional frames, plane stress, plate bending and mixed finite element problems are discussed. These results indicate that both implementations of the algorithm yielded a well-conditioned, banded, basis matrixB whenA is well-conditioned. However, the orthogonal implementation yielded a better conditionedB for large, illconditioned problems.The research by these authors was supported by the U.S. Air Force under grant No. AFOSR-83-0255 and by the National Science Foundation under grant No. MCS-82-19500The research by these authors was supported by the Applied Mathematical Sciences Program of the U.S. Department of Energy, under contract to Martin Marietta Energy Systems, Inc.     

Spectral Tau approximation of the two-dimensional stokes problem

Summary We analyse the Spectral Tau method for the approximation of the Stokes system on a square with Dirichlet boundary conditions. We provide an error estimate, in the norm of the Sobolev spaceH ^s, for the approximation of a divergence free vector field with polynomial divergence free vector fields. We apply this result to prove some convergence estimates for the solution of the discrete Stokes problem.This work has been partially supported by the U.S. Army through its European Research Office under contract No. DAJA-84-C 0035     

On the solution of large,structured linear complementarity problems: The block partitioned case

This paper delineates the underlying theory of an efficient method for solving a class of specially-structured linear complementarity problems of potentially very large size. We propose a block-iterative method in which the subproblems are linear complementarity problems. Problems of the type considered here arise in the process of making discrete approximations to differential equations in the presence of special side conditions. This problem source is exemplified by the free boundary problem for finite-length journal bearings. Some of the authors' computational experience with the method is presented.Research partially supported by the Office of Naval Research under contract N-00014-67-A-0112-0011; U.S. Atomic Energy Commission Contract AT (04-3)-326 PA No. 18.Research partially supported by U.S. Atomic Energy Commission Contract AT(04-3)-326 PA No. 30; and National Science Foundation Grant GJ 35135X.     

On a computational method for the second fundamental tensor and its application to bifurcation problems

Summary An algorithm is presented for the computation of the second fundamental tensorV of a Riemannian submanifoldM ofR ⁿ . FromV the riemann curvature tensor ofM is easily obtained. Moreover,V has a close relation to the second derivative of certain functionals onM which, in turn, provides a powerful new tool for the computational determination of multiple bifurcation directions. Frequently, in applications, thed-dimensional manifoldM is defined implicitly as the zero set of a submersionF onR ⁿ . In this case, the principal cost of the algorithm for computingV(p) at a given pointpM involves only the decomposition of the JacobianDF(p) ofF atp and the projection ofd(d+1) neighboring points ontoM by means of a local iterative process usingDF(p). Several numerical examples are given which show the efficiency and dependability of the method.Dedicated to R. S. Varga on the occasion of his sixtieth birthdayThis work was in part supported by the National Science Foundation (DCR-8309926) and the Office of Naval Research (N-00014-80-C09455). The second author began some of the work while visiting the University of Heidelberg/Germany as an Alexander von Humboldt Senior U.S. Scientist     

On the zeros and poles of Padé approximants toe z . III

Summary In this paper, we continue our study of the location of the zeros and poles of general Padé approximants toe ^z . We state and prove here new results for the asymptotic location of the normalized zeros and poles for sequences of Padé approximants toe ^z , and for the asymptotic location of the normalized zeros for the associated Padé remainders toe ^z . In so doing, we obtain new results for nontrivial zeros of Whittaker functions, and also generalize earlier results of Szegö and Olver.Research supported in part by the Air Force Office of Scientific Research under Grant AFOSR-74-2688Research supported in part by the Air Force Office of Scientific Research under Grant AFOSR-74-2729, and by the Energy Research and Development Administration (ERDA) under Grant EY-76-S-02-2075