Showing papers in "SIAM Journal on Matrix Analysis and Applications in 2001"

A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling

Abstract: In this paper, we analyze the main features and discuss the tuning of the algorithms for the direct solution of sparse linear systems on distributed memory computers developed in the context of a long term European research project. The algorithms use a multifrontal approach and are especially designed to cover a large class of problems. The problems can be symmetric positive definite, general symmetric, or unsymmetric matrices, both possibly rank deficient, and they can be provided by the user in several formats. The algorithms achieve high performance by exploiting parallelism coming from the sparsity in the problem and that available for dense matrices. The algorithms use a dynamic distributed task scheduling technique to accommodate numerical pivoting and to allow the migration of computational tasks to lightly loaded processors. Large computational tasks are divided into subtasks to enhance parallelism. Asynchronous communication is used throughout the solution process to efficiently overlap communication with computation. We illustrate our design choices by experimental results obtained on an SGI Origin 2000 and an IBM SP2 for test matrices provided by industrial partners in the PARASOL project.

A Krylov--Schur Algorithm for Large Eigenproblems

Abstract: Sorensen's implicitly restarted Arnoldi algorithm is one of the most successful and flexible methods for finding a few eigenpairs of a large matrix. However, the need to preserve the structure of the Arnoldi decomposition on which the algorithm is based restricts the range of transformations that can be performed on the decomposition. In consequence, it is difficult to deflate converged Ritz vectors from the decomposition. Moreover, the potential forward instability of the implicit QR algorithm can cause unwanted Ritz vectors to persist in the computation. In this paper we introduce a general Krylov decomposition that solves both problems in a natural and efficient manner.

Rank-One Approximation to High Order Tensors

Abstract: The singular value decomposition (SVD) has been extensively used in engineering and statistical applications. This method was originally discovered by Eckart and Young in [Psychometrika, 1 (1936), pp. 211--218], where they considered the problem of low-rank approximation to a matrix. A natural generalization of the SVD is the problem of low-rank approximation to high order tensors, which we call the multidimensional SVD. In this paper, we investigate certain properties of this decomposition as well as numerical algorithms.

Orthogonal Tensor Decompositions

Abstract: We explore the orthogonal decomposition of tensors (also known as multidimensional arrays or n-way arrays) using two different definitions of orthogonality. We present numerous examples to illustrate the difficulties in understanding such decompositions. We conclude with a counterexample to a tensor extension of the Eckart--Young SVD approximation theorem by Leibovici and Sabatier [Linear Algebra Appl., 269 (1998), pp. 307--329].

On the Best Rank-1 Approximation of Higher-Order Supersymmetric Tensors

Abstract: Recently the problem of determining the best, in the least-squares sense, rank-1 approximation to a higher-order tensor was studied and an iterative method that extends the well-known power method for matrices was proposed for its solution. This higher-order power method is also proposed for the special but important class of supersymmetric tensors, with no change. A simplified version, adapted to the special structure of the supersymmetric problem, is deemed unreliable, as its convergence is not guaranteed. The aim of this paper is to show that a symmetric version of the above method converges under assumptions of convexity (or concavity) for the functional induced by the tensor in question, assumptions that are very often satisfied in practical applications. The use of this version entails significant savings in computational complexity as compared to the unconstrained higher-order power method. Furthermore, a novel method for initializing the iterative process is developed which has been observed to yield an estimate that lies closer to the global optimum than the initialization suggested before. Moreover, its proximity to the global optimum is a priori quantifiable. In the course of the analysis, some important properties that the supersymmetry of a tensor implies for its square matrix unfolding are also studied.

Structured Pseudospectra for Polynomial Eigenvalue Problems, with Applications

Abstract: Pseudospectra associated with the standard and generalized eigenvalue problems have been widely investigated in recent years. We extend the usual definitions in two respects, by treating the polynomial eigenvalue problem and by allowing structured perturbations of a type arising in control theory. We explore connections between structured pseudospectra, structured backward errors, and structured stability radii. Two main approaches for computing pseudospectra are described. One is based on a transfer function and employs a generalized Schur decomposition of the companion form pencil. The other, specific to quadratic polynomials, finds a solvent of the associated quadratic matrix equation and thereby factorizes the quadratic $\lambda$-matrix. Possible approaches for large, sparse problems are also outlined. A collection of examples from vibrating systems, control theory, acoustics, and fluid mechanics is given to illustrate the techniques.

Nonsymmetric Algebraic Riccati Equations and Wiener--Hopf Factorization for M -Matrices

Abstract: We consider the nonsymmetric algebraic Riccati equation for which the four coefficient matrices form an M-matrix. Nonsymmetric algebraic Riccati equations of this type appear in applied probability and transport theory. The minimal nonnegative solution of these equations can be found by Newton's method and basic fixed-point iterations. The study of these equations is also closely related to the so-called Wiener--Hopf factorization for M-matrices. We explain how the minimal nonnegative solution can be found by the Schur method and compare the Schur method with Newton's method and some basic fixed-point iterations. The development in this paper parallels that for symmetric algebraic Riccati equations arising in linear quadratic control.

On an Iteration Method for Solving a Class of Nonlinear Matrix Equations

Abstract: This paper treats a set of equations of the form $X+A^{\star}{\cal F}(X)A =Q$, where ${\cal F}$ maps positive definite matrices either into positive definite matrices or into negative definite matrices, and satisfies some monotonicity property. Here A is arbitrary and Q is a positive definite matrix. It is shown that under some conditions an iteration method converges to a positive definite solution. An estimate for the rate of convergence is given under additional conditions, and some numerical results are given. Special cases are considered, which cover also particular cases of the discrete algebraic Riccati equation.

Solving a Quadratic Matrix Equation by Newton's Method with Exact Line Searches

Abstract: We show how to incorporate exact line searches into Newton's method for solving the quadratic matrix equation AX2 + BX + C = 0, where A, B and C are square matrices. The line searches are relatively inexpensive and improve the global convergence properties of Newton's method in theory and in practice. We also derive a condition number for the problem and show how to compute the backward error of an approximate solution.

The Multishift QR Algorithm. Part I: Maintaining Well-Focused Shifts and Level 3 Performance

Abstract: This paper presents a small-bulge multishift variation of the multishift QR algorithm that avoids the phenomenon of shift blurring, which retards convergence and limits the number of simultaneous shifts. It replaces the large diagonal bulge in the multishift QR sweep with a chain of many small bulges. The small-bulge multishift QR sweep admits nearly any number of simultaneous shifts---even hundreds---without adverse effects on the convergence rate. With enough simultaneous shifts, the small-bulge multishift QR algorithm takes advantage of the level 3 BLAS, which is a special advantage for computers with advanced architectures.

The Multishift QR Algorithm. Part II: Aggressive Early Deflation

Abstract: Aggressive early deflation is a QR algorithm deflation strategy that takes advantage of matrix perturbations outside of the subdiagonal entries of the Hessenberg QR iterate. It identifies and deflates converged eigenvalues long before the classic small-subdiagonal strategy would. The new deflation strategy enhances the performance of conventional large-bulge multishift QR algorithms, but it is particularly effective in combination with the small-bulge multishift QR algorithm. The small-bulge multishift QR sweep with aggressive early deflation maintains a high rate of execution of floating point operations while significantly reducing the number of operations required.

Twice Differentiable Spectral Functions

Abstract: A function F on the space of n × n real symmetric matrices is called spectral if it depends only on the eigenvalues of its argument. Spectral functions are just symmetric functions of the eigenvalues. We show that a spectral function is twice (continuously) differentiable at a matrix if and only if the corresponding symmetric function is twice (continuously) differentiable at the vector of eigenvalues. We give a concise and usable formula for the Hessian.

Implicitly Restarted Arnoldi Methods and Subspace Iteration

Abstract: This goal of this paper is to present an elegant relationship between an implicitly restarted Arnoldi method (IRAM) and nonstationary (subspace) simultaneous iteration. This relationship allows the geometric convergence theory developed for nonstationary simultaneous iteration due to Watkins and Elsner [Linear Algebra Appl., 143 (1991), pp. 19--47] to be used for analyzing the rate of convergence of an IRAM. We also comment on the relationship with other restarting schemes. A set of experiments demonstrates that implicit restarted methods can converge at a much faster rate than simultaneous iteration when iterating on a subspace of equal dimension.

Robust eigenstructure assignment in quadratic matrix polynomials: Nonsingular case

Abstract: Feedback design for a second-order control system leads to an eigenstructure assignment problem for a quadratic matrix polynomial. It is desirable that the feedback controller not only assigns specified eigenvalues to the second-order closed loop system but also that the system is robust, or insensitive to perturbations. We derive here new sensitivity measures, or condition numbers, for the eigenvalues of the quadratic matrix polynomial and define a measure of the robustness of the corresponding system. We then show that the robustness of the quadratic inverse eigenvalue problem can be achieved by solving a generalized linear eigenvalue assignment problem subject to structured perturbations. Numerically reliable methods for solving the structured generalized linear problem are developed that take advantage of the special properties of the system in order to minimize the computational work required. In this part of the work we treat the case where the leading coefficient matrix in the quadratic polynomial is nonsingular, which ensures that the polynomial is regular. In a second part, we will examine the case where the open loop matrix polynomial is not necessarily regular.

A Regularized Robust Design Criterion for Uncertain Data

Abstract: This paper formulates and solves a robust criterion for least-squares designs in the presence of uncertain data. Compared with earlier studies, the proposed criterion incorporates simultaneously both regularization and weighting and applies to a large class of uncertainties. The solution method is based on reducing a vector optimization problem to an equivalent scalar minimization problem of a provably unimodal cost function, thus achieving considerable reduction in computational complexity.

Low-Rank Approximations with Sparse Factors I: Basic Algorithms and Error Analysis

Abstract: We consider the problem of computing low-rank approximations of matrices. The novel aspects of our approach are that we require the low-rank approximations to be written in a factorized form with sparse factors, and the degree of sparsity of the factors can be traded off for reduced reconstruction error by certain user-determined parameters. We give a detailed error analysis of our proposed algorithms and compare the computed sparse low-rank approximations with those obtained from singular value decomposition. We present numerical examples arising from some application areas to illustrate the efficiency and accuracy of our algorithms.

A Stabilized Superfast Solver for Nonsymmetric Toeplitz Systems

Abstract: We present a stabilized superfast solver for nonsymmetric Toeplitz systems Tx=b. An explicit formula for T-1 is expressed in such a way that the matrix-vector product T^-1b can be calculated via FFTs and Hadamard products. This inversion formula involves certain polynomials that can be computed by solving two linearized rational interpolation problems on the unit circle. The heart of our Toeplitz solver is a superfast algorithm to solve these interpolation problems. To stabilize the algorithm, i.e., to improve the accuracy, several techniques are used: pivoting, iterative improvement, downdating, and giving "difficult" interpolation points an adequate treatment. We have implemented our algorithm in Fortran 90. Numerical examples illustrate the effectiveness of our approach.

An Iterative Method with Variable Relaxation Parameters for Saddle-Point Problems

Abstract: In this paper, we propose an inexact Uzawa method with variable relaxation parameters for iteratively solving linear saddle-point problems. The method involves two variable relaxation parameters, which can be updated easily in each iteration, similar to the evaluation of the two iteration parameters in the conjugate gradient method. This new algorithm has an advantage over most existing Uzawa-type algorithms: it is always convergent without any a priori estimates on the spectrum of the preconditioned Schur complement matrix, which may not be easy to achieve in applications. The rate of the convergence of the inexact Uzawa method is analyzed. Numerical results of the algorithm applied for the Stokes problem and a purely linear system of algebraic equations are presented.

The Singularity-Induced Bifurcation and its Kronecker Normal Form

Abstract: It is shown that the singularity-induced bifurcation theorem due to Venkatasubramanian, Schattler, and Zaborszky [ Proceedings of the IEEE, 83 (1995), pp. 1530--1558] can be expressed as the perturbation of an infinite eigenvalue of a particular class of parameterized index-1 matrix pencil, denoted $(M,L(\lambda))$. It is shown that the matrix pencil at the singularity-induced bifurcation point, $(M,L(\lambda_{0}))$, has Kronecker index-2. It is also shown that a two-parameter unfolding of a singularity-induced bifurcation point results in a locus of index-0 pencils, denoted $(M(\epsilon),L(\lambda(\epsilon)))$, which has two purely imaginary eigenvalues near infinity.

Existence, Uniqueness, and Parametrization of Lagrangian Invariant Subspaces

Abstract: The existence, uniqueness, and parametrization of Lagrangian invariant subspaces for Hamiltonian matrices is studied. Necessary and sufficient conditions and a complete parametrization are given. Some necessary and sufficient conditions for the existence of Hermitian solutions of algebraic Riccati equations follow as simple corollaries.

A Spectral Characterization of Generalized Real Symmetric Centrosymmetric and Generalized Real Symmetric Skew-Centrosymmetric Matrices

Abstract: We show that the only real symmetric matrices whose spectrum is invariant modulo sign changes after either row or column reversal are the centrosymmetric matrices; moreover, we prove that the class of real symmetric centrosymmetric matrices can be completely characterized by this property. We also show that the only real symmetric matrices whose spectrum changes by multiplication by i after either row or column reversal are the skew-centrosymmetric matrices; here, too, we show that the class of real symmetric skew-centrosymmetric matrices can be completely characterized by this property of their eigenvalues. We prove both of these spectral characterizations as special cases of results for what we've called generalized centrosymmetric K-matrices and generalized skew-centrosymmetric K-matrices. Some results illustrating the application of the generalized centrosymmetric spectral characterization to other classes of real symmetric matrices are also given.

Generalized Polar Decompositions for the Approximation of the Matrix Exponential

Abstract: In this paper we describe the use of the theory of generalized polar decompositions [H. Munthe-Kaas, G. R. W. Quispel, and A. Zanna, Found. Comput. Math., 1 (2001), pp. 297--324] to approximate a matrix exponential. The algorithms presented have the property that, if $Z \in {\frak{g}}$, a Lie algebra of matrices, then the approximation for exp(Z) resides in G, the matrix Lie group of ${\frak{g}}$. This property is very relevant when solving Lie-group ODEs and is not usually fulfilled by standard approximations to the matrix exponential. We propose algorithms based on a splitting of Z into matrices having a very simple structure, usually one row and one column (or a few rows and a few columns), whose exponential is computed very cheaply to machine accuracy. The proposed methods have a complexity of ${\cal O}(\kappa n^{3})$, with constant $\kappa$ small, depending on the order and the Lie algebra ${\frak{g}}$. % The algorithms are recommended in cases where it is of fundamental importance that the approximation for the exponential resides in G, and when the order of approximation needed is not too high. We present in detail algorithms up to fourth order.

A Shifted Cyclic Reduction Algorithm for Quasi-Birth-Death Problems

Abstract: The problem of the computation of the stochastic matrix G associated with discrete-time quasi-birth-death (QBD) Markov chains is analyzed. We present a shifted cyclic reduction algorithm and show that the speed of convergence of the latter modified algorithm is always faster than that of the original cyclic reduction.

Spectral Factorizations and Sums of Squares Representations via Semidefinite Programming

Abstract: In this paper we find a characterization for when a multivariable trigonometric polynomial can be written as a sum of squares. In addition, the truncated moment problem is addressed. A numerical algorithm for finding a sum of squares representation is presented as well. In the one-variable case, the algorithm finds a spectral factorization. The latter may also be used to find inner-outer factorizations.

Adaptive Eigenvalue Computations Using Newton's Method on the Grassmann Manifold

Abstract: We consider the problem of updating an invariant subspace of a large and structured Hermitian matrix when the matrix is modified slightly. The problem can be formulated as that of computing stationary values of a certain function with orthogonality constraints. The constraint is formulated as the requirement that the solution must be on the Grassmann manifold, and Newton's method on the manifold is used. In each Newton iteration a Sylvester equation is to be solved. We discuss the properties of the Sylvester equation and conclude that for large problems preconditioned iterative methods can be used. Preconditioning techniques are discussed. Numerical examples from signal subspace computations are given in which the matrix is Toeplitz and we compute a partial singular value decomposition corresponding to the largest singular values. Further we solve numerically the problem of computing the smallest eigenvalues and corresponding eigenvectors of a large sparse matrix that has been slightly modified.

Nested-Dissection Orderings for Sparse LU with Partial Pivoting

Abstract: We describe the implementation and performance of a novel fill-minimization ordering technique for sparse LU factorization with partial pivoting. The technique was proposed by Gilbert and Schreiber in 1980 but never implemented and tested. Like other techniques for ordering sparse matrices for LU with partial pivoting, our new method preorders the columns of the matrix (the row permutation is chosen by the pivoting sequence during the numerical factorization). Also like other methods, the column permutation Q that we select is a permutation that attempts to reduce the fill in the Cholesky factor of QT ATAQ. Unlike existing column-ordering techniques, which all rely on minimum-degree heuristics, our new method is based on a nested-dissection ordering of A T A. Our algorithm, however, never computes a representation of AT A, which can be expensive. We only work with a representation of A itself. Our experiments demonstrate that the method is efficient and that it can reduce fill significantly relative to the best existing methods. The method reduces the LU running time on some very large matrices (tens of millions of nonzeros in the factors) by more than a factor of 2.

Convergence Rate of an Iterative Method for a Nonlinear Matrix Equation

Abstract: We prove a convergence result for an iterative method, proposed recently by Meini, for finding the maximal Hermitian positive definite solution of the matrix equation X+A*X-1A=Q, where Q is Hermitian positive definite.

Analysis of the Cholesky Method with Iterative Refinement for Solving the Symmetric Definite Generalized Eigenproblem

Abstract: A standard method for solving the symmetric definite generalized eigenvalue problem $Ax = \lambda Bx$, where A is symmetric and B is symmetric positive definite, is to compute a Cholesky factorization B = LLT (optionally with complete pivoting) and solve the equivalent standard symmetric eigenvalue problem $C y = \lambda y$, where C = L-1 AL-T. Provided that a stable eigensolver is used, standard error analysis says that the computed eigenvalues are exact for $A+\Delta A$ and $B+\Delta B$ with $\max({\|\Delta A\|_2}/{\|A\|_2}, {\|\Delta B\|_2}/{\|B\|_2})$ bounded by a multiple of $\kappa_2(B)u$, where $u$ is the unit roundoff. We take the Jacobi method as the eigensolver and give a detailed error analysis that yields backward error bounds potentially much smaller than $\kappa_2(B)u$. To show the practical utility of our bounds we describe a vibration problem from structural engineering in which $B$ is ill conditioned yet the error bounds are small. We show how, in cases of instability, iterative refinement based on Newton's method can be used to produce eigenpairs with small backward errors. Our analysis and experiments also give insight into the popular Cholesky--QR method, in which the QR method is used as the eigensolver. We argue that it is desirable to augment current implementations of this method with pivoting in the Cholesky factorization.

On a Question Concerning Condition Numbers for Markov Chains

Abstract: Let S be an irreducible stochastic matrix of order n with left stationary vector $\pi^T,$ and let S(i) denote the principal submatrix of S formed by deleting the ith row and column. We prove that $\max_{1 \le i \le n}\pi_i||(I- S_{(i)})^{-1}||_{\infty} \le \min_{1 \le j \le n}||(I- S_{(j)})^{-1}||_{\infty},$ thus answering a question posed by Cho and Meyer. We provide an attainable lower bound on $\max_{1 \le i \le n}\pi_i||(I- S_{(i)})^{-1}||_{\infty},$ and discuss the case that equality holds in that bound.

Perturbation Bounds of Unitary and Subunitary Polar Factors

Abstract: In this paper, we present some new perturbation bounds for (generalized) polar decompositions under the Frobenius norm for both complex and real matrices. For subunitary polar factors, we show that our bounds always improve the existing bounds. Based on some interesting properties of the matrix equation W+W*=W* W, some new bounds involving both the Frobenius norm and the spectral norm of the perturbation are given. The optimality of bounds is discussed.