Search

The roundoff-error-free (REF) LU and Cholesky factorizations, combined with the REF substitution algorithms, allow rational systems of linear equations to be solved exactly and efficiently by working entirely in integer arithmetic. These REF computational tools share two key properties: their constituent divisions are exact, and their matrix entries' bit-lengths (i.e., required number of storage bits) are bounded polynomially. However, as this paper explains, updating the REF factorizations via the common delete-insert-reduce approach turns out to be inefficient in terms of operand bit-length growth and increased computational effort. In fact, we prove that this inefficiency can eliminate the computational savings expected of the factorization update process. Hence, the current work develops REF update algorithms that differ significantly from their traditional counterparts. To achieve this, it introduces the frame matrix methodology for representing and applying the REF factorizations, and then it devises new REF operations and shortcuts for transitioning between "nearby" frame matrices. The presented push-and-swap update approach utilizes these transitions to preserve the special structure of the REF factorizations and to avoid additional growth in their matrix entries' bit-lengths; this update process requires the expected O(n²) operations. The featured updates are column addition, deletion, and replacement with respect to the REF LU factorization. We also prove that analogous row updates can be performed via the column updates and discuss special considerations for updating the REF Cholesky factorization. An accompanying set of computational tests attests that for fully dense random instances, the push-and-swap column replacement update for the REF LU factorization is on average an order of magnitude faster than the corresponding delete-insert-reduce update for the exact rational arithmetic LU factorization. Altogether, the REF update algorithms represent a promising preliminary step toward enhancing the basic solution validation subroutines employed by state-of-the-art mixed-precision methodologies for exact linear and mixed-integer programming. [ABSTRACT FROM AUTHOR]

The development of strong linearizations preserving whatever structure a matrix polynomial might possess has been a very active area of research in the last years, since such linearizations are the starting point of numerical algorithms for computing eigenvalues of structured matrix polynomials with the properties imposed by the considered structure. In this context, Hermitian matrix polynomials are one of the most important classes of matrix polynomials arising in applications and their real eigenvalues are of great interest. The sign characteristic is a set of signs attached to these real eigenvalues which is crucial for determining the behavior of systems described by Hermitian matrix polynomials and, therefore, it is desirable to develop linearizations that preserve the sign characteristic of these polynomials, but, at present, only one such linearization is known. In this paper, we present a complete characterization of all the Hermitian strong linearizations that preserve the sign characteristic of a given Hermitian matrix polynomial and identify several families of such linearizations that can be constructed very easily from the coefficients of the polynomial. [ABSTRACT FROM AUTHOR]

Fast algorithms for matrix multiplication, namely those that perform asymptotically fewer scalar operations than the classical algorithm, have been considered primarily of theoretical interest. Apart from Strassen's original algorithm, few fast algorithms have been efficiently im- plemented or used in practical applications. However, there exist many practical alternatives to Strassen's algorithm with varying performance and numerical properties. Fast algorithms are known to be numerically stable, but because their error bounds are slightly weaker than the classical al- gorithm, they are not used even in cases where they provide a performance benefit. We argue in this paper that the numerical sacrifice of fast algorithms, particularly for the typical use cases of practical algorithms, is not prohibitive, and we explore ways to improve the accuracy both theoret- ically and empirically. The numerical accuracy of fast matrix multiplication depends on properties of the algorithm and of the input matrices, and we consider both contributions independently. We generalize and tighten previous error analyses of fast algorithms and compare their properties. We discuss algorithmic techniques for improving the error guarantees from two perspectives: manipulat- ing the algorithms, and reducing input anomalies by various forms of diagonal scaling. Finally, we benchmark performance and demonstrate our improved numerical accuracy. [ABSTRACT FROM AUTHOR]

Randomized sampling has recently been proven a highly efficient technique for computing approximate factorizations of matrices that have low numerical rank. This paper describes an extension of such techniques to a wider class of matrices that are not themselves rank-deficient but have off-diagonal blocks that are; specifically, the class of so-called hierarchically semiseparable (HSS) matrices. HSS matrices arise frequently in numerical analysis and signal processing, particularly in the construction of fast methods for solving differential and integral equations numerically. The HSS structure admits algebraic operations (matrix-vector multiplications, matrix factorizations, matrix inversion, etc.) to be performed very rapidly, but only once the HSS representation of the matrix has been constructed. How to rapidly compute this representation in the first place is much less well understood. The present paper demonstrates that if an N x N matrix can be applied to a vector in O(N) time, and if individual entries of the matrix can be computed rapidly, then provided that an HSS representation of the matrix exists, it can be constructed in O(N k²) operations, where k is an upper bound for the numerical rank of the off-diagonal blocks. The point is that when legacy codes (based on, e.g., the fast multipole method) can be used for the fast matrix-vector multiply, the proposed algorithm can be used to obtain the HSS representation of the matrix, and then well-established techniques for HSS matrices can be used to invert or factor the matrix. [ABSTRACT FROM AUTHOR]

In this paper, we show how to perform the Hessenberg reduction of a rank structured matrix under unitary similarity operations in an efficient way, using the Givens-weight representation. This reduction can be used as a first step for eigenvalue computation. We also show how the algorithm can be modified to compute the bidiagonal reduction of a rank structured matrix, this latter method being a preprocessing step for computing the singular values of the matrix. For the main cases of interest, the algorithms we describe in this paper are of complexity O((ar + bs)n2), where n is the matrix size, r is some measure for the average rank index of the rank structure, s is some measure for the bandwidth of the unstructured matrix part around the main diagonal, and a, b ? R are certain weighting parameters. Numerical experiments demonstrate the stability of this approach. [ABSTRACT FROM AUTHOR]

In this paper we use a displacement structure approach to design a class of new preconditioners for the conjugate gradient method applied to the solution of large Toeplitz linear equations. Explicit formulas are suggested for the G.Strang-type and for the T.Chan-type preconditioners belonging to any of eight classes of matrices diagonalized by the corresponding discrete cosine or sine transforms. Under the standard Wiener class assumption the clustering property is established for all of these preconditioners, guaranteeing rapid convergence of the preconditioned con- jugate gradient method. All the computations related to the new preconditioners can be done in real arithmetic, and to fully exploit this advantageous property one has to suggest a fast real-arithmetic algorithm for multiplication of a Toeplitz matrix by a vector. It turns out that the obtained formulas for the Strang-type preconditioners allow a number of representations for Toeplitz matrices leading to a wide variety of real-arithmetic multiplication algorithms based on any of eight discrete cosine or sine transforms. Recently, transformations of Toeplitz matrices to Vandermonde-like or Cauchy-like matrices have been found to be useful in developing accurate direct methods for Toeplitz linear equations. In this paper we suggest further extending the range of the transformation approach by exploring it for iterative methods; this technique allowed us to reduce the complexity of each iteration of the preconditioned conjugate gradient method. The results of this paper were announced in [T. Kailath and V. Olshevsky, Calcolo, 33 (1996), pp. 191-208]. [ABSTRACT FROM AUTHOR]

In (SIAM J. Matrix Anal. Appl., 23 (2002), pp. 706–727, we developed numerical algorithms for computing sparse low-rank approximations of matrices, and we also provided a detailed error analysis of the proposed algorithms together with some numerical experiments. The low-rank approximations are constructed in a certain factored form with the degree of sparsity of the factors controlled by some user-specified parameters. In this paper, we cast the sparse low-rank approximation problem in the framework of penalized optimization problems. We discuss various approximation schemes for the penalized optimization problem which are more amenable to numerical computations. We also include some analysis to show the relations between the original optimization problem and the reduced one. We then develop a globally convergent discrete Newton-like iterative method for solving the approximate penalized optimization problems. We also compare the reconstruction errors of the sparse low-rank approximations computed by our new methods with those obtained using the methods in the earlier paper and several other existing methods for computing sparse low-rank approximations. Numerical examples show that the penalized methods are more robust and produce approximations with factors which have fewer columns and are sparser. [ABSTRACT FROM AUTHOR]

This paper presents a new algorithm for accurate floating-point computation of the singular value decomposition (SVD) of the product A = BτSC, where B ∈ Rp×m, C ∈ Rq×n, S ∈ Rp×q, and p ≤ m, q ≤ n. The new algorithm uses diagonal scalings, the LU factorization with complete pivoting, the QR factorization with column pivoting, and matrix multiplication to replace A by A′ = B′τS′C′, where A and A′ have the same singular values and the matrix A′ is computed explicitly. The singular values of A′ are computed using the Jacobi SVD algorithm. It is shown that the accuracy of the new algorithm is determined by (i) the accuracy of the QR factorizations of Bτ and Cτ; (ii) the accuracy of the LU factorization with complete pivoting of S; and (iii) the accuracy of the computation of the SVD of a matrix A′ with moderate minD=diagκ2(A′D). Theoretical analysis and numerical evidence show that, in the case of rank(B) = rank(C) = p and full rank S, the accuracy of the new algorithm is unaffected by replacing B, S, C with, respectively, D1B, D2SD3, D4C, where Di, i = 1,…,4, are arbitrary diagonal matrices. As an application, the paper proposes new accurate algorithms for computing the (H,K)­SVD and (H-1,K)­SVD of S. [ABSTRACT FROM AUTHOR]

This paper considers an implicitly restarted Krylov subspace method that approximates a stable, linear transfer function f(s) of order n by one of order m, where n >> m. It is well known that oblique projections onto a Krylov subspace may generate unstable partial realizations. To remedy this situation, the oblique projectors obtained via classical Krylov subspace methods are supplemented with further projectors which enable the formation of stable partial realizations directly from f(s). A key feature of this process is that it may be incorporated into an implicit restart scheme. A second difficulty arises from the fact that Krylov subspace methods often generate partial realizations that contain nonessential modes. To this end, balanced truncation may be employed to discard the unwanted part of the reduced-order model. This paper proposes oblique projection methods for large-scale model reduction that simultaneously compute stable reduced-order models while discarding all nonessential modes. It is shown that both of these tasks may be effected by a single oblique projection process. Furthermore, the process is shown to naturally fit into an implicit restart framework. The theoretical properties of these methods are thoroughly investigated, and exact low-dimensional expressions for the L∞-norm of the residual errors are derived. Finally, the behavior of the algorithm is illustrated on two large-scale examples. [ABSTRACT FROM AUTHOR]

The generalized null space decomposition (GNSD) is a unitary reduction of a general matrix A of order n to a block upper triangular form that reveals the structure of the Jordan blocks of A corresponding to a zero eigenvalue. The reduction was introduced by Kublanovskaya. It was extended first by Ruhe and then by Golub and Wilkinson, who based the reduction on the singular value decomposition. Unfortunately, if A has large Jordan blocks, the complexity of these algorithms can approach the order of n4. This paper presents an alternative algorithm, based on repeated updates of a QR decomposition of A, that is guaranteed to be of order n³. Numerical experiments confirm the stability of this algorithm, which turns out to produce essentially the same form as that of Golub and Wilkinson. The effect of errors in A on the ability to recover the original structure is investigated empirically. Several applications are discussed, including the computation of the Drazin inverse. [ABSTRACT FROM AUTHOR]

In this paper we discuss the deflation criterion used in the extended QR algorithm based on the chasing of rotations. We provide absolute and relative perturbation bounds for this deflation criterion. Further, we present a generalization of aggressive early deflation to the extended QR algorithms. Aggressive early deflation is the key technique for the identification and deflation of already converged, but hidden, eigenvalues. Often these possibilities for deflation are not detected by the standard technique. We present numerical results underpinning the power of aggressive early deflation also in the context of extended QR algorithms. We further generalize these ideas by the transcription of middle deflations. [ABSTRACT FROM AUTHOR]

In this paper we study the problem of finding a weighted geometric mean sufficiently close to the weighted Karcher mean of positive definite matrices. This problem arises in various numerical algorithms for computing the weighted Karcher mean and in a physical problem of averaging on the Riemannian manifold of positive definite matrices. We prove that this is possible for any "convex" geometric means satisfying the Jensen-type inequality for geodesically convex functions varying over an open set of weights in the simplex of positive probability vectors. This follows from a useful estimation of the Riemannian distance between the Karcher mean and any convex geometric means. A better upper bound between the Karcher mean and the Riemannian convex combination admitting an explicit formula is presented. [ABSTRACT FROM AUTHOR]

In this paper a new O(N log³ N) solver for N x N Toeplitz-like systems, based on a divide and conquer technique, is presented. Similarly to the superfast algorithm MBA for the inversion of a Toeplitz-like matrix [R. R. Bitmead and B. D. O. Anderson, Linear Algebra Appl., 34 (1980), pp. 103-116; M. Morf, Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing, 1980, pp. 954-959], it exploits the displacement properties. In order to avoid the well-known numerical instability of the explicit inversion, the new algorithm relies on the triangular factorization and back-substitution formula for the system seen as a 2x2 block system with blocks of half size. This idea is the one used in [M. Stewart, SIAM J. Matrix Anal. Appl., 25 (2003), pp. 669-693] to improve the numerical stability of superfast methods based on the generalized Schur algorithm for positive definite Toeplitz matrices, but the algorithm we propose can be applied also to nonsymmetric Toeplitz-like systems. The stability of the algorithm is examined through numerical experiments. [ABSTRACT FROM AUTHOR]

Recently a generalization of Francis's implicitly shifted QR algorithm was proposed, notably widening the class of matrices admitting low-cost implicit QR steps. This unifying framework covered the methods and theory for Hessenberg and inverse Hessenberg matrices and furnished also new, single-shifted, QR-type methods for, e.g., CMV matrices. Convergence of this approach was only suggested by numerical experiments. No theoretical claims supporting the results were presented. In this paper we present multishift variants of these new algorithms. We also provide a convergence theory that shows that the new algorithm performs nested subspace iterations on rational Krylov subspaces. Numerical experiments confirm the validity of the theory. [ABSTRACT FROM AUTHOR]

The dqds algorithm computes all the singular values of an n x n bidiagonal matrix to high relative accuracy in O(n²) cost. Its efficient implementation is now available as a LAPACK subroutine and is the preferred algorithm for this purpose. In this paper we incorporate into dqds a technique called aggressive early deflation, which has been applied successfully to the Hessenberg QR algorithm. Extensive numerical experiments show that aggressive early deflation often reduces the dqds runtime significantly. In addition, our theoretical analysis suggests that with aggressive early deflation, the performance of dqds is largely independent of the shift strategy. We confirm through experiments that the zero-shift version is often as fast as the shifted version. We give a detailed error analysis to prove that with our proposed deflation strategy, dqds computes all the singular values to high relative accuracy. [ABSTRACT FROM AUTHOR]

In this paper, based on Patel's algorithm (1993), we propose a structure-preserving algorithm for solving palindromic quadratic eigenvalue problems (QEPs). We also show the relationship between the structure-preserving algorithm and the URV-based structure-preserving algorithm by Schröder (2007). For large sparse palindromic QEPs, we develop a generalized T-skew-Hamiltonian implicitly restarted shift-and-invert Arnoldi algorithm for solving the resulting T-skew-Hamiltonian pencils. Numerical experiments show that our proposed structure-preserving algorithms perform well on the palindromic QEP arising from a finite element model of high-speed trains and rails. [ABSTRACT FROM AUTHOR]

We consider the problem of performing accurate computations with rectangular (m × n) totally nonnegative matrices. The matrices under consideration have the property of having a unique representation as products of nonnegative bidiagonal matrices. Given that representation, one can compute the inverse, LDU decomposition, eigenvalues, and SVD of a totally nonnegative matrix to high relative accuracy in O(max(m3, n3)) time-much more accurately than conventional algorithms that ignore that structure. The contribution of this paper is to show that the high relative accuracy is preserved by operations that preserve the total nonnegativity-taking a product, re-signed inverse (when m = n), converse, Schur complement, or submatrix of a totally nonnegative matrix, any of which costs at most O(max(m3, n3)). In other words, the class of totally nonnegative matrices for which we can do numerical linear algebra very accurately in O(max(m3, n3)) time (namely, those for which we have a product representation via nonnegative bidiagonals) is closed under the operations listed above. [ABSTRACT FROM AUTHOR]

This paper proves the stability of a variant of Gragg's unitary Hessenberg QR algorithm (UHQR). It is shown that a single UHQR iteration with numerically unimodular shift applied to the Schur parameters and complementary parameters of a unitary Hessenberg matrix H gives computed parameters for QHHQ that are close to those obtained from a perturbed set of parameters using a perturbed shift with no numerical error. The perturbations of the parameters and the shift are bounded. [ABSTRACT FROM AUTHOR]

Algebraic multigrid (AMG) is a very efficient iterative solver and preconditioner for large unstructured sparse linear systems. Traditional coarsening schemes for AMG can, however, lead to computational complexity growth as problem size increases, resulting in increased memory use and execution time, and diminished scalability. Two new parallel AMG coarsening schemes are proposed that are based solely on enforcing a maximum independent set property, resulting in sparser coarse grids. The new coarsening techniques remedy memory and execution time complexity growth for various large three-dimensional (3D) problems. If used within AMG as a preconditioner for Krylov subspace methods, the resulting iterative methods tend to converge fast. This paper discusses complexity issues that can arise in AMG, describes the new coarsening schemes, and examines the performance of the new preconditioners for various large 3D problems. [ABSTRACT FROM AUTHOR]

In this paper, we present new algorithms that can replace the diagonal entries of a Hermitian matrix by any set of diagonal entries that majorize the original set without altering the eigenvalues of the matrix. They perform this feat by applying a sequence of (N-1) or fewer plane rotations, where N is the dimension of the matrix. Both the Bendel--Mickey and the Chan--Li algorithms are special cases of the proposed procedures. Using the fact that a positive semidefinite matrix can always be factored as $\mtx{X^\adj X}$, we also provide more efficient versions of the algorithms that can directly construct factors with specified singular values and column norms. We conclude with some open problems related to the construction of Hermitian matrices with joint diagonal and spectral properties. [ABSTRACT FROM AUTHOR]

In this paper we show how to compute recursively an approximation of the left and right dominant singular subspaces of a given matrix. In order to perform as few as possible operations on each column of the matrix, we use a variant of the classical Gram­Schmidt algorithm to estimate this subspace. The method is shown to be particularly suited for matrices with many more rows than columns. Bounds for the accuracy of the computed subspace are provided. Moreover, the analysis of error propagation in this algorithm provides new insights in the loss of orthogonality typically observed in the classical Gram­Schmidt method. [ABSTRACT FROM AUTHOR]

In the convergence theory of multisplittings for symmetric positive definite (s.p.d.) matrices it is usually assumed that the weighting matrices are scalar matrices, i.e., multiples of the identity. In this paper, this restrictive condition is eliminated. In its place it is assumed that more than one (inner) iteration is performed in each processor (or block). The theory developed here is applied to nonstationary multisplittings for s.p.d. matrices, as well as to two-stage multisplittings for symmetric positive semidefinite matrices. [ABSTRACT FROM AUTHOR]

Various forms of sparse approximate inverses (SAI) have been shown to be useful for preconditioning. Their potential usefulness in a parallel environment has motivated much interest in recent years. However, the capability of an approximate inverse in eliminating the local error has not yet been fully exploited in multigrid algorithms. A careful examination of the iteration matrices of these approximate inverses indicates their superiority in smoothing the high-frequency error in addition to their inherent parallelism. We propose a new class of SAI smoothers in this paper and present their analytic smoothing factors for constant coefficient PDEs. The following are several distinctive features that make this technique special: By adjusting the quality of the approximate inverse, the smoothing factor can be improved accordingly. For hard problems, this is useful. In contrast to the ordering sensitivity of other smoothing techniques, this technique is ordering independent. In general, the sequential performance of many superior parallel algorithms is not very competitive. This technique is useful in both parallel and sequential computations. Our theoretical and numerical results have demonstrated the effectiveness of this new technique. [ABSTRACT FROM AUTHOR]

The purpose of this paper is to derive new computable convergence bounds for GMRES. The new bounds depend on the initial guess and are thus conceptually different from standard “worst-case” bounds. Most importantly, approximations to the new bounds can be computed from information generated during the run of a certain GMRES implementation. The approximations allow predictions of how the algorithm will perform. Heuristics for such predictions are given. Numerical experiments illustrate the behavior of the new bounds as well as the use of the heuristics. [ABSTRACT FROM AUTHOR]

The BR algorithm, a new method for calculating the eigenvalues of an upper Hes- senberg matrix, is introduced. It is a bulge-chasing algorithm like the QR algorithm, but, unlike the QR algorithm, it is well adapted to computing the eigenvalues of the narrow-band, nearly tridiagonal matrices generated by the look-ahead Lanczos process. This paper describes the BR algorithm and gives numerical evidence that it works well in conjunction with the Lanczos process. On the biggest problems run so far, the BR algorithm beats the QR algorithm by a factor of 30–60 in computing time and a factor of over 100 in matrix storage space. [ABSTRACT FROM AUTHOR]