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     

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.     

计算部分奇异值分解的隐式重新启动的双对角化Lanczos方法和精化的双对角化Lanczos方法

The singular value decomposition problem is mathematically equivalent to the eigenproblem of an argumented matrix. Golub et al. give a bidiagonalization Lanczos method for computing a number of largest or smallest singular values and corresponding singular vertors, but the method may encounter some convergence problems. In this paper we analyse the convergence of the method and show why it may fail to converge. To correct this possible nonconvergence, we propose a refined bidiagonalization Lanczos method and apply the implicitly restarting technique to it, and we then present an implicitly restarted bidiagonalization Lanczos algorithm(IRBL) and an implicitly restarted refined bidiagonalization Lanczos algorithm (IRRBL). A new implicitly restarting scheme and a reliable and efficient algorithm for computing refined shifts are developed for this special structure eigenproblem.Theoretical analysis and numerical experiments show that IRRBL performs much better than IRBL.     

On Kogbetliantz's SVD Algorithm in the Presence of Clusters

We consider matrices with off-diagonal blocks of small norm and derive tight bounds for the approximation of their singular values by those of their diagonal blocks. These results are used to show that triangular matrices with clusters of singular values must possess a principal submatrix of “nearly” diagonal form. From the latter we then derive results pertaining to the quadratic convergence of Kogbetliantz's algorithm for computing the SVD, in the presence of clusters.     

Low-rank revealing QR factorizations

Rank revealing factorizations are used extensively in signal processing in connection with, for example, linear prediction and signal subspace algorithms. We present an algorithm for computing rank revealing QR factorizations of low-rank matrices. The algorithm produces tight upper and lower bounds for all the largest singular values, thus making it particularly useful for treating rank deficient problems by means of subset selection, truncated QR, etc. The algorithm is similar in spirit to an algorithm suggested earlier by Chan for matrices with a small nullity, and it can also be considered as an extension of ordinary column pivoting.     

A comparison between the complex symmetric based and classical computation of the singular value decomposition of normal matrices

An algorithm for computing the singular value decomposition of normal matrices using intermediate complex symmetric matrices is proposed. This algorithm, as most eigenvalue and singular value algorithms, consists of two steps. It is based on combining the unitarily equivalence of normal matrices to complex symmetric tridiagonal form with the symmetric singular value decomposition of complex symmetric matrices. Numerical experiments are included comparing several algorithms, with respect to speed and accuracy, for computing the symmetric singular value decomposition (also known as the Takagi factorization). Next we compare the novel approach with the classical Golub-Kahan method for computing the singular value decomposition of normal matrices: it is faster, consumes less memory, but on the other hand the results are significantly less accurate.     

An efficient and reliable algorithm for computing the singular subspace of a matrix,associated with its smallest singular values

In this paper, an improved algorithm PSVD for computing the singular subspace of a matrix corresponding to its smallest singular values is presented. As only a basis of the desired singular subspace is needed, the classical Singular Value Decomposition (SVD) algorithm can be modified in three ways. First, the Householder transformations of the bidiagonalization need only to be applied on the base vectors of the desired singular subspace. Second, the bidiagonal must only be partially diagonalized and third, the convergence rate of the iterative diagonalization can be improved by an appropriate choice between QR and QL iteration steps. An analysis of the operation counts, as well as computational results, show the relative efficiency of PSVD with respect to the classical SVD algorithm. Depending on the gap, the desired numerical accuracy and the dimension of the desired subspace, PSVD can be three times faster than the classical SVD algorithm while the same accuracy can be maintained. The new algorithm can be successfully used in total least squares applications, in the computation of the null space of a matrix and in solving (non) homogeneous linear equations. Based on PSVD a very efficient and reliable algorithm is also derived for solving nonhomogeneous equations.     

Verified Computation of Eigenpairs in the Generalized Eigenvalue Problem for Nonsquare Matrix Pencils

Consider an optimization problem arising from the generalized eigenvalue problem Ax=λBx,where A,B∈Cm×n and m>n.Ito et al.showed that the optimization problem can be solved by utilizing right singular vectors of C:=[B,A].In this paper,we focus on computing intervals containing the solution.When some singular values of C are multiple or nearly multiple,we can enclose bases of corresponding invariant subspaces of CHC,where CH denotes the conjugate transpose of C,but cannot enclose the corresponding right singular vectors.The purpose of this paper is to prove that the solution can be obtained even when we utilize the bases instead of the right singular vectors.Based on the proved result,we propose an algorithm for computing the intervals.Numerical results show property of the algorithm.     

On the spectral problem for trivariate functions

Using a variational approach applied to generalized Rayleigh functionals, we extend the concepts of singular values and singular functions to trivariate functions defined on a rectangular parallelepiped. We also consider eigenvalues and eigenfunctions for trivariate functions whose domain is a cube. For a general finite-rank trivariate function, we describe an algorithm for computing the canonical polyadic (CP) decomposition, provided that the CP factors are linearly independent in two variables. All these notions are computed using Chebfun3; a part of Chebfun for numerical computing with 3D functions. Application in finding the best rank-1 approximation of trivariate functions is investigated. We also prove that if the function is analytic and two-way orthogonally decomposable (odeco), then the CP values decay geometrically, and optimal finite-rank approximants converge at the same rate.     

Hopf bifurcation for a class of three-dimensional nonlinear dynamic systems

In this paper, Hopf bifurcation for a class of three-dimensional nonlinear dynamic systems is studied, a new algorithm of the formal series for the flow on center manifold is discussed, from this, a recursion formula for computation of the singular point quantities is obtained for the corresponding bifurcation equation, which is linear and then avoids complex integrating operations, therefore the calculation can be readily done with using computer symbol operation system such as Mathematica, and more the algebraic equivalence of the singular point quantities and corresponding focal values is proved, thus Hopf bifurcation can be considered easily. Finally an example is studied, by computing the singular point quantities and constructing a bifurcation function, the existence of 5 limit cycles bifurcated from the origin for the flow on center manifold is proved.     

A sort-Jacobi algorithm for semisimple lie algebras

A structure preserving sort-Jacobi algorithm for computing eigenvalues or singular values is presented. It applies to an arbitrary semisimple Lie algebra on its (−1)-eigenspace of the Cartan involution. Local quadratic convergence for arbitrary cyclic schemes is shown for the regular case. The proposed method is independent of the representation of the underlying Lie algebra and generalizes well-known normal form problems such as e.g. the symmetric, Hermitian, skew-symmetric, symmetric and skew-symmetric R-Hamiltonian eigenvalue problem and the singular value decomposition.     

A block algorithm for computing rank-revealing QR factorizations

We present a block algorithm for computing rank-revealing QR factorizations (RRQR factorizations) of rank-deficient matrices. The algorithm is a block generalization of the RRQR-algorithm of Foster and Chan. While the unblocked algorithm reveals the rank by peeling off small singular values one by one, our algorithm identifies groups of small singular values. In our block algorithm, we use incremental condition estimation to compute approximations to the nullvectors of the matrix. By applying another (in essence also rank-revealing) orthogonal factorization to the nullspace matrix thus created, we can then generate triangular blocks with small norm in the lower right part ofR. This scheme is applied in an iterative fashion until the rank has been revealed in the (updated) QR factorization. We show that the algorithm produces the correct solution, under very weak assumptions for the orthogonal factorization used for the nullspace matrix. We then discuss issues concerning an efficient implementation of the algorithm and present some numerical experiments. Our experiments show that the block algorithm is reliable and successfully captures several small singular values, in particular in the initial block steps. Our experiments confirm the reliability of our algorithm and show that the block algorithm greatly reduces the number of triangular solves and increases the computational granularity of the RRQR computation.This work was supported by the Applied Mathematical Sciences subprogram of the Office of Energy Research, US Department of Energy, under Contract W-31-109-Eng-38. The second author was also sponsored by a travel grant from the Knud Højgaards Fond.This work was partially completed while the author was visiting the IBM Scientific Center in Heidelberg, Germany.     

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 <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.     

An implicitly restarted block Lanczos bidiagonalization method using Leja shifts

In this paper, we propose an implicitly restarted block Lanczos bidiagonalization (IRBLB) method for computing a few extreme or interior singular values and associated right and left singular vectors of a large matrix A. Our method combines the advantages of a block routine, implicit shifting, and the application of Leja points as shifts in the accelerating polynomial. The method neither requires factorization of A nor of a matrix that contains A. This makes the method well suited for very large matrices. For many matrices, the method presented in this paper requires less computational work and computer storage than other available approaches. Computed examples illustrate the performance of the method described.     

An Optimization algorithm computing the maximal singular number of a real matrix

The problem of calculating the maximal singular number of a given real matrix is considered. The existing solution methods are briefly surveyed. A new optimization-type algorithm for computing the maximal singular number is suggested and substantiated. Its rate of convergence is proved to be linear. A relationship between the row sums of the matrix and one of its singular numbers is established, and new localization theorems are proved. It is shown how the suggested algorithm is related to the Relay relation relaxation method. Exceptional situations in which the algorithm converges to a non-maximal singular number are described. A computational trick for avoiding such situations with fairly high reliability is suggested.     

The partial total least squares algorithm

In this paper, an improved algorithm PTLS for solving total least squares (TLS) problems AX ≈ B is presented. As only a basis of the right singular subspace associated with the smallest singular values of the data [A; B] is needed, the computational cost can be reduced considerably by using the partial SVD algorithm. This algorithm computes in an efficient way a basis for the left and/or right singular subspace of a matrix associated with its smallest singular values.An analysis of the operation counts, as well as computational results, show the relative efficiency of PTLS with respect to the classical TLS algorithm. Typically, PTLS reduces the computation time with a factor 2.     

用奇异值分解方法计算具有重特征值矩阵的特征矢量

若当(Jordan)形是矩阵在相似条件下的一个标准形,在代数理论及其工程应用中都具有十分重要的意义．针对具有重特征值的矩阵,提出了一种运用奇异值分解方法计算它的特征矢量及若当形的算法．大量数值例子的计算结果表明,该算法在求解具有重特征值的矩阵的特征矢量及若当形上效果良好,优于商用软件MATLAB和MATHEMATICA．     

On the equivalence of singular point quantities and the integrability of a fine critical singular point

In this paper, the equivalence of singular point quantities and the integrability of a fine critical singular point are discussed. For an extensive class of complex autonomous differential systems, the theoretical basis for the methods of using formal series and integrating factors to calculate the singular point quantities for the origin (a fine critical singular point) and judge its integrability is obtained. Two recursion formulas, which are much simpler than focal values for computation because the operation is rational, for computing singular point quantities are introduced.     

An algorithm for computing wavefront amplitudes and inverse problems (tsunami,electrodynamics, acoustics,and viscoelasticity)

An algorithm for computing the amplitude of the leading wavefront generated by an impulse source of oscillations is proposed. According to the algorithm, the fundamental solution is represented in the form of the sum of singular and regular components. As a result, the time required for the amplitude computation is reduced by one order of magnitude. Examples of wavefront amplitudes of tsunami, electromagnetic, acoustic, and viscoelastic waves are given.