Showing papers in "Mathematics of Computation in 1970"

Conditioning of Quasi-Newton Methods for Function Minimization

Abstract: Quasi-Newton methods accelerate the steepest-descent technique for function minimization by using computational history to generate a sequence of approximations to the inverse of the Hessian matrix. This paper presents a class of approximating matrices as a function of a scalar parameter. The problem of optimal conditioning of these matrices under an appropriate norm as a function of the scalar parameter is investigated. A set of computational results verifies the superiority of the new methods arising from conditioning considerations to known methods.

A family of variable-metric methods derived by variational means

Abstract: A new rank-two variable-metric method is derived using Greenstadt's variational approach [Math. Comp., this issue]. Like the Davidon-Fletcher-Powell (DFP) variable-metric method, the new method preserves the positive-definiteness of the approximating matrix. Together with Greenstadt's method, the new method gives rise to a one-parameter family of variable-metric methods that includes the DFP and rank-one methods as special cases. It is equivalent to Broyden's one-parameter family [Math. Comp., v. 21, 1967, pp. 368-381]. Choices for the inverse of the weighting matrix in the variational approach are given that lead to the derivation of the DFP and rank-one methods directly. In the preceding paper [6], Greenstadt derives two variable-metric methods, using a classical variational approach. Specifically, two iterative formulas are developed for updating the matrix Hk, (i.e., the inverse of the variable metric), where Hk is an approximation to the inverse Hessian G-'(Xk) of the function being minimized.* Using the iteration formula Hk+1 = Hk + Ek to provide revised estimates to the inverse Hessian at each step, Greenstadt solves for the correction term Ek that minimizes the norm N(Ek) = Tr (WEkWEkJ) subject to the conditions

Factoring polynomials over large finite fields

Abstract: This paper reviews some of the known algorithms for factoring polynomials over finite fields and presents a new deterministic procedure for reducing the problem of factoring an arbitrary polynomial over the Galois field GF(p m) to the problem of finding the roots in GF(p) of certain other polynomials over GF(p). The amount of computation and the storage space required by these algorithms are algebraic in both the degree of the polynomial to be factored and the logarithm of the order of the finite field. Certain observations on the application of these methods to the factorization of polynomials over the rational integers are also included.

Solution of Vandermonde Systems of Equations

Abstract: We obtain in this paper a considerable improvement over a method developed earlier by Ballester and Pereyra for the solution of systems of linear equations with Vandermonde matrices of coefficients. This is achieved by observing that a part of the earlier algorithm is equivalent to Newton's interpolation method. This allows also to produce a progressive algorithm which is significantly more efficient than previous available methods. Algol-60 programs and numerical results are included. Confluent Vandermonde systems are also briefly discussed.

Triangular elements in the finite element method

Abstract: For a plane polygonal domain a and a corresponding (general) triangulation we define classes of functions pm(x, y) which are polynomials on each triangle and which are in Cm)(Q) and also belong to the Sobolev space Wn'"'1(Q). Approximation theoretic properties are proved concerning these functions. These results are then applied to the approximate solution of arbitrary-order elliptic boundary value problems by the Galerkin method. Estimates for the error are given. The case of second-order problems is discussed in conjunction with special choices of approximating polynomials.

Modification of a quasi-Newton method for nonlinear equations with a sparse Jacobian

Abstract: For solving large systems of nonlinear equations by quasi-Newton methods it may often be preferable to store an approximation to the Jacobian rather than an approximation to the inverse Jacobian The main reason is that when the Jacobian is sparse and the locations of the zeroes are known, the updating procedure can be made more efficient for the approximate Jacobian than for the approximate inverse Jacobian

On the construction of Gaussian quadrature rules from modified moments.

Abstract: : Given a weight function omega(x) on (alpha, beta), and a system of polynomials (p sub k)(x), k = 0 to infinity, with degree p sub k (x) = k, we consider the problem of constructing Gaussian quadrature rules from 'modified moments'. Classical procedures take p sub k (x) = x, but suffer from progressive ill-conditioning as n increases. A more recent procedure, due to Sack and Donovan, takes for p sub k (x) a system of (classical) orthogonal polynomials. The problem is then remarkably well-conditioned, at least for finite intervals (alpha, beta). In support of this observation, we obtain upper bounds for the respective asymptotic condition number. In special cases, these bounds grow like a fixed power of n. We also derive an algorithm for solving the problem considered which generalizes one due to Golub and Welsch. Finally, some numerical examples are presented.

Optimal conditioning of quasi-Newton methods

Abstract: Quasi-Newton methods accelerate gradient methods for minimizing a function by approximating the inverse Hessian matrix of the function. Several papers in recent literature have dealt with the generation of classes of approximating matrices as a function of a scalar parameter. This paper derives necessary and sufficient conditions on the range of one such parameter to guarantee stability of the method. It further shows that the parameter effects only the length, not the direction, of the search vector at each step, and uses this result to derive several computational algorithms. The algorithms are evaluated on a series of test problems.

The convergence of single-rank quasi-Newton methods

Abstract: Analyses of the convergence properties of general quasi-Newton methods are presented, particular attention being paid to how the approximate solutions and the itera- tion matrices approach their final values. It is further shown that when Broyden's algorithm is applied to linear systems, the error norms are majorised by a superlinearly convergent sequence of an unusual kind.

Variations on Variable-Metric Methods

Abstract: In unconstrained minimization of a functions, the method of Davidon-Fletcher- Powell (a "variable-metric" method) enables the inverse of the Hessian H of f to be ap- proximated stepwise, using only values of the gradient of f. It is shown here that, by solv- ing a certain variational problem, formulas for the successive corrections to H can be de- rived which closely resemble Davidon's. A symmetric correction matrix is sought which minimizes a weighted Euclidean norm, and also satisfies the "DFP condition." Numerical tests are described, comparing the performance (on four "standard" test functions) of two variationally-derived formulas with Davidon's. A proof by Y. Bard, modelled on Fletcher and Powell's, showing that the new formulas give the exact H after N steps, is included in an appendix. 1. The DFP Method. The class of gradient methods for finding the uncon- strained minimum of a function f(x)* in which the direction Sk of the next iterative step from Xk to Xk+l is computed from a formula such as:

Symmetric quadrature formulae for simplexes

Abstract: Symmetric interpolation polynomials are defined for N-dimensional simplexes with the aid of a symmetric coordinate notation. These polynomials are used to produce symmetric interpolatory quadrature formulae of arbitrary degree of precision over simplexes of arbitrary dimensionality. Tabulated values of weight coefficients are given for triangles and tetrahedra. 1. Introduction. Interpolative quadrature formulae for N-dimensional simplexes have been given by various authors, e.g., Stroud (1) or Hammer and Stroud (2). Their principal attraction in applications lies in the fact that multidimensional regions of integration can often be closely approximated by unions of simplexes. The object of this paper is to show that, for any N, it is possible to define quadrature formulae of any degree of precision n, symmetric in the sense of Hammer, Marlowe, and Stroud (3); and to give a straightforward procedure for finding the weights and node locations. Although the resulting quadrature formulae are not efficient in the sense of (3), they possess the advantage of being very convenient computationally, and can be generated easily for any reasonable values of N and n. They represent a natural generalization of the Newton-Cotes formulae to the N-dimensional case, and include the latter for N = 1.

The simscript ii programming language

Abstract: : A user's and programmer's manual for SIMSCRIPT II that requires only a basic knowledge of computers and programming language translators (compilers). Sections that are unusually difficult or contain features of limited use are marked with an asterisk. The manual is divided into five chapters, corresponding to five language 'levels.' Level 1 is a teaching language designed to introduce programming concepts to nonprogrammers. Level 2 is a language roughly comparable in power to FORTRAN, but departs from it in specific features. Level 3 is comparable in power to ALGOL or PL/I, but with specific differences, and contains information on the new ALPHA mode for alpha-numeric manipulations, on writing formatted reports, and on internal editing. Level 4 contains the entity-attribute-set features of SIMSCRIPT, which have been up- dated and augmented to provide a more powerful list-processing capability. Level 5, the simulation-oriented part of SIMSCRIPT II, contains statements for time advance, event and activity processing, generation of statistical variates, and accumulation and analysis of simulation-generated data. Two new debugging routines, BEFORE and AFTER, enable the monitoring of six complex processes.

Chebyshev approximations for Dawson’s integral

Abstract: Rational Chebyshev approximations to Dawson's integral are presented in well-conditioned forms for 1xj < 2.5, 2.5 ? jxj ? 3.5, 3.5 < jxj ? 5.0 and 5.0 < lxi. Maximal relative errors range down to between 2 X 10-20 and 7 X 10-22.

Some results on sparse matrices

Abstract: A comparison in the context of sparse matrices is made between the Product Form of the Inverse PFI (a form of Gauss-Jordan elimination) and the Elimination Form of the Inverse EFI (a form of Gaussian elimination). The precise relation of the elements of these two forms of the inverse is given in terms of the nontrivial elements of the three matrices L, U, U-1 associated with the triangular factorization of the coefficient matrix A; i.e., A = L. U, where L is lower triangular and U is unit upper triangular. It is shown that the zero- nonzero structure of the PFI always has more nonzeros than the EFI. It is proved that Gaussian elimination is a minimal algorithm with respect to preserving sparseness if the diagonal elements of the matrix A are nonzero. However, Gaussian elimination is not nec- essarily minimal if A has some zero diagonal elements. The same statements hold for the PFI as well. A probabilistic study of fill-in and computing times for the PFI and EF sparse ma- trix algorithms is presented. This study suggests quantitatively how rapidly sparse matrices fill up for increasing densities, and emphasizes the necessity for reordering to minimize fill-in. I. Introduction. A sparse matrix is a matrix with very few nonzero elements. In many applications, a rough rule seems to be that there are O(N) nonzero entries; typically, say 2 to 10 nonzero entries per row. If the dimension N of the matrix is not large, then there is no compelling reason to treat sparse matrices differently from full matrices. It is when N becomes large and one attempts computations with the sparse matrix that it becomes necessary to take advantage of the zeros. The reason for this is obvious: there is a storage requirement of the order of N2 and arithmetic operations count of the order of N3 for many matrix algorithms using the full matrix. On the other hand, by storing only nonzero quantities and using logical operations to decide when an arithmetic operation is necessary, the storage requirement and arithmetic operations count can be reduced by a factor of N in many instances. Of course, this not only becomes a sizable savings of computer time, but also dictates whether or not some problems can be attempted. Computations with sparse matrices are not new. Iterative techniques for these matrices, especially those related to the solution of partial differential equations have been extensively developed (1). Sparse matrix methods for solving linear equa- tions by direct methods have been used for a long time in linear programming and there is a large body of literature, computational experience, programs, and artfulness, which has been built up in this area (2)-(15). In most linear programming codes, the product form of the inverse (PFI) is the method used to solve linear equations (16)-(19), although there are exceptions (20)-(21). Methods for scaling, pivoting for accuracy and sparseness, structuring data, and handling input-output have been extensively developed (22)-(72). However, there do not seem to exist rigorous results

Optimal approximation in Hilbert spaces with reproducing kernel functions

Abstract: Characterisations of optimal linear estimation rules are given in terms of the reproducing kernel function of a suitable Hilbert space The results are illustrated by means of three different, useful function spaces, showing, among other things, how Gaussian quadrature rules, and the Whittaker Cardinal Function, relate to optimal linear estimation rules in particular spaces

Roundoff error analysis of the fast Fourier transform

Abstract: This paper presents an analysis of roundoff errors occurring in the floating-point computation of the fast Fourier transform. Upper bounds are derived for the ratios of the root-mean-square (RMS) and maximum roundoff errors in the output data to the RMS value of the input data for both single and multidimensional transformations. These bounds are compared experimentally with actual roundoff errors.

Inequalities on the elements of the inverse of a certain tridiagonal matrix

Abstract: During an investigation into the convergence properties of natural splines it was found useful to have bounds on the inverse of a tridiagonal matrix with positive off-diagonal elements. Matrices of this type arise in other branches of numerical analysis, in particular in the discrete analogue of certain second-order differential operators, and so it may be useful to record these results. The matrix is

Convergence Properties of the Symmetric and Unsymmetric Successive Overrelaxation Methods and Related Methods

Abstract: The paper is concerned with variants of the successive overrelaxation method (SOR method) for solving the linear system Au = b. Necessary and sufficient conditions are given for the convergence of the symmetric and unsymmetric SOR methods when A is symmetric. The modified SOR, symmetric SOR, and unsymmetric SOR methods are also considered for systems of the form Diu, - CuU2 = bi, - CLU1 + D2u2 = b2 where DI and D2 are square diagonal matrices. Different values of the relaxation factor are used on each set of equations. It is shown that if the matrix corresponding to the Jacobi method of iteration has real eigenvalues and has spectral radius ji < 1, then the spectral radius of the matrix G associated with any of the methods is not less than that of the ordinary SOR method with w = 2(1 + (1 - TA2)112)-l. Moreover, if the eigenvalues of G are real then no improvement is possible by the use of semi-iterative methods. Introduction. In this paper we study convergence properties of several iterative methods for solving the linear system (1.1) Au- b, where A is a given real nonsingular N X N matrix with nonvanishing diagonal elements, b is a given column vector, and u is a column vector to be determined. Each method can be characterized by an equation

Computing irreducible representations of groups

Abstract: How can you find a complete set of inequivalent irreducible (ordinary) representations of a finite group? The theory is classical but, except when the group was very small or had a rather special structure, the actual computations were prohibitive before the advent of high-speed computers; and there remain practical difficulties even for groups of relatively small orders (< 100). The present paper describes three techniques to help solve this problem. These are: the reduction of a reducible unitary representation into its irreducible components; the construction of a complete set of irreducible unitary representations from a single faithful representation; and the calculation of the precise values of a group character from values which have only been computed approximately.

Error bounds for polynomial spline interpolation

Abstract: New upper and lower bounds for the L2 and L- norms of derivatives of the error in polynomial spline interpolation are derived. These results improve corresponding results of Ahlberg, Nilson, and Walsh, cf. (1), and Schultz and Varga, cf. (5). 1. Introduction. In this paper, we derive new bounds for the L2 and LX norms of derivatives of the error in polynomial spline interpolation. These bounds improve and generalize the known error bounds, cf. (1) and (5), in the following important ways: (1) these bounds can be explicitly calculated and are not merely asymptotic error bounds such as those given in (1) and (5); (2) explicit lower bounds are given for the error for a class of functions; (3) the degree of regularity required of the func- tion, f, being interpolated is extended, i.e., in L1) and (5) we demand that the mth or 2mth derivative of f be in L2, if we are interpolating by splines of degree 2m - 1, while here we demand only that some pth derivative of f, where m ? p ? 2m, be in L2; and (4) bounds are given for high-order derivatives of the interpolation errors. 2. Notations. Let - o < a < b < o and for each positive integer, m, let Km(a, b) denote the collection of all real-valued functions u(x) defined on (a, b) such that u E Cm'-(a, b) and such that Dm-lu is absolutely continuous, with Dmu E L2 (a, b), where Du _ du/dx denotes the derivative of u. For each nonnegative integer, M,

A note on trapezoidal methods for the solution of initial value problems

Abstract: The trapezoidal rule for the numerical integration of first-order ordinary differential equations is shown to possess, for a certain type of problem, an undesirable property. The removal of this difficulty is shown to be straightforward, resulting in a modified trapezoidal rule. Whilst this latent difficulty is slight (and probably rare in practice), the fact that the proposed modification involves negligible additional programming effort would suggest that it is worthwhile. A corresponding modification for the trapezoidal rule for the Goursat problem is also included. 1. The Trapezoidal Rule. We consider the numerical solution of the initial value problem (1.1) y f(x, y), y(xO) = yo in the region xo < x ? X. The trapezoidal rule for the numerical solution of (1.1) is given by (1.2) Yn+l Yn = h [f(Xn+l, Y.+,) + f(x, YA,)], where x,+1 = xn + h, h being the mesh length in the x direction. Equation (1.2) is a one-step implicit finite-difference method which is frequently employed for the numerical solution of (I.1). In fact, it is well known that it is the most accurate A-stable multistep method [1]. A method M is said to be A-stable if all solutions of M tend to zero as n --+ co when the method is applied with fixed positive h to any differential equation of the form (1.3) y= -Xy, where X is a complex constant with positive real part. For such an equation, (1.2) can be seen to reduce to