A method for computing the Perron root for primitive matrices

Following the Perron theorem, the spectral radius of a primitive matrix is a simple eigenvalue. It is shown that for a primitive matrix A, there is a positive rank one matrix X such that B = A ° X , where ° denotes the Hadamard product of matrices, and such that the row (column) sums of matrix B are the same and equal to the Perron root. An iterative algorithm is presented to obtain matrix B without an explicit knowledge of X. The convergence rate of this algorithm is similar to that of the power method but it uses less computational load. A byproduct of the proposed algorithm is a new method for calculating the first eigenvector.     

A New Algorithm for Simulating a Correlation Matrix Based on Parameter Expansion and Reparameterization

The correlation matrix (denoted by R) plays an important role in many statistical models. Unfortunately, sampling the correlation matrix in Markov chain Monte Carlo (MCMC) algorithms can be problematic. In addition to the positive definite constraint of covariance matrices, correlation matrices have diagonal elements fixed at one. In this article, we propose an efficient two-stage parameter expanded reparameterization and Metropolis-Hastings (PX-RPMH) algorithm for simulating R. Using this algorithm, we draw all elements of R simultaneously by first drawing a covariance matrix from an inverse Wishart distribution, and then translating it back to a correlation matrix through a reduction function and accepting it based on a Metropolis-Hastings acceptance probability. This algorithm is illustrated using multivariate probit (MVP) models and multivariate regression (MVR) models with a common correlation matrix across groups. Via both a simulation study and a real data example, the performance of the PX-RPMH algorithm is compared with those of other common algorithms. The results show that the PX-RPMH algorithm is more efficient than other methods for sampling a correlation matrix.     

A note on the solution of not balanced banded Toeplitz systems

A direct algorithm is presented for the solution of linear systems having banded Toeplitz coefficient matrix with unbalanced bandwidths. It is derived from the cyclic reduction algorithm, it makes use of techniques based on the displacement rank and it relies on the Morrison–Sherman–Woodbury formula. The algorithm always equals and sometimes outperforms the already known direct ones in terms of asymptotic computational cost. The case where the coefficient matrix is a block banded block Toeplitz matrix in block Hessenberg form is analyzed as well. The algorithm is numerically stable if applied to M‐matrices that are point diagonally dominant by columns. Copyright © 2007 John Wiley & Sons, Ltd.     

A binary powering Schur algorithm for computing primary matrix roots

An algorithm for computing primary roots of a nonsingular matrix A is presented. In particular, it computes the principal root of a real matrix having no nonpositive real eigenvalues, using real arithmetic. The algorithm is based on the Schur decomposition of A and has an order of complexity lower than the customary Schur based algorithm, namely the Smith algorithm.     

Rounding Error Analysis in Solving M-Matrix Linear Systems of Block Hessenberg Form

The problem of solving large M-matrix linear systems with sparse coefficient matrix in block Hessenberg form is here addressed. In previous work of the authors a divide-and-conquer strategy was proposed and a backward error analysis of the resulting algorithm was presented showing its effectiveness for the solution of computational problems of queueing theory and Markov chains. In particular, it was shown that for block Hessenberg M-matrices the algorithm is weakly backward stable in the sense that the computed solution is the exact solution of a nearby linear system, where the norm of the perturbation is proportional to the condition number of the coefficient matrix. In this note a better error estimate is given by showing that for block Hessenberg M-matrices the algorithm is even backward stable.     

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.     

Computing the inertia from sign patterns

A symmetric matrix A is said to be sign-nonsingular if every symmetric matrix with the same sign pattern as A is nonsingular. Hall, Li and Wang showed that the inertia of a sign-nonsingular symmetric matrix is determined uniquely by its sign pattern. The purpose of this paper is to present an efficient algorithm for computing the inertia of such symmetric matrices. The algorithm runs in time for a symmetric matrix of order n with m nonzero entries. In addition, it is shown to be NP-complete to decide whether the inertia of a given symmetric matrix is not determined by its sign pattern.     

Efficient eigenvalue computation for quasiseparable Hermitian matrices under low rank perturbations

In this paper we address the problem of efficiently computing all the eigenvalues of a large N×N Hermitian matrix modified by a possibly non Hermitian perturbation of low rank. Previously proposed fast adaptations of the QR algorithm are considerably simplified by performing a preliminary transformation of the matrix by similarity into an upper Hessenberg form. The transformed matrix can be specified by a small set of parameters which are easily updated during the QR process. The resulting structured QR iteration can be carried out in linear time using linear memory storage. Moreover, it is proved to be backward stable. Numerical experiments show that the novel algorithm outperforms available implementations of the Hessenberg QR algorithm already for small values of N.      

A pivoting strategy for symmetric tridiagonal matrices

The LBL^T factorization of Bunch for solving linear systems involving a symmetric indefinite tridiagonal matrix T is a stable, efficient method. It computes a unit lower triangular matrix L and a block 1 × 1 and 2 × 2 matrix B such that T=LBL^T. Choosing the pivot size requires knowing a priori the largest element σ of T in magnitude. In some applications, it is required to factor T as it is formed without necessarily knowing σ. In this paper, we present a modification of the Bunch algorithm that can satisfy this requirement. We demonstrate that this modification exhibits the same bound on the growth factor as the Bunch algorithm and is likewise normwise backward stable. Copyright © 2005 John Wiley & Sons, Ltd.     

A fast minimal residual algorithm for shifted unitary matrices

A new iterative scheme is described for the solution of large linear systems of equations with a matrix of the form A = ρU + ζI, where ρ and ζ are constants, U is a unitary matrix and I is the identity matrix. We show that for such matrices a Krylov subspace basis can be generated by recursion formulas with few terms. This leads to a minimal residual algorithm that requires little storage and makes it possible to determine each iterate with fairly little arithmetic work. This algorithm provides a model for iterative methods for non-Hermitian linear systems of equations, in a similar way to the conjugate gradient and conjugate residual algorithms. Our iterative scheme illustrates that results by Faber and Manteuffel [3,4] on the existence of conjugate gradient algorithms with short recurrence relations, and related results by Joubert and Young [13], can be extended.     

A Quantum Version of the Désarménien Matrix

We use elements in the quantum hyperalgebra to define a quantum version of the Désarménien matrix. We prove that our matrix is upper triangular with ones on the diagonal and that, as in the classical case, it gives a quantum straightening algorithm for quantum bideterminants. We use our matrix to give a new proof of the standard basis theorem for the q-Weyl module. As well, we show that the standard basis for the q-Weyl module and the basis dual to the standard basis for the q-Schur module are related by the quantum Désarménien matrix.     

Weighted and deflated global GMRES algorithms for solving large Sylvester matrix equations

The solution of a large-scale Sylvester matrix equation plays an important role in control and large scientific computations. In this paper, we are interested in the large Sylvester matrix equation with large dimensionA and small dimension B, and a popular approach is to use the global Krylov subspace method. In this paper, we propose three new algorithms for this problem. We first consider the global GMRES algorithm with weighting strategy, which can be viewed as a precondition method. We present three new schemes to update the weighting matrix during iterations. Due to the growth of memory requirements and computational cost, it is necessary to restart the algorithm effectively. The deflation strategy is efficient for the solution of large linear systems and large eigenvalue problems; to the best of our knowledge, little work is done on applying deflation to the (weighted) global GMRES algorithm for large Sylvester matrix equations. We then consider how to combine the weighting strategy with deflated restarting, and propose a weighted global GMRES algorithm with deflation for solving large Sylvester matrix equations. In particular, we are interested in the global GMRES algorithm with deflation, which can be viewed as a special case when the weighted matrix is chosen as the identity. Theoretical analysis is given to show rationality of the new algorithms. Numerical experiments illustrate the numerical behavior of the proposed algorithms.

     

An integer arithmetic method to compute generalized matrix inverse and solve linear equations exactly

An algorithm that uses integer arithmetic is suggested. It transforms anm ×n matrix to a diagonal form (of the structure of Smith Normal Form). Then it computes a reflexive generalized inverse of the matrix exactly and hence solves a system of linear equations error-free.     

An alternating direction algorithm for matrix completion with nonnegative factors

This paper introduces an algorithm for the nonnegative matrix factorization-and-completion problem, which aims to find nonnegative low-rank matrices X and Y so that the product XY approximates a nonnegative data matrix M whose elements are partially known (to a certain accuracy). This problem aggregates two existing problems: (i) nonnegative matrix factorization where all entries of M are given, and (ii) low-rank matrix completion where nonnegativity is not required. By taking the advantages of both nonnegativity and low-rankness, one can generally obtain superior results than those of just using one of the two properties. We propose to solve the non-convex constrained least-squares problem using an algorithm based on the classical alternating direction augmented Lagrangian method. Preliminary convergence properties of the algorithm and numerical simulation results are presented. Compared to a recent algorithm for nonnegative matrix factorization, the proposed algorithm produces factorizations of similar quality using only about half of the matrix entries. On tasks of recovering incomplete grayscale and hyperspectral images, the proposed algorithm yields overall better qualities than those produced by two recent matrix-completion algorithms that do not exploit nonnegativity.     

An extended fast algorithm for constructing the Dixon resultant matrix

In recent years,the Dixon resultant matrix has been used widely in the re-sultant elimination to solve nonlinear polynomial equations and many researchers havestudied its efficient algorithms.The recursive algorithm is a very efficient algorithm,butwhich deals with the case of three polynomial equations with two variables at most.Inthis paper,we extend the algorithm to the general case of n 1 polynomial equations in nvariables.The algorithm has been implemented in Maple 9.By testing the random polyno-mial equations,the results demonstrate that the efficiency of our program is much betterthan the previous methods,and it is exciting that the necessary condition for the existenceof common intersection points on four general surfaces in which the degree with respectto every variable is not greater than 2 is given out in 48×48 Dixon matrix firstly by ourprogram.     

An algorithm for determining if the inverse of a strictly diagonally dominant Stieltjes matrix is strictly ultrametric

Summary It was recently shown that the inverse of a strictly ultrametric matrix is a strictly diagonally dominant Stieltjes matrix. On the other hand, as it is well-known that the inverse of a strictly diagonally dominant Stieltjes matrix is a real symmetric matrix with nonnegative entries, it is natural to ask, conversely, if every strictly diagonally dominant Stieltjes matrix has a strictly ultrametric inverse. Examples show, however, that the converse is not true in general, i.e., there are strictly diagonally dominant Stieltjes matrices in ^n×n (for everyn3) whose inverses are not strictly ultrametric matrices. Then, the question naturally arises if one can determine which strictly diagonally dominant Stieltjes matrices, in ^n×n (n3), have inverses which are strictly ultrametric. Here, we develop an algorithm, based on graph theory, which determines if a given strictly diagonally dominant Stieltjes matrixA has a strictly ultrametric inverse, where the algorithm is applied toA and requires no computation of inverse. Moreover, if this given strictly diagonally dominant Stieltjes matrix has a strictly ultrametric inverse, our algorithm uniquely determines this inverse as a special sum of rank-one matrices.Research supported by the National Science FoundationResearch supported by the Deutsche Forschungsgemeinschaft     

QR factorization of Toeplitz matrices

Summary This paper presents a new algorithm for computing theQR factorization of anm×n Toeplitz matrix inO(mn) operations. The algorithm exploits the procedure for the rank-1 modification and the fact that both principal (m–1)×(n–1) submatrices of the Toeplitz matrix are identical. An efficient parallel implementation of the algorithm is possible.     

A Finite Smoothing Algorithm for Quantile Regression

This article introduces a new method for computing regression quantile functions. This method applies a finite smoothing algorithm based on smoothing the nondifferentiable quantile regression objective function ρ_τ. The smoothing can be done for all τ ∈ (0, 1), and the convergence is finite for any finite number of τ_i ∈ (0, 1), i = 1,…,N. Numerical comparison shows that the finite smoothing algorithm outperforms the simplex algorithm in computing speed. Compared with the powerful interior point algorithm, which was introduced in an earlier article, it is competitive overall; however, it is significantly faster than the interior point algorithm when the design matrix in quantile regression has a large number of covariates. Additionally, the new algorithm provides the same accuracy as the simplex algorithm. In contrast, the interior point algorithm gives only the approximate solutions in theory, and rounding may be necessary to improve the accuracy of these solutions in practice.     

Computing the Matrix Cosine

An algorithm is developed for computing the matrix cosine, building on a proposal of Serbin and Blalock. The algorithm scales the matrix by a power of 2 to make the -norm less than or equal to 1, evaluates a Padé approximant, and then uses the double angle formula cos(2A)=2cos(A)²–I to recover the cosine of the original matrix. In addition, argument reduction and balancing is used initially to decrease the norm. We give truncation and rounding error analyses to show that an [8,8] Padé approximant produces the cosine of the scaled matrix correct to machine accuracy in IEEE double precision arithmetic, and we show that this Padé approximant can be more efficiently evaluated than a corresponding Taylor series approximation. We also provide error analysis to bound the propagation of errors in the double angle recurrence. Numerical experiments show that our algorithm is competitive in accuracy with the Schur–Parlett method of Davies and Higham, which is designed for general matrix functions, and it is substantially less expensive than that method for matrices of -norm of order 1. The dominant computational kernels in the algorithm are matrix multiplication and solution of a linear system with multiple right-hand sides, so the algorithm is well suited to modern computer architectures.     

A new algorithm for computing the Geronimus transformation with large shifts

A monic Jacobi matrix is a tridiagonal matrix which contains the parameters of the three-term recurrence relation satisfied by the sequence of monic polynomials orthogonal with respect to a measure. The basic Geronimus transformation with shift α transforms the monic Jacobi matrix associated with a measure dμ into the monic Jacobi matrix associated with dμ/(x − α) + Cδ(x − α), for some constant C. In this paper we examine the algorithms available to compute this transformation and we propose a more accurate algorithm, estimate its forward errors, and prove that it is forward stable. In particular, we show that for C = 0 the problem is very ill-conditioned, and we present a new algorithm that uses extended precision.