SIAM Digital Library
 
 
 

SIAM J. on Matrix Analysis and Applications

BROWSE VOLUMES

Year Range: 

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

Search Issue | RSS Feeds RSS
Previous Issue Next Issue

2009

Volume 31, Issue 3, pp. 865-1539


Polynomial-Time Computation of the Joint Spectral Radius for Some Sets of Nonnegative Matrices

Vincent D. Blondel and Yurii Nesterov

SIAM. J. Matrix Anal. & Appl. 31, pp. 865-876 (12 pages) | Cited 1 time

Online Publication Date: August 05, 2009

Full Text: | Download PDF

Show Abstract
We propose two simple upper bounds for the joint spectral radius of sets of nonnegative matrices. These bounds, the joint column radius and the joint row radius, can be computed in polynomial time as solutions of convex optimization problems. We show that these bounds are within a factor $1/n$ of the exact value, where $n$ is the size of the matrices. Moreover, for sets of matrices with independent column uncertainties or with independent row uncertainties, the corresponding bounds coincide with the joint spectral radius. In these cases, the joint spectral radius is also given by the largest spectral radius of the matrices in the set. As a by-product of these results, we propose a polynomial-time technique for solving Boolean optimization problems related to the spectral radius. We also describe economics and engineering applications of our results.

Convergence Analysis of Iterative Solvers in Inexact Rayleigh Quotient Iteration

Fei Xue and Howard C. Elman

SIAM. J. Matrix Anal. & Appl. 31, pp. 877-899 (23 pages) | Cited 1 time

Online Publication Date: August 05, 2009

Full Text: | Download PDF

Show Abstract
We present a detailed convergence analysis of preconditioned MINRES for approximately solving the linear systems that arise when Rayleigh quotient iteration is used to compute the lowest eigenpair of a symmetric positive definite matrix. We provide insight into the initial stagnation of MINRES iteration in both a qualitative and quantitative way and show that the convergence of MINRES mainly depends on how quickly the unique negative eigenvalue of the preconditioned shifted coefficient matrix is approximated by its corresponding harmonic Ritz value. By exploring when the negative Ritz value appears in MINRES iteration, we obtain a better understanding of the limitation of preconditioned MINRES in this context and the virtue of a new type of preconditioner with “tuning.” A comparison of MINRES with SYMMLQ in this context is also given. Finally, we show that tuning based on a rank-2 modification can be applied with little additional cost to guarantee positive definiteness of the tuned preconditioner.

Structured Eigenvalue Condition Number and Backward Error of a Class of Polynomial Eigenvalue Problems

Shreemayee Bora

SIAM. J. Matrix Anal. & Appl. 31, pp. 900-917 (18 pages) | Cited 1 time

Online Publication Date: August 05, 2009

Full Text: | Download PDF

Show Abstract
Necessary and sufficient conditions are obtained for simple eigenvalues of complex matrix polynomials with $\star$-even/odd and $\star$-palindromic/antipalindromic structures to have the same normwise condition number with respect to structure preserving and arbitrary perturbations. Here, $\star$ denotes either the transpose $T$ or the conjugate transpose $*$. Exact expressions for the normwise structured condition number of simple eigenvalues of $*$-palindromic/antipalindromic and $T$-even/odd polynomials and tight bounds that localize the structured condition number of simple eigenvalues of $T$-palindromic/antipalindromic and $*$-even/odd polynomials are also obtained. The backward error of complex numbers $z$ as approximate eigenvalues of these polynomials is also considered, and necessary and sufficient conditions are obtained for a given complex number $z$ to have the same structured and unstructured backward errors.

Regularized Total Least Squares: Computational Aspects and Error Bounds

Shuai Lu, Sergei V. Pereverzev, and Ulrich Tautenhahn

SIAM. J. Matrix Anal. & Appl. 31, pp. 918-941 (24 pages)

Online Publication Date: August 05, 2009

Full Text: | Download PDF

Show Abstract
For solving linear ill-posed problems, regularization methods are required when the right-hand side and/or the operator are corrupted by some noise. In the present paper, regularized solutions are constructed using regularized total least squares (RTLS) and dual regularized total least squares (DRTLS). We discuss computational aspects and provide order optimal error bounds that characterize the accuracy of the regularized solutions. The results extend earlier results where the operator is exactly given. We also present some numerical experiments, which shed light on the relationship between RTLS, DRTLS, and the standard Tikhonov regularization.

Shift-Invert Arnoldi's Method with Preconditioned Iterative Solves

Melina A. Freitag and Alastair Spence

SIAM. J. Matrix Anal. & Appl. 31, pp. 942-969 (28 pages) | Cited 1 time

Online Publication Date: August 12, 2009

Full Text: | Download PDF

Show Abstract
We consider the computation of a few eigenvectors and corresponding eigenvalues of a large sparse nonsymmetric matrix using shift-invert Arnoldi's method with and without implicit restarts. For the inner iterations we use preconditioned GMRES as the inexact iterative solver. The costs of the solves are measured by the number of inner iterations needed by the iterative solver at each outer step of the algorithm. We first extend the relaxation strategy developed by Simoncini [SIAM J. Numer. Anal., 43 (2005), pp. 1155–1174] to implicitly restarted Arnoldi's method, which yields an improvement in the overall costs of the method. Secondly, we apply a new preconditioning strategy to the inner solver. We show that small rank changes to the preconditioner can produce significant savings in the total number of iterations. The combination of the new preconditioner with the relaxation strategy in implicitly restarted Arnoldi produces enhancement in the overall costs of around 50 percent in the examples considered here. Numerical experiments illustrate the theory throughout the paper.

A New Scaling and Squaring Algorithm for the Matrix Exponential

Awad H. Al-Mohy and Nicholas J. Higham

SIAM. J. Matrix Anal. & Appl. 31, pp. 970-989 (20 pages) | Cited 3 times

Online Publication Date: August 12, 2009

Full Text: | Download PDF

Show Abstract
The scaling and squaring method for the matrix exponential is based on the approximation $e^A \approx (r_m(2^{-s}A))^{2^s}$, where $r_m(x)$ is the $[m/m]$ Padé approximant to $e^x$ and the integers $m$ and $s$ are to be chosen. Several authors have identified a weakness of existing scaling and squaring algorithms termed overscaling, in which a value of $s$ much larger than necessary is chosen, causing a loss of accuracy in floating point arithmetic. Building on the scaling and squaring algorithm of Higham [SIAM J. Matrix Anal. Appl., 26 (2005), pp. 1179–1193], which is used by the MATLAB function expm, we derive a new algorithm that alleviates the overscaling problem. Two key ideas are employed. The first, specific to triangular matrices, is to compute the diagonal elements in the squaring phase as exponentials instead of from powers of $r_m$. The second idea is to base the backward error analysis that underlies the algorithm on members of the sequence $\{\|A^k\|^{1/k}\}$ instead of $\|A\|$, since for nonnormal matrices it is possible that $\|A^k\|^{1/k}$ is much smaller than $\|A\|$, and indeed this is likely when overscaling occurs in existing algorithms. The terms $\|A^k\|^{1/k}$ are estimated without computing powers of $A$ by using a matrix 1-norm estimator in conjunction with a bound of the form $\|A^k\|^{1/k} \le \max\bigl( \|A^p\|^{1/p}, \|A^q\|^{1/q} \bigr)$ that holds for certain fixed $p$ and $q$ less than $k$. The improvements to the truncation error bounds have to be balanced by the potential for a large $\|A\|$ to cause inaccurate evaluation of $r_m$ in floating point arithmetic. We employ rigorous error bounds along with some heuristics to ensure that rounding errors are kept under control. Our numerical experiments show that the new algorithm generally provides accuracy at least as good as the existing algorithm of Higham at no higher cost, while for matrices that are triangular or cause overscaling it usually yields significant improvements in accuracy, cost, or both.

Optimal Conditioning of Bernstein Collocation Matrices

Jorge Delgado and J. M. Peña

SIAM. J. Matrix Anal. & Appl. 31, pp. 990-996 (7 pages)

Online Publication Date: August 12, 2009

Full Text: | Download PDF

Show Abstract
Some results on the optimal conditioning of the collocation matrices of the Bernstein basis are provided. Numerical examples illustrate the results.

Preconditioned Hermitian and Skew-Hermitian Splitting Method for Finite Element Approximations of Convection-Diffusion Equations

Alessandro Russo and Cristina Tablino Possio

SIAM. J. Matrix Anal. & Appl. 31, pp. 997-1018 (22 pages)

Online Publication Date: August 21, 2009

Full Text: | Download PDF

Show Abstract
A two-step preconditioned iterative method based on the Hermitian and skew-Hermitian splitting is applied to the solution of nonsymmetric linear systems arising from the finite element approximation of diffusion-dominated convection-diffusion equations. The theoretical spectral analysis focuses on the case of matrix sequences related to finite element approximations on uniform structured meshes, by referring to spectral tools derived from Toeplitz theory. In such a setting, if the problem is coercive and the diffusive and convective coefficients are regular enough, then the proposed preconditioned matrix sequence shows a strong clustering at unity, i.e., a superlinear preconditioning sequence is obtained. Under the same assumptions, the optimality of the preconditioned Hermitian and skew-Hermitian splitting (PHSS) method is proved, and some numerical experiments confirm the theoretical results. Tests on unstructured meshes are also presented, showing the same convergence behavior.

Best Nonspherical Symmetric Low Rank Approximation

Mili I. Shah and Danny C. Sorensen

SIAM. J. Matrix Anal. & Appl. 31, pp. 1019-1039 (21 pages)

Online Publication Date: August 21, 2009

Full Text: | Download PDF

Show Abstract
The symmetry preserving singular value decomposition (SPSVD) produces the best symmetric (low rank) approximation to a set of data. These symmetric approximations are characterized via an invariance under the action of a symmetry group on the set of data. The symmetry groups of interest consist of all the nonspherical symmetry groups in three dimensions. This set includes the rotational, reflectional, dihedral, and inversion symmetry groups. In order to calculate the best symmetric (low rank) approximation, the symmetry of the data set must be determined. Therefore, matrix representations for each of the nonspherical symmetry groups have been formulated. These new matrix representations lead directly to a novel reweighting iterative method to determine the symmetry of a given data set by solving a series of minimization problems. Once the symmetry of the data set is found, the best symmetric (low rank) approximation in the Frobenius norm and matrix 2-norm can be established by using the SPSVD.

Computing Expected Transition Events in Reducible Markov Chains

Brian D. Ewald, Jeffrey Humpherys, and Jeremy M. West

SIAM. J. Matrix Anal. & Appl. 31, pp. 1040-1054 (15 pages)

Online Publication Date: August 21, 2009

Full Text: | Download PDF

Show Abstract
We present a closed-form, computable expression for the expected number of times any transition event occurs during the transient phase of a reducible Markov chain. Examples of events include time to absorption, number of visits to a state, traversals of a particular transition, loops from a state to itself, and arrivals to a state from a particular subset of states. We give an analogous expression for time-average events, which describe the steady-state behavior of reducible chains as well as the long-term behavior of irreducible chains.

Riemannian Metric and Geometric Mean for Positive Semidefinite Matrices of Fixed Rank

Silvère Bonnabel and Rodolphe Sepulchre

SIAM. J. Matrix Anal. & Appl. 31, pp. 1055-1070 (16 pages)

Online Publication Date: August 28, 2009

Full Text: | Download PDF

Show Abstract
This paper introduces a new metric and mean on the set of positive semidefinite matrices of fixed-rank. The proposed metric is derived from a well-chosen Riemannian quotient geometry that generalizes the reductive geometry of the positive cone and the associated natural metric. The resulting Riemannian space has strong geometrical properties: it is geodesically complete, and the metric is invariant with respect to all transformations that preserve angles (orthogonal transformations, scalings, and pseudoinversion). A meaningful approximation of the associated Riemannian distance is proposed, that can be efficiently numerically computed via a simple algorithm based on SVD. The induced mean preserves the rank, possesses the most desirable characteristics of a geometric mean, and is easy to compute.

Solving Ellipsoid-Constrained Integer Least Squares Problems

Xiao-Wen Chang and Gene H. Golub

SIAM. J. Matrix Anal. & Appl. 31, pp. 1071-1089 (19 pages)

Online Publication Date: August 28, 2009

Full Text: | Download PDF

Show Abstract
A new method is proposed to solve an ellipsoid-constrained integer least squares (EILS) problem arising in communications. In this method, the LLL reduction, which is cast as a QRZ factorization of a matrix, is used to transform the original EILS problem to a reduced EILS problem, and then a search algorithm is proposed to solve the reduced EILS problem. Simulation results indicate the new method can be much more computationally efficient than the existing method. The method is extended to solve a more general EILS problem.

Finding the Largest Eigenvalue of a Nonnegative Tensor

Michael Ng, Liqun Qi, and Guanglu Zhou

SIAM. J. Matrix Anal. & Appl. 31, pp. 1090-1099 (10 pages) | Cited 4 times

Online Publication Date: August 28, 2009

Full Text: | Download PDF

Show Abstract
In this paper we propose an iterative method for calculating the largest eigenvalue of an irreducible nonnegative tensor. This method is an extension of a method of Collatz (1942) for calculating the spectral radius of an irreducible nonnegative matrix. Numerical results show that our proposed method is promising. We also apply the method to studying higher-order Markov chains.

A Randomized Algorithm for Principal Component Analysis

Vladimir Rokhlin, Arthur Szlam, and Mark Tygert

SIAM. J. Matrix Anal. & Appl. 31, pp. 1100-1124 (25 pages) | Cited 4 times

Online Publication Date: August 28, 2009

Full Text: | Download PDF

Show Abstract
Principal component analysis (PCA) requires the computation of a low-rank approximation to a matrix containing the data being analyzed. In many applications of PCA, the best possible accuracy of any rank-deficient approximation is at most a few digits (measured in the spectral norm, relative to the spectral norm of the matrix being approximated). In such circumstances, efficient algorithms have not come with guarantees of good accuracy, unless one or both dimensions of the matrix being approximated are small. We describe an efficient algorithm for the low-rank approximation of matrices that produces accuracy that is very close to the best possible accuracy, for matrices of arbitrary sizes. We illustrate our theoretical results via several numerical examples.

The Exact Distribution of the Condition Number of a Gaussian Matrix

William Anderson and Martin T. Wells

SIAM. J. Matrix Anal. & Appl. 31, pp. 1125-1130 (6 pages)

Online Publication Date: September 18, 2009

Full Text: | Download PDF

Show Abstract
Suppose $G_{p\times n}$ is a real random matrix whose elements are independent and identically distributed standard normal random variables. Let $W_{p\times p}=G_{p\times n}^{}G_{n\times p}^{{\scriptscriptstyle\mathsf{T}}}$, which is the usual Wishart matrix. In addition, let $\lambda_{1}>\lambda_{2}>\cdots>\lambda_{p}>0$ and $\sigma_{1}>\sigma_{2}>\cdots>\sigma_{p}>0$ denote the distinct eigenvalues of the matrix $W_{p\times p}$ and singular values of $G_{p\times n}$, respectively. The 2-norm condition number of $G_{p\times n}$ is $\kappa_{2}(G_{p\times n})=\sqrt{\lambda_{1}/\lambda_{p}}=\sigma_{1}/\sigma_{p}$, the square root of the ratio of largest to smallest eigenvalues of the Wishart matrix. In this article we derive an exact expression, albeit somewhat complex, for the probability density function of $\kappa_{2}(G_{p\times n})$.

An Improved Arc Algorithm for Detecting Definite Hermitian Pairs

Chun-Hua Guo, Nicholas J. Higham, and Françoise Tisseur

SIAM. J. Matrix Anal. & Appl. 31, pp. 1131-1151 (21 pages)

Online Publication Date: September 18, 2009

Full Text: | Download PDF

Show Abstract
A 25-year-old and somewhat neglected algorithm of Crawford and Moon attempts to determine whether a given Hermitian matrix pair $(A,B)$ is definite by exploring the range of the function $f(x) = x^*(A+iB)x / | x^*(A+iB)x |$, which is a subset of the unit circle. We revisit the algorithm and show that with suitable modifications and careful attention to implementation details it provides a reliable and efficient means of testing definiteness. A clearer derivation of the basic algorithm is given that emphasizes an arc expansion viewpoint and makes no assumptions about the definiteness of the pair. Convergence of the algorithm is proved for all $(A,B$), definite or not. It is shown that proper handling of three details of the algorithm is crucial to the efficiency and reliability: how the midpoint of an arc is computed, whether shrinkage of an arc is permitted, and how directions of negative curvature are computed. For the latter, several variants of Cholesky factorization with complete pivoting are explored and the benefits of pivoting demonstrated. The overall cost of our improved algorithm is typically just a few Cholesky factorizations. Three applications of the algorithm are described: testing the hyperbolicity of a Hermitian quadratic matrix polynomial, constructing conjugate gradient methods for sparse linear systems in saddle point form, and computing the Crawford number of the pair $(A,B)$ via a quasiconvex univariate minimization problem.

Spectral Analysis of Saddle Point Matrices with Indefinite Leading Blocks

N. I. M. Gould and V. Simoncini

SIAM. J. Matrix Anal. & Appl. 31, pp. 1152-1171 (20 pages) | Cited 1 time

Online Publication Date: September 18, 2009

Full Text: | Download PDF

Show Abstract
We provide eigenvalue intervals for symmetric saddle point and regularized saddle point matrices in the case where the (1,1) block may be indefinite. These generalize known results for the definite (1,1) case. We also study the spectral properties of the equivalent augmented formulation, which is an alternative to explicitly dealing with the indefinite (1,1) block. Such an analysis may be used to assess the convergence of suitable Krylov subspace methods. We conclude with spectral analyses of the effects of common block-diagonal preconditioners.

Proximity of Weighted and Layered Least Squares Solutions

Tomonari Kitahara and Takashi Tsuchiya

SIAM. J. Matrix Anal. & Appl. 31, pp. 1172-1186 (15 pages)

Online Publication Date: September 18, 2009

Full Text: | Download PDF

Show Abstract
In this paper, we analyze the limiting behavior of the weighted least squares (WLS) problem $\min_{x\in\Re^n}\sum_{i=1}^p\|D_i(A_ix-b_i)\|^2$, where each $D_i$ is a positive definite diagonal matrix. We consider the situation where the magnitude of the weights differs drastically from one block to the next so that $\max(D_1)\geq\min(D_1)\gg\max(D_2)\geq\min(D_2)\gg\max(D_3)\geq\cdots\gg\max(D_{p-1})\geq\min(D_{p-1})\gg\max(D_p)$. Here $\max(\cdot)$ and $\min(\cdot)$ represent the maximum and minimum entries of diagonal elements, respectively. Specifically, we consider the case when the gap $g\equiv\min_i1/(\|D_i^{-1}\|\|D_{i+1}\|)$ is very large or tends to infinity. Vavasis and Ye proved that the limiting solution exists (when the proportion of diagonal elements within each block $D_i$ is unchanged and only the gap $g$ tends to $\infty$), and showed that the limit is characterized as the solution of a variant of the least squares problem called the layered least squares (LLS) problem. We analyze the difference between the solutions of WLS and LLS quantitatively and show that the norm of the difference of the two solutions is bounded above by $O(\chi_A\bar{\chi}_A^{2(p+1)}g^{-2}\|b\|)$ and $O(\bar{\chi}_A^{2p+3}g^{-2}\|b\|)$ in the variable and the residual spaces, respectively, using the two condition numbers $\chi_A\equiv\max_{B\in{\cal B}}\|B^{-1}\|$ and $\bar{\chi}_A\equiv\max_{B\in{\cal B}}\|AB^{-1}\|$ of $A$, where ${\cal B}$ is the set of all nonsingular $n\times n$ submatrices of $A$, $A=[A_1;\ldots;A_p]$ and $b = [b_1;\ldots;b_p]$. A remarkable feature of this result is that the error bound is represented in terms of $A$, $g$ (and $b$) and independent of the weights $D_i$, $i=1,\ldots,p$. The analysis is carried out by making the change of variables to convert the matrix $A$ into a basis lower-triangular form and then by applying the Sherman–Morrison–Woodbury formula.

Inertia and Rank Characterizations of Some Matrix Expressions

Delin Chu, Y. S. Hung, and Hugo J. Woerdeman

SIAM. J. Matrix Anal. & Appl. 31, pp. 1187-1226 (40 pages)

Online Publication Date: October 29, 2009

Full Text: | Download PDF

Show Abstract
In this paper we consider the admissible inertias and ranks of the expressions $A-BXB^*-CYC^*$ and $A-BXC^*\pm CX^*B^*$ with unknowns $X$ and $Y$ in the four cases when these expressions are: (i) complex self-adjoint, (ii) complex skew-adjoint, (iii) real symmetric, (iv) real skew symmetric. We also provide a construction for $X$ and $Y$ to achieve the desired inertia/rank that uses only unitary/orthogonal transformation, thus leading to a numerically reliable construction. In addition, we look at related block matrix completion problems $\left[\begin{array}{@{}ccc@{}} {\cal A} & {\cal B} & {\cal C} \ \pm {\cal B}^\star & {\cal X} & {\cal E} \\\pm {\cal C}^\star & \pm {\cal E}^\star & {\cal Y} \end{array}\right]$ with either two diagonal unknown blocks and $\left[\begin{array}{@{}ccc@{}} {\cal A} & {\cal B} & {\cal X} \ \pm {\cal B}^\star & {\cal D} & {\cal C} \ \pm {\cal X}^\star & \pm {\cal C}^\star & {\cal E} \end{array}\right]$ with an unknown off-diagonal block. Finally, we also provide all admissible ranks in the case when we drop any adjointness/symmetry constraint.

On a Generalization of Soules Bases

S. D. Eubanks and J. J. McDonald

SIAM. J. Matrix Anal. & Appl. 31, pp. 1227-1234 (8 pages)

Online Publication Date: October 29, 2009

Full Text: | Download PDF

Show Abstract
We characterize ordered orthonormal bases $\{s_1,\dots,s_n\} \subset \mathbb{R}^n$ for which the matrix $s_1s_1^T+\dots+s_ts_t^T$ is nonnegative for each $t=1,\dots,n$. In particular we show that the associated orthogonal matrix $S=[s_1,s_2,\dots,s_n]$ has certain submatrices that are Soules matrices. Based on the construction of Soules bases developed by Elsner, Nabben, and Neumann [Linear Algebra Appl., 271 (1998), pp. 323–343] we are able to construct all possible such bases. Additionally, we show that this set is the closure of the set of Soules bases. We also include some results regarding matrix functions of $S \Lambda S^T$, where $\Lambda$ is a diagonal matrix with nonincreasing diagonal entries.

Interior-Point Method for Nuclear Norm Approximation with Application to System Identification

Zhang Liu and Lieven Vandenberghe

SIAM. J. Matrix Anal. & Appl. 31, pp. 1235-1256 (22 pages) | Cited 3 times

Online Publication Date: November 04, 2009

Full Text: | Download PDF

Show Abstract
The nuclear norm (sum of singular values) of a matrix is often used in convex heuristics for rank minimization problems in control, signal processing, and statistics. Such heuristics can be viewed as extensions of $\ell_1$-norm minimization techniques for cardinality minimization and sparse signal estimation. In this paper we consider the problem of minimizing the nuclear norm of an affine matrix-valued function. This problem can be formulated as a semidefinite program, but the reformulation requires large auxiliary matrix variables, and is expensive to solve by general-purpose interior-point solvers. We show that problem structure in the semidefinite programming formulation can be exploited to develop more efficient implementations of interior-point methods. In the fast implementation, the cost per iteration is reduced to a quartic function of the problem dimensions and is comparable to the cost of solving the approximation problem in the Frobenius norm. In the second part of the paper, the nuclear norm approximation algorithm is applied to system identification. A variant of a simple subspace algorithm is presented in which low-rank matrix approximations are computed via nuclear norm minimization instead of the singular value decomposition. This has the important advantage of preserving linear matrix structure in the low-rank approximation. The method is shown to perform well on publicly available benchmark data.

General Matrix Pencil Techniques for Solving Discrete-Time Nonsymmetric Algebraic Riccati Equations

Marc Jungers, Cristian Oară, Hisham Abou-Kandil, and Radu Ştefan

SIAM. J. Matrix Anal. & Appl. 31, pp. 1257-1278 (22 pages)

Online Publication Date: November 04, 2009

Full Text: | Download PDF

Show Abstract
A discrete-time nonsymmetric algebraic Riccati system which incorporates as special cases of various discrete-time nonsymmetric algebraic Riccati equations is introduced and studied without any restrictive assumptions on the matrix coefficients. Necessary and sufficient existence conditions together with computable formulas for the stabilizing solution are given in terms of proper deflating subspaces of an associated matrix pencil. The theory is applied in the framework of game theory with an open-loop information structure to design Nash strategy without the classical assumptions on the invertibility of some matrix coefficients.

Verified Computation of Square Roots of a Matrix

Andreas Frommer and Behnam Hashemi

SIAM. J. Matrix Anal. & Appl. 31, pp. 1279-1302 (24 pages)

Online Publication Date: November 04, 2009

Full Text: | Download PDF

Show Abstract
We present methods to compute verified square roots of a square matrix $A$. Given an approximation $X$ to the square root, obtained by a classical floating point algorithm, we use interval arithmetic to find an interval matrix which is guaranteed to contain the error of $X$. Our approach is based on the Krawczyk method, which we modify in two different ways in such a manner that the computational complexity for an $n\times n$ matrix is reduced to $n^3$. The methods are based on the spectral decomposition or, in the case that the eigenvector matrix is ill conditioned, on a similarity transformation to block diagonal form. Numerical experiments prove that our methods are computationally efficient and that they yield narrow enclosures provided $X$ is a good approximation. This is particularly true for symmetric matrices, since their eigenvector matrix is perfectly conditioned.

Why Lagrange Multipliers with Extreme Magnitudes Give Extrema of Definite Hermitian Forms on Quadric Surfaces

Yves Nievergelt

SIAM. J. Matrix Anal. & Appl. 31, pp. 1303-1328 (26 pages)

Online Publication Date: November 13, 2009

Full Text: | Download PDF

Show Abstract
This work explains the geometry of analytic theorems by Forsythe and Golub, Gander, and Moré to the effect that the constrained minimum (respectively, maximum) of a positive definite Hermitian form on a level set of any nonnecessarily definite Hermitian polynomial corresponds to the Lagrange multiplier with the smallest (respectively, largest) absolute value. Locally, the law of sines for the triangle with vertices at the center and at two stationary points reveals that the objective values are in the same order as the magnitudes of the Lagrange multipiers. Global geometry explains the same results for global minima and global maxima by showing that the constraining quadric surface and the Lagrange multiplier form a set of geodetic coordinates for the entire ambient space. Duality and the symmetry of the constraining quadratic hypersurface also explain why the difference between two stationary points with the same objective value is a generalized eigenvector.

A Global Convergence Proof for Cyclic Jacobi Methods with Block Rotations

Zlatko Drmač

SIAM. J. Matrix Anal. & Appl. 31, pp. 1329-1350 (22 pages)

Online Publication Date: November 13, 2009

Full Text: | Download PDF

Show Abstract
This paper introduces a globally convergent block (column- and row-) cyclic Jacobi method for diagonalization of Hermitian matrices and for computation of the singular value decomposition of general matrices. It is shown that a block rotation (a generalization of the Jacobi $2\times2$ rotation) can be computed and implemented in a particular way to guarantee global convergence. The proof includes the convergence of the eigenspaces in the general case of multiple eigenvalues. This solves a long standing open problem of convergence of block cyclic Jacobi methods.

A Matrix Computation View of FastMap and RobustMap Dimension Reduction Algorithms

George Ostrouchov

SIAM. J. Matrix Anal. & Appl. 31, pp. 1351-1360 (10 pages)

Online Publication Date: December 04, 2009

Full Text: | Download PDF

Show Abstract
Given a set of pairwise object distances and a dimension $k$, FastMap and RobustMap algorithms compute a set of $k$-dimensional coordinates for the objects. These metric space embedding methods implicitly assume a higher-dimensional coordinate representation and are a sequence of translations and orthogonal projections based on a sequence of object pair selections (called pivot pairs). We develop a matrix computation viewpoint of these algorithms that operates on the coordinate representation explicitly using Householder reflections. The resulting coordinate mapping algorithm (CMA) is a fast approximate alternative to truncated principal component analysis (PCA), and it brings the FastMap and RobustMap algorithms into the mainstream of numerical computation where standard BLAS building blocks are used. Motivated by the geometric nature of the embedding methods, we further show that truncated PCA can be computed with CMA by specific pivot pair selections. Describing FastMap, RobustMap, and PCA as CMA computations with different pivot pair choices unifies the methods along a pivot pair selection spectrum. We also sketch connections to the semidiscrete decomposition and the QLP decomposition.

Continuation of Invariant Subspaces for Parameterized Quadratic Eigenvalue Problems

Wolf Jürgen Beyn and Vera Thümmler

SIAM. J. Matrix Anal. & Appl. 31, pp. 1361-1381 (21 pages) | Cited 1 time

Online Publication Date: December 04, 2009

Full Text: | Download PDF

Show Abstract
We consider quadratic eigenvalue problems with large and sparse matrices depending on a parameter. Problems of this type occur, for example, in the stability analysis of spatially discretized and parameterized nonlinear wave equations. The aim of the paper is to present and analyze a continuation method for invariant subspaces that belong to a group of eigenvalues, the number of which is much smaller than the dimension of the system. The continuation method is of predictor-corrector type, similar to the approach for the linear eigenvalue problem in [Beyn, Kleß, and Thümmler, Ergodic Theory, Analysis, and Efficient Simulation of Dynamical Systems, Springer, Berlin, 2001], but we avoid linearizing the problem, which will double the dimension and change the sparsity pattern. The matrix equations that occur in the predictor and the corrector step are solved by a bordered version of the Bartels–Stewart algorithm. Furthermore, we set up an update procedure that handles the transition from real to complex conjugate eigenvalues, which occurs when eigenvalues from inside the continued cluster collide with eigenvalues from outside. The method is demonstrated on several numerical examples: a homotopy between random matrices, a fluid conveying pipe problem, and a traveling wave of a damped wave equation.

Superfast Multifrontal Method for Large Structured Linear Systems of Equations

Jianlin Xia, Shivkumar Chandrasekaran, Ming Gu, and Xiaoye S. Li

SIAM. J. Matrix Anal. & Appl. 31, pp. 1382-1411 (30 pages) | Cited 8 times

Online Publication Date: December 04, 2009

Full Text: | Download PDF

Show Abstract
In this paper we develop a fast direct solver for large discretized linear systems using the supernodal multifrontal method together with low-rank approximations. For linear systems arising from certain partial differential equations such as elliptic equations, during the Gaussian elimination of the matrices with proper ordering, the fill-in has a low-rank property: all off-diagonal blocks have small numerical ranks with proper definition of off-diagonal blocks. Matrices with this low-rank property can be efficiently approximated with semiseparable structures called hierarchically semiseparable (HSS) representations. We reveal the above low-rank property by ordering the variables with nested dissection and eliminating them with the multifrontal method. All matrix operations in the multifrontal method are performed in HSS forms. We present efficient ways to organize the HSS structured operations along the elimination. Some fast HSS matrix operations using tree structures are proposed. This new structured multifrontal method has nearly linear complexity and a linear storage requirement. Thus, we call it a superfast multifrontal method. It is especially suitable for large sparse problems and also has natural adaptability to parallel computations and great potential to provide effective preconditioners. Numerical results demonstrate the efficiency.

Reducible Spectral Theory with Applications to the Robustness of Matrices in Max-Algebra

P. Butkovič, R. A. Cuninghame-Green, and S. Gaubert

SIAM. J. Matrix Anal. & Appl. 31, pp. 1412-1431 (20 pages)

Online Publication Date: December 04, 2009

Full Text: | Download PDF

Show Abstract
Let $a\oplus b=\max(a,b)$ and $a\otimes b=a+b$ for $a,b\in\overline{\mathbb{R}}:=\mathbb{R}\cup\{-\infty\}$. By max-algebra we understand the analogue of linear algebra developed for the pair of operations $(\oplus,\otimes)$, extended to matrices and vectors. The symbol $A^{k}$ stands for the $k$th max-algebraic power of a square matrix $A$. Let us denote by $\varepsilon$ the max-algebraic “zero” vector, all the components of which are $-\infty$. The max-algebraic eigenvalue-eigenvector problem is the following: Given $A\in\overline{\mathbb{R}}^{n\times n}$, find all $\lambda\in\overline{\mathbb{R}}$ and $x\in\overline{\mathbb{R}}^{n}$, $x\neq\varepsilon$, such that $A\otimes x=\lambda\otimes x$. Certain problems of scheduling lead to the following question: Given $A\in\overline{\mathbb{R}}^{n\times n}$, is there a $k$ such that $A^{k}\otimes x$ is a max-algebraic eigenvector of $A$? If the answer is affirmative for every $x\neq\varepsilon$, then $A$ is called robust. First, we give a complete account of the reducible max-algebraic spectral theory, and then we apply it to characterize robust matrices.

Calibrating Least Squares Semidefinite Programming with Equality and Inequality Constraints

Yan Gao and Defeng Sun

SIAM. J. Matrix Anal. & Appl. 31, pp. 1432-1457 (26 pages) | Cited 3 times

Online Publication Date: December 18, 2009

Full Text: | Download PDF

Show Abstract
In this paper, we consider the least squares semidefinite programming with a large number of equality and inequality constraints. One difficulty in finding an efficient method for solving this problem is due to the presence of the inequality constraints. In this paper, we propose to overcome this difficulty by reformulating the problem as a system of semismooth equations with two level metric projection operators. We then design an inexact smoothing Newton method to solve the resulting semismooth system. At each iteration, we use the BiCGStab iterative solver to obtain an approximate solution to the generated smoothing Newton linear system. Our numerical experiments confirm the high efficiency of the proposed method.

On the Semidefinite B-Arnoldi Method

G. W. Stewart

SIAM. J. Matrix Anal. & Appl. 31, pp. 1458-1468 (11 pages)

Online Publication Date: December 18, 2009

Full Text: | Download PDF

Show Abstract
The B-Arnoldi method is a variant of the ordinary Arnoldi method in which orthogonalization is done with respect to the inner product generated by a positive definite matrix $B$. It arises in connection with the generalized eigenvalue problem $Ax = \lambda Bx$. When $B$ is semidefinite, the algorithm can proceed formally, with “orthogonalization” taking place in the semi-inner product generated by $B$. However, it has been observed that components of the Arnoldi vectors lying in the null space of $B$ can grow rapidly. In this paper we examine the source and consequences of this growth.

Uniqueness Conditions for Constrained Three-Way Factor Decompositions with Linearly Dependent Loadings

Alwin Stegeman and André L. F. de Almeida

SIAM. J. Matrix Anal. & Appl. 31, pp. 1469-1490 (22 pages) | Cited 2 times

Online Publication Date: December 23, 2009

Full Text: | Download PDF

Show Abstract
In this paper, we derive uniqueness conditions for a constrained version of the parallel factor (Parafac) decomposition, also known as canonical decomposition (Candecomp). Candecomp/Parafac (CP) decomposes a three-way array into a prespecified number of outer product arrays. The constraint is that some vectors forming the outer product arrays are linearly dependent according to a prespecified pattern. This is known as the PARALIND family of models. An important subclass is where some vectors forming the outer product arrays are repeated according to a prespecified pattern. These are known as CONFAC decompositions. We discuss the relation between PARALIND, CONFAC, and the three-way decompositions CP, Tucker3, and the decomposition in block terms. We provide both essential uniqueness conditions and partial uniqueness conditions for PARALIND and CONFAC and discuss the relation with uniqueness of constrained Tucker3 models and the block decomposition in rank-$(L,L,1)$ terms. Our results are demonstrated by means of examples.

Convexity Properties of the Condition Number

Carlos Beltrán, Jean-Pierre Dedieu, Gregorio Malajovich, and Mike Shub

SIAM. J. Matrix Anal. & Appl. 31, pp. 1491-1506 (16 pages) | Cited 1 time

Online Publication Date: January 06, 2010

Full Text: | Download PDF

Show Abstract
We define in the space of $n\times m$ matrices of rank $n$, $n\leq m$, the condition Riemannian structure as follows: For a given matrix $A$ the tangent space at $A$ is equipped with the Hermitian inner product obtained by multiplying the usual Frobenius inner product by the inverse of the square of the smallest singular value of $A$ denoted $\sigma_n(A)$. When this smallest singular value has multiplicity 1, the function $A\rightarrow\log(\sigma_n(A)^{-2})$ is a convex function with respect to the condition Riemannian structure that is $t\rightarrow\log(\sigma_n(A(t))^{-2})$ is convex, in the usual sense for any geodesic $A(t)$. In a more abstract setting, a function $\alpha$ defined on a Riemannian manifold $(\mathcal{M},\langle,\rangle)$ is said to be self-convex when $\log\alpha(\gamma(t))$ is convex for any geodesic in $(\mathcal{M},\alpha\,\langle,\rangle)$. Necessary and sufficient conditions for self-convexity are given when $\alpha$ is $C^2$. When $\alpha(x)=d(x,\mathcal{N})^{-2}$, where $d(x,\mathcal{N})$ is the distance from $x$ to a $C^2$ submanifold $\mathcal{N}\subset\mathbb{R}^j$, we prove that $\alpha$ is self-convex when restricted to the largest open set of points $x$ where there is a unique closest point in $\mathcal{N}$ to $x$. We also show, using this more general notion, that the square of the condition number $\|A\|_F\,/\sigma_n(A)$ is self-convex in projective space and the solution variety.

The Stable Perturbation of the Drazin Inverse of the Square Matrices

Qingxiang Xu, Chuanning Song, and Yimin Wei

SIAM. J. Matrix Anal. & Appl. 31, pp. 1507-1520 (14 pages)

Online Publication Date: January 08, 2010

Full Text: | Download PDF

Show Abstract
For any $n\times n$ complex matrix $A$, let $A^D$ and $A^\pi$ be the Drazin inverse and the spectral projector of $A$, respectively, where $A^\pi=I-AA^D$. When $A$ is singular, an $n\times n$ complex matrix $B$ is said to be a stable perturbation of $A$ if $I-A^\pi-B^\pi$ is nonsingular or, equivalently, if the matrix $B$ satisfies condition (${\cal C}_s$) recently introduced by Castro-González, Robles, and Vélez-Cerrada [SIAM J. Matrix Anal. Appl., 30 (2008), pp. 882–897]. In the perturbation analysis of the Drazin inverse, the condition of $\Vert B-A\Vert$ being small is usually implicitly assumed in the literature. In this case, the condition of $B$ being a stable perturbation of $A$ is necessary in order to ensure the continuity of the Drazin inverse. In this paper, only under the condition that $B$ is a stable perturbation of $A$, two explicit formulas for the Drazin inverse $B^D$ and the spectral projector $B^\pi$ are provided, respectively, and some upper bounds for $\Vert B^D-A^D\Vert/\Vert A^D\Vert$ and $\Vert B^\pi-A^\pi\Vert$ are derived from these formulas under certain conditions. In the case where $\Vert B-A\Vert$ is small, numerical examples are given indicating the sharpness of these norm upper bounds. Furthermore, a numerical example is provided to illustrate that a perturbation analysis of the Drazin inverse can also be conducted, even if $\Vert B-A\Vert$ is not small.

Rayleigh–Ritz Majorization Error Bounds with Applications to FEM

Andrew V. Knyazev and Merico E. Argentati

SIAM. J. Matrix Anal. & Appl. 31, pp. 1521-1537 (17 pages) | Cited 1 time

Online Publication Date: January 08, 2010

Full Text: | Download PDF

Show Abstract
The Rayleigh–Ritz (RR) method finds the stationary values, called Ritz values, of the Rayleigh quotient on a given trial subspace as approximations to eigenvalues of a Hermitian operator $A$. If the trial subspace is $A$-invariant, the Ritz values are exactly some of the eigenvalues of $A$. Given two subspaces $\mathcal{X}$ and $\mathcal{Y}$ of the same finite dimension, such that $\mathcal{X}$ is $A$-invariant, the absolute changes in the Ritz values of $A$ with respect to $\mathcal{X}$ compared to the Ritz values with respect to $\mathcal{Y}$ represent the RR absolute eigenvalue approximation error. Our first main result is a sharp majorization-type RR error bound in terms of the principal angles between $\mathcal{X}$ and $\mathcal{Y}$ for an arbitrary $A$-invariant $\mathcal{X}$, which was a conjecture in [SIAM J. Matrix Anal. Appl., 30 (2008), pp. 548–559]. Second, we prove a novel type of RR error bound that deals with the products of the errors, rather than the sums. Third, we establish majorization bounds for the relative errors. We extend our bounds to the case $\dim\mathcal{X}\leq\dim\mathcal{Y}<\infty$ in Hilbert spaces and apply them in the context of the finite element method.

Note on “A New Scaling for Newton's Iteration for the Polar Decomposition and Its Backward Stability” by R. Byers and H. Xu

Andrzej KieŁbasiński and Krystyna Ziętak

SIAM. J. Matrix Anal. & Appl. 31, pp. 1538-1539 (2 pages) | Cited 1 time

Online Publication Date: January 08, 2010

Full Text: | Download PDF

Show Abstract
We show that the proof of Theorem 4.5 in the paper by Byers and Xu is not complete for the considered class of matrices. The first order analysis is not applicable to the considered problem of rounding error propagation.
Close

Your Recent History Recent History Help

Recently Viewed

  1. Splitting a Matrix of Laurent Polynomials with Symmetry and itsApplication to Symmetric Framelet Fil...
  2. Distribution of Entries in a Substochastic Matrix Having Eigenvalues Near 1
  3. Logarithmic Norms for Matrix Pencils
  4. Radiative Decay of Bubble Oscillations in a Compressible Fluid
  5. Dynamics of Compressible Non-isothermal Fluids of Non-Newtonian Korteweg Type
© SIAM By using SIAM Publications Online you agree to abide by the Terms and Conditions of Use. Banner art adapted from a figure by Hinke M. Osinga and Bernd Krauskopf (University of Bristol, UK).
close