SIAM Digital Library
 
 
 

SIAM J. on Matrix Analysis and Applications

Year Range: 

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

Search Issue | RSS Feeds RSS
Previous Issue

2010

Volume 31, Issue 5, pp. 2239-2996


Quasi-Birth-and-Death Processes and Matrix-Valued Orthogonal Polynomials

Aubrey Clayton

SIAM. J. Matrix Anal. & Appl. 31, pp. 2239-2260 (22 pages) | Cited 1 time

Online Publication Date: July 01, 2010

Full Text: | Download PDF

Show Abstract
We consider a matrix-valued spectral decomposition of a family of block-tridiagonal matrices arising as the transition matrices of so-called quasi-birth-and-death processes. This representation is a generalization of the Karlin–McGregor representation for the $n$-step transition probabilities of a birth-and-death process via a system of orthogonal polynomials. At the heart of the representation is a self-adjoint matrix-valued measure associated to the process. We make use of a previously known formula relating the Stieltjes transform of this measure to that of the measure associated to the “0th associated process,” generalizing a theorem of Karlin and McGregor, to compute the Stieltjes transform of the spectral measure for several examples. In addition, we apply matrix-valued orthogonal polynomial techniques to the study of “sin-graphs” and higher-dimensional birth-and-death processes, for which the relevant polynomials are multivariate.

On the Numerical Rank of the Off-Diagonal Blocks of Schur Complements of Discretized Elliptic PDEs

S. Chandrasekaran, P. Dewilde, M. Gu, and N. Somasunderam

SIAM. J. Matrix Anal. & Appl. 31, pp. 2261-2290 (30 pages) | Cited 1 time

Online Publication Date: July 01, 2010

Full Text: | Download PDF

Show Abstract
It is shown that the numerical rank of the off-diagonal blocks of certain Schur complements of matrices that arise from the finite-difference discretization of constant coefficient, elliptic PDEs in two spatial dimensions is bounded by a constant independent of the grid size. Moreover, in three-dimensional problems the Schur complements are shown to have off-diagonal blocks whose numerical rank is a slowly growing function.

Preconditioned Conjugate Gradient Method for Optimal Control Problems with Control and State Constraints

Roland Herzog and Ekkehard Sachs

SIAM. J. Matrix Anal. & Appl. 31, pp. 2291-2317 (27 pages) | Cited 1 time

Online Publication Date: July 01, 2010

Full Text: | Download PDF

Show Abstract
Optimality systems and their linearizations arising in optimal control of partial differential equations with pointwise control and (regularized) state constraints are considered. The preconditioned conjugate gradient (PCG) method in a nonstandard inner product is employed for their efficient solution. Preconditioned condition numbers are estimated for problems with pointwise control constraints, mixed control-state constraints, and of Moreau–Yosida penalty type. Numerical results for elliptic problems demonstrate the performance of the PCG iteration. Regularized state-constrained problems in three dimensions with more than 750,000 variables are solved.

Combinatorics of Column Minimal Indices and Matrix Pencil Completion Problems

Marija Dodig and Marko Stošić

SIAM. J. Matrix Anal. & Appl. 31, pp. 2318-2346 (29 pages)

Online Publication Date: July 13, 2010

Full Text: | Download PDF

Show Abstract
In this paper we describe the possible Kronecker invariants of a matrix pencil with a prescribed subpencil, when the prescribed subpencil has only column (row) minimal indices as Kronecker invariants. This is done by exploiting combinatorics of column minimal indices and by introducing some novel labels. The result is given over arbitrary fields.

An Augmented Stability Result for the Lanczos Hermitian Matrix Tridiagonalization Process

Christopher C. Paige

SIAM. J. Matrix Anal. & Appl. 31, pp. 2347-2359 (13 pages)

Online Publication Date: July 13, 2010

Full Text: | Download PDF

Show Abstract
It is shown that a good implementation of the Hermitian matrix tridiagonalization process of Lanczos [J. Research Nat. Bur. Standards, 45 (1950), pp. 255–282] produces a tridiagonal matrix that is, at each step, the exact result for the process applied to a strange augmented problem. Since the process is not stable in the standard sense, this augmented stability result cannot be transformed to prove standard stability. The intent is to obtain an increased understanding of the Lanczos tridiagonalization process, and this result could later be used to analyze the many applications of the process to large sparse matrix problems, such as the solution of the eigenproblem, compatible linear systems, least squares, and the singular value decomposition.

Dynamical Tensor Approximation

Othmar Koch and Christian Lubich

SIAM. J. Matrix Anal. & Appl. 31, pp. 2360-2375 (16 pages) | Cited 2 times

Online Publication Date: July 15, 2010

Full Text: | Download PDF

Show Abstract
For the approximation of time-dependent data tensors and of solutions to tensor differential equations by tensors of low Tucker rank, we study a computational approach that can be viewed as a continuous-time updating procedure. This approach works with the increments rather than the full tensor and avoids the computation of decompositions of large matrices. In this method, the derivative is projected onto the tangent space of the manifold of tensors of Tucker rank $(r_1,\dots,r_N)$ at the current approximation. This yields nonlinear differential equations for the factors in a Tucker decomposition, suitable for numerical integration. Approximation properties of this approach are analyzed.

A Fast Algorithm for Updating and Downsizing the Dominant Kernel Principal Components

Nicola Mastronardi, Eugene E. Tyrtyshnikov, and Paul Van Dooren

SIAM. J. Matrix Anal. & Appl. 31, pp. 2376-2399 (24 pages)

Online Publication Date: July 22, 2010

Full Text: | Download PDF

Show Abstract
Many important kernel methods in the machine learning area, such as kernel principal component analysis, feature approximation, denoising, compression, and prediction require the computation of the dominant set of eigenvectors of the symmetric kernel Gram matrix. Recently, an efficient incremental approach was presented for the fast calculation of the dominant kernel eigenbasis. In this paper we propose faster algorithms for incrementally updating and downsizing the dominant kernel eigenbasis. These methods are well-suited for large scale problems since they are efficient in terms of both complexity and data management.

GMRES Methods for Least Squares Problems

Ken Hayami, Jun-Feng Yin, and Tokushi Ito

SIAM. J. Matrix Anal. & Appl. 31, pp. 2400-2430 (31 pages) | Cited 1 time

Online Publication Date: August 12, 2010

Full Text: | Download PDF

Show Abstract
The standard iterative method for solving large sparse least squares problems $\min\|\mbox{\boldmath$b$}-A\mbox{\boldmath$x$}\|_2$, $A\in\mathbf{R}^{m\times n}$, is the CGLS method, or its stabilized version, LSQR, which is mathematically equivalent to applying the conjugate gradient method to the normal equation $A^{\mbox{\tiny T}}A\mbox{\boldmath$x$}=A^{\mbox{\tiny T}}\mbox{\boldmath$b$}$. We consider alternative methods using a matrix $B\in\mathbf{R}^{n\times m}$ and applying the generalized minimal residual (GMRES) method to $\min\|\mbox{\boldmath$b$}-AB\mbox{\boldmath$z$}\|_2$ or $\min\|B\mbox{\boldmath$b$}-BA\mbox{\boldmath$x$}\|_2$. We give a sufficient condition concerning $B$ for the GMRES methods to give a least squares solution without breakdown for arbitrary $\mbox{\boldmath$b$}$, for overdetermined, underdetermined, and possibly rank-deficient problems. We then give a convergence analysis of the GMRES methods as well as the CGLS method. Then, we propose using the robust incomplete factorization (RIF) for $B$. Finally, we show by numerical experiments on overdetermined and underdetermined problems that, for ill-conditioned problems, the GMRES methods with RIF give least squares solutions faster than the CGLS and LSQR methods with RIF, and are similar in performance to the reorthogonalized CGLS with RIF.

Improved Balanced Incomplete Factorization

Rafael Bru, José Marín, José Mas, and Miroslav Tůma

SIAM. J. Matrix Anal. & Appl. 31, pp. 2431-2452 (22 pages)

Online Publication Date: August 12, 2010

Full Text: | Download PDF

Show Abstract
In this paper we improve the BIF algorithm which computes simultaneously the LU factors (direct factors) of a given matrix and their inverses (inverse factors). This algorithm was introduced in [R. Bru, J. Marín, J. Mas, and M. Tůma, SIAM J. Sci. Comput., 30 (2008), pp. 2302–2318]. The improvements are based on a deeper understanding of the inverse Sherman–Morrison (ISM) decomposition, and they provide a new insight into the BIF decomposition. In particular, it is shown that a slight algorithmic reformulation of the basic algorithm implies that the direct and inverse factors numerically influence each other even without any dropping for incompleteness. Algorithmically, the nonsymmetric version of the improved BIF algorithm is formulated. Numerical experiments show very high robustness of the incomplete implementation of the algorithm used for preconditioning nonsymmetric linear systems.

Eigenvalue Estimates for Preconditioned Nonsymmetric Saddle Point Matrices

Shu-Qian Shen, Ting-Zhu Huang, and Juan Yu

SIAM. J. Matrix Anal. & Appl. 31, pp. 2453-2476 (24 pages)

Online Publication Date: August 12, 2010

Full Text: | Download PDF

Show Abstract
We present a valid lower bound for the real eigenvalues of a nonsymmetric saddle point matrix even if its $(2,2)$ block is singular. This bound, together with bounds introduced by Benzi and Simoncini in [Numer. Math., 103 (2006), pp. 173–196], is applied for the analysis of symmetric indefinite preconditioners and primal-based penalty preconditioners. Eigenvalue estimates for the Hermitian and skew-Hermitian splitting preconditioned saddle point matrices are also derived. The model problem of Stokes equations is used to illustrate the presented theoretical results.

Accurate Matrix Factorization: Inverse LU and Inverse QR Factorizations

Takeshi Ogita

SIAM. J. Matrix Anal. & Appl. 31, pp. 2477-2497 (21 pages)

Online Publication Date: August 17, 2010

Full Text: | Download PDF

Show Abstract
In this paper, algorithms for accurate matrix factorizations named inverse LU and inverse QR factorizations for extremely ill-conditioned matrices are proposed. The proposed algorithms are based on standard numerical algorithms using pure floating-point arithmetic and accurate dot product. Detailed analysis of the algorithms is presented. As an application of the proposed algorithms, a method of computing accurate solutions of linear systems is also proposed. Numerical results are presented for illustrating the performance of the proposed algorithms. Computing times for the algorithms adaptively change according to the difficulty of given problems.

On Uniqueness of the $n$th Order Tensor Decomposition into Rank-1 Terms with Linear Independence in One Mode

Alwin Stegeman

SIAM. J. Matrix Anal. & Appl. 31, pp. 2498-2516 (19 pages) | Cited 3 times

Online Publication Date: August 17, 2010

Full Text: | Download PDF

Show Abstract
We study uniqueness of the decomposition of an $n$th order tensor (also called $n$-way array) into a sum of $R$ rank-1 terms (where each term is the outer product of $n$ vectors). This decomposition is also known as Parafac or Candecomp, and a general uniqueness condition for $n=3$ has been obtained by Kruskal in 1977 [Linear Algebra Appl., 18 (1977), pp. 95–138]. More recently, Kruskal's uniqueness condition has been generalized to $n\geq3$, and less restrictive uniqueness conditions have been obtained for the case where the vectors of the rank-1 terms are linearly independent in (at least) one of the $n$ modes. For this case, only $n=3$ and $n=4$ have been studied. We generalize these results by providing a framework of analysis for arbitrary $n\geq3$. Our results include necessary, sufficient, necessary and sufficient, and generic uniqueness conditions. For the sufficient uniqueness conditions, the rank of a matrix needs to be checked. The generic uniqueness conditions have the form of a bound on $R$ in terms of the dimensions of the tensor to be decomposed.

Further Results for Perron–Frobenius Theorem for Nonnegative Tensors

Yuning Yang and Qingzhi Yang

SIAM. J. Matrix Anal. & Appl. 31, pp. 2517-2530 (14 pages) | Cited 1 time

Online Publication Date: August 31, 2010

Full Text: | Download PDF

Show Abstract
We give further results on the Perron–Frobenius theorem for tensors, generalize other theorems from matrices to tensors, and give an equivalent condition for nonnegative irreducible tensors.

Stability of the Levinson Algorithm for Toeplitz-Like Systems

P. Favati, G. Lotti, and O. Menchi

SIAM. J. Matrix Anal. & Appl. 31, pp. 2531-2552 (22 pages)

Online Publication Date: August 31, 2010

Full Text: | Download PDF

Show Abstract
Numerical stability of the Levinson algorithm, generalized for Toeplitz-like systems, is studied. Arguments based on the analytic results of an error analysis for floating point arithmetic produce an upper bound on the norm of the residual vector, which grows exponentially with respect to the size of the problem. The base of such an exponential function can be small for diagonally dominant Toeplitz-like matrices. Numerical experiments show that, for these matrices, Gaussian elimination by row and the Levinson algorithm have residuals of the same order of magnitude. As expected, the empirical results point out that the theoretical bound is too pessimistic.

A Riemannian Optimization Approach for Computing Low-Rank Solutions of Lyapunov Equations

Bart Vandereycken and Stefan Vandewalle

SIAM. J. Matrix Anal. & Appl. 31, pp. 2553-2579 (27 pages)

Online Publication Date: September 09, 2010

Full Text: | Download PDF

Show Abstract
We propose a new framework based on optimization on manifolds to approximate the solution of a Lyapunov matrix equation by a low-rank matrix. The method minimizes the error on the Riemannian manifold of symmetric positive semidefinite matrices of fixed rank. We detail how objects from differential geometry, like the Riemannian gradient and Hessian, can be efficiently computed for this manifold. As a minimization algorithm we use the Riemannian trust-region method of [P.-A. Absil, C. Baker, and K. Gallivan, Found. Comput. Math., 7 (2007), pp. 303–330] based on a second-order model of the objective function on the manifold. Together with an efficient preconditioner, this method can find low-rank solutions with very little memory. We illustrate our results with numerical examples.

The Condition Metric in the Space of Rectangular Full Rank Matrices

Paola Boito and Jean-Pierre Dedieu

SIAM. J. Matrix Anal. & Appl. 31, pp. 2580-2602 (23 pages)

Online Publication Date: September 14, 2010

Full Text: | Download PDF

Show Abstract
The condition metric in spaces of polynomial systems has been introduced and studied in a series of papers by Beltrán, Dedieu, Malajovich, and Shub. The interest of this metric comes from the fact that the associated geodesics avoid ill-conditioned problems and are a useful tool to improve classical complexity bounds for Bézout's theorem. The linear case is examined here: using nonsmooth nonconvex analysis techniques, we study the properties of condition geodesics in the space of full rank, real, or complex rectangular matrices. Our main results include an existence theorem for the boundary problem, a differential inclusion for such geodesics based on Clarke's generalized gradients, regularity properties, and a detailed description of a few particular cases: diagonal and unitary matrices. Moreover, we study condition geodesics from a numerical viewpoint, and we develop an effective algorithm that allows us to compute geodesics with given endpoints and helps to illustrate theoretical results and formulate new conjectures.

Computing a Nearest Correlation Matrix with Factor Structure

Rüdiger Borsdorf, Nicholas J. Higham, and Marcos Raydan

SIAM. J. Matrix Anal. & Appl. 31, pp. 2603-2622 (20 pages)

Online Publication Date: September 14, 2010

Full Text: | Download PDF

Show Abstract
An $n\times n$ correlation matrix has $k$ factor structure if its off-diagonal agrees with that of a rank $k$ matrix. Such correlation matrices arise, for example, in factor models of collateralized debt obligations (CDOs) and multivariate time series. We analyze the properties of these matrices and, in particular, obtain an explicit formula for the rank in the one factor case. Our main focus is on the nearness problem of finding the nearest $k$ factor correlation matrix $C(X) = \diag(I-XX^T) + XX^T$ to a given symmetric matrix, subject to natural nonlinear constraints on the elements of the $n\times k$ matrix $X$, where distance is measured in the Frobenius norm. For a special one parameter case we obtain an explicit solution. For the general $k$ factor case we obtain the gradient and Hessian of the objective function and derive an instructive result on the positive definiteness of the Hessian when $k=1$. We investigate several numerical methods for solving the nearness problem: the alternating directions method; a principal factors method used by Anderson, Sidenius, and Basu in the CDO application, which we show is equivalent to the alternating projections method and lacks convergence results; the spectral projected gradient method of Birgin, Martínez, and Raydan; and Newton and sequential quadratic programming methods. The methods differ in whether or not they can take account of the nonlinear constraints and in their convergence properties. Our numerical experiments show that the performance of the methods depends strongly on the problem, but that the spectral projected gradient method is the clear winner.

Structured Total Maximum Likelihood: An Alternative to Structured Total Least Squares

Amir Beck and Yonina C. Eldar

SIAM. J. Matrix Anal. & Appl. 31, pp. 2623-2649 (27 pages)

Online Publication Date: September 16, 2010

Full Text: | Download PDF

Show Abstract
Linear inverse problems with uncertain measurement matrices appear in many different applications. One of the standard techniques for solving such problems is the total least squares (TLS) method. Recently, an alternative approach has been suggested, based on maximizing an appropriate likelihood function assuming that the measurement matrix consists of random Gaussian variables. We refer to this technique as the total maximum likelihood (TML) method. Here we extend this strategy to the case in which the measurement matrix is structured so that the perturbations are not arbitrary but rather follow a fixed pattern. The resulting estimate is referred to as the structured TML (STML). As we show, the STML can be viewed as a regularized version of the structured TLS (STLS) approach in which the regularization consists of a logarithmic penalty. In contrast to the STLS solution, the STML always exists. Furthermore, its performance in practice tends to be superior to that of the STLS and competitive to other regularized solvers, as we illustrate via several examples. We also consider a few interesting special cases in which the STML can be computed efficiently either by reducing it into a one-dimensional problem regardless of the problem size or by a decomposition via a discrete Fourier transform.

Direction-Preserving and Schur-Monotonic Semiseparable Approximations of Symmetric Positive Definite Matrices

Ming Gu, Xiaoye S. Li, and Panayot S. Vassilevski

SIAM. J. Matrix Anal. & Appl. 31, pp. 2650-2664 (15 pages)

Online Publication Date: September 16, 2010

Full Text: | Download PDF

Show Abstract
For a given symmetric positive definite matrix $A\in\mathbf{R}^{N\times N}$, we develop a fast and backward stable algorithm to approximate $A$ by a symmetric positive definite semiseparable matrix, accurate to a constant multiple of any prescribed tolerance. In addition, this algorithm preserves the product, $AZ$, for a given matrix $Z\in\mathbf{R}^{N\times d}$, where $d\ll N$. Our algorithm guarantees the positive-definiteness of the semiseparable matrix by embedding an approximation strategy inside a Cholesky factorization procedure to ensure that the Schur complements during the Cholesky factorization all remain positive definite after approximation. It uses a robust direction-preserving approximation scheme to ensure the preservation of $AZ$. We present numerical experiments and discuss the potential implications of our work.

On a Parametrization of Positive Semidefinite Matrices with Zeros

Mathias Drton and Josephine Yu

SIAM. J. Matrix Anal. & Appl. 31, pp. 2665-2680 (16 pages)

Online Publication Date: September 23, 2010

Full Text: | Download PDF

Show Abstract
We study a class of parametrizations of convex cones of positive semidefinite matrices with prescribed zeros. Each such cone corresponds to a graph whose nonedges determine the prescribed zeros. Each parametrization in this class is a polynomial map associated with a simplicial complex supported on cliques of the graph. The images of the maps are convex cones, and the maps can only be surjective onto the cone of zero-constrained positive semidefinite matrices when the associated graph is chordal and the simplicial complex is the clique complex of the graph. Our main result gives a semialgebraic description of the images of the parametrizations for chordless cycles. The work is motivated by the fact that the considered maps correspond to Gaussian statistical models with hidden variables.

Spectral Methods for Parameterized Matrix Equations

Paul G. Constantine, David F. Gleich, and Gianluca Iaccarino

SIAM. J. Matrix Anal. & Appl. 31, pp. 2681-2699 (19 pages) | Cited 3 times

Online Publication Date: September 23, 2010

Full Text: | Download PDF

Show Abstract
We apply polynomial approximation methods—known in the numerical PDEs context as spectral methods—to approximate the vector-valued function that satisfies a linear system of equations where the matrix and the right-hand side depend on a parameter. We derive both an interpolatory pseudospectral method and a residual-minimizing Galerkin method, and we show how each can be interpreted as solving a truncated infinite system of equations; the difference between the two methods lies in where the truncation occurs. Using classical theory, we derive asymptotic error estimates related to the region of analyticity of the solution, and we present a practical residual error estimate. We verify the results with two numerical examples.

Optimizing Halley's Iteration for Computing the Matrix Polar Decomposition

Yuji Nakatsukasa, Zhaojun Bai, and François Gygi

SIAM. J. Matrix Anal. & Appl. 31, pp. 2700-2720 (21 pages)

Online Publication Date: September 29, 2010

Full Text: | Download PDF

Show Abstract
We introduce a dynamically weighted Halley (DWH) iteration for computing the polar decomposition of a matrix, and we prove that the new method is globally and asymptotically cubically convergent. For matrices with condition number no greater than $10^{16}$, the DWH method needs at most six iterations for convergence with the tolerance $10^{-16}$. The Halley iteration can be implemented via QR decompositions without explicit matrix inversions. Therefore, it is an inverse free communication friendly algorithm for the emerging multicore and hybrid high performance computing systems.

Decay in Functions of Multiband Matrices

N. Mastronardi, M. Ng, and E. E. Tyrtyshnikov

SIAM. J. Matrix Anal. & Appl. 31, pp. 2721-2737 (17 pages)

Online Publication Date: October 07, 2010

Full Text: | Download PDF

Show Abstract
The Benzi–Golub result on decay properties for matrix functions of a banded Hermitian matrix [BIT, 39 (1999), pp. 417–438] is extended to the case of multiband matrices. It is shown how the simple diagonal dominance technique applies to the general non-Hermitian case. We also present $O(1)$ algorithms computing matrix functions of multiband and multi-Toeplitz (multilevel Toeplitz) matrices in time that depends on the bandwidth and prescribed approximation accuracy but does not depend on the size of matrices.

$\mathcal{H}_2$-Optimal Model Reduction with Higher-Order Poles

Paul Van Dooren, Kyle A. Gallivan, and P.-A. Absil

SIAM. J. Matrix Anal. & Appl. 31, pp. 2738-2753 (16 pages)

Online Publication Date: October 14, 2010

Full Text: | Download PDF

Show Abstract
We revisit the problem of approximating a multiple-input multiple-output $p\times m$ rational transfer function $H(s)$ of high degree by another $p\times m$ rational transfer function $\widehat{H}(s)$ of much smaller degree, so that the $\mathcal{H}_2$-norm of the approximation error is minimized. We show that in the general case of higher-order poles in the reduced-order model, called the defective case, the stationary points of the $\mathcal{H}_2$-norm of the approximation error can still be characterized by tangential interpolation conditions. We also indicate that the sensitivity of the solution of this problem depends on the parameterization used.

Acquired Clustering Properties and Solution of Certain Saddle Point Systems

M. A. Olshanskii and V. Simoncini

SIAM. J. Matrix Anal. & Appl. 31, pp. 2754-2768 (15 pages)

Online Publication Date: October 14, 2010

Full Text: | Download PDF

Show Abstract
Many mathematical models involve flow equations characterized by nonconstant viscosity, and a Stokes-type problem with variable viscosity coefficient arises. Appropriate block diagonal preconditioners for the resulting algebraic saddle point linear system produce well-clustered spectra, except for a few interior isolated eigenvalues which may tend to approach zero. These outliers affect the convergence of Krylov subspace system solvers, causing a possibly long stagnation phase. In this paper we characterize the influence of the spectral properties of the preconditioner on the final spectrum of the saddle point matrix by providing accurate spectral intervals depending on the involved operators. Moreover, we suggest that the stagnation phase may be completely eliminated by means of an augmentation procedure, where approximate spectral eigenspace information can be injected. We show that the modifications to the original code are minimal and can be easily implemented. Numerical experiments confirm our findings.

Smoothed Analysis of Moore–Penrose Inversion

Peter Bürgisser and Felipe Cucker

SIAM. J. Matrix Anal. & Appl. 31, pp. 2769-2783 (15 pages)

Online Publication Date: October 14, 2010

Full Text: | Download PDF

Show Abstract
We perform a smoothed analysis of the condition number of rectangular matrices. We prove that, asymptotically, the expected value of this condition number depends only on the elongation of the matrix and not on the center and variance of the underlying probability distribution.

Solving a Structured Quadratic Eigenvalue Problem by a Structure-Preserving Doubling Algorithm

Chun-Hua Guo and Wen-Wei Lin

SIAM. J. Matrix Anal. & Appl. 31, pp. 2784-2801 (18 pages)

Online Publication Date: October 14, 2010

Full Text: | Download PDF

Show Abstract
In studying the vibration of fast trains, we encounter a palindromic quadratic eigenvalue problem (QEP) $(\lambda^{2}A^{T}+\lambda Q+A)z=0$, where $A,Q\in\mathbb{C}^{n\times n}$ and $Q^{T}=Q$. Moreover, the matrix $Q$ is block tridiagonal and block Toeplitz, and the matrix $A$ has only one nonzero block in the upper-right corner. So most of the eigenvalues of the QEP are zero or infinity. In a linearization approach, one typically starts with deflating these known eigenvalues for the sake of efficiency. However, this initial deflation process involves the inverses of two potentially ill-conditioned matrices. As a result, large error might be introduced into the data for the reduced problem. In this paper we propose using the solvent approach directly on the original QEP, without any deflation process. We apply a structure-preserving doubling algorithm to compute the stabilizing solution of the matrix equation $X+A^{T}X^{-1}A=Q$, whose existence is guaranteed by a result on the Wiener–Hopf factorization of rational matrix functions associated with semi-infinite block Toeplitz matrices and a generalization of Bendixson's theorem to bounded linear operators on Hilbert spaces. The doubling algorithm is shown to be well defined and quadratically convergent. The complexity of the doubling algorithm is drastically reduced by using the Sherman–Morrison–Woodbury formula and the special structures of the problem. Once the stabilizing solution is obtained, all nonzero finite eigenvalues of the QEP can be found efficiently and with the automatic reciprocal relationship, while the known eigenvalues at zero or infinity remain intact.

Matrix $p$-Norms Are NP-Hard to Approximate If $p\neq1,2,\infty$

Julien M. Hendrickx and Alex Olshevsky

SIAM. J. Matrix Anal. & Appl. 31, pp. 2802-2812 (11 pages)

Online Publication Date: October 28, 2010

Full Text: | Download PDF

Show Abstract
We show that, for any rational $p\in[1,\infty)$ except $p=1,2$, unless $P=NP$, there is no polynomial time algorithm which approximates the matrix $p$-norm to arbitrary relative precision. We also show that, for any rational $p\in[1,\infty)$ including $p=1,2$, unless $P=NP$, there is no polynomial-time algorithm which approximates the $\infty,p$ mixed norm to some fixed relative precision.

Preconditioning Stochastic Galerkin Saddle Point Systems

Catherine E. Powell and Elisabeth Ullmann

SIAM. J. Matrix Anal. & Appl. 31, pp. 2813-2840 (28 pages) | Cited 2 times

Online Publication Date: November 04, 2010

Full Text: | Download PDF

Show Abstract
Mixed finite element discretizations of deterministic second-order elliptic PDEs lead to saddle point systems for which the study of iterative solvers and preconditioners is mature. Galerkin approximation of solutions of stochastic second-order elliptic PDEs, which couple standard mixed finite element discretizations in physical space with global polynomial approximation on a probability space, also give rise to linear systems with familiar saddle point structure. For stochastically nonlinear problems, the solution of such systems presents a serious computational challenge. The blocks are sums of Kronecker products of pairs of matrices associated with two distinct discretizations, and the systems are large, reflecting the curse of dimensionality inherent in most stochastic approximation schemes. Moreover, for the problems considered herein, the leading blocks of the saddle point matrices are block-dense, and the cost of a matrix vector product is nontrivial. We implement a stochastic Galerkin discretization for the steady-state diffusion problem written as a mixed first-order system. The diffusion coefficient is assumed to be a lognormal random field, approximated via a nonlinear function of a finite number of Gaussian random variables. We study the resulting saddle point systems and investigate the efficiency of block-diagonal preconditioners of Schur complement and augmented type for use with the minimal residual method (MINRES). By introducing so-called Kronecker product preconditioners, we improve the robustness of cheap, mean-based preconditioners with respect to the statistical properties of the stochastically nonlinear diffusion coefficients.

Rigorous Perturbation Bounds of Some Matrix Factorizations

X.-W. Chang and D. Stehlé

SIAM. J. Matrix Anal. & Appl. 31, pp. 2841-2859 (19 pages)

Online Publication Date: November 04, 2010

Full Text: | Download PDF

Show Abstract
This article presents rigorous normwise perturbation bounds for the Cholesky, LU, and QR factorizations with normwise or componentwise perturbations in the given matrix. The considered componentwise perturbations have the form of backward rounding errors for the standard factorization algorithms. The used approach is a combination of the classic and refined matrix equation approaches. Each of the new rigorous perturbation bounds is a small constant multiple of the corresponding first-order perturbation bound obtained by the refined matrix equation approach in the literature and can be estimated efficiently. These new bounds can be much tighter than the existing rigorous bounds obtained by the classic matrix equation approach, while the conditions for the former to hold are almost as moderate as the conditions for the latter to hold.

Structured Pseudospectra and the Condition of a Nonderogatory Eigenvalue

Michael Karow

SIAM. J. Matrix Anal. & Appl. 31, pp. 2860-2881 (22 pages) | Cited 1 time

Online Publication Date: November 30, 2010

Full Text: | Download PDF

Show Abstract
Let $\lambda$ be a nonderogatory eigenvalue of $A\in\mathbb{C}^{n\times n}$ of algebraic multiplicity $m$. The sensitivity of $\lambda$ with respect to matrix perturbations of the form $A\leadsto A+\Delta$, $\Delta\in\boldsymbol{\Delta}$, is measured by the structured condition number $\kappa_{\boldsymbol{\Delta}}(A,\lambda)$. Here $\boldsymbol{\Delta}$ denotes the set of admissible perturbations. However, if $\boldsymbol{\Delta}$ is not a vector space over $\mathbb{C}$, then $\kappa_{\boldsymbol{\Delta}}(A,\lambda)$ provides only incomplete information about the mobility of $\lambda$ under small perturbations from $\boldsymbol{\Delta}$. The full information is then given by the set $K_{\boldsymbol{\Delta}}(x,y)=\{y^*\Delta x;$ $\Delta\in\boldsymbol{\Delta},$ $\|\Delta\|\leq1\}\subset\mathbb{C}$ that depends on $\boldsymbol{\Delta}$, a pair of normalized right and left eigenvectors $x,y$, and the norm $\|\cdot\|$ that measures the size of the perturbations. We always have $\kappa_{\boldsymbol{\Delta}}(A,\lambda)=\max\{|z|^{1/m};$ $z\in K_{\boldsymbol{\Delta}}(x,y)\}$. Furthermore, $K_{\boldsymbol{\Delta}}(x,y)$ determines the shape and growth of the $\boldsymbol{\Delta}$-structured pseudospectrum in a neighborhood of $\lambda$. In this paper we study the sets $K_{\boldsymbol{\Delta}}(x,y)$ and obtain methods for computing them. In doing so we obtain explicit formulae for structured eigenvalue condition numbers with respect to many important perturbation classes.

Improved Bounds on Restricted Isometry Constants for Gaussian Matrices

Bubacarr Bah and Jared Tanner

SIAM. J. Matrix Anal. & Appl. 31, pp. 2882-2898 (17 pages)

Online Publication Date: November 30, 2010

Full Text: | Download PDF

Show Abstract
The restricted isometry constant (RIC) of a matrix $A$ measures how close to an isometry is the action of $A$ on vectors with few nonzero entries, measured in the $\ell^2$ norm. Specifically, the upper and lower RICs of a matrix $A$ of size $n\times N$ are the maximum and the minimum deviation from unity (one) of the largest and smallest, respectively, square of singular values of all ${N\choose k}$ matrices formed by taking $k$ columns from $A$. Calculation of the RIC is intractable for most matrices due to its combinatorial nature; however, many random matrices typically have bounded RIC in some range of problem sizes $(k,n,N)$. We provide the best known bound on the RIC for Gaussian matrices, which is also the smallest known bound on the RIC for any large rectangular matrix. Our results are built on the prior bounds of Blanchard, Cartis, and Tanner [SIAM Rev., to appear], with improvements achieved by grouping submatrices that share a substantial number of columns.

Robust Approximate Cholesky Factorization of Rank-Structured Symmetric Positive Definite Matrices

Jianlin Xia and Ming Gu

SIAM. J. Matrix Anal. & Appl. 31, pp. 2899-2920 (22 pages) | Cited 1 time

Online Publication Date: November 30, 2010

Full Text: | Download PDF

Show Abstract
Given a symmetric positive definite matrix $A$, we compute a structured approximate Cholesky factorization $A\approx\mathbf{R}^{T}\mathbf{R}$ up to any desired accuracy, where $\mathbf{R}$ is an upper triangular hierarchically semiseparable (HSS) matrix. The factorization is stable, robust, and efficient. The method compresses off-diagonal blocks with rank-revealing orthogonal decompositions. In the meantime, positive semidefinite terms are automatically and implicitly added to Schur complements in the factorization so that the approximation $\mathbf{R}^{T}\mathbf{R}$ is guaranteed to exist and be positive definite. The approximate factorization can be used as a structured preconditioner which does not break down. No extra stabilization step is needed. When $A$ has an off-diagonal low-rank property, or when the off-diagonal blocks of $A$ have small numerical ranks, the preconditioner is data sparse and is especially efficient. Furthermore, the method has a good potential to give satisfactory preconditioning bounds even if this low-rank property is not obvious. Numerical experiments are used to demonstrate the performance of the method. The method can be used to provide effective structured preconditioners for large sparse problems when combined with some sparse matrix techniques. The hierarchical compression scheme in this work is also useful in the development of more HSS algorithms.

Schur Decompositions of a Matrix and the Boundary of Its Pseudospectrum

Lyonell Boulton and Peter Lancaster

SIAM. J. Matrix Anal. & Appl. 31, pp. 2921-2933 (13 pages)

Online Publication Date: November 30, 2010

Full Text: | Download PDF

Show Abstract
We discuss the notion of irreducible block Schur decomposition of a complex square matrix and show how such a decomposition provides information about singularities in the boundaries of pseudospectra of the matrix. Detailed examples are provided to illustrate the theory.

Spectral Condition Numbers of Orthogonal Projections and Full Rank Linear Least Squares Residuals

Joseph F. Grcar

SIAM. J. Matrix Anal. & Appl. 31, pp. 2934-2949 (16 pages)

Online Publication Date: December 08, 2010

Full Text: | Download PDF

Show Abstract
A simple formula is proved to be a tight estimate for the condition number of the full rank linear least squares residual with respect to the matrix of least squares coefficients and scaled 2-norms. The tight estimate reveals that the condition number depends on three quantities, two of which can cause ill-conditioning. The numerical linear algebra literature presents several estimates of various instances of these condition numbers. All the prior values exceed the formula introduced here, sometimes by large factors.

The Trace Ratio Optimization Problem for Dimensionality Reduction

T. T. Ngo, M. Bellalij, and Y. Saad

SIAM. J. Matrix Anal. & Appl. 31, pp. 2950-2971 (22 pages)

Online Publication Date: December 16, 2010

Full Text: | Download PDF

Show Abstract
This paper considers the problem of optimizing the ratio $\mathrm{Tr}[V^{T}AV]/\mathrm{Tr}[V^{T}BV]$ over all unitary matrices $V$ with $p$ columns, where $A,B$ are two positive definite matrices. This problem is common in supervised learning techniques. However, because its numerical solution is typically expensive it is often replaced by the simpler optimization problem which consists of optimizing $\mathrm{Tr}[V^{T}AV]$ under the constraint that $V^{T}BV=I$, the identity matrix. The goal of this paper is to examine this trace ratio optimization problem in detail, to consider different algorithms for solving it, and to illustrate the use of these algorithms for dimensionality reduction.

Riemannian Newton Method for the Multivariate Eigenvalue Problem

Lei-Hong Zhang

SIAM. J. Matrix Anal. & Appl. 31, pp. 2972-2996 (25 pages)

Online Publication Date: December 16, 2010

Full Text: | Download PDF

Show Abstract
The multivariate eigenvalue problem (MEP) which originally arises from the canonical correlation analysis is an important generalization of the classical eigenvalue problem. Recently, the MEP also finds applications in many other areas and continues to receive interest. However, the existing algorithms for the MEP are the generalization of the power iteration for the classical eigenvalue problem and converge slowly. In this paper, we propose a Riemannian Newton method for the MEP, which is a generalization of the classical Rayleigh quotient iteration (RQI). Under a mild condition, the local quadratic convergence can be guaranteed. We also develop the inexact implementation by employing some Krylov subspace method and establishing the preconditioning technique to obtain an inexact Riemannian Newton step efficiently. Preliminary but promising numerical experiments are reported which show a good convergence performance in terms of the proposed method's speed and global convergence.
Close

close