SIAM Digital Library
 
 
 

SIAM J. on Matrix Analysis and Applications

All Journal content published prior to 1997 is part of LOCUS.

Search Issue | RSS Feeds RSS
Previous Issue

2007

Volume 29, Issue 4, pp. 1051-1410


Approximate Diagonalization

E. B. Davies

SIAM. J. Matrix Anal. & Appl. 29, pp. 1051-1064 (14 pages)

Online Publication Date: November 02, 2007

Full Text: | Download PDF

Show Abstract
We describe a new method of computing functions of highly nonnormal matrices by using the concept of approximate diagonalization. We formulate a conjecture about its efficiency and provide both theoretical and numerical evidence in support of the conjecture. We apply the method to compute arbitrary real powers of highly nonnormal matrices.

Analysis and Exploitation of Matrix Structure Arising in Linearized Optical Tomographic Imaging

Damon Hyde, Misha Kilmer, Dana H. Brooks, and Eric Miller

SIAM. J. Matrix Anal. & Appl. 29, pp. 1065-1082 (18 pages)

Online Publication Date: November 02, 2007

Full Text: | Download PDF

Show Abstract
We present a novel method by which the large dense forward matrix $\mathbf{A}$ involved in a linear inverse diffusion problem can be decomposed into a number of sparse easily computed matrices. We begin by introducing an errorless decomposition which is applicable to a wide array of such imaging problems. Next, we incorporate interpolation into the construction of the matrices to reduce the computational complexity involved in the matrix-vector multiplications necessary to obtain an inverse solution. Error and computational complexity analysis are provided to support these developments. We then present numerical results that illustrate the gain in computational efficiency when the approximation is used in the Tikhonov regularized inverse problem, and show that the use of the approximation has virtually no negative effect on the quality of the reconstructed images. Finally, we discuss applicability to other imaging problems.

On the Doubling Algorithm for a (Shifted) Nonsymmetric Algebraic Riccati Equation

Chun-Hua Guo, Bruno Iannazzo, and Beatrice Meini

SIAM. J. Matrix Anal. & Appl. 29, pp. 1083-1100 (18 pages) | Cited 4 times

Online Publication Date: November 09, 2007

Full Text: | Download PDF

Show Abstract
Nonsymmetric algebraic Riccati equations for which the four coefficient matrices form an irreducible $M$-matrix $M$ are considered. The emphasis is on the case where $M$ is an irreducible singular $M$-matrix, which arises in the study of Markov models. The doubling algorithm is considered for finding the minimal nonnegative solution, the one of practical interest. The algorithm has been recently studied by others for the case where $M$ is a nonsingular $M$-matrix. A shift technique is proposed to transform the original Riccati equation into a new Riccati equation for which the four coefficient matrices form a nonsingular matrix. The convergence of the doubling algorithm is accelerated when it is applied to the shifted Riccati equation.

Block Diagonal and Schur Complement Preconditioners for Block-Toeplitz Systems with Small Size Blocks

Wai-Ki Ching, Michael K. Ng, and You-Wei Wen

SIAM. J. Matrix Anal. & Appl. 29, pp. 1101-1119 (19 pages)

Online Publication Date: November 09, 2007

Full Text: | Download PDF

Show Abstract
In this paper we consider the solution of Hermitian positive definite block-Toeplitz systems with small size blocks. We propose and study block diagonal and Schur complement preconditioners for such block-Toeplitz matrices. We show that for some block-Toeplitz matrices, the spectra of the preconditioned matrices are uniformly bounded except for a fixed number of outliers where this fixed number depends only on the size of the block. Hence, conjugate gradient type methods, when applied to solving these preconditioned block-Toeplitz systems with small size blocks, converge very fast. Recursive computation of such block diagonal and Schur complement preconditioners is considered by using the nice matrix representation of the inverse of a block-Toeplitz matrix. Applications to block-Toeplitz systems arising from least squares filtering problems and queueing networks are presented. Numerical examples are given to demonstrate the effectiveness of the proposed method.

Matrix Nearness Problems with Bregman Divergences

Inderjit S. Dhillon and Joel A. Tropp

SIAM. J. Matrix Anal. & Appl. 29, pp. 1120-1146 (27 pages) | Cited 4 times

Online Publication Date: November 21, 2007

Full Text: | Download PDF

Show Abstract
This paper discusses a new class of matrix nearness problems that measure approximation error using a directed distance measure called a Bregman divergence. Bregman divergences offer an important generalization of the squared Frobenius norm and relative entropy, and they all share fundamental geometric properties. In addition, these divergences are intimately connected with exponential families of probability distributions. Therefore, it is natural to study matrix approximation problems with respect to Bregman divergences. This article proposes a framework for studying these problems, discusses some specific matrix nearness problems, and provides algorithms for solving them numerically. These algorithms apply to many classical and novel problems, and they admit a striking geometric interpretation.

A Givens-Weight Representation for Rank Structured Matrices

Steven Delvaux and Marc Van Barel

SIAM. J. Matrix Anal. & Appl. 29, pp. 1147-1170 (24 pages) | Cited 3 times

Online Publication Date: November 21, 2007

Full Text: | Download PDF

Show Abstract
In this paper we introduce a Givens-weight representation for rank structured matrices, where the rank structure is defined by certain submatrices starting from the bottom left or upper right matrix corner being of low rank. We proceed in two steps. First, we introduce a unitary-weight representation. This representation will be compared to the (block) quasiseparable representations introduced by P. Dewilde and A.-J. van der Veen [Time-varying Systems and Computations, Kluwer Academic Publishers, Boston, 1998]. More specifically, we show that our unitary-weight representations are theoretically equivalent to the so-called block quasiseparable representations in input or output normal form introduced by Dewilde and van der Veen [Time varying Systems and Computations, Kluwer Academic Publishers, Boston, 1998]. Next, we move from the unitary-weight to the Givens-weight representation. We then provide some basic algorithms for the unitary/Givens-weight representation, showing how to obtain such a representation for a dense matrix by means of numerical approximation. We also show how to “swap” the representation and how to reduce the number of parameters of the representation, whenever appropriate. As such, these results will become the basis for algorithms on unitary/Givens-weight representations to be described in subsequent papers.

When is the Adjoint of a Matrix a Low Degree Rational Function in the Matrix?

Jörg Liesen

SIAM. J. Matrix Anal. & Appl. 29, pp. 1171-1180 (10 pages) | Cited 2 times

Online Publication Date: December 07, 2007

Full Text: | Download PDF

Show Abstract
We show that the adjoint $A^+$ of a matrix $A$ with respect to a given inner product is a rational function in $A$, if and only if $A$ is normal with respect to the inner product. We consider such matrices and analyze the McMillan degrees of the rational functions $r$ such that $A^+=r(A)$. We introduce the McMillan degree of $A$ as the smallest among these degrees, characterize this degree in terms of the number and distribution of the eigenvalues of $A$, and compare the McMillan degree with the normal degree of $A$, which is defined as the smallest degree of a polynomial $p$ for which $A^+=p(A)$. We show that unless the eigenvalues of $A$ lie on a single circle in the complex plane, the ratio of the normal degree and the McMillan degree of $A$ is bounded by a small constant that depends neither on the number nor on the distribution of the eigenvalues of $A$. Our analysis is motivated by applications in the area of short recurrence Krylov subspace methods.

On the Index of Conditional Stability of Stable Invariant Lagrangian Subspaces

André C. M. Ran and Leiba Rodman

SIAM. J. Matrix Anal. & Appl. 29, pp. 1181-1190 (10 pages)

Online Publication Date: December 07, 2007

Full Text: | Download PDF

Show Abstract
Given a nondegenerate sesquilinear inner product on a finite dimensional complex vector space, or a nondegenerate symmetric or skewsymmetric inner product on finite dimensional real vector space, subspaces that are simultaneously Lagrangian and invariant for a selfadjoint or a skewadjoint matrix with respect to the inner product are considered. The rate of conditional stability of such subspaces is studied, under small perturbations of both the inner product and the matrix. The concept of conditional stability (in contrast with unconditional stability) presupposes that one considers only those perturbed matrix and inner product for which the existence of invariant Lagrangian subspaces can be guaranteed a priori. Open problems regarding the index (= exact rate) of conditional stability are stated. Several inaccurate statements in the authors' previous works concerning the index are made precise. Finally, an application is given to conditional stability of hermitian solutions of continuous type algebraic Riccati equations.

Recursive Solution of Certain Structured Linear Systems

André Klein and Peter Spreij

SIAM. J. Matrix Anal. & Appl. 29, pp. 1191-1217 (27 pages)

Online Publication Date: December 07, 2007

Full Text: | Download PDF

Show Abstract
We provide explicit representations of the null space $\mathcal{S}$ of adjoints of companion-related matrices and of certain rectangular generalized Vandermonde matrices of block Toeplitz type which are encountered in the Fisher information matrix of time series processes. A formula for the right-inverse of this class of matrices $A$ is provided which allows one to express the solution of the system $Ax=b$ as $x=A^{-}b+$ $\mathcal{S}$. The formulas can be easily turned into solution algorithms.

Backward Error of Polynomial Eigenproblems Solved by Linearization

Nicholas J. Higham, Ren-Cang Li, and Françoise Tisseur

SIAM. J. Matrix Anal. & Appl. 29, pp. 1218-1241 (24 pages) | Cited 10 times

Online Publication Date: December 07, 2007

Full Text: | Download PDF

Show Abstract
The most widely used approach for solving the polynomial eigenvalue problem $P(\lambda)x = (\sum_{i=0}^m \l^i A_i) x = 0$ in $n\times n$ matrices $A_i$ is to linearize to produce a larger order pencil $L(\lambda) = \lambda X + Y$, whose eigensystem is then found by any method for generalized eigenproblems. For a given polynomial $P$, infinitely many linearizations $L$ exist and approximate eigenpairs of $P$ computed via linearization can have widely varying backward errors. We show that if a certain one-sided factorization relating $L$ to $P$ can be found then a simple formula permits recovery of right eigenvectors of $P$ from those of $L$, and the backward error of an approximate eigenpair of $P$ can be bounded in terms of the backward error for the corresponding approximate eigenpair of $L$. A similar factorization has the same implications for left eigenvectors. We use this technique to derive backward error bounds depending only on the norms of the $A_i$ for the companion pencils and for the vector space $\mathbb{DL}(P)$ of pencils recently identified by Mackey, Mackey, Mehl, and Mehrmann. In all cases, sufficient conditions are identified for an optimal backward error for $P$. These results are shown to be entirely consistent with those of Higham, Mackey, and Tisseur on the conditioning of linearizations of $P$. Other contributions of this work are a block scaling of the companion pencils that yields improved backward error bounds; a demonstration that the bounds are applicable to certain structured linearizations of structured polynomials; and backward error bounds specialized to the quadratic case, including analysis of the benefits of a scaling recently proposed by Fan, Lin, and Van Dooren. The results herein make no assumptions on the stability of the method applied to $L$ or whether the method is direct or iterative.

Further Results on the Reverse Order Law for Generalized Inverses

Dragan S. Djordjević

SIAM. J. Matrix Anal. & Appl. 29, pp. 1242-1246 (5 pages)

Online Publication Date: December 13, 2007

Full Text: | Download PDF

Show Abstract
The reverse order rule $(AB)^\dag=B^\dag A^\dag$ for the Moore–Penrose inverse is established in several equivalent forms. Results related to other generalized inverses are also proved.

A Superfast Algorithm for Toeplitz Systems of Linear Equations

S. Chandrasekaran, M. Gu, X. Sun, J. Xia, and J. Zhu

SIAM. J. Matrix Anal. & Appl. 29, pp. 1247-1266 (20 pages) | Cited 2 times

Online Publication Date: December 13, 2007

Full Text: | Download PDF

Show Abstract
In this paper we develop a new superfast solver for Toeplitz systems of linear equations. To solve Toeplitz systems many people use displacement equation methods. With displacement structures, Toeplitz matrices can be transformed into Cauchy-like matrices using the FFT or other trigonometric transformations. These Cauchy-like matrices have a special property, that is, their off-diagonal blocks have small numerical ranks. This low-rank property plays a central role in our superfast Toeplitz solver. It enables us to quickly approximate the Cauchy-like matrices by structured matrices called sequentially semiseparable (SSS) matrices. The major work of the constructions of these SSS forms can be done in precomputations (independent of the Toeplitz matrix entries). These SSS representations are compact because of the low-rank property. The SSS Cauchy-like systems can be solved in linear time with linear storage. Excluding precomputations the main operations are the FFT and SSS system solve, which are both very efficient. Our new Toeplitz solver is stable in practice. Numerical examples are presented to illustrate the efficiency and the practical stability.

Steepest Descent and Conjugate Gradient Methods with Variable Preconditioning

Andrew V. Knyazev and Ilya Lashuk

SIAM. J. Matrix Anal. & Appl. 29, pp. 1267-1280 (14 pages)

Online Publication Date: December 19, 2007

Full Text: | Download PDF

Show Abstract
We analyze the conjugate gradient (CG) method with variable preconditioning for solving a linear system with a real symmetric positive definite (SPD) matrix of coefficients $A$. We assume that the preconditioner is SPD on each step, and that the condition number of the preconditioned system matrix is bounded above by a constant independent of the step number. We show that the CG method with variable preconditioning under this assumption may not give improvement, compared to the steepest descent (SD) method. We describe the basic theory of CG methods with variable preconditioning with the emphasis on “worst case” scenarios, and provide complete proofs of all facts not available in the literature. We give a new elegant geometric proof of the SD convergence rate bound. Our numerical experiments, comparing the preconditioned SD and CG methods, not only support and illustrate our theoretical findings, but also reveal two surprising and potentially practically important effects. First, we analyze variable preconditioning in the form of inner-outer iterations. In previous such tests, the unpreconditioned CG inner iterations are applied to an artificial system with some fixed preconditioner as a matrix of coefficients. We test a different scenario, where the unpreconditioned CG inner iterations solve linear systems with the original system matrix $A$. We demonstrate that the CG-SD inner-outer iterations perform as well as the CG-CG inner-outer iterations in these tests. Second, we compare the CG methods using a two-grid preconditioning with fixed and randomly chosen coarse grids, and observe that the fixed preconditioner method is twice as slow as the method with random preconditioning.

PageRank Computation, with Special Attention to Dangling Nodes

Ilse C. F. Ipsen and Teresa M. Selee

SIAM. J. Matrix Anal. & Appl. 29, pp. 1281-1296 (16 pages) | Cited 3 times

Online Publication Date: December 19, 2007

Full Text: | Download PDF

Show Abstract
We present a simple algorithm for computing the PageRank (stationary distribution) of the stochastic Google matrix $G$. The algorithm lumps all dangling nodes into a single node. We express lumping as a similarity transformation of $G$ and show that the PageRank of the nondangling nodes can be computed separately from that of the dangling nodes. The algorithm applies the power method only to the smaller lumped matrix, but the convergence rate is the same as that of the power method applied to the full matrix $G$. The efficiency of the algorithm increases as the number of dangling nodes increases. We also extend the expression for PageRank and the algorithm to more general Google matrices that have several different dangling node vectors, when it is required to distinguish among different classes of dangling nodes. We also analyze the effect of the dangling node vector on the PageRank and show that the PageRank of the dangling nodes depends strongly on that of the nondangling nodes but not vice versa. Last we present a Jordan decomposition of the Google matrix for the (theoretical) extreme case when all Web pages are dangling nodes.

Maximum Attainable Accuracy of Inexact Saddle Point Solvers

Pavel Jiránek and Miroslav Rozložník

SIAM. J. Matrix Anal. & Appl. 29, pp. 1297-1321 (25 pages)

Online Publication Date: January 04, 2008

Full Text: | Download PDF

Show Abstract
In this paper we study numerical behavior of several iterative Krylov subspace solvers applied to the solution of large-scale saddle point problems. Two main representatives of segregated solution approach are analyzed: the Schur complement reduction method based on the elimination of primary unknowns and the null-space projection method, which relies on a basis for the subspace described by the constraints. We show that the choice of the back-substitution formula may considerably influence the maximum attainable accuracy of approximate solutions computed in finite precision arithmetic.

New Fast and Accurate Jacobi SVD Algorithm. I

Zlatko Drmač and Krešimir Veselić

SIAM. J. Matrix Anal. & Appl. 29, pp. 1322-1342 (21 pages) | Cited 2 times

Online Publication Date: January 04, 2008

Full Text: | Download PDF

Show Abstract
This paper is the result of concerted efforts to break the barrier between numerical accuracy and run-time efficiency in computing the fundamental decomposition of numerical linear algebra—the singular value decomposition (SVD) of general dense matrices. It is an unfortunate fact that the numerically most accurate one-sided Jacobi SVD algorithm is several times slower than generally less accurate bidiagonalization-based methods such as the QR or the divide-and-conquer algorithm. Our quest for a highly accurate and efficient SVD algorithm has led us to a new, superior variant of the Jacobi algorithm. The new algorithm has inherited all good high accuracy properties of the Jacobi algorithm, and it can outperform the QR algorithm.

New Fast and Accurate Jacobi SVD Algorithm. II

Zlatko Drmač and Krešimir Veselić

SIAM. J. Matrix Anal. & Appl. 29, pp. 1343-1362 (20 pages) | Cited 2 times

Online Publication Date: January 04, 2008

Full Text: | Download PDF

Show Abstract
This paper presents a new one-sided Jacobi SVD algorithm for triangular matrices computed by revealing QR factorizations. If used in the preconditioned Jacobi SVD algorithm, described in part one of this paper, it delivers superior performance leading to the currently fastest method for computing SVD decomposition with high relative accuracy. Furthermore, the efficiency of the new algorithm is comparable to the less accurate bidiagonalization-based methods. The paper also discusses underflow issues in floating point implementation and shows how to use perturbation theory to fix the imperfectness of the machine arithmetic.

Algorithmic Aspects of Elimination Trees for Sparse Unsymmetric Matrices

Stanley C. Eisenstat and Joseph W. H. Liu

SIAM. J. Matrix Anal. & Appl. 29, pp. 1363-1381 (19 pages)

Online Publication Date: January 04, 2008

Full Text: | Download PDF

Show Abstract
The elimination tree of a symmetric matrix plays an important role in sparse elimination. We recently defined a generalization of this structure to the unsymmetric case that retains many of its properties. Here we present an algorithm for constructing the elimination tree of an unsymmetric matrix and show how it can be used to find a symmetric reordering of the matrix into a recursive, bordered block triangular form. We also present two symbolic factorization algorithms that use the elimination tree to determine the nonzero structures of the triangular factors of such matrices. Numerical experiments demonstrate that these algorithms are efficient and compare the new symbolic factorization schemes with existing ones.

On the Convergence of General Stationary Linear Iterative Methods for Singular Linear Systems

Zhi-Hao Cao

SIAM. J. Matrix Anal. & Appl. 29, pp. 1382-1388 (7 pages) | Cited 2 times

Online Publication Date: January 04, 2008

Full Text: | Download PDF

Show Abstract
Recently, Lee et al. have published an interesting paper [SIAM J. Matrix Anal. Appl., 28 (2006), pp. 634–641] concerning the energy norm convergence of general stationary linear iterative methods for semidefinite linear systems. In this paper, we first consider the convergence of general stationary linear iterative methods for general singular consistent linear systems and show that the convergence and the quotient convergence are equivalent. Then we consider the convergence of general stationary iterative methods for the semidefinite systems and clarify some issues in Lee et al.'s paper.

Structured Mapping Problems for Matrices Associated with Scalar Products. Part I: Lie and Jordan Algebras

D. Steven Mackey, Niloufer Mackey, and Françoise Tisseur

SIAM. J. Matrix Anal. & Appl. 29, pp. 1389-1410 (22 pages) | Cited 3 times

Online Publication Date: January 09, 2008

Full Text: | Download PDF

Show Abstract
Given a class of structured matrices $\mathbb{S}$, we identify pairs of vectors $x,b$ for which there exists a matrix $A\in\mathbb{S}$ such that $Ax=b$, and we also characterize the set of all matrices $A\in\mathbb{S}$ mapping $x$ to $b$. The structured classes we consider are the Lie and Jordan algebras associated with orthosymmetric scalar products. These include (skew-)symmetric, (skew-)Hamiltonian, pseudo(skew-)Hermitian, persymmetric, and perskew-symmetric matrices. Structured mappings with extremal properties are also investigated. In particular, structured mappings of minimal rank are identified and shown to be unique when rank one is achieved. The structured mapping of minimal Frobenius norm is always unique, and explicit formulas for it and its norm are obtained. Finally the set of all structured mappings of minimal 2-norm is characterized. Our results generalize and unify existing work, answer a number of open questions, and provide useful tools for structured backward error investigations.
Close

close