Showing papers in "Siam Journal on Scientific and Statistical Computing in 1984"
TL;DR: In this article, the use of Partial Least Squares (PLS) for handling collinearities among the independent variables X in multiple regression is discussed, and successive estimates are obtained using the residuals from previous rank as a new dependent variable y.
Abstract: The use of partial least squares (PLS) for handling collinearities among the independent variables X in multiple regression is discussed. Consecutive estimates $({\text{rank }}1,2,\cdots )$ are obtained using the residuals from previous rank as a new dependent variable y. The PLS method is equivalent to the conjugate gradient method used in Numerical Analysis for related problems.To estimate the “optimal” rank, cross validation is used. Jackknife estimates of the standard errors are thereby obtained with no extra computation.The PLS method is compared with ridge regression and principal components regression on a chemical example of modelling the relation between the measured biological activity and variables describing the chemical structure of a set of substituted phenethylamines.
2,290 citations
TL;DR: In this article, the upwind-differencing first-order schemes of Godunov, Engquist-Osher and Roe are discussed on the basis of the inviscid Burgers equations.
Abstract: The upwind-differencing first-order schemes of Godunov, Engquist–Osher and Roe are discussed on the basis of the inviscid Burgers equations. The differences between the schemes are interpreted as differences between the approximate Riemann solutions on which their numerical flux-functions are based. Special attention is given to the proper formulation of these schemes when a source term is present. Second-order two-step schemes, based on the numerical flux-functions of the first-order schemes are also described. The schemes are compared in a numerical experiment and recommendations on their use are included.
444 citations
TL;DR: In this paper, a method for producing monotone piecewise cubic interpolants to monotonous data is described, which is completely local and which is extremely simple to implement.
Abstract: A method is described for producing monotone piecewise cubic interpolants to monotone data which is completely local and which is extremely simple to implement.
294 citations
TL;DR: It is shown that general sparse sets of linear equations whose pattern is symmetric (or nearly so) can be solved efficiently by a multifrontal technique and there is scope for extra performance during factorization and solution on a vector or parallel machine.
Abstract: We show that general sparse sets of linear equations whose pattern is symmetric (or nearly so) can be solved efficiently by a multifrontal technique. The main advantages are that the analysis time is small compared to the factorization time and that analysis can be performed in a predictable amount of storage. Additionally, there is scope for extra performance during factorization and solution on a vector or parallel machine. We show performance figures for examples run on the IBM 3081K and CRAY-1 computers.
209 citations
TL;DR: In this paper, the authors give a necessary and sufficient condition for equality in (1) and discuss nonuniqueness of the representation of the nonlinear functions of linear combinations in projection pursuit algorithms.
Abstract: Projection pursuit algorithms approximate a function of p variables by a sum of nonlinear functions of linear combinations: \[ (1)\qquad f\left( {x_1 , \cdots ,x_p } \right) \doteq \sum_{i = 1}^n {g_i \left( {a_{i1} x_1 + \cdots + a_{ip} x_p } \right)} . \] We develop some approximation theory, give a necessary and sufficient condition for equality in (1), and discuss nonuniqueness of the representation.
159 citations
TL;DR: In this article, a numerical method is given for the solution of a system of ordinary differential equations and algebraic, unilateral constraints, where contacts between the bodies are created and disappear in the time interval of interest.
Abstract: A numerical method is given for the solution of a system of ordinary differential equations and algebraic, unilateral constraints. The equations govern the motion of a mechanical system of rigid bodies, where contacts between the bodies are created and disappear in the time interval of interest. The ordinary differential equations are discretized by linear multistep methods. In order to satisfy the constraints, a quadratic programming problem is solved at each time step. The fact that the variation of the objective function is small from step to step is utilized to save computing time. A discrete friction model, based on Coulomb’s law of friction and suitable for efficient computation, is proposed for planar problems where dry friction cannot be neglected. The normal forces and the friction forces are the optimal solution to a quadratic programming problem. The methods are tested on four model problems. A data structure and possible generalizations are discussed.
156 citations
TL;DR: A new method is described, suitable for any decreasing or symmetric unimodal density, that is faster and more easily implemented, thereby providing a standard procedure for developing both the fast and the slow part for many given densities.
Abstract: The fastest computer methods for sampling from a given density are those based on a mixture of a fast and slow part. This paper describes a new method of this type, suitable for any decreasing or symmetric unimodal density. It differs from others in that it is faster and more easily implemented, thereby providing a standard procedure for developing both the fast and the slow part for many given densities. It is called the ziggurat method, after the shape of a single, convenient density that provides for both the fast and the slow parts of the generating process. Examples are given for REXP and RNOR, subroutines that generate exponential and normal variates that, as assembler routines, are nearly twice as fast as the previous best assembler routines, and that, as Fortran subroutines, approach the limiting possible speed: the time for Fortran subroutine linkage conventions plus the time to generate one uniform variate.
132 citations
TL;DR: A nonlinear SSOR type preconditioning is derived which numerical experiments show to be as effective as the linear SSOR preconditionsing that uses the Jacobian explicitly.
Abstract: We propose an algorithm for implementing Newton's method for a general nonlinear system $f(x) = 0$ where the linear systems that arise at each step of Newton's method are solved by a preconditioned Krylov subspace iterative method. The algorithm requires only function evaluations and does not require the evaluation or storage of the Jacobian matrix. Matrix-vector products involving the Jacobian matrix are approximated by directional differences. We develop a framework for constructing preconditionings for this inner iterative method which do not reference the Jacobian matrix explicitly. We derive a nonlinear SSOR type preconditioning which numerical experiments show to be as effective as the linear SSOR preconditioning that uses the Jacobian explicitly.
123 citations
TL;DR: The authors describes a variety of methods for generating random correlation matrices, with emphasis on choice of random variables and distributions so as to provide matrices with given structure, expected values or eigenvalues.
Abstract: This paper describes a variety of methods for generating random correlation matrices, with emphasis on choice of random variables and distributions so as to provide matrices with given structure, expected values or eigenvalues.
116 citations
TL;DR: A multiprocessor structure for solving a dense system of n linear equations that is well suited for VLSI implementation as identical processors with a simple and regular interconnection pattern are required.
Abstract: We propose a multiprocessor structure for solving a dense system of n linear equations. The solution is obtained in two stages. First, the matrix of coefficients is reduced to upper triangular form via Givens rotations. Second, a back substitution process is applied to the triangular system. A two-dimensional array of $\theta (n^2 )$ processors is employed to implement the first step, and (using a previously known scheme) a one-dimensional array of $\theta (n)$ processors is employed to implement the second step. These processor arrays allow both stages to be carried out in time $O(n)$, and they are well suited for VLSI implementation as identical processors with a simple and regular interconnection pattern are required.
99 citations
TL;DR: It is shown how an algorithm similar to the SYMMLQ can be derived for nonsymmetric problems and a more economical algorithm based upon the $LU$ factorization with partial pivoting is described.
Abstract: The main purpose of this paper is to develop stable versions of some Krylov subspace methods for solving linear systems of equations $Ax = b$ As in the case of Paige and Saunders's SYMMLQ [SIAM J Numer Anal, 12 (1975), pp 617–624], our algorithms are based on stable factorizations of the banded Hessenberg matrix representing the restriction of the linear application A to a Krylov subspace We will show how an algorithm similar to the SYMMLQ can be derived for nonsymmetric problems and we will describe a more economical algorithm based upon the $LU$ factorization with partial pivoting In the particular case where A is symmetric indefinite the new algorithm is theoretically equivalent to SYMMLQ but slightly more economical As a consequence, an advantage of the new approach is that nonsymmetric or symmetric indefinite or both nonsymmetric and indefinite systems of linear equations can be handled by a single algorithm
TL;DR: In this paper, the existence of multiple steady states with the same farfield behavior is discussed for simple 1-D transonic model problems, and numerical experiments show that implicit schemes can converge to the physically unstable steady states and this phenomenon is analysed.
Abstract: The existence of multiple steady states with the same farfield behavior is discussed for simple 1-D transonic model problems. These multiple solutions all have only entropy satisfying compressive steady shock waves. Only some of these solutions are stable in the time-dependent system and are accessible through physical time-dependent perturbations. This is demonstrated by some elementary explicit solutions in a scalar model problem. However, for a large class of initial data and large C.F.L. numbers, numerical experiments show that implicit schemes can converge to the physically unstable steady states and this phenomenon is analysed. The scalar model is also discussed as a very simple numerical test problem for implicit schemes with rich structure in both the steady state and time-dependent regimes.
TL;DR: In this article, a new finite difference scheme for the Stokes equations and incompressible Navier-Stokes equations for low Reynolds number is presented. The scheme uses the primitive variable formulation of the equations and is applicable with nonuniform grids and nonrectangular geometries.
Abstract: This paper presents a new finite difference scheme for the Stokes equations and incompressible Navier–Stokes equations for low Reynolds number. The scheme uses the primitive variable formulation of the equations and is applicable with nonuniform grids and nonrectangular geometries. Several other methods of solving the Navier–Stokes equations are also examined in this paper and some of their strengths and weaknesses are described. Computational results using the new scheme are presented for the Stokes equations for a region with curved boundaries and for a disk with polar coordinates. The results show the method to be second-order accurate.
TL;DR: Several incomplete factorization methods (ILU) for strongly nonsymmetric blockbanded systems are developed in this paper, which result from the coupled partial differential equations which occur in the simulation of enhanced oil recovery.
Abstract: Several incomplete factorization methods (ILU) for strongly nonsymmetric block-banded systems are developed. These systems result from the coupled partial differential equations which occur in the simulation of enhanced oil recovery. Several approximations using natural, diagonal (D2) and alternating diagonal (D4) orderings are used with ORTHOMIN acceleration. These methods can all be used in conjunction with the COMBINATIVE method for multi-phase problems. Use of the modified first order factorization is also investigated. Test results are given for single- and multi-phase problems. Timings for both scalar and vector mode on the CRAY are presented.
TL;DR: In this paper, the Boussinesq equation is studied numerically using Galerkin methods and soliton solutions are shown to exist, and an analytic formula for the two-soliton interaction is presented and verified by numerical experiment.
Abstract: The “good” Boussinesq equation is studied numerically using Galerkin methods and soliton solutions are shown to exist. An analytic formula for the two-soliton interaction is presented and verified by numerical experiment.
TL;DR: It is hoped that future improvements to linear-equation software will be oriented more specifically to the case of related matrices, as well as to avoid refactorization, or to speed the convergence of an iterative method.
Abstract: Optimization algorithms typically require the solution of many systems of linear equations $B_k y_k = b_k $. When large numbers of variables or constraints are present, these linear systems could account for much of the total computation time.Both direct and iterative equation solvers are needed in practice. Unfortunately, most of the off-the-shelf solvers are designed for single systems, whereas optimization problems give rise to hundreds or thousands of systems. To avoid refactorization, or to speed the convergence of an iterative method, it is essential to note that $B_k $ is related to $B_{k - 1} $.We review various sparse matrices that arise in optimization, and discuss compromises that are currently being made in dealing with them. Since significant advances continue to be made with single-system solvers, we give special attention to methods that allow such solvers to be used repeatedly on a sequence of modified systems (e.g., the product-form update; use of the Schur complement). The speed of factorizing a matrix then becomes relatively less important than the efficiency of subsequent solves with very many right-hand sides.At the same time, we hope that future improvements to linear-equation software will be oriented more specifically to the case of related matrices $B_k $.
TL;DR: This paper shows that the block-elimination algorithm may become unstable and inaccurate when A is nearly singular, and proposes a stable variant which employs deflation techniques for solving the two systems with A.
Abstract: In numerical continuation methods for solving parametrized nonlinear systems, one often has to solve linear systems with matrices of the following form: \[ M = \left[ {\begin{array}{*{20}c} A & b \\ {c^T } & d \\ \end{array} } \right] \] where A may become singular but M is well conditioned. If A has special structures (e.g. sparseness, special data structure, special solver), then direct Gaussian elimination on M with pivoting will destroy the structures in A. An often used method that does exploit structures in A is the block-elimination (BE) algorithm which involves solving two systems with A for each system with M. In this paper, we show that the BE algorithm may become unstable and inaccurate when A is nearly singular. We then propose a stable variant which employs deflation techniques for solving the two systems with A. The deflation techniques can be viewed as working in coordinate systems orthogonal to the approximate null vectors of A, enabling an accurate representation of the solution to be computed. The extra work amounts to a few (e.g. 2) more backsolves with A. Backward error bounds and numerical results are presented.
TL;DR: A survey of numerical methods for large sparse linear least squares problems, focusing mainly on developments since the last comprehensive surveys of the subject published in 1976, considers direct methods based on elimination and on orthogonalization, as well as various iterative methods.
Abstract: Large sparse least squares problems arise in many applications, including geodetic network adjustments and finite element structural analysis. Although geodesists and engineers have been solving such problems for years, it is only relatively recently that numerical analysts have turned attention to them. In this paper we present a survey of numerical methods for large sparse linear least squares problems, focusing mainly on developments since the last comprehensive surveys of the subject published in 1976. We consider direct methods based on elimination and on orthogonalization, as well as various iterative methods. The ramifications of rank deficiency, constraints, and updating are also discussed.
TL;DR: It is concluded that gradual underflow makes it significantly easier to write good numerical codes than store zero, and that this remains true even if extra range and precision are available for intermediate calculations.
Abstract: We examine the effects of different underflow mechanisms on the reliability of numerical software. Software is considered reliable in the face of underflow if the effects of underflow are no worse than the uncertainty due to roundoff alone. The two primary underflow mechanisms discussed are store zero and gradual underflow, although we consider other mechanisms as well. By examining a variety of codes, including Gaussian elimination, polynomial evaluation, and eigenvalue calculation, we conclude that gradual underflow makes it significantly easier to write good numerical codes than store zero, and that this remains true even if extra range and precision are available for intermediate calculations.
TL;DR: In this article, a new open boundary condition applicable to rotating stratified fluid flows is proposed, which allows the fluid to be forced everywhere, including the open boundaries, and some numerical experiments with a fluid contained within an infinitely long channel are considered.
Abstract: Previous studies on open boundary conditions in unbounded rotating fluid flows have concentrated on how to accurately simulate the outflow of free waves through an open boundary. In most limited area integrations of rotating fluids, however, the generated waves are forced rather than free, in which case the Sommerfeld radiation condition, applied in previous studies, is not valid. A new open boundary condition applicable to rotating stratified fluid flows is suggested below, which allows the fluid to be forced everywhere, including the open boundaries. To prove the applicability of the suggested open boundary conditions, some numerical experiments with a fluid contained within an infinitely long channel are considered. The new open boundary conditions are applied at two cross sections along the channel.
TL;DR: This paper gives an implementation of the force method which is numerically stable and preserves sparsity, and in its approach each of the two main phases is carried out using orthogonal factorizaton techniques recently developed for linear least squares problems.
Abstract: Historically there are two principal methods of matrix structural analysis, the displacement (or stiffness) method and the force (or flexibility) method. In recent times the force method has been used relatively little because the displacement method has been deemed easier to implement on digital computers, especially for large sparse systems. The force method has theoretical advantages, however, for multiple redesign problems or nonlinear elastic analysis because it allows the solution of modified problems without restarting the computation from the beginning. In this paper we give an implementation of the force method which is numerically stable and preserves sparsity. Although it is motivated by earlier elimination schemes, in our approach each of the two main phases of the force method is carried out using orthogonal factorizaton techniques recently developed for linear least squares problems.
TL;DR: In this article, a method for solving a form of the inverse Sturm-Liouville problem is presented, which is based on modifying the given differential eigenvalues so that the tridiagonal matrix eigenvalue problem recovered from the first N eigen values can be identified with the finite difference approximation of the required differential Eigen value problem.
Abstract: In this paper we present a method for solving a form of the inverse Sturm–Liouville problem. The basis of the method is to modify the given differential eigenvalues so that the $N \times N$ tridiagonal matrix eigenvalue problem recovered from the first N eigenvalues can (after a suitable transformation) be identified with the finite difference approximation of the required differential eigenvalue problem. Numerical results are presented to illustrate the effectiveness of this method.
TL;DR: A spline approximation method for computation of the gain operators in feedback controls is proposed and tested numerically and compares with a method based on “averaging” approximations.
Abstract: We consider the infinite interval regulator problem for systems with delays. A spline approximation method for computation of the gain operators in feedback controls is proposed and tested numerically. Comparison with a method based on “averaging” approximations is made.
TL;DR: In this article, incomplete orthogonal decompositions of A obtained by modifying either the Gram-Schmidt or the Givens rotation method were applied to four small structural analyses and the numerical efficiencies are compared.
Abstract: Incomplete methods are considered for the solution of $A^T Ax = b$. One possibility is to construct the product $A^T A$ and proceed with an ICCG solution. However, it is also possible to utilize incomplete orthogonal decompositions of A obtained by modifying either the Gram–Schmidt or the Givens rotation method. These three methods are applied to four small structural analyses and the numerical efficiencies are compared.
TL;DR: The design features of a code which solves sparse unsymmetric systems of linear equations using a frontal method, in particular the user interface, the internal data structures, the pivoting strategy, and the isolation of machine dependencies are considered.
Abstract: We discuss design features of a code which solves sparse unsymmetric systems of linear equations using a frontal method. We consider in particular the user interface, the internal data structures, the pivoting strategy, and the isolation of machine dependencies. We illustrate the performance of our code on a variety of test problems on both the IBM 3033 and the CRAY-1.
TL;DR: A comparison of an implementation of a simple direct $LU$ factorization method, suggested by Funderlic and Mankin, with other direct methods recommended by Paige, Styan and Wachter, for computing stationary distributions of Markov chains is reported.
Abstract: The purpose of this paper is to report on a comparison of an implementation of a simple direct $LU$ factorization method, suggested by Funderlic and Mankin [SIAM J. Sci. Stat. Comput., 2 (1981), pp. 375–383], with other direct methods recommended by Paige, Styan and Wachter [J. Stat. Comput.Simul., 4 (1975), pp. 173–186], for computing stationary distributions of Markov chains. A backward error analysis is developed and conditioning problems are addressed. The method is stable without pivoting, it is numerically efficient and it lends itself to symmetric pivoting to preserve sparsity.
TL;DR: A new method is described to compute the solutions of linear BVP in an efficient and stable way by decoupling the multiple shooting recursion; this means that the choice of output points can be made virtually without regard to restrictions.
Abstract: A new method is described to compute the solutions of linear BVP in an efficient and stable way. The stability is achieved by decoupling the multiple shooting recursion; this means that the choice of output points can be made virtually without regard to restrictions. By fixing the number of integration steps per “shooting” interval and assembling as many of them as is needed to fit the user's requirements, high efficiency is gained. Apart from a mathematical description, we also give a stability analysis of the method. A large number of numerical examples confirm this analysis and illustrate the possibilities of the algorithm.
TL;DR: This paper develops an updating algorithm for the solution of linear least squares problems which are sparse except for a small subset of dense equations, which can be divided into sparse and dense subsets.
Abstract: Linear least squares problems which are sparse except for a small subset of dense equations can be efficiently solved by an updating method. Often the least squares solution is also required to satisfy a set of linear constraints, which again can be divided into sparse and dense subsets. This paper develops an updating algorithm for the solution of such problems. The algorithm is completely general in that no restrictive assumption on the rank of any subset of equations or constraints is made.
TL;DR: In this paper, the inverse problem of electrocardiography is treated by a formulation of a suitable model and a powerful transformation of the time domain; regularization techniques using the singular value decomposition.
Abstract: The inverse problem of electrocardiography is treated by a) the formulation of a suitable model and a powerful transformation of the time domain; b) regularization techniques using the singular value decomposition.The solution calculated from the potentials at the body surface is expressed in terms of the activation times on the heart surface. Promising results of numerical experiments with simulated measurements errors are presented.
TL;DR: A new method for computing simple turning points of nonlinear equations of the form G(u,\lambda ) = 0 which is based on applying Newton's method to the characterization of digma, where $\sigma $ is a pseudo-arclength parameter used in a continuation method for following the solution paths.
Abstract: We present a new method for computing simple turning points of nonlinear equations of the form $G(u,\lambda ) = 0$ which is based on applying Newton's method to the characterization ${{d\lambda (\sigma )} / {d\sigma }} = 0$, where $\sigma $ is a pseudo-arclength parameter used in a continuation method for following the solution paths. The method is quadratically convergent and needs only one starting point on the solution path. Second derivatives of G (or difference approximations of them) have to be computed but the method is relatively insensitive to their values and they also give rise to a more accurate second order predictor in the continuation method. We present a chord-Newton variant for improving the efficiency of the algorithm which requires only one factorization of a Jacobian matrix. We also present a damped-Newton variant for improving the robustness and the global convergence of the algorithm. Results of numerical experiments on two standard nonlinear elliptic problems of Simpson's [SIAM J. Numer. Anal., 12 (1975), pp. 439–451] show that the new algorithm compares favorably with the best of the existing methods in terms of efficiency and robustness.