A <Emphasis Type="Italic">QR</Emphasis>–method for computing the singular values via semiseparable matrices

Summary. The standard procedure to compute the singular value decomposition of a dense matrix, first reduces it into a bidiagonal one by means of orthogonal transformations. Once the bidiagonal matrix has been computed, the QR–method is applied to reduce the latter matrix into a diagonal one. In this paper we propose a new method for computing the singular value decomposition of a real matrix. In a first phase, an algorithm for reducing the matrix A into an upper triangular semiseparable matrix by means of orthogonal transformations is described. A remarkable feature of this phase is that, depending on the distribution of the singular values, after few steps of the reduction, the largest singular values are already computed with a precision depending on the gaps between the singular values. An implicit QR–method for upper triangular semiseparable matrices is derived and applied to the latter matrix for computing its singular values. The numerical tests show that the proposed method can compete with the standard method (using an intermediate bidiagonal matrix) for computing the singular values of a matrix.Mathematics Subject Classification (2000): 65F15, 15A18The research of the first two authors was partially supported by the Research Council K.U.Leuven, project OT/00/16 (SLAP: Structured Linear Algebra Package), by the Fund for Scientific Research–Flanders (Belgium), projects G.0078.01 (SMA: Structured Matrices and their Applications), G.0176.02 (ANCILA: Asymptotic aNalysis of the Convergence behavior of Iterative methods in numerical Linear Algebra), G.0184.02 (CORFU: Constructive study of Orthogonal Functions) and G.0455.0 (RHPH: Riemann-Hilbert problems, random matrices and Padé-Hermite approximation), and by the Belgian Programme on Interuniversity Poles of Attraction, initiated by the Belgian State, Prime Ministers Office for Science, Technology and Culture, project IUAP V-22 (Dynamical Systems and Control: Computation, Identification & Modelling). The work of the third author was partially supported by MIUR, grant number 2002014121. The scientific responsibility rests with the authors.Acknowledgments.We thank the referees for their suggestions which increased the readability of the paper.     

Accurate singular values and differential qd algorithms

Summary. We have discovered a new implementation of the qd algorithm that has a far wider domain of stability than Rutishauser's version. Our algorithm was developed from an examination of the {Cholesky~LR} transformation and can be adapted to parallel computation in stark contrast to traditional qd. Our algorithm also yields useful a posteriori upper and lower bounds on the smallest singular value of a bidiagonal matrix. The zero-shift bidiagonal QR of Demmel and Kahan computes the smallest singular values to maximal relative accuracy and the others to maximal absolute accuracy with little or no degradation in efficiency when compared with the LINPACK code. Our algorithm obtains maximal relative accuracy for all the singular values and runs at least four times faster than the LINPACK code. Received August 8, 1993/Revised version received May 26, 1993     

Results on the relative perturbation of the singular values of a matrix

A highly accurate computation of the singular values of a matrix is a topic of current interest in the literature. In this paper we develop general bounds on relative perturbation of singular values. These bounds permit slight improvements in a unified derivation of some previous inequalities. The main result is a better criterion to neglect off-diagonal elements in the bidiagonal singular values decomposition.The present paper has been developed under the M. U. R. S. T. 40% National Program and the Interuniversity Center of Numerical Analysis and Computational Mathematics.     

Computing the complete CS decomposition

An algorithm for computing the complete CS decomposition of a partitioned unitary matrix is developed. Although the existence of the CS decomposition (CSD) has been recognized since 1977, prior algorithms compute only a reduced version. This reduced version, which might be called a 2-by-1 CSD, is equivalent to two simultaneous singular value decompositions. The algorithm presented in this article computes the complete 2-by-2 CSD, which requires the simultaneous diagonalization of all four blocks of a unitary matrix partitioned into a 2-by-2 block structure. The algorithm appears to be the only fully specified algorithm available. The computation occurs in two phases. In the first phase, the unitary matrix is reduced to bidiagonal block form, as described by Sutton and Edelman. In the second phase, the blocks are simultaneously diagonalized using techniques from bidiagonal SVD algorithms of Golub, Kahan, Reinsch, and Demmel. The algorithm has a number of desirable numerical features.      

Estimating the largest singular values of large sparse matrices via modified moments

We describe a procedure for determining a few of the largest singular values of a large sparse matrix. The method by Golub and Kent which uses the method of modified moments for estimating the eigenvalues of operators used in iterative methods for the solution of linear systems of equations is appropriately modified in order to generate a sequence of bidiagonal matrices whose singular values approximate those of the original sparse matrix. A simple Lanczos recursion is proposed for determining the corresponding left and right singular vectors. The potential asynchronous computation of the bidiagonal matrices using modified moments with the iterations of an adapted Chebyshev semi-iterative (CSI) method is an attractive feature for parallel computers. Comparisons in efficiency and accuracy with an appropriate Lanczos algorithm (with selective re-orthogonalization) are presented on large sparse (rectangular) matrices arising from applications such as information retrieval and seismic reflection tomography. This procedure is essentially motivated by the theory of moments and Gauss quadrature.This author's work was supported by the National Science Foundation under grants NSF CCR-8717492 and CCR-910000N (NCSA), the U.S. Department of Energy under grant DOE DE-FG02-85ER25001, and the Air Force Office of Scientific Research under grant AFOSR-90-0044 while at the University of Illinois at Urbana-Champaign Center for Supercomputing Research and Development.This author's work was supported by the U.S. Army Research Office under grant DAAL03-90-G-0105, and the National Science Foundation under grant NSF DCR-8412314.     

An implicit QR algorithm for symmetric semiseparable matrices

The QR algorithm is one of the classical methods to compute the eigendecomposition of a matrix. If it is applied on a dense n × n matrix, this algorithm requires O(n³) operations per iteration step. To reduce this complexity for a symmetric matrix to O(n), the original matrix is first reduced to tridiagonal form using orthogonal similarity transformations. In the report (Report TW360, May 2003) a reduction from a symmetric matrix into a similar semiseparable one is described. In this paper a QR algorithm to compute the eigenvalues of semiseparable matrices is designed where each iteration step requires O(n) operations. Hence, combined with the reduction to semiseparable form, the eigenvalues of symmetric matrices can be computed via intermediate semiseparable matrices, instead of tridiagonal ones. The eigenvectors of the intermediate semiseparable matrix will be computed by applying inverse iteration to this matrix. This will be achieved by using an O(n) system solver, for semiseparable matrices. A combination of the previous steps leads to an algorithm for computing the eigenvalue decompositions of semiseparable matrices. Combined with the reduction of a symmetric matrix towards semiseparable form, this algorithm can also be used to calculate the eigenvalue decomposition of symmetric matrices. The presented algorithm has the same order of complexity as the tridiagonal approach, but has larger lower order terms. Numerical experiments illustrate the complexity and the numerical accuracy of the proposed method. Copyright © 2005 John Wiley & Sons, Ltd.     

Superquadratic convergence of DLASQ for computing matrix singular values

DLASQ is a routine in LAPACK for computing the singular values of a real upper bidiagonal matrix with high accuracy. The basic algorithm, the so-called dqds algorithm, was first presented by Fernando-Parlett, and implemented as the DLASQ routine by Parlett-Marques. DLASQ is now recognized as one of the most efficient routines for computing singular values. In this paper, we prove the asymptotic superquadratic convergence of DLASQ in exact arithmetic.     

Computation of exact inertia and inclusions of eigenvalues (singular values) of tridiagonal (bidiagonal) matrices

This report may be considered as a non-trivial extension of an unpublished report by William Kahan (Accurate Eigenvalues of a symmetric tri-diagonal matrix, Technical Report CS 41, Computer Science Department, Stanford University, 1966). His interplay between matrix theory and computer arithmetic led to the development of algorithms for computing accurate eigenvalues and singular values. His report is generally considered as the precursor for the development of IEEE standard 754 for binary arithmetic. This standard has been universally adopted by virtually all PC, workstation and midrange hardware manufactures and tens of billions of such machines have been produced. Now we use the features in this standard to improve the original algorithm.In this paper, we describe an algorithm in floating-point arithmetic to compute the exact inertia of a real symmetric (shifted) tridiagonal matrix. The inertia, denoted by the integer triplet (π, ν, ζ), is defined as the number of positive, negative and zero eigenvalues of a real symmetric (or complex Hermitian) matrix and the adjective exact refers to the eigenvalues computed in exact arithmetic. This requires the floating-point computation of the diagonal matrix D of the LDL^t factorization of the shifted tridiagonal matrix T − τI with +∞ and −∞ rounding modes defined in IEEE 754 standard. We are not aware of any other algorithm which gives the exact answer to a numerical problem when implemented in floating-point arithmetic in standard working precisions. The guaranteed intervals for eigenvalues are obtained by bisection or multisection with this exact inertia information. Similarly, using the Golub-Kahan form, guaranteed intervals for singular values of bidiagonal matrices can be computed. The diameter of the eigenvalue (singular value) intervals depends on the number of shifts with inconsistent inertia in two rounding modes. Our algorithm not only guarantees the accuracy of the solutions but is also consistent across different IEEE 754 standard compliant architectures. The unprecedented accuracy provided by our algorithms could be also used to debug and validate standard floating-point algorithms for computation of eigenvalues (singular values). Accurate eigenvalues (singular values) are also required by certain algorithms to compute accurate eigenvectors (singular vectors).We demonstrate the accuracy of our algorithms by using standard matrix examples. For the Wilkinson matrix, the eigenvalues (in IEEE double precision) are very accurate with an (open) interval diameter of 6 ulps (units of the last place held of the mantissa) for one of the eigenvalues and lesser (down to 2 ulps) for others. These results are consistent across many architectures including Intel, AMD, SGI and DEC Alpha. However, by enabling IEEE double extended precision arithmetic in Intel/AMD 32-bit architectures at no extra computational cost, the (open) interval diameters were reduced to one ulp, which is the best possible solution for this problem. We have also computed the eigenvalues of a tridiagonal matrix which manifests in Gauss-Laguerre quadrature and the results are extremely good in double extended precision but less so in double precision. To demonstrate the accuracy of computed singular values, we have also computed the eigenvalues of the Kac₃₀ matrix, which is the Golub-Kahan form of a bidiagonal matrix. The tridiagonal matrix has known integer eigenvalues. The bidiagonal Cholesky factor of the Gauss-Laguerre tridiagonal is also included in the singular value study.     

Orthogonal similarity transformation into block‐semiseparable matrices of semiseparability rank k

Very recently, an algorithm, which reduces any symmetric matrix into a semiseparable one of semi‐ separability rank 1 by similar orthogonality transformations, has been proposed by Vandebril, Van Barel and Mastronardi. Partial execution of this algorithm computes a semiseparable matrix whose eigenvalues are the Ritz‐values obtained by the Lanczos' process applied to the original matrix. Also a kind of nested subspace iteration is performed at each step. In this paper, we generalize the above results and propose an algorithm to reduce any symmetric matrix into a similar block‐semiseparable one of semiseparability rank k, with k ∈ ?, by orthogonal similarity transformations. Also in this case partial execution of the algorithm computes a block‐semiseparable matrix whose eigenvalues are the Ritz‐values obtained by the block‐Lanczos' process with k starting vectors, applied to the original matrix. Subspace iteration is performed at each step as well. Copyright © 2005 John Wiley & Sons, Ltd.     

Structural and computational properties of possibly singular semiseparable matrices

A classical result of structured numerical linear algebra states that the inverse of a nonsingular semiseparable matrix is a tridiagonal matrix. Such a property of a semiseparable matrix has been proved to be useful for devising linear complexity solvers, for establishing recurrence relations among its columns or rows and, moreover, for efficiently evaluating its characteristic polynomial. In this paper, we provide sparse structured representations of a semiseparable matrix A which hold independently of the fact that A is singular or not. These relations are found by pointing out the band structure of the inverse of the sum of A plus a certain sparse perturbation of minimal rank. Further, they can be used to determine in a computationally efficient way both a reflexive generalized inverse of A and its characteristic polynomial.     

Equality conditions for lower bounds on the smallest singular value of a bidiagonal matrix

Several lower bounds have been proposed for the smallest singular value of a square matrix, such as Johnson’s bound, Brauer-type bound, Li’s bound and Ostrowski-type bound. In this paper, we focus on a bidiagonal matrix and investigate the equality conditions for these bounds. We show that the former three bounds give strict lower bounds if all the bidiagonal elements are non-zero. For the Ostrowski-type bound, we present an easily verifiable necessary and sufficient condition for the equality to hold.     

基于反Krylov矩阵正交分解的半可分矩阵

在本文中，我们证明了对一个反Krylov矩阵作QR分解后，利用得到的正交矩阵可以将一个具有互异特征值的对称矩阵转化为一个半可分矩阵的形式，这个结果表明了反Krylov矩阵与半可分矩阵之间的联系．另外，我们还证明了这类对称半可分矩阵在QR达代下矩阵结构保持不变性．     

Orthogonal similarity transformation of a symmetric matrix into a diagonal-plus-semiseparable one with free choice of the diagonal

In this paper we describe an orthogonal similarity transformation for transforming arbitrary symmetric matrices into a diagonal-plus-semiseparable matrix, where we can freely choose the diagonal. Very recently an algorithm was proposed for transforming arbitrary symmetric matrices into similar semiseparable ones. This reduction is strongly connected to the reduction to tridiagonal form. The class of semiseparable matrices can be considered as a subclass of the diagonalplus- semiseparable matrices. Therefore we can interpret the proposed algorithm here as an extension of the reduction to semiseparable form. A numerical experiment is performed comparing thereby the accuracy of this reduction algorithm with respect to the accuracy of the traditional reduction to tridiagonal form, and the reduction to semiseparable form. The experiment indicates that all three reduction algorithms are equally accurate. Moreover it is shown in the experiments that asymptotically all the three approaches have the same complexity, i.e. that they have the same factor preceding the n³ term in the computational complexity. Finally we illustrate that special choices of the diagonal create a specific convergence behavior. The research was partially supported by the Research Council K.U.Leuven, project OT/05/40 (Large rank structured matrix computations), by the Fund for Scientific Research–Flanders (Belgium), projects G.0078.01 (SMA: Structured Matrices and their Applications), G.0176.02 (ANCILA: Asymptotic aNalysis of the Convergence behavior of Iterative methods in numerical Linear Algebra), G.0184.02 (CORFU: Constructive study of Orthogonal Functions) and G.0455.0 (RHPH: Riemann-Hilbert problems, random matrices and Padé-Hermite approximation), and by the Belgian Programme on Interuniversity Poles of Attraction, initiated by the Belgian State, Prime Minister's Office for Science, Technology and Culture, project IUAP V-22 (Dynamical Systems and Control: Computation, Identification & Modelling). The scientific responsibility rests with the authors.     

On the convergence properties of the orthogonal similarity transformations to tridiagonal and semiseparable (plus diagonal) form

In this paper, we will compare the convergence properties of three basic reduction methods, by placing them in a general framework. It covers the reduction to tridiagonal, semiseparable and semiseparable plus diagonal form. These reductions are often used as the first step in the computation of the eigenvalues and/or eigenvectors of arbitrary matrices. In this way, the calculation of the eigenvalues using, for example, the QR-algorithm reduces in complexity. First we will investigate the convergence properties of these three reduction algorithms. It will be shown that for the partially reduced matrices at step k of any of these reduction algorithms, the lower right k × k (already reduced) sub-block will have the Lanczos–Ritz values, w.r.t. a certain starting vector. It will also be shown that the reductions to semiseparable and to semiseparable plus diagonal form have an extra convergence behavior a special type of subspace iteration is performed on the lower right k × k submatrix, which contains these Ritz-values. Secondly we look in more detail at the behavior of the involved subspace iteration. It will be shown that the reduction method can be interpreted as a nested type of multi-shift iteration. Theoretical results will be presented, making it possible to predict the convergence behavior of these reduction algorithms. Also a theoretical bound on the convergence rate is presented. Finally we illustrate by means of numerical examples, how it is possible to tune the convergence behavior such that it can become a powerful tool for certain applications.     

An efficient and accurate parallel algorithm for the singular value problem of bidiagonal matrices

Summary. In this paper we propose an algorithm based on Laguerre's iteration, rank two divide-and-conquer technique and a hybrid strategy for computing singular values of bidiagonal matrices. The algorithm is fully parallel in nature and evaluates singular values to tiny relative error if necessary. It is competitive with QR algorithm in serial mode in speed and advantageous in computing partial singular values. Error analysis and numerical results are presented. Received March 15, 1993 / Revised version received June 7, 1994     

Accurate bidiagonal decomposition of collocation matrices of weighted ϕ‐transformed systems

Given a system of functions, we introduce the concept of weighted φ‐transformed system, which will include a very large class of useful representations in Statistics and Computer Aided Geometric Design. An accurate bidiagonal decomposition of the collocation matrices of these systems is obtained. This decomposition is used to present computational methods with high relative accuracy for solving algebraic problems with collocation matrices of weighted φ‐transformed systems such as the computation of eigenvalues, singular values, and the solution of some linear systems. Numerical examples illustrate the accuracy of the performed computations.     

Error analysis of the mdLVs algorithm for computing bidiagonal singular values

Some of the authors show in Iwasaki and Nakamura (Inverse Probl 20:553?C563, 2004) that the integrable discrete Lotka?CVolterra (dLV) system is applicable for computing singular values of bidiagonal matrix. The resulting numerical algorithm is referred to as the dLV algorithm. They also observe in Iwasaki and Nakamura (Electron Trans Numer Anal 38:184?C201, 2011) that the singular values are numerically computed with high relative accuracy by using the mdLVs algorithm, which is an acceleration version by introducing a shift of origin. In this paper, we investigate the perturbations on singular values and the forward errors of the mdLVs variables, which occur in the mdLVs algorithm, through two kinds of error analysis in floating point arithmetic. Therefore the forward stability of the mdLVs algorithm in the sense of Bueno?CMarcellan?CDopico is proved.     

Fast structured LU factorization for nonsymmetric matrices

In this paper, an approximate LU factorization algorithm is developed for nonsymmetric matrices based on the hierarchically semiseparable matrix techniques. It utilizes a technique involving orthogonal transformations and approximations to avoid the explicit computation of the Schur complement in each factorization step. A modified compression method is further developed for the case when some diagonal blocks have small singular values. The complexity of the methods proposed in this paper is analyzed and shown to be $O(N^2k)$ , where $N$ is the dimension of matrix and $k$ is the maximum off-diagonal (numerical) rank. Depending on the accuracy and efficiency requirements in the approximation, this factorization can be used either as a direct solver or a preconditioner. Numerical results from applications are included to show the efficiency of our method.     

Accurate computations with Laguerre matrices

This paper provides an accurate method to obtain the bidiagonal factorization of collocation matrices of generalized Laguerre polynomials and of Lah matrices, which in turn can be used to compute with high relative accuracy the eigenvalues, singular values, and inverses of these matrices. Numerical examples are included.