Simultaneous multidiagonalization for the CS decomposition

When an orthogonal matrix is partitioned into a two-by-two block structure, its four blocks can be simultaneously bidiagonalized. This observation underlies numerically stable algorithms for the CS decomposition and the existence of CMV matrices for orthogonal polynomial recurrences. We discover a new matrix decomposition for simultaneous multidiagonalization, which reduces the blocks to any desired bandwidth. Its existence is proved, and a backward stable algorithm is developed. The resulting matrix with banded blocks is parameterized by a product of Givens rotations, guaranteeing orthogonality even on a finite-precision computer. The algorithm relies heavily on Level 3 BLAS routines and supports parallel computation.     

Computing the CS and the generalized singular value decompositions

Summary If the columns of a matrix are orthonormal and it is partitioned into a 2-by-1 block matrix, then the singular value decompositions of the blocks are related. This is the essence of the CS decomposition. The computation of these related SVD's requires some care. Stewart has given an algorithm that uses the LINPACK SVD algorithm together with a Jacobitype clean-up operation on a cross-product matrix. Our technique is equally stable and fast but avoids the cross product matrix. The simplicity of our technique makes it more amenable to parallel computation on systolic-type computer architectures. These developments are of interest because a good way to compute the generalized singular value decomposition of a matrix pair (A, B) is to compute the CS decomposition of a certain orthogonal column matrix related toA andB.The research associated with this paper was partially supported by the Office of Naval Research contract N00014-83-K-0640, USA     

Algorithms for Finding the Inverses of Factor Block Circulant Matrices

In this paper,algorithms for finding the inverse of a factor block circulant matrix, a factor block retrocirculant matrix and partitioned matrix with factor block circulant blocks over the complex field are presented respectively.In addition,two algorithms for the inverse of a factor block circulant matrix over the quaternion division algebra are proposed.     

A real-valued block conjugate gradient type method for solving complex symmetric linear systems with multiple right-hand sides

We consider solving complex symmetric linear systems with multiple right-hand sides. We assume that the coefficient matrix has indefinite real part and positive definite imaginary part. We propose a new block conjugate gradient type method based on the Schur complement of a certain 2-by-2 real block form. The algorithm of the proposed method consists of building blocks that involve only real arithmetic with real symmetric matrices of the original size. We also present the convergence property of the proposed method and an efficient algorithmic implementation. In numerical experiments, we compare our method to a complex-valued direct solver, and a preconditioned and nonpreconditioned block Krylov method that uses complex arithmetic.     

Bidiagonal decomposition of totally positive Bernstein-Vandermonde matrices

The class of Bernstein-Vandermonde matrices (a generalization of Vandermonde matrices arising when the monomial basis is replaced by the Bernstein basis) is considered. A convenient ordering of their rows makes these matrices strictly totally positive. By using results related to total positivity and Neville elimination, an algorithm for computing the bidiagonal decomposition of a Bernstein-Vandermonde matrix is constructed. The use of explicit expressions for the determinants involved in the process serves to make the algorithm both fast and accurate. One of the applications of our algorithm is the design of fast and accurate algorithms for solving Lagrange interpolation problems when using the Bernstein basis, an approach useful for the field of Computer Aided Geometric Design since it avoids the stability problems involved with basis transformations between the Bernstein and the monomial bases. A different application consists of the use of the bidiagonal decomposition as an intermediate step of the computation of the eigenvalues and the singular value decomposition of a totally positive Bernstein-Vandermonde matrix. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Non‐recursive solution of sparse block Hessenberg systems

A double‐phase algorithm, based on the block recursive LU decomposition, has been recently proposed to solve block Hessenberg systems with sparsity properties. In the first phase the information related to the factorization of A and required to solve the system, is computed and stored. The solution of the system is then computed in the second phase. In the present paper the algorithm is modified: the two phases are merged into a one‐phase version having the same computational cost and allowing a saving of storage. Moreover, the corresponding non‐recursive version of the new algorithm is presented, which is especially suitable to solve infinite systems where the coefficient matrix dimension is not a priori fixed and a subsequent size enlargement technique is used. A special version of the algorithm, well suited to deal with block Hessenberg matrices having also a block band structure, is presented. Theoretical asymptotic bounds on the computational costs are proved. They indicate that, under suitable sparsity conditions, the proposed algorithms outperform Gaussian elimination. Numerical experiments have been carried out, showing the effectiveness of the algorithms when the size of the system is of practical interest. Copyright © 2004 John Wiley & Sons, Ltd.     

Accurate computations with Said-Ball-Vandermonde matrices

A generalization of the Vandermonde matrices which arise when the power basis is replaced by the Said-Ball basis is considered. When the nodes are inside the interval (0,1), then those matrices are strictly totally positive. An algorithm for computing the bidiagonal decomposition of those Said-Ball-Vandermonde matrices is presented, which allows us to use known algorithms for totally positive matrices represented by their bidiagonal decomposition. The algorithm is shown to be fast and to guarantee high relative accuracy. Some numerical experiments which illustrate the good behaviour of the algorithm are included.     

An inverse free parallel spectral divide and conquer algorithm for nonsymmetric eigenproblems

Summary. We discuss an inverse-free, highly parallel, spectral divide and conquer algorithm. It can compute either an invariant subspace of a nonsymmetric matrix , or a pair of left and right deflating subspaces of a regular matrix pencil . This algorithm is based on earlier ones of Bulgakov, Godunov and Malyshev, but improves on them in several ways. This algorithm only uses easily parallelizable linear algebra building blocks: matrix multiplication and QR decomposition, but not matrix inversion. Similar parallel algorithms for the nonsymmetric eigenproblem use the matrix sign function, which requires matrix inversion and is faster but can be less stable than the new algorithm. Received September 20, 1994 / Revised version received February 5, 1996     

The design of a parallel dense linear algebra software library: Reduction to Hessenberg,tridiagonal, and bidiagonal form

This paper discusses issues in the design of ScaLAPACK, a software library for performing dense linear algebra computations on distributed memory concurrent computers. These issues are illustrated using the ScaLAPACK routines for reducing matrices to Hessenberg, tridiagonal, and bidiagonal forms. These routines are important in the solution of eigenproblems. The paper focuses on how building blocks are used to create higher-level library routines. Results are presented that demonstrate the scalability of the reduction routines. The most commonly-used building blocks used in ScaLAPACK are the sequencing BLAS, the parallel BLAS (PBLAS) and the Basic Linear Algebra Communication Subprograms (BLACS). Each of the matrix reduction algorithms consists of a series of steps in each of which one block column (orpanel), and/or block row, of the matrix is reduced, followed by an update of the portion of the matrix that has not been factorized so far. This latter phase is performed using Level 3 PBLAS operations and contains the bulk of the computation. However, the panel reduction phase involves a significant amount of communication, and is important in determining the scalability of the algorithm. The simplest way to parallelize the panel reduction phase is to replace the BLAS routines appearing in the LAPACK routine (mostly matrix-vector and matrix-matrix multiplications) with the corresponding PBLAS routines. However, in some cases it is possible to reduce communication startup costs by performing the communication necessary for consecutive BLAS operations in a single communication using a BLACS call. Thus, there is a tradeoff between efficiency and software engineering considerations, such as ease of programming and simplicity of code.Research was supported in part by the Applied Mathematical Sciences Research Program of the Office of Energy Research, U.S. Department of Energy, by the Defense Advanced Research Projects Agency under contract DAAL03-91-C-0047, administered by the Army Research Office, and in part by the Center for Research on Parallel Computing.     

酉延拓矩阵的奇异值分解及其广义逆

从普通奇异值分解出发,导出了酉延拓矩阵的奇异值和奇异向量与母矩阵的奇异值和奇异向量间的定量关系,同时对酉延拓矩阵的满秩分解及g逆,反射g逆,最小二乘g逆,最小范数g逆作了定量分析,得到了酉延拓矩阵的满秩分解矩阵F*和G*与母矩阵A的分解矩阵F和G之间的关系.最后给出了相应的快速求解算法,并举例说明该算法大大降低了分解的计算量和存储量,提高了计算效率.     

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.     

More on generalized inverses of partitioned matrices with Banachiewicz–Schur forms

Necessary and sufficient conditions are derived for a 2-by-2 partitioned matrix to have {1}-, {1,2}-, {1,3}-, {1,4}-inverses and the Moore–Penrose inverse with Banachiewicz–Schur forms. As applications, the Banachiewicz–Schur forms of {1}-, {1,2}-, {1,3}-, {1,4}-inverses and the Moore–Penrose inverse of a 2-by-2 partitioned Hermitian matrix are also given.     

A parallel balance scheme for banded linear systems

A parallel algorithm is proposed for the solution of narrow banded non‐symmetric linear systems. The linear system is partitioned into blocks of rows with a small number of unknowns common to multiple blocks. Our technique yields a reduced system defined only on these common unknowns which can then be solved by a direct or iterative method. A projection based extension to this approach is also proposed for computing the reduced system implicitly, which gives rise to an inner–outer iteration method. In addition, the product of a vector with the reduced system matrix can be computed efficiently on a multiprocessor by concurrent projections onto subspaces of block rows. Scalable implementations of the algorithm can be devized for hierarchical parallel architectures by exploiting the two‐level parallelism inherent in the method. Our experiments indicate that the proposed algorithm is a robust and competitive alternative to existing methods, particularly for difficult problems with strong indefinite symmetric part. Copyright © 2001 John Wiley & Sons, Ltd.     

Solving the eigenvalue problem for matrices

One presents some algorithms related among themselves for solving the partial and the complete eigenvalue problem for an arbitrary matrix. Algorithm 1 allows us to construct the invariant subspaces and to obtain with their aid a matrix whose eigenvalues coincide with the eigenvalues of the initial matrix and belong to a given semiplane. Algorithm 2 solves the same problem for a given strip. The algorithms 3 and 4 reduces the complete eigenvalue problem of an arbitrary matrix to some problem for a quasitriangular matrix whose diagonal blocks have eigenvalues with identical real parts. Algorithm 4 finds also the unitary matrix which realizes this transformation. One gives Algol programs which realize the algorithms 1–3 for real matrices and testing examples.     

High‐performance direct algorithms for computing the sign function of triangular matrices

Algorithms and implementations for computing the sign function of a triangular matrix are fundamental building blocks for computing the sign of arbitrary square real or complex matrices. We present novel recursive and cache‐efficient algorithms that are based on Higham's stabilized specialization of Parlett's substitution algorithm for computing the sign of a triangular matrix. We show that the new recursive algorithms are asymptotically optimal in terms of the number of cache misses that they generate. One algorithm that we present performs more arithmetic than the nonrecursive version, but this allows it to benefit from calling highly optimized matrix multiplication routines; the other performs the same number of operations as the nonrecursive version, suing custom computational kernels instead. We present implementations of both, as well as a cache‐efficient implementation of a block version of Parlett's algorithm. Our experiments demonstrate that the blocked and recursive versions are much faster than the previous algorithms and that the inertia strongly influences their relative performance, as predicted by our analysis.     

Fast and Stable Reduction of Diagonal Plus Semi-Separable Matrices to Tridiagonal and Bidiagonal Form

A two-way chasing algorithm to reduce a diagonal plus a symmetric semi-separable matrix to a symmetric tridiagonal one and an algorithm to reduce a diagonal plus an unsymmetric semi-separable matrix to a bidiagonal one are considered. Both algorithms are fast and stable, requiring a computational cost of N ², where N is the order of the considered matrix.     

Computing the generalized singular values/vectors of large sparse or structured matrix pairs

Summary. We present a numerical algorithm for computing a few extreme generalized singular values and corresponding vectors of a sparse or structured matrix pair . The algorithm is based on the CS decomposition and the Lanczos bidiagonalization process. At each iteration step of the Lanczos process, the solution to a linear least squares problem with as the coefficient matrix is approximately computed, and this consists the only interface of the algorithm with the matrix pair . Numerical results are also given to demonstrate the feasibility and efficiency of the algorithm. Received April 1, 1994 / Revised version received December 15, 1994     

A unitary joint diagonalization algorithm for nonsymmetric higher‐order tensors based on Givens‐like rotations

Based on Givens‐like rotations, we present a unitary joint diagonalization algorithm for a set of nonsymmetric higher‐order tensors. Each unitary rotation matrix only depends on one unknown parameter which can be analytically obtained in an independent way following a reasonable assumption and a complex derivative technique. It can serve for the canonical polyadic decomposition of a higher‐order tensor with orthogonal factors. Furthermore, based on cross‐high‐order cumulants of observed signals, we show that the proposed algorithm can be applied to solve the joint blind source separation problem. The simulation results reveal that the proposed algorithm has a competitive performance compared with those of several existing related methods.     

A quaternion QR algorithm

Summary This paper extends the Francis QR algorithm to quaternion and antiquaternion matrices. It calculates a quaternion version of the Schur decomposition using quaternion unitary similarity transformations. Following a finite step reduction to a Hessenberg-like condensed form, a sequence of implicit QR steps reduces the matrix to triangular form. Eigenvalues may be read off the diagonal. Eigenvectors may be obtained from simple back substitutions. For serial computation, the algorithm uses only half the work and storage of the unstructured Francis QR iteration. By preserving quaternion structure, the algorithm calculates the eigenvalues of a nearby quaternion matrix despite rounding errors.     

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