A backward stability analysis of diagonal pivoting methods for solving unsymmetric tridiagonal systems without interchanges

This paper concerns the LBM ^T factorization of unsymmetric tridiagonal matrices, where L and M are unit lower triangular matrices and B is block diagonal with 1×1 and 2×2 blocks. In some applications, it is necessary to form this factorization without row or column interchanges while the tridiagonal matrix is formed. Bunch and Kaufman proposed a pivoting strategy without interchanges specifically for symmetric tridiagonal matrices, and more recently, Bunch and Marcia proposed pivoting strategies that are normwise backward stable for linear systems involving such matrices. In this paper, we extend these strategies to the unsymmetric tridiagonal case and demonstrate that the proposed methods both exhibit bounded growth factors and are normwise backward stable. Copyright © 2009 John Wiley & Sons, Ltd.     

Pivoting strategy for rank-one modification ofLDM t-like factorization

The rank-one modification algorithm of theLDM ^t factorization was given by Bennett [1]. His method, however, could break down even when the matrix is nonsingular and well-conditioned. We introduce a pivoting strategy for avoiding possible break-down as well as for suppressing error growth in the modification process. The method is based on a symbolic formula of the rank-one modification of the factorization of a possibly singular nonsymmetric matrix. A new symbolic formula is also obtained for the inverses of the factor matrices. Repeated application of our method produces theLDM ^t-like product form factorization of a matrix. A numerical example is given to illustrate our pivoting method. An incomplete factorization algorithm is also introduced for updating positive definite matrix useful in quasi-Newton methods, in which the Fletcher and Powell algorithm [2] and the Gill, Murray and Saunders algorithm [4] are usually used.This paper is presented at the Japan SIAM Annual Meeting held at University of Tokyo, Japan, October 7–9, 1991.     

Partial pivoting strategies for symmetric gaussian elimination

In an earlier paper [6] it is shown how the use of symmetric additions of rows and columns enables a stableLDL ^T factorization of symmetric indefinite matrices. In this paper we describe partial pivoting strategies. These strategies are faster than the complete pivoting strategies that were introduced in the first paper. Numerical experiments are included. The results show that some of the new strategies share the stable behaviour of complete pivoting while running almost as fast as the well-known Cholesky decomposition.     

Detecting “dense” columns in interior point methods for linear programs

During the iterations of interior point methods symmetric indefinite systems are decomposed by LD̂L ^T factorization. This step can be performed in a special way where the symmetric indefinite system is transformed to a positive definite one, called the normal equations system. This approach proved to be efficient in most of the cases and numerically reliable, due to the positive definite property. It has been recognized, however, that in case the linear program contains “dense” columns, this approach results in an undesirable fill–in during the computations and the direct factorization of the symmetric indefinite system is more advantageous. The paper describes a new approach to detect cases where the system of normal equations is not preferable for interior point methods and presents a new algorithm for detecting the set of columns which is responsible for the excessive fill–in in the matrix AA ^T . By solving large–scale linear programming problems we demonstrate that our heuristic is reliable in practice. This work was supported in part by the Hungarian Scientific Research Fund OTKA K60480.     

On the eigenvalue distribution of a class of preconditioning methods

Summary A class of preconditioning methods depending on a relaxation parameter is presented for the solution of large linear systems of equationAx=b, whereA is a symmetric positive definite matrix. The methods are based on an incomplete factorization of the matrixA and include both pointwise and blockwise factorization. We study the dependence of the rate of convergence of the preconditioned conjugate gradient method on the distribution of eigenvalues ofC ^–1 A, whereC is the preconditioning matrix. We also show graphic representations of the eigenvalues and present numerical tests of the methods.     

Sparse quadratic programming in chemical process optimization

The quadratic programming aspects of a full space successive quadratic programming (SQP) method are described. In particular, fill-in, matrix factor and active set updating, numerical stability, and indefiniteness of the Hessian matrix are discussed in conjunction with a sparse modification of Bunch and Parlett factorization of symmetric indefinite (Kuhn-Tucker) matrices of the type often encountered in optimization. A new pivoting strategy, called constrained pivoting, is proposed to reduce fill-in and compared with complete, partial and threshold pivoting. It is shown that constrained pivoting often significantly reduces fill-in and thus the iterative computational burdens associated with the factorization and solution of Kuhn-Tucker conditions within the QP subproblem. A numerical algorithm for updating the lower triangular and diagonal factors is presented and shown to be very fast, usually requiring only about 5% of the cost of refactorization. Two active set strategies are also presented. These include the options of adding inequalities either one or several at a time. In either case, the effects on matrix factor updating is shown to be small. Finally, a simple test is used to maintain iterative descent directions in the quadratic program. Our sparse symmetric indefinite QP algorithm is tested in the context of a family of SQP algorithms that include a full space Newton method with analytical derivatives, a full space BFGS method and a Range and Null space Decomposition (RND) method in which the projected Hessian is calculated from either analytical second derivatives or the BFGS update. Several chemical process optimization problems, with small and large degrees of freedom, are used as test problems. These include minimum work calculations for multistage isothermal compression, minimum area targeting for heat exchanger networks, and distillation optimizations involving some azeotropic and extractive distillations. Numerical results show uniformly that both the proposed QP and SQP algorithms, particularly the full space Newton method, are reliable and efficient. No failures were experienced at either level.     

A twisted factorization method for symmetric SVD of a complex symmetric tridiagonal matrix

This paper presents an O(n²) method based on the twisted factorization for computing the Takagi vectors of an n‐by‐n complex symmetric tridiagonal matrix with known singular values. Since the singular values can be obtained in O(n²) flops, the total cost of symmetric singular value decomposition or the Takagi factorization is O(n²) flops. An analysis shows the accuracy and orthogonality of Takagi vectors. Also, techniques for a practical implementation of our method are proposed. Our preliminary numerical experiments have verified our analysis and demonstrated that the twisted factorization method is much more efficient than the implicit QR method, divide‐and‐conquer method and Matlab singular value decomposition subroutine with comparable accuracy. Copyright © 2009 John Wiley & Sons, Ltd.     

Computing Row and Column Counts for Sparse QR and LU Factorization

We present algorithms to determine the number of nonzeros in each row and column of the factors of a sparse matrix, for both the QR factorization and the LU factorization with partial pivoting. The algorithms use only the nonzero structure of the input matrix, and run in time nearly linear in the number of nonzeros in that matrix. They may be used to set up data structures or schedule parallel operations in advance of the numerical factorization.The row and column counts we compute are upper bounds on the actual counts. If the input matrix is strong Hall and there is no coincidental numerical cancellation, the counts are exact for QR factorization and are the tightest bounds possible for LU factorization.These algorithms are based on our earlier work on computing row and column counts for sparse Cholesky factorization, plus an efficient method to compute the column elimination tree of a sparse matrix without explicitly forming the product of the matrix and its transpose.This revised version was published online in October 2005 with corrections to the Cover Date.     

Using a hybrid preconditioner for solving large-scale linear systems arising from interior point methods

We devise a hybrid approach for solving linear systems arising from interior point methods applied to linear programming problems. These systems are solved by preconditioned conjugate gradient method that works in two phases. During phase I it uses a kind of incomplete Cholesky preconditioner such that fill-in can be controlled in terms of available memory. As the optimal solution of the problem is approached, the linear systems becomes highly ill-conditioned and the method changes to phase II. In this phase a preconditioner based on the LU factorization is found to work better near a solution of the LP problem. The numerical experiments reveal that the iterative hybrid approach works better than Cholesky factorization on some classes of large-scale problems.     

A factored implicit scheme for numerical weather prediction

The nonlinear partial differential equations of atmospheric dynamics govern motion on two time scales, a fast one and a slow one. Only the slow-scale motions are relevant in predicting the evolution of large weather patterns. Implicit numerical methods are therefore attractive for weather prediction, since they permit a large time step chosen to resolve only the slow motions. To develop an implicit method which is efficient for problems in more than one spatial dimension, one must approximate the problem by smaller, usually one-dimensional problems. A popular way to do so is to approximately factor the multidimensional implicit operator into one-dimensional operators. The factorization error incurred in such methods, however, is often unacceptably large for problems with multiple time scales. We propose a new factorization method for numerical weather prediction which is based on factoring separately the fast and slow parts of the implicit operator. We show analytically that the new method has small factorization error, which is comparable to other discretization errors of the overall scheme. The analysis is based on properties of the shallow water equations, a simple two-dimensional version of the fully three-dimensional equations of atmospheric dynamics.     

An Inversion-Free Estimating Equations Approach for Gaussian Process Models

One of the scalability bottlenecks for the large-scale usage of Gaussian processes is the computation of the maximum likelihood estimates of the parameters of the covariance matrix. The classical approach requires a Cholesky factorization of the dense covariance matrix for each optimization iteration. In this work, we present an estimating equations approach for the parameters of zero-mean Gaussian processes. The distinguishing feature of this approach is that no linear system needs to be solved with the covariance matrix. Our approach requires solving an optimization problem for which the main computational expense for the calculation of its objective and gradient is the evaluation of traces of products of the covariance matrix with itself and with its derivatives. For many problems, this is an O(nlog?n) effort, and it is always no larger than O(n²). We prove that when the covariance matrix has a bounded condition number, our approach has the same convergence rate as does maximum likelihood in that the Godambe information matrix of the resulting estimator is at least as large as a fixed fraction of the Fisher information matrix. We demonstrate the effectiveness of the proposed approach on two synthetic examples, one of which involves more than 1 million data points.     

Numerical solution of nonlinear partial quadratic integro‐differential equations of fractional order via hybrid of block‐pulse and parabolic functions

In this paper, an effective numerical approach based on a new two‐dimensional hybrid of parabolic and block‐pulse functions (2D‐PBPFs) is presented for solving nonlinear partial quadratic integro‐differential equations of fractional order. Our approach is based on 2D‐PBPFs operational matrix method together with the fractional integral operator, described in the Riemann–Liouville sense. The main characteristic behind this approach is to reduce such problems to those of solving systems of algebraic equations, which greatly simplifies the problem. By using Newton's iterative method, this system is solved, and the solution of fractional nonlinear partial quadratic integro‐differential equations is achieved. Convergence analysis and an error estimate associated with the proposed method is obtained, and it is proved that the numerical convergence order of the suggested numerical method is O(h³) . The validity and applicability of the method are demonstrated by solving three numerical examples. Numerical examples are presented in the form of tables and graphs to make comparisons with the exact solutions much easier.     

Comparative computer results of a new complementary pivot algorithm for solving equality and inequality constrained optimization problems

This paper reports the development of a new algorithm for solving the general constrained optimization problem (that of optimizing an objective function subject to both equality and inequality constraints). The approach is based on the complementary pivoting algorithms which have been developed to solve certain classes of fixed point problems. The specific approach is to use the equality constraints to solve for some variables in terms of the remaining ones thus enabling one to eliminate the equality constraints altogether. The result, under certain circumstances, is an optimization problem which may be transformed into a fixed point problem in such a way that a complementary pivoting code may be used to search for a solution.Seventeen test problems have been solved by this method and the results are compared against those obtained from GRG (Generalized Reduced Gradient method). The results of the tests indicate that the fixed point approach is robust (all 17 problems were solved by this method where as GRG solved 16). As to the computer times, the fixed point code proved to be as fast or faster than GRG on the lower dimensional problems; however, as the dimension increased, the trend reversed and on a 40 dimensional problem GRG was approximately 11 times faster. The conclusion from these tests is that when the dimension of the original problem can be reduced sufficiently by the equality constraints, the fixed point approach appears to be more effective than GRG.     

Time accurate fast wavelet‐Taylor Galerkin method for partial differential equations

We introduce the concept of fast wavelet‐Taylor Galerkin methods for the numerical solution of partial differential equations. In wavelet‐Taylor Galerkin method discretization in time is performed before the wavelet based spatial approximation by introducing accurate generalizations of the standard Euler, θ and leap‐frog time‐stepping scheme with the help of Taylor series expansions in the time step. We will present two different time‐accurate wavelet schemes to solve the PDEs. First, numerical schemes taking advantage of the wavelet bases capabilities to compress the operators and sparse representation of functions which are smooth, except for in localized regions, up to any given accuracy are presented. Here numerical experiments deal with advection equation with the spiky solution in one dimension, two dimensions, and nonlinear equation with a shock in solution in two dimensions. Second, our schemes deal with more regular class of problems where wavelets are not efficient procedure for data compression but we can use the good approximation properties of wavelet. Here time‐accurate schemes lead to consistent mass matrix in an explicit time stepping, which can be solved by approximate factorization techniques. Numerical experiment deals with more regular class of problems like heat equation as well as coupled linear system in two dimensions. © 2005 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2006     

Product-form Cholesky factorization in interior point methods for second-order cone programming

Second-order cone programming (SOCP) problems are typically solved by interior point methods. As in linear programming (LP), interior point methods can, in theory, solve SOCPs in polynomial time and can, in practice, exploit sparsity in the problem data. Specifically, when cones of large dimension are present, the density that results in the normal equations that are solved at each iteration can be remedied in a manner similar to the treatment of dense columns in an LP. Here we propose a product-form Cholesky factorization (PFCF) approach, and show that it is more numerically stable than the alternative Sherman-Morrison-Woodbury approach. We derive several PFCF variants and compare their theoretical perfomance. Finally, we prove that the elements of L in the Cholesky factorizations LDL^T that arise in interior point methods for SOCP are uniformly bounded as the duality gap tends to zero as long as the iterates remain is some conic neighborhood of the cental path.Mathematics Subject Classification (1991): 90C25, 90C51, 15A23Research supported in part by NSF Grants CDA 97-26385, DMS 01-04282, ONR Grant N000140310514 and DOE Grant GE-FG01-92ER-25126     

The importance of structure in incomplete factorization preconditioners

In this paper, we consider level-based preconditioning, which is one of the basic approaches to incomplete factorization preconditioning of iterative methods. It is well-known that while structure-based preconditioners can be very useful, excessive memory demands can limit their usefulness. Here we present an improved strategy that considers the individual entries of the system matrix and restricts small entries to contributing to fewer levels of fill than the largest entries. Using symmetric positive-definite problems arising from a wide range of practical applications, we show that the use of variable levels of fill can yield incomplete Cholesky factorization preconditioners that are more efficient than those resulting from the standard level-based approach. The concept of level-based preconditioning, which is based on the structural properties of the system matrix, is then transferred to the numerical incomplete decomposition. In particular, the structure of the incomplete factorization determined in the symbolic factorization phase is explicitly used in the numerical factorization phase. Further numerical results demonstrate that our level-based approach can lead to much sparser but efficient incomplete factorization preconditioners.     

ILUCP: a Crout ILU preconditioner with pivoting

We introduce a new preconditioner, ILUCP, to be used with an iterative method for solving sparse linear systems. It is based on an incomplete LU factorization combining Crout's formulation of Gaussian elimination with pivoting by columns. It is usually faster than ILUTP, which is based on a delayed update version of Gaussian elimination with pivoting, but requires more memory. For applications where memory is not a primary concern, ILUCP can be an attractive alternative to ILUTP. Copyright © 2005 John Wiley & Sons, Ltd.     

On large-scale diagonalization techniques for the Anderson model of localization

We propose efficient preconditioning algorithms for an eigenvalue problem arising in quantum physics, namely the computation of a few interior eigenvalues and their associated eigenvectors for large-scale sparse real and symmetric indefinite matrices of the Anderson model of localization. Our preconditioning approaches for the shift-and-invert symmetric indefinite linear system are based on maximum weighted matchings and algebraic multi-level incomplete LDL^T factorizations. These techniques can be seen as a complement to the alternative idea of using more complete pivoting techniques for the highly ill-conditioned symmetric indefinite Anderson matrices. Our numerical examples reveal that recent algebraic multi-level preconditioning solvers can accelerate the computation of a large-scale eigenvalue problem corresponding to the Anderson model of localization by several orders of magnitude. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Fourth-Order Symmetric DIRK Methods for Periodic Stiff Problems

New symmetric DIRK methods specially adapted to the numerical integration of first-order stiff ODE systems with periodic solutions are obtained. Our interest is focused on the dispersion (phase errors) of the dominant components in the numerical oscillations when these methods are applied to the homogeneous linear test model. Based on this homogeneous test model we derive the dispersion conditions for symmetric DIRK methods as well as symmetric stability functions with real poles and maximal dispersion order. Two new fourth-order symmetric methods with four and five stages are obtained. One of the methods is fourth-order dispersive whereas the other method is symplectic and sixth-order dispersive. These methods have been applied to a number of test problems (linear as well as nonlinear) and some numerical results are presented to show their efficiency when they are compared with the symplectic DIRK method derived by Sanz-Serna and Abia (SIAM J. Numer. Anal. 28 (1991) 1081–1096).     

Methods and algorithms of solving spectral problems for polynomial and rational matrices

In the present paper, methods and algorithms for numerical solution of spectral problems and some problems in algebra related to them for one- and two-parameter polynomial and rational matrices are considered. A survey of known methods of solving spectral problems for polynomial matrices that are based on the rank factorization of constant matrices, i.e., that apply the singular value decomposition (SVD) and the normalized decomposition (the QR factorization), is given. The approach to the construction of methods that makes use of rank factorization is extended to one- and two-parameter polynomial and rational matrices. Methods and algorithms for solving some parametric problems in algebra based on ideas of rank factorization are presented. Bibliography: 326titles.Dedicated to the memory of my son AlexanderTranslated fromZapiski Nauchnykh Seminarov POMI, Vol. 238, 1997, pp. 7–328.Translated by V. N. Kublanovskaya.