Commuting Hermitian matrices may be simultaneously diagonalized by a common unitary matrix. However, the numerical aspects are delicate. We revisit a previously rejected numerical approach in a new algorithm called 'do-one-then-do-the-other'. One of two input matrices is diagonalized by a unitary similarity, and then the computed eigenvectors are applied to the other input matrix. Additional passes are applied as necessary to resolve invariant subspaces associated with repeated eigenvalues and eigenvalue clusters. The algorithm is derived by first developing a spectral divide-and-conquer method and then allowing the method to break the spectrum into, not just two invariant subspaces, but as many as safely possible. Most computational work is delegated to a black-box eigenvalue solver, which can be tailored to specific computer architectures. The overall running time is a small multiple of a single eigenvalue-eigenvector computation, even on difficult problems with tightly clustered eigenvalues. The article concludes with applications to a structured eigenvalue problem and a highly sensitive eigenvector computation.

The minimum semidefinite rank (msr) of a graph is defined to be the minimum rank among all positive semidefinite matrices whose zero/nonzero pattern corresponds to that graph. We recall some known facts and present new results, including results concerning the effects of vertex or edge removal from a graph on msr.

We consider a structured inverse eigenvalue problem in which the eigenvalues of a real symmetric matrix are specified and selected entries may be constrained to take specific numerical values or to be nonzero. This includes the problem of specifying the graph of the matrix, which is determined by the locations of zero and nonzero entries. In this article, we develop a numerical method for constructing a solution to the structured inverse eigenvalue problem. The problem is recast as a constrained optimization problem over the orthogonal manifold, and a numerical optimization routine seeks its solution.

My new book Numerical Analysis: Theory and Experiments presents rapidly converging numerical methods and develops techniques, both theoretical and experimental, for analyzing them. The book is based on my teaching experiences and important recent developments in the field. In my own numerical analysis course, I emphasize a multifaceted approach to problem solving. Students are expected to assess properties such as smoothness and conditioning, to predict performance of numerical methods, to implement solutions in computer code, to measure rates of convergence with experiments, and to interpret results using graphics. The book reflects this structure, interleaving mathematical analysis, computer code, graphics, and discussion. Lab problems from my courses appear as exercises throughout the book, and hundreds of additional exercises support a variety of student needs and course structures. The outside development that motivated a new book was the recent embrace of Chebyshev technology by the numerical analysis community. Compared with older methods, Chebyshev methods can converge blazingly fast. However, they also require a more subtle analysis. Whereas many older methods achieve their full potential on functions possessing a few derivatives, Chebyshev methods achieve their full potential on functions that are analytic in the complex plane. My book presents both kinds of methods in parallel from the earliest chapters. This requires additional investment in understanding function smoothness, but the benefit is immense; we can achieve geometric or even supergeometric rates of convergence and exhaust the precision of the computer in a fraction of a second. Numerical Analysis: Theory and Experiments provides hands-on experience with, and careful analysis of, some of the best numerical methods available today.

We introduce a backward stable algorithm for computing the CS decomposition of a partitioned $2n \times n$ matrix with orthonormal columns, or a rank-deficient partial isometry. The algorithm computes two $n \times n$ polar decompositions (which can be carried out in parallel) followed by an eigendecomposition of a judiciously crafted $n \times n$ Hermitian matrix. We prove that the algorithm is backward stable whenever the aforementioned decompositions are computed in a backward stable way. Since the polar decomposition and the symmetric eigendecomposition are highly amenable to parallelization, the algorithm inherits this feature. We illustrate this fact by invoking recently developed algorithms for the polar decomposition and symmetric eigendecomposition that leverage Zolotarev's best rational approximations of the sign function. Numerical examples demonstrate that the resulting algorithm for computing the CS decomposition enjoys excellent numerical stability.

“Low temperature” random matrix theory is the study of random eigenvalues as energy is removed. In standard notation, β is identified with inverse temperature, and low temperatures are achieved through the limit β → ∞. In this paper, we derive statistics for low-temperature random matrices at the “soft edge,” which describes the extreme eigenvalues for many random matrix distributions. Specifically, new asymptotics are found for the expected value and standard deviation of the general-β Tracy-Widom distribution. The new techniques utilize beta ensembles, stochastic differential operators, and Riccati diffusions. The asymptotics fit known high-temperature statistics curiously well and contribute to the larger program of general-β random matrix theory.

When an orthogonal matrix is partitioned into a two-by-two block structure, its four blocks can be simultaneously bidiagonalized. This observation underlies numerically stable algorithms for the CS decomposition and the existence of CMV matrices for orthogonal polynomial recurrences. We discover a new matrix decomposition for simultaneous multidiagonalization, which reduces the blocks to any desired bandwidth. Its existence is proved, and a backward stable algorithm is developed. The resulting matrix with banded blocks is parameterized by a product of Givens rotations, guaranteeing orthogonality even on a finite-precision computer. The algorithm relies heavily on Level 3 BLAS routines and supports parallel computation.

This paper serves to prove the thesis that a computational trick can open entirely new approaches to theory. We illustrate this by describing such random matrix techniques as the stochastic operator approach, the method of ghosts and shadows, and the method of "Riccatti Diffusion/Sturm Sequences." We thereby provide new insights into the deeper mathematics underlying random matrix theory.

This paper serves to prove the thesis that a computational trick can open entirely new approaches to theory. We illustrate this by describing such random matrix techniques as the stochastic operator approach, the method of ghosts and shadows, and the method of “Riccatti Diffusion/Sturm Sequences.” We thereby provide new insights into the deeper mathematics underlying random matrix theory.

We develop a divide-and-conquer algorithm for the bidiagonal CS decomposition (CSD). This complements an earlier algorithm based on simultaneous QR iteration. The new algorithm is designed to provide the efficiency gains of familiar divide-and-conquer algorithms on both serial and parallel architectures. The solution uses many components of existing algorithms, particularly the bidiagonal SVD algorithm of Gu and Eisenstat, but extra steps and reparameterizations are required to maintain orthogonality and consistent singular vectors, especially when the singular vectors are ill conditioned. The algorithm supports the stable computation of the generalized singular value decomposition (GSVD) in addition to the CSD.

Since its discovery in 1977, the CS decomposition (CSD) has resisted computation, even though it is a sibling of the well-understood eigenvalue and singular value decompositions. Several algorithms have been developed for the reduced 2-by-1 form of the decomposition, but none have been extended to the complete 2-by-2 form of the decomposition in Stewart's original paper. In this article, we present an algorithm for simultaneously bidiagonalizing the four blocks of a unitary matrix partitioned into a 2-by-2 block structure. This serves as the first, direct phase of a two-stage algorithm for the CSD, much as Golub-Kahan-Reinsch bidiagonalization serves as the first stage in computing the singular value decomposition. Backward stability is proved.

An algorithm for computing the complete CS decomposition of a partitioned unitary matrix is developed. Although the existence of the CS decomposition (CSD) has been recognized since 1977, prior algorithms compute only a reduced version. This reduced version, which might be called a 2-by-1 CSD, is equivalent to two simultaneous singular value decompositions. The algorithm presented in this article computes the complete 2-by-2 CSD, which requires the simultaneous diagonalization of all four blocks of a unitary matrix partitioned into a 2-by-2 block structure. The algorithm appears to be the only fully specified algorithm available. The computation occurs in two phases. In the first phase, the unitary matrix is reduced to bidiagonal block form, as described by Sutton and Edelman. In the second phase, the blocks are simultaneously diagonalized using techniques from bidiagonal SVD algorithms of Golub, Kahan, Reinsch, and Demmel. The algorithm has a number of desirable numerical features.

We are generally concerned with the possible lists of multiplicities for the eigenvalues of a real symmetric matrix with a given graph. Many restrictions are known, but it is often problematic to construct a matrix with desired multiplicities, even if a matrix with such multiplicities exists. Here, we develop a technique for construction using the implicit function theorem in a certain way. We show that the technique works for a large variety of trees, give examples and determine all possible multiplicities for a large class of trees for which this was not previously known.

Let $\mathcal{P}(G)$ be the set of all positive semidefinite matrices whose graph is $G$, and $\operatorname{msr}(G)$ be the minimum rank of all matrices in $\mathcal{P}(G)$. Upper and lower bounds for $\operatorname{msr}(G)$ are given and used to determine $\operatorname{msr}(G)$ for some well-known graphs, including chordal graphs, and for all simple graphs on less than seven vertices.

We provide a solution to the β-Jacobi matrix model problem posed by Dumitriu and the first author. The random matrix distribution introduced here, called a matrix model, is related to the model of Killip and Nenciu, but the development is quite different. We start by introducing a new matrix decomposition and an algorithm for computing this decomposition. Then we run the algorithm on a Haar-distributed random matrix to produce the β-Jacobi matrix model. The Jacobi ensemble on ${\Bbb R}^{n}$ , parametrized by β > 0, a > -1,and b > -1, is the probability distribution whose density is proportional to $\prod_{i}\lambda_{i}^{({\beta}/{2})(a+1)-1}(1-\lambda_{i})^{({\beta}/{2})(b+1)-1}\prod_{i 0. Observing a connection between Haar measure on the orthogonal (resp., unitary) group and pairs of real (resp., complex) Gaussian matrices, we find a direct connection between multivariate analysis of variance (MANOVA) and the new matrix model.

We propose that classical random matrix models are properly viewed as finite difference schemes for stochastic differential operators. Three particular stochastic operators commonly arise, each associated with a familiar class of local eigenvalue behavior. The stochastic Airy operator displays soft edge behavior, associated with the Airy kernel. The stochastic Bessel operator displays hard edge behavior, associated with the Bessel kernel. The article concludes with suggestions for a stochastic sine operator, which would display bulk behavior, associated with the sine kernel.

Let $\kappa$ be the condition number of an $m$-by-$n$ matrix with independent standard Gaussian entries, either real ($\beta = 1$) or complex ($\beta = 2$). The major result is the existence of a constant $C$ (depending on $m$, $n$, and $\beta$) such that $P[\kappa > x] < C \, x^{-\beta}$ for all $x$. As $x \rightarrow \infty$, the bound is asymptotically tight. An analytic expression is given for the constant $C$, and simple estimates are given, one involving a Tracy--Widom largest eigenvalue distribution. All of the results extend beyond real and complex entries to general $\beta$.

Given an n-by-n Hermitian matrix A and a real number $\lambda$, index $i$ is said to be Parter (resp., neutral, downer) if the multiplicity of $\lambda$ as an eigenvalue of the principal submatrix $A(i)$ is one more (resp., the same, one less) than that in $A$. In case the multiplicity of $\lambda$ in $A$ is at least 2 and the graph of $A$ is a tree, there are always Parter vertices. Our purpose here is to advance the classification of vertices and, in particular, to relate classification to the combinatorial structure of eigenspaces. Some general results are given and then used to deduce some rather specific facts not otherwise easily observed. Examples are given.