Search

This thesis presents new results on probabilistic error analysis. By freeing ourselves of the burden of having to account for every possible worst-case scenario, we provide error bounds that are more informative than the long standing worst-case bounds. We also see in practice how these new bounds can guide the development of new algorithms. We begin by studying stochastic rounding, a rounding mode which rounds a computed result to the next larger or smaller floating-point number randomly. We compare stochastic rounding with round-to-nearest, finding some similarities but also a large number of properties that hold for round-to-nearest but fail to hold for stochastic rounding. Importantly, we show that the rounding errors produced by stochastic rounding are mean zero, mean independent random variables. Using this fact, we build on earlier probabilistic error analysis to show that for a wide range of inner product based linear algebra computations, stochastic rounding provides backward error bounds where we can take the square root of the dimensional constants in the worst-case bounds. This has been a well known rule of thumb for some time, but we show it to be a rule for stochastic rounding. Expanding on this work for inner-product based algorithms, we investigate orthogonal transformations. Using the same model of rounding errors as that which stochastic rounding satisfies, we show that for Householder QR factorization, and other related algorithms, we see the same square rooting of the dimensional constant in the backward error bound. The analysis makes use of a matrix concentration inequality, where previous probabilistic error analyses employed scalar concentration inequalities. We also derive a new backward error formula for QR factorization, used in our numerical experiments which validate our bounds. Finally, we consider the randomized SVD, a well known method for computing low rank matrix approximations. We perform a complete rounding error analysis of the fixed-rank problem, the most straightforward version of the randomized SVD where we have a priori specified a target rank. This is complementary to the existing literature providing error bounds on this procedure in exact arithmetic. While our initial analysis is worst-case, we use our previous probabilistic results to refine these error bounds. Using these refined bounds, we propose a mixed-precision version of the algorithm that offers potential speedups by gradually reducing the precision during the execution of the algorithm.

Functions of matrices play an important role in many applications in science and engineering. Their reliable computation has been a topic of interest in numerical linear algebra over the decades, and a wide variety of methods for computing different functions have been studied. Meanwhile, the interest in multiple precision computing environments has been growing in recent years and nowadays there is an explosion of floating-point arithmetics beyond the most widely used standard IEEE binary32 and binary64 formats. Under such background, there are several works dedicated to the computation of matrix functions in arbitrary precision arithmetic, but overall this research topic has not yet attracted much attention. In this thesis, we study methods for computing functions of matrices in arbitrary precision arithmetic. Unlike many existing algorithms that are tightly coupled to a specific precision of floating-point arithmetic, the algorithms we develop take the unit roundoff of the working precision as an input argument since this is known only at runtime, and so works in an arbitrary precision. First, we provide a version of the Schur--Parlett algorithm that requires only function values and not derivatives. The algorithm requires access to arithmetic of a matrix-dependent precision at least double the working precision, which is used to evaluate the function on the diagonal blocks of order greater than 2 (if there are any) of the reordered and blocked Schur form. The key idea is to compute by diagonalization the function of a small random diagonal perturbation of each diagonal block, where the perturbation ensures that diagonalization will succeed. The algorithm is inspired by Davies's randomized approximate diagonalization method, but we explain why that is not a reliable numerical method for computing matrix functions. This multiprecision Schur--Parlett algorithm is applicable to arbitrary analytic functions f and, like the original Schur--Parlett algorithm, it generally behaves in a numerically stable fashion. The algorithm is especially useful when the derivatives of f are not readily available or accurately computable. We apply our algorithm to the matrix Mittag--Leffler function and show that it yields results of accuracy similar to, and in some cases much greater than, the state-of-the-art algorithm for this function. Second, we develop an algorithm for computing the matrix cosine in arbitrary precision. The algorithm employs a Taylor approximation with scaling and recovering and it can be used with a Schur decomposition or in a decomposition-free manner. We also derive a framework for computing the Frechet derivative, construct an efficient evaluation scheme for computing the cosine and its Frechet derivative simultaneously in arbitrary precision, and show how this scheme can be extended to compute the matrix sine, cosine, and their Frechet derivatives all together. Numerical experiments show that the new algorithms behave in a forward stable way over a wide range of precisions. The transformation-free version of the algorithm for computing the cosine is competitive in accuracy with the state-of-the-art algorithms in double precision and surpasses existing alternatives in both speed and accuracy in working precisions higher than double. Finally, we consider the problem of computing the square root of a perturbation of the scaled identity matrix, A = al + UV*, where U and V are n x k matrices with k ≤ n. This problem arises in various applications, including computer vision and optimization methods for machine learning. We derive a new formula for the pth root of A that involves a weighted sum of powers of the pth root of the k ˆ k matrix aIk ` V °U. This formula is particularly attractive for the square root, since the sum has just one term when p " 2. We derive a new class of Newton iterations for computing the square root that exploit the low-rank structure. Theses methods can be employed in arbitrary precision by simply executing all elementary scalar operations in arbitrary precision, and, additionally, for the iterative algorithms, by properly adjusting the internal tolerance that is used as stopping criterion. We also proposed several Schur-based methods that can utilize the structure of A. In particular, we established a scheme to obtain the Schur decomposition of the n ˆ n matrix UV° from the Schur decomposition of the k ˆ k matrix V °U. We test these new methods on random matrices and on positive definite matrices arising in applications. Numerical experiments show that the new approaches can yield much smaller residual than existing alternatives and can be significantly faster when the perturbation UV° has low rank.

In this thesis we focus on the theoretical and computational aspects of nonlinear eigenvalue problems (NEPs), which arise in several fields of computational science and engineering, such as fluid dynamics, optics, and structural engineering. In the last twenty years several researchers devoted their time in studying efficient and precise ways to solve NEPs, which cemented their importance in numerical linear algebra. The most successful algorithms developed towards this goal are either based on contour integrals, or on rational approximations and linearizations. The first part of the thesis is dedicated to contour integral algorithms. In this framework, one computes specific integrals of a holomorphic function G(z) over the contour of a target region and exploits results of complex analysis to retrieve the eigenvalues of G(z) inside it. Our main contribution consists in having expanded the theory to include meromorphic functions, i.e., functions with poles inside the target region. We showed that under some circumstances, these algorithms can mistake a pole for an eigenvalue, but these situations are easily recognised and the main results from the holomorphic case can be extended. Furthermore, we proposed several heuristics to automatically choose the parameters that are needed to precisely retrieve the eigenpairs. In the second part of the thesis, we focus on rational approximations. Our goal was developing algorithms that construct robust, i.e., reliable for a given tolerance and scaling independent, rational approximations for functions given in split form or in black-box form. In the first case, we proposed a variant of the set-valued AAA, named weighted AAA, which guarantees to return an approximation with the chosen accuracy. In the second one, we built a two-phase approach, where an initial step performed by the surrogate AAA is followed by a cyclic Leja-Bagby refinement. We concluded the section with numerous numerical experiments based on the NLEVP library, comparing contour integral and rational approximation algorithms. The third and final part of the thesis is about tropical linear algebra. Our research on this topic started has a way to set the parameters of the aforementioned contour integral algorithms: in order to do that, we extended the theory of tropical roots from tropical polynomials to tropical Laurent series. Unlike in the polynomial case, a tropical Laurent series may have infinite tropical roots, but they are still in bijection with the slopes of the associated Newton polygon and they still provide annuli of exclusion for the eigenvalues of the Laurent series.

In the first part of this thesis, methods for the partial solution of generalized eigenvalue problems arising from structural dynamics are studied. A natural choice for partially solving the generalized eigenvalue problem (GEP) is the Lanczos iteration, or the shift-and-invert Lanczos (SIL) algorithm if a large number of eigenpairs is required. When the external loading of the associated structural dynamics problem is specific to that of earthquake engineering, the contribution of an eigenvector of a generalized eigenvalue problem can be quantified by a simple expression. Based on this property, a novel shifting strategy for the SIL algorithm which arises from the theory of orthogonal polynomials is presented. We discuss the Lanczos algorithm and its variant suited for the partial solution of the GEP arising from structural dynamics, followed by the shift-and-invert Lanczos algorithm. Various theoretical aspects and numerical issues of these algorithms are examined and addressed. The Ritz vector and the Lanczos vector methods are also introduced. These methods use vectors other than eigenvectors for the solution of structural dynamics problems. The Ritz vector method is widely used by engineers and is often implemented in engineering software. Orthogonal polynomials and the theory required to devise a new shifting strategy for the SIL algorithm are then introduced. With a specific choice of the starting vector for the Lanczos process it becomes possible to identify the intervals of the spectrum of the associated GEP in which the eigenvalues corresponding to the eigenvectors of interest lie. The shifts for the SIL algorithm are then chosen in the middle of these intervals, and a stopping criterion is provided. The algorithm is termed MASIL. The numerical experiments are performed on real engineering problems with our implementations of SIL, MASIL, Ritz vector, and Lanczos vector methods. While providing a comparable approximation to the solution of the structural dynamics problem, the MASIL approach computes up to 70% fewer eigenvectors and requires fewer shifts, on average, when compared with the standard shifting strategy used for this problem. The Ritz and Lanczos vector methods often provide an accurate approximation to the solution of the structural dynamics problem, but we show that there are cases where these methods do not provide a satisfactory approximation. In the second part of this thesis we explore max-plus algebra. It has been recently observed that an order of magnitude approximations to the absolute values of the eigenvalues of linear and polynomial eigenvalue problems can be obtained by the use of max-plus algebra, when the underlying matrices are unstructured and have large variations between the magnitudes of their entries. We extend the existing algorithms to return only some of the eigenvalues of the standard and generalized max-plus eigenvalue problems, and provide a result which allows for an estimation of the number of eigenvalues in modulus in a given interval. We then present numerical experiments assessing the quality of the max-plus approximations, and find that an order of magnitude approximations to the eigenvalues of symmetric definite eigenvalue problems can also be obtained.

Numerical methods based on rational Krylov spaces have become an indispensable tool of scientific computing. In this thesis we study rational Krylov spaces by considering rational Krylov decompositions; matrix relations which, under certain conditions, are associated with these spaces. We investigate the algebraic properties of such decompositions and present an implicit Q theorem for rational Krylov spaces. We derive standard and harmonic Ritz extraction strategies for approximating the eigenpairs of a matrix and for approximating the action of a matrix function onto a vector. While these topics have been considered previously, our approach does not require the last pole to be infinite, which makes the extraction procedure computationally more efficient. Typically, the computationally most expensive component of the rational Arnoldi algorithm for computing a rational Krylov basis is the solution of a large linear system of equations at each iteration. We explore the option of solving several linear systems simultaneously, thus constructing the rational Krylov basis in parallel. If this is not done carefully, the basis being orthogonalized may become poorly conditioned, leading to numerical instabilities in the orthogonalization process. We introduce the new concept of continuation pairs which gives rise to a near-optimal parallelization strategy that allows to control the growth of the condition number of this non orthogonal basis. As a consequence we obtain a more accurate and reliable parallel rational Arnoldi algorithm. The computational benefits are illustrated using our high performance C++ implementation. We develop an iterative algorithm for solving nonlinear rational least squares problems. The difficulty is in finding the poles of a rational function. For this purpose, at each iteration a rational Krylov decomposition is constructed and a modified linear problem is solved in order to relocate the poles to new ones. Our numerical results indicate that the algorithm, called RKFIT, is well suited for model order reduction of linear time invariant dynamical systems and for optimisation problems related to exponential integration. Furthermore, we derive a strategy for the degree reduction of the approximant obtained by RKFIT. The rational function obtained by RKFIT is represented with the aid of a scalar rational Krylov decomposition and an additional coefficient vector. A function represented in this form is called an RKFUN. We develop efficient methods for the evaluation, pole and root finding, and for performing basic arithmetic operations with RKFUNs. Lastly, we discuss RKToolbox, a rational Krylov toolbox for MATLAB, which implements all our algorithms and is freely available from http://rktoolbox.org. RKToolbox also features an extensive guide and a growing number of examples. In particular, most of our numerical experiments are easily reproducible by downloading the toolbox and running the corresponding example files in MATLAB.

A range of numerical linear algebra problems that arise in finite element-based structural analysis are considered. These problems were encountered when implementing the finite element method in the software package Oasys GSA. We present novel solutions to these problems in the form of a new method for error detection, algorithms with superior numerical effeciency and algorithms with scalable performance on parallel computers. The solutions and their corresponding software implementations have been integrated into GSA's program code and we present results that demonstrate the use of these implementations by engineers to solve real-world structural analysis problems.

In this thesis we develop new theoretical and numerical results for matrix polynomials and polynomial eigenproblems. This includes the cases of standard and generalized eigenproblems. Two chapters concern quadratic eigenproblems $(M\lambda^2+D\lambda+K)x=0$, where $M$, $D$ and $K$ enjoy special properties that are commonly encountered in modal analysis. We discuss this application in some detail, in particular the mathematics behind discrete dampers. We show how the physical intuition of a damper that gets stronger and stronger can be mathematically proved using matrix analysis. We then develop an algorithm for quadratic eigenvalue problems with low rank damping, which outperforms existing algorithm both in terms of speed and accuracy. The first part of our algorithm requires the solution of a generalized eigenproblem with semidefinite coefficient matrices. To solve this problem we develop a new algorithm based on an algorithm proposed by Wang and Zhao [SIAM J. Matrix Anal. Appl. 12-4 (1991), pp. 654--660]. The new algorithm computes all eigenvalues in a backward stable and symmetry preserving manner. The next two chapters are about equivalences of matrix polynomials. We show, for an algebraically closed field $\mathbb{F}$, that any matrix polynomial $P(\lambda)\in\mathbb{F}[\lambda]^{n\times m}$, $n\leq m$, can be reduced to triangular form, while preserving the degree and the finite and infinite elementary divisors. We then show that the same result holds for real matrix polynomials if we replace "triangular" with "quasi-triangular", that is, block-triangular with diagonal blocks of size $1\times 1$ and $2 \times 2$. The proofs are constructive in the sense that we build up triangular and quasi-triangular matrix polynomials starting the Smith form. In this sense we are solving structured inverse problems. In particular, our results imply that the necessary constraints that make list of elementary divisors admissible for a real square matrix polynomial of degree $\ell$ are also sufficient conditions. For the case of matrix polynomials with invertible leading coefficients, we show how triangular/quasi-triangular forms, as well as diagonal and Hessenberg forms, can be computed numerically. Finally, we present a backward error analysis of the shift-and-invert Arnoldi algorithm for matrices. This algorithm is also of interest to polynomial eigenproblems with easily constructible monic linearizations. The analysis shows how errors from the linear system solves and orthonormalization process affect the Arnoldi recurrence. Residual bounds for linear systems and columnwise backward error bounds for QR factorizations come to play, so we discuss these in some detail. The main result is a set of backward error bounds that can be estimated cheaply. We also use our error analysis to define a sensible condition for "breakdown", that is, a condition for when the Arnoldi iteration should be stopped.

This thesis considers Hermitian/symmetric, alternating and palindromic matrix polynomials which all arise frequently in a variety of applications, such as vibration analysis of dynamical systems and optimal control problems. A classification of Hermitian matrix polynomials whose eigenvalues belong to the extended real line, with each eigenvalue being of definite type, is provided first. We call such polynomials quasidefinite. Definite pencils, definitizable pencils, overdamped quadratics, gyroscopically stabilized quadratics, (quasi)hyperbolic and definite matrix polynomials are all quasidefinite. We show, using homogeneous rotations, special Hermitian linearizations and a new characterization of hyperbolic matrix polynomials, that the main common thread between these many subclasses is the distribution of their eigenvalue types. We also identify, amongst all quasihyperbolic matrix polynomials, those that can be diagonalized by a congruence transformation applied to a Hermitian linearization of the matrix polynomial while maintaining the structure of the linearization. Secondly, we generalize the notion of self-adjoint standard triples associated with Hermitian matrix polynomials in Gohberg, Lancaster and Rodman's theory of matrix polynomials to present spectral decompositions of structured matrix polynomials in terms of standard pairs (X,T), which are either real or complex, plus a parameter matrix S that acquires particular properties depending on the structure under investigation. These decompositions are mainly an extension of the Jordan canonical form for a matrix over the real or complex field so we investigate the important special case of structured Jordan triples. Finally, we use the concept of structured Jordan triples to solve a structured inverse polynomial eigenvalue problem. As a consequence, we can enlarge the collection of nonlinear eigenvalue problems [NLEVP, 2010] by generating quadratic and cubic quasidefinite matrix polynomials in different subclasses from some given spectral data by solving an appropriate inverse eigenvalue problem. For the quadratic case, we employ available algorithms to provide tridiagonal definite matrix polynomials.