A NUMERICALLY STABLE BLOCK MODIFIED GRAM-SCHMIDT ALGORITHM FOR SOLVING STIFF WEIGHTED LEAST SQUARES PROBLEMS

Recently, Wei in proved that perturbed stiff weighted pseudoinverses and stiff weighted least squares problems are stable, if and only if the original and perturbed coefficient matrices A and A^- satisfy several row rank preservation conditions. According to these conditions, in this paper we show that in general, ordinary modified Gram-Schmidt with column pivoting is not numerically stable for solving the stiff weighted least squares problem. We then propose a row block modified Gram-Schmidt algorithm with column pivoting, and show that with appropriately chosen tolerance, this algorithm can correctly determine the numerical ranks of these row partitioned sub-matrices, and the computed QR factor R^- contains small roundoff error which is row stable. Several numerical experiments are also provided to compare the results of the ordinary Modified Gram-Schmidt algorithm with column pivoting and the row block Modified Gram-Schmidt algorithm with column pivoting.     

Compressions of Stable Contractions

The stability of compressions of stable contractions is studied and a sufficient orbit condition is given. On the other hand, it is shown that there are non-stable compressions of the 1-dimensional backward shift and a complete characterization of weighted unilateral shifts with this property is provided. Dilations of bilateral weighted shifts to backward shifts are also considered.     

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.     

Algorithm engineering for optimal alignment of protein structure distance matrices

Protein structural alignment is an important problem in computational biology. In this paper, we present first successes on provably optimal pairwise alignment of protein inter-residue distance matrices, using the popular dali scoring function. We introduce the structural alignment problem formally, which enables us to express a variety of scoring functions used in previous work as special cases in a unified framework. Further, we propose the first mathematical model for computing optimal structural alignments based on dense inter-residue distance matrices. We therefore reformulate the problem as a special graph problem and give a tight integer linear programming model. We then present algorithm engineering techniques to handle the huge integer linear programs of real-life distance matrix alignment problems. Applying these techniques, we can compute provably optimal dali alignments for the very first time.     

A backward stability analysis of diagonal pivoting methods for solving unsymmetric tridiagonal systems without interchanges

This paper concerns the LBM ^T factorization of unsymmetric tridiagonal matrices, where L and M are unit lower triangular matrices and B is block diagonal with 1×1 and 2×2 blocks. In some applications, it is necessary to form this factorization without row or column interchanges while the tridiagonal matrix is formed. Bunch and Kaufman proposed a pivoting strategy without interchanges specifically for symmetric tridiagonal matrices, and more recently, Bunch and Marcia proposed pivoting strategies that are normwise backward stable for linear systems involving such matrices. In this paper, we extend these strategies to the unsymmetric tridiagonal case and demonstrate that the proposed methods both exhibit bounded growth factors and are normwise backward stable. Copyright © 2009 John Wiley & Sons, Ltd.     

Stability and Sensitivity of Tridiagonal LU Factorization without Pivoting

In this paper the accuracy of LU factorization of tridiagonal matrices without pivoting is considered. Two types of componentwise condition numbers for the L and U factors of tridiadonal matrices are presented and compared. One type is a condition number with respect to small relative perturbations of each entry of the matrix. The other type is a condition number with respect to small componentwise perturbations of the kind appearing in the backward error analysis of the usual algorithm for the LU factorization. We show that both condition numbers are of similar magnitude. This means that the algorithm is componentwise forward stable, i.e., the forward errors are of similar magnitude to those produced by a componentwise backward stable method. Moreover the presented condition numbers can be computed in O(n) flops, which allows to estimate with low cost the forward errors. AMS subject classification (2000) 65F35, 65F50, 15A12, 15A23, 65G50.Received October 2003. Accepted August 2004. Communicated by Per Christian Hansen.Froilán M. Dopico: This research has been partially supported by the Ministerio de Ciencia y Tecnología of Spain through grants BFM2003-06335-C03-02 (M. I. Bueno) and BFM2000-0008 (F. M. Dopico).     

Backward error of polynomial eigenvalue problems solved by linearization

It is commonplace in many application domains to utilize polynomial eigenvalue problems to model the behaviour of physical systems. Many techniques exist to compute solutions of these polynomial eigenvalue problems. One of the most frequently used techniques is linearization, in which the polynomial eigenvalue problem is turned into an equivalent linear eigenvalue problem with the same eigenvalues, and with easily recoverable eigenvectors. The eigenvalues and eigenvectors of the linearization are usually computed using a backward stable solver such as the QZ algorithm. Such backward stable algorithms ensure that the computed eigenvalues and eigenvectors of the linearization are exactly those of a nearby linear pencil, where the perturbations are bounded in terms of the machine precision and the norms of the matrices defining the linearization. Although we have solved a nearby linear eigenvalue problem, we are not certain that our computed solution is in fact the exact solution of a nearby polynomial eigenvalue problem. Here, we perform a backward error analysis for the solution of a specific linearization for polynomials expressed in the monomial basis. We use a suitable one-sided factorization of the linearization that allows us to map generic perturbations of the linearization onto structured perturbations of the polynomial coefficients. (© 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

On solution of large systems of linear equations with block-Toeplitz banded matrices

A novel numerical algorithm is proposed for solving systems of linear equations with block-Toeplitz narrow-banded matrices. Its arithmetical cost is about double of that of the block cyclic reduction. The backward roundoff error stability of the method guarantees its reliability for the matrices that are not symmetric positive definite or diagonally dominant.     

ON STABLE PERTURBATIONS OF THE STIFFLY WEIGHTED PSEUDOINVERSE AND WEIGHTED LEAST SQUARES PROBLEM

In this paper we study perturbations of the stiffly weighted pseudoinverse （W^1/2 A）^＋W^1/2 and the related stiffly weighted least squares problem, where both the matrices A and W are given with W positive diagonal and severely stiff. We show that the perturbations to the stiffly weighted pseudoinverse and the related stiffly weighted least squares problem are stable, if and only if the perturbed matrices A = A ＋ δA satisfy several row rank preserving conditions.     

STABLE DOUBLE LR ALGORITHM AND ITS ERROR ANALYSIS

In this paper, the normative matrices and their double LR transformation with origin shifts are defined, and the essential relationship between the double LR transformation of a normative matrix and the QR transformation of the related symmetric tridiagonal matrix is proved. We obtain a stable double LR algorithm for double LR transformation of normative matrices and give the error analysis of our algorithm. The operation number of the stable double LR algorithm for normative matrices is only four sevenths of the rational QR algorithm for reed symmetric tridiagonal matrices.     

Rounding Errors in Solving Block Hessenberg Systems

A rounding error analysis is presented for a divide-and-conquer algorithm to solve linear systems with block Hessenberg matrices. Conditions are derived under which the algorithm computes a stable solution. The algorithm is shown to be stable for block diagonally dominant matrices and for M-matrices.

     

Factoring symmetric totally nonpositive matrices and inverses with a diagonal pivoting method

In this paper, we consider how to factor symmetric totally nonpositive matrices and their inverses by taking advantage of the symmetric property. It is well-known that the Bunch-Kaufman algorithm is the most commonly used pivoting strategy which can, however, produce arbitrarily large entries in the lower triangular factor for such matrices as illustrated by our example. Therefore, it is interesting to show that when the Bunch-Parlett algorithm is simplified for these matrices, it only requires O(n ²) comparisons with the growth factor being nicely bounded by 4. These facts, together with a nicely bounded lower triangular factor and a pleasantly small relative backward error, show that the Bunch-Parlett algorithm is more preferable than the Bunch-Kaufman algorithm when dealing with these matrices.     

Block‐row Hankel weighted low rank approximation

This paper extends the weighted low rank approximation (WLRA) approach to linearly structured matrices. In the case of Hankel matrices with a special block structure, an equivalent unconstrained optimization problem is derived and an algorithm for solving it is proposed. Copyright © 2005 John Wiley & Sons, Ltd.     

Finding stationary points on bounded-rank matrices: a geometric hurdle and a smooth remedy

Mathematical Programming - We consider the problem of provably finding a stationary point of a smooth function to be minimized on the variety of bounded-rank matrices. This turns out to be...     

Existence and uniqueness of weighted pseudoinverse matrices and weighted normal pseudosolutions with singular weights

For one definition of weighted pseudoinversion with singular weights, we establish necessary and sufficient conditions for the existence and uniqueness of a solution of a system of matrix equations. Expansions of weighted pseudoinverse matrices in matrix power series and matrix power products are obtained. A relationship between weighted pseudoinverse matrices the weighted normal pseudosolutions is established. Iterative methods for the calculation of weighted pseudoinverse matrices and weighted normal pseudosolutions are constructed.     

Finite element analysis of time-dependent viscous flows utilizing an implicit mixed interpolation algorithm with a frontal solution technique

A finite element algorithm is described which implements the Galerkin approximation to the Navier-Stokes equations and incorporates five predominate features. Although none of these five features is unique to this algorithm, their orchestration, as described in this paper, results in an algorithm which is not only easy to implement but also stable, accurate, and robust, as well as computationally efficient. The zero stress natural boundary condition is implemented which permits calculation of the outflow velocity distribution. A nine-node, Lagrangian, isoparametric, quadrilateral element is used to represent the velocity while the pressure uses a four-node, Lagrangian, superparametric element coincident with the velocity element. The easily implemented, computationally efficient frontal solution technique is used to assemble the element coefficient matrices, impose the boundary conditions, and solve the resulting linear system of equations. An implicit backward Euler time integration rule provides a very stable solution method for time-dependent problems. A Picard scheme with a relatively large radius of convergence is used for iteration of the nonlinear equations at each time step. Results are given for the calculation of two-dimensional, steady- state and time-dependent convection dominated flows of viscous incompressible fluids.     

A finite algorithm for ϵ-optimal solutions of adaptive queueing control

In this paper, by introducing a weighted supremum norm, then using non-adaptive approach and backward induction, we suggest an ε-optimal algorithm for the problem of adaptive queueing control. In case of Bernoulli queues, we obtain an ε-optimal cost and a 2ε-optimal policy. Furthermore, a simple inequality is given as an upper bound of the error, which can be used for determining the number of backward induction steps.     

A new parallel chasing algorithm for transforming arrowhead matrices to tridiagonal form

Rutishauser, Gragg and Harrod and finally H.Y. Zha used the same class of chasing algorithms for transforming arrowhead matrices to tridiagonal form. Using a graphical theoretical approach, we propose a new chasing algorithm. Although this algorithm has the same sequential computational complexity and backward error properties as the old algorithms, it is better suited for a pipelined approach. The parallel algorithm for this new chasing method is described, with performance results on the Paragon and nCUBE. Comparison results between the old and the new algorithms are also presented.

     

The application of an implicit finite element algorithm with a frontal solution technique to the analysis of time dependent viscous flows

A finite element algorithm is described which implements the Galerkin approximation to the Navier-Stokes equations and incorporates five predominate features. Although none of these five features is unique to this algorithm, their orchestration, as described in this paper, results in an algorithm which is not only easy to implement but also stable, accurate, and robust, as well as computationally efficient. The zero stress natural boundary condition is implemented which permits calculation of the outflow velocity distribution. A nine-node, Lagrangian, isoparametric, quadrilateral elementis used to represent the velocity while the pressure uses a four-node, Lagrangian, superparametric element coincident with the velocity element. The easily implemented, computationally efficient frontal solution technique is used to assemble the element coefficient matrices, impose the boundary conditions, and solve the resulting linear system of equations. An implicit backward Euler time integration rule provides a very stable solution method for time dependent problems. A Picard scheme with a relatively large radius of convergence is used for iteration of the non-linear equations at each time step. Results are given from the calculation of two dimensional, steady-state and time-dependent convection dominated flows of viscous incompressible fluids.     

On the modified Gram-Schmidt algorithm for weighted and constrained linear least squares problems

A framework and an algorithm for using modified Gram-Schmidt for constrained and weighted linear least squares problems is presented. It is shown that a direct implementation of a weighted modified Gram-Schmidt algorithm is unstable for heavily weighted problems. It is shown that, in most cases it is possible to get a stable algorithm by a simple modification free from any extra computational costs. In particular, it is not necessary to perform reorthogonalization.Solving the weighted and constrained linear least squares problem with the presented weighted modified Gram-Schmidt algorithm is seen to be numerically equivalent to an algorithm based on a weighted Householder-likeQR factorization applied to a slightly larger problem. This equivalence is used to explain the instability of the weighted modified Gram-Schmidt algorithm. If orthogonality, with respect to a weighted inner product, of the columns inQ is important then reorthogonalization can be used. One way of performing such reorthogonalization is described.Computational tests are given to show the main features of the algorithm.