On the sensitivity of the LU factorization

This paper gives sensitivity analyses by two approaches forL andU in the factorizationA=LU for general perturbations inA which are sufficiently small in norm. By the matrix-vector equation approach, we derive the condition numbers for theL andU factors. By the matrix equation approach we derive corresponding condition estimates. We show how partial pivoting and complete pivoting affect the sensitivity of the LU factorization. The material presented here is a part of the first author's PhD thesis under the supervision of the second author. This research was supported by NSERC of Canada Grant OGP0009236.     

Rounding error analysis of the classical Gram-Schmidt orthogonalization process

This paper provides two results on the numerical behavior of the classical Gram-Schmidt algorithm. The first result states that, provided the normal equations associated with the initial vectors are numerically nonsingular, the loss of orthogonality of the vectors computed by the classical Gram-Schmidt algorithm depends quadratically on the condition number of the initial vectors. The second result states that, provided the initial set of vectors has numerical full rank, two iterations of the classical Gram-Schmidt algorithm are enough for ensuring the orthogonality of the computed vectors to be close to the unit roundoff level.The work of the second author was supported in part by the US Department of Energy, Office of Basic Energy Science under LAB03-17 initiative, DOE contract No. DE-FG02-03ER25584, and in part by the TeraScale Optimal PDE Simulations (TOPS) SciDAC, DoE Contract No. DE-FC02-01ER25480The work of the third author was supported by the project 1ET400300415 within the National Program of Research ‘‘Information Society’’ and by the GA AS CR under grant No. IAA1030405.     

Multiple LU factorizations of a singular matrix

A singular matrix A may have more than one LU factorizations. In this work the set of all LU factorizations of A is explicitly described when the lower triangular matrix L is nonsingular. To this purpose, a canonical form of A under left multiplication by unit lower triangular matrices is introduced. This canonical form allows us to characterize the matrices that have an LU factorization and to parametrize all possible LU factorizations. Formulae in terms of quotient of minors of A are presented for the entries of this canonical form.     

Componentwise error analysis for the block LU factorization of totally nonnegative matrices

For the first time, perturbation bounds including componentwise perturbation bounds for the block LU factorization have been provided by Dopico and Molera (2005) [5]. In this paper, componentwise error analysis is presented for computing the block LU factorization of nonsingular totally nonnegative matrices. We present a componentwise bound on the equivalent perturbation for the computed block LU factorization. Consequently, combining with the componentwise perturbation results we derive componentwise forward error bounds for the computed block factors.     

Saving flops in LU based shift-and-invert strategy

The shift-and-invert method is very efficient in eigenvalue computations, in particular when interior eigenvalues are sought. This method involves solving linear systems of the form (A−σI)z=b. The shift σ is variable, hence when a direct method is used to solve the linear system, the LU factorization of (A−σI) needs to be computed for every shift change. We present two strategies that reduce the number of floating point operations performed in the LU factorization when the shift changes. Both methods perform first a preprocessing step that aims at eliminating parts of the matrix that are not affected by the diagonal change. This leads to about 43% and 50% flops savings respectively for the dense matrices.     

Numerical enclosure for each eigenvalue in generalized eigenvalue problem

An algorithm for enclosing all eigenvalues in generalized eigenvalue problem Ax=λBx is proposed. This algorithm is applicable even if A∈C^n×n is not Hermitian and/or B∈C^n×n is not Hermitian positive definite, and supplies nerror bounds while the algorithm previously developed by the author supplies a single error bound. It is proved that the error bounds obtained by the proposed algorithm are equal or smaller than that by the previous algorithm. Computational cost for the proposed algorithm is similar to that for the previous algorithm. Numerical results show the property of the proposed algorithm.     

Factoring matrices with a tree-structured sparsity pattern

Let A be a matrix whose sparsity pattern is a tree with maximal degree d_max. We show that if the columns of A are ordered using minimum degree on |A|+|A|^∗, then factoring A using a sparse LU with partial pivoting algorithm generates only O(d_maxn) fill, requires only O(d_maxn) operations, and is much more stable than LU with partial pivoting on a general matrix. We also propose an even more efficient and just-as-stable algorithm called sibling-dominant pivoting. This algorithm is a strict partial pivoting algorithm that modifies the column preordering locally to minimize fill and work. It leads to only O(n) work and fill. More conventional column pre-ordering methods that are based (usually implicitly) on the sparsity pattern of |A|^∗|A| are not as efficient as the approaches that we propose in this paper.     

Numerical stability of fast cosine transforms

In this paper we consider the numerical stability of fast algorithms for discrete cosine transform (DCT) of type III and II, respectively. We show that various fast DCTs can possess a very different behaviour of numerical stability. By matrix factorizations we find that a complex fast DCT which is based mainly on a fast Fouier transform has a better numerical stability than a real fast DCT despite its larger arithmetical complexity. Numerical tests illustrate our theoretical results.     

Polynomial least squares fitting in the Bernstein basis

The problem of polynomial least squares fitting in which the usual monomial basis is replaced by the Bernstein basis is considered. The coefficient matrix of the overdetermined system to be solved in the least squares sense is then a rectangular Bernstein-Vandermonde matrix. In order to use the method based on the QR decomposition of A, the first stage consists of computing the bidiagonal decomposition of the coefficient matrix A. Starting from that bidiagonal decomposition, an algorithm for obtaining the QR decomposition of A is the applied. Finally, a triangular system is solved by using the bidiagonal decomposition of the R-factor of A. Some numerical experiments showing the behavior of this approach are included.     

Relative perturbation bounds for the unitary polar factor

LetB be anm×n (m≥n) complex (or real) matrix. It is known that there is a uniquepolar decomposition B=QH, whereQ*Q=I, then×n identity matrix, andH is positive definite, providedB has full column rank. Existing perturbation bounds suggest that in the worst case, for complex matrices the change inQ be proportional to the reciprocal ofB's least singular value, or the reciprocal of the sum ofB's least and second least singular values if matrices are real. However, there are situations where this unitary polar factor is much more accurately determined by the data than the existing perturbation bounds would indicate. In this paper the following question is addressed: how much mayQ change ifB is perturbed to $\tilde B = D_1^* BD_2 $ , whereD ₁ andD ₂ are nonsingular and close to the identity matrices of suitable dimensions? It is shown that for a such kind of perturbation, the change inQ is bounded only by the distances fromD ₁ andD ₂ to identity matrices and thus is independent ofB's singular values. Such perturbation is restrictive, but not unrealistic. We show how a frequently used scaling technique yields such a perturbation and thus scaling may result in better-conditioned polar decompositions.     

Adaptive data distribution for concurrent continuation

Summary Continuation methods compute paths of solutions of nonlinear equations that depend on a parameter. This paper examines some aspects of the multicomputer implementation of such methods. The computations are done on a mesh connected multicomputer with 64 nodes.One of the main issues in the development of concurrent programs is load balancing, achieved here by using appropriate data distributions. In the continuation process, many linear systems have to be solved. For nearby points along the solution path, the corresponding system matrices are closely related to each other. Therefore, pivots which are good for theLU-decomposition of one matrix are likely to be acceptable for a whole segment of the solution path. This suggests to choose certain data distributions that achieve good load balancing. In addition, if these distributions are used, the resulting code is easily vectorized.To test this technique, the invariant manifold of a system of two identical nonlinear oscillators is computed as a function of the coupling between them. This invariant manifold is determined by the solution of a system of nonlinear partial differential equations that depends on the coupling parameter. A symmetry in the problem reduces this system to one single equation, which is discretized by finite differences. The solution of the discrete nonlinear system is followed as the coupling parameter is changed.This material is based upon work supported by the NSF under Cooperative Agreement No. CCR-8809615. The government has certain rights in this material.     

Numerical experiments with two approximate inverse preconditioners

We present the results of numerical experiments aimed at comparing two recently proposed sparse approximate inverse preconditioners from the point of view of robustness, cost, and effectiveness. Results for a standard ILU preconditioner are also included. The numerical experiments were carried out on a Cray C98 vector processor. This work was partially supported by the GA AS CR under grant 2030706 and by the grant GA CR 205/96/0921.     

Relative perturbation results for eigenvalues and eigenvectors of diagonalisable matrices

Let and be a perturbed eigenpair of a diagonalisable matrixA. The problem is to bound the error in and . We present one absolute perturbation bound and two relative perturbation bounds. The absolute perturbation bound is an extension of Davis and Kahan's sin θ Theorem from Hermitian to diagonalisable matrices. The two relative perturbation bounds assume that and are an exact eigenpair of a perturbed matrixD ₁ AD ₂ , whereD ₁ andD ₂ are non-singular, butD ₁ AD ₂ is not necessarily diagonalisable. We derive a bound on the relative error in and a sin θ theorem based on a relative eigenvalue separation. The perturbation bounds contain both the deviation ofD ₁ andD ₂ from similarity and the deviation ofD ₂ from identity. This work was partially supported by NSF grant CCR-9400921.     

A modification to the linpack downdating algorithm

Alinpack downdating algorithm is being modified by interleaving its two different phases, the forward solving a triangular system and the backward sweep of Givens rotations, to yield a new forward method for finding the Cholesky decomposition ofR ^T R –zz ^T . By showing that the new algorithm saves forty percent purely redundant operations of the original, better stability properties are expected. In addition, various other downdating algorithms are rederived and analyzed under a uniform framework.     

The {\cal {D Q R}} algorithm, basic theory, convergence, and conditional stability

Summary. We describe a fast matrix eigenvalue algorithm that uses a matrix factorization and reverse order multiply technique involving three factors and that is based on the symmetric matrix factorization as well as on –orthogonal reduction techniques where is computed from the given matrix . It operates on a similarity reduction of a real matrix to general tridiagonal form and computes all of 's eigenvalues in operations, where the part of the operations is possibly performed over , instead of the 7–8 real flops required by the eigenvalue algorithm. Potential breakdo wn of the algorithm can occur in the reduction to tridiagonal form and in the –orthogonal reductions. Both, however, can be monitored during the computations. The former occurs rather rarely for dimensions and can essentially be bypassed, while the latter is extremely rare and can be bypassed as well in our conditionally stable implementation of the steps. We prove an implicit theorem which allows implicit shifts, give a convergence proof for the algorithm and show that is conditionally stable for general balanced tridiagonal matrices . Received April 25, 1995 / Revised version received February 9, 1996     

A posteriori error estimates for linear equations

Summary This paper describes upper and lowerp-norm error bounds for approximate solutions of the linear system of equationsAx=b. These bounds imply that the error is proportional to the quantity wherer is the residual andq is the conjugate index top. The constant of proportionality is larger than 1 and lies in a specified range. Similar results are obtained for approximations toA ^–1 and solutions of nonsingular linear equations on general spaces.Research was partially supported by NSF Grant DMS8901477     

Circular trichotomy of the spectrum of regular matrix pencils

We introduce a method for approximating the right and left deflating subspaces of a regular matrix pencil corresponding to the eigenvalues inside, on and outside the unit circle. The method extends the iteration used in the context of spectral dichotomy, where the assumption on the absence of eigenvalues on the unit circle is removed. It constructs two matrix sequences whose null spaces and the null space of their sum lead to approximations of the deflating subspaces corresponding to the eigenvalues of modulus less than or equal to 1, equal to 1 and larger than or equal to 1. An orthogonalization process is then used to extract the desired delating subspaces. The resulting algorithm is an inverse free, easy to implement, and sufficiently fast. The derived convergence estimates reveal the key parameters, which determine the rate of convergence. The method is tested on several numerical examples.     

Stable solutions of linear systems involving long chain of matrix multiplications

This paper is concerned with solving linear system (I_n+B_L?B₂B₁)x=b arising from the Green’s function calculation in the quantum Monte Carlo simulation of interacting electrons. The order of the system and integer L are adjustable. Also adjustable is the conditioning of the coefficient matrix to give rise an extreme ill-conditioned system. Two numerical methods based on the QR decomposition with column pivoting and the singular value decomposition, respectively, are studied in this paper. It is proved that the computed solution by each of the methods is weakly backward stable in the sense that the computed is close to the exact solution of a nearby linear system

     

Computing permanents via determinants for some classes of sparse matrices

Starting from recent formulas for calculating the permanents of some sparse circulant matrices, we obtain more general formulas expressing the permanents of a wider class of matrices as a linear combination of appropriate determinants.     

Characterizations of the polynomial numerical hull of degree k

Six characterizations of the polynomial numerical hull of degree k are established for bounded linear operators on a Hilbert space. It is shown how these characterizations provide a natural distinction between interior and boundary points. One of the characterizations is used to prove that the polynomial numerical hull of any fixed degree k for a Toeplitz matrix whose symbol is piecewise continuous approaches all or most of that of the infinite-dimensional Toeplitz operator, as the matrix size goes to infinity.