Search

Small- to medium-sized polynomial eigenvalue problems can be solved by linearizing the matrix polynomial and solving the resulting generalized eigenvalue problem using the QZ algorithm. The QZ algorithm, in turn, requires an initial reduction of a matrix pair to Hessenberg-triangular (HT) form. In this paper, we discuss the design and evaluation of high-performance parallel algorithms and software for HT reduction of a specific linearization of matrix polynomials of arbitrary degree. The proposed algorithm exploits the sparsity structure of the linearization to reduce the number of operations and improve the cache reuse compared to existing algorithms for unstructured inputs. Experiments on both a workstation and a high-performance computing system demonstrate that our structure-exploiting parallel implementation can outperform both the general LAPACK routine DGGHRD and the prototype implementation DGGHR3 of a general blocked algorithm.

This paper investigates the effect of structure-preserving perturbations on the eigenvalues of linearly and nonlinearly structured eigenvalue problems. Particular attention is paid to structures that form Jordan algebras, Lie algebras, and automorphism groups of a scalar product. Bounds and computable expressions for structured eigenvalue condition numbers are derived for these classes of matrices, which include complex symmetric, pseudo-symmetric, persymmetric, skew-symmetric, Hamiltonian, symplectic, and orthogonal matrices. In particular we show that under reasonable assumptions on the scalar product, the structured and unstructured eigenvalue condition numbers are equal for structures in Jordan algebras. For Lie algebras, the effect on the condition number of incorporating structure varies greatly with the structure. We identify Lie algebras for which structure does not affect the eigenvalue condition number.

We develop algorithms that construct robust (i.e., reliable for a given tolerance, and scaling independent) rational approximants of matrix-valued functions on a given subset of the complex plane. We consider matrix-valued functions provided in both split form (i.e., as a sum of scalar functions times constant coefficient matrices) and in black box form. We develop a new error analysis and use it for the construction of stopping criteria, one for each form. Our criterion for split forms adds weights chosen relative to the importance of each scalar function, leading to the weighted AAA algorithm, a variant of the set-valued AAA algorithm that can guarantee to return a rational approximant with a user-chosen accuracy. We propose two-phase approaches for black box matrix-valued functions that construct a surrogate AAA approximation in phase one and refine it in phase two, leading to the surrogate AAA algorithm with exact search and the surrogate AAA algorithm with cyclic Leja--Bagby refinement. The stopping criterion for black box matrix-valued functions is updated at each step of phase two to include information from the previous step. When convergence occurs, our two-phase approaches return rational approximants with a user-chosen accuracy. We select problems from the NLEVP collection that represent a variety of matrix-valued functions of different sizes and properties and use them to benchmark our algorithms.

The shift-and-invert Lanczos algorithm is a commonly used solution procedure to compute the eigenpairs of large, sparse eigenvalue problems that arise when approximating the elastic dynamic response of large structures under the influence of seismic forces. Not all eigenvectors are equally important to the response when the structure is exposed to a mass-dependent external force of the form $g(t) Mb$, where $M$ is the mass matrix of the system and $b$ the rigid body vector. Structural engineers select eigenvectors $x_j$, $j=1,\ldots,\ell$, such that their cumulative mass participation, measured as $\sum_{j=1}^\ell (x_j^T Mb)^2/(b^TMb)$,is above a target threshold $\xi$. We show that when the starting vector for the unshifted Lanczos algorithm is the spatial distribution vector $b$, the Lanczos procedure can be used to provide an estimate of the cumulative mass participation. This allows us to identify intervals containing eigenvalues whose eigenvectors have a large contribution to the cumulative mass participation and filter out intervals containing eigenvalues whose eigenvectors have a negligible contribution. We use this information to devise a sequence of shifts $\sigma_1, \ldots,\sigma_p$ for the shift-and-invert Lanczos algorithm as well as a stopping criterion for the iteration with shift $\sigma_i$ so that the cumulative mass participation of the computed eigenvectors reaches the required level $\xi$. Numerical experiments on real engineering problems show that our approach computes up to $70\%$ fewer eigenvectors and requires fewer shifts, on average, than the more general shifting strategy proposed by Ericsson and Ruhe (Math. Comp., 35 (1980)) together with its modification presented in Grimes et al. (SIAM J. Matrix Anal. and Appl. 40(4), 1994).

Given a matrix polynomial $A(\lambda)$ of degree $d$ and the associated vector space of pencils $\DLP(A)$ described in Mackey, Mackey, Mehl, and Mehrmann [SIAM J. Matrix Anal. Appl., 28 (2006), pp. 971-1004], we construct a parametrization for the set of left and right transformations that preserve the block structure of such pencils, and hence produce a new matrix polynomial $\At(\lambda)$ that is still of degree $d$ and is unimodularly equivalent to $A(\lambda)$. We refer to such left and right transformations as structure-preserving transformations (SPTs). Unlike previous work on SPTs, we do not require the leading matrix coefficient of $A(\lambda)$ to be nonsingular. We show that additional constraints on the parametrization lead to SPTs that also preserve extra structures in $A(\lambda)$ such as symmetric, alternating, and $T$-palindromic structures. Our parametrization allows easy construction of SPTs that are low-rank modifications of the identity matrix. The latter transform $A(\lambda)$ into an equivalent matrix polynomial $\At(\lambda)$ whose $j$th matrix coefficient $\At_j$ is a low-rank modification of $A_j$. We expect such SPTs to be one of the key tools for developing algorithms that reduce a matrix polynomial to Hessenberg form or tridiagonal form in a finite number of steps and without the use of a linearization.

We develop algorithms that construct robust (i.e., reliable for a given tolerance and scaling independent) rational approximations of matrix-valued functions on a given subset of the complex plane. We consider matrix-valued functions provided in both split form (i.e., as a sum of scalar functions times constant coefficient matrices) and in black box form. We use an error analysis to construct stopping criteria, one for each form. The criterion for split forms adds weights chosen relative to the importance of each scalar function, which we use in our weighted AAA algorithm, a variant of the set-valued AAA algorithm that guarantees to return a rational approximant with a user chosen accuracy. We propose two-phase approaches for black box matrix-valued functions that construct a surrogate AAA approximation in phase one and refine it in phase two, leading to the surrogate AAA algorithm with exact search and the surrogate AAA algorithm with cyclic Leja--Bagby refinement. The stopping criterion for black box matrix-valued functions is updated at each step of phase two to include information from the previous step. When convergence occurs, our two-phase approaches return rational approximants with a user chosen accuracy. We select problems from the NLEVP collection that represent a variety of matrix-valued functions of different sizes and properties and use them to benchmark our algorithms.

This is the description of the new changes made in version 4.0 of the NLEVP MATLAB toolbox. The collection and its organization are described in a separate paper. A user's guide describes how to install and use the NLEVP MATLAB toolbox.

Bilateral cruise control (BCC) suppresses traffic flow instabilities. Previously, for simplicity of analysis, vehicles in BCC traffic flow were assumed to be identical, i.e., using the same gains for control. In this study, we analyze the stability of an inhomogeneous vehicular chain in which the gains used by different vehicles are not the same. Not unexpectedly, mathematical analysis becomes more difficult, and leads to a quadratic eigenvalue problem. We study several different cases, and shows that a chain of vehicles under bilateral cruise control is stable even when the vehicles do not all have the same control system properties. Numerical simulations validate the analysis.

The shift-and-invert Lanczos algorithm is a commonly used solution procedure to compute the eigenpairs of large, sparse eigenvalue problems that arise when approximating the elastic dynamic response of large structures under the influence of seismic forces. Not all eigenvectors are equally important to the response when the structure is exposed to a mass-dependent external force of the form $g(t) Mb$, where $M$ is the mass matrix of the system and $b$ the rigid body vector. Structural engineers select eigenvectors $x_j$, $j=1,\ldots,\ell$, such that their cumulative mass participation, measured as $\sum_{j=1}^\ell (x_j^T Mb)^2/(b^TMb)$, is above a target threshold $\xi$. We show that when the starting vector for the unshifted Lanczos algorithm is the spatial distribution vector $b$, the Lanczos procedure can be used to provide an estimate of the cumulative mass participation. This allows us to identify intervals containing eigenvalues whose eigenvectors have a large contribution to the cumulative mass participation and filter out intervals containing eigenvalues whose eigenvectors have a negligible contribution. We use this information to devise a sequence of shifts $\sigma_1, \ldots,\sigma_p$ for the shift-and-invert Lanczos algorithm as well as a stopping criterion for the iteration with shift $\sigma_i$ so that the cumulative mass participation of the computed eigenvectors reaches the required level $\xi$. Numerical experiments on real engineering problems show that our approach computes up to $80\%$ fewer eigenvectors and requires fewer shifts, on average, than the more general shifting strategy proposed by Ericsson and Ruhe (Math. Comp., 35 (1980)).

A Hungarian scaling is a diagonal scaling of a matrix that is typically applied along with a permutation to a sparse symmetric or nonsymmetric indefinite linear system before calling a direct or iterative solver. A Hungarian scaled and reordered matrix has all its entries of modulus less than or equal to 1 and entries of modulus 1 on the diagonal. We use max-plus algebra to characterize the set of all Hungarian scalings for a given matrix and show that max-balancing a Hungarian scaled matrix yields the most ``diagonally dominant" Hungarian scaled matrix possible with respect to some ordering. We also propose a new scaling, called centre of mass scaling, which can be seen as an approximate max-balancing Hungarian scaling and whose computation is embarrassingly parallel. Numerical experiments illustrate the increased diagonal dominance produced by max-balancing and centre of mass scaling of Hungarian scaled matrices as well as the reduced need for pivoting in Gaussian elimination with partial pivoting and the improved stability of LU factorizations without pivoting.

We present a new method for constructing incomplete Cholesky factorization preconditioners for use in solving large sparse symmetric positive-definite linear systems. This method uses max-plus algebra to predict the positions of the largest entries in the Cholesky factor and then uses these positions as the sparsity pattern for the preconditioner. Our method builds on the max-plus incomplete LU factorization preconditioner recently proposed in [J. Hook and F. Tisseur, Incomplete LU preconditioner based on max-plus approximation of LU factorization, MIMS Eprint 2016.47, Manchester, 2016] but applied to symmetric positive-definite matrices, which comprise an important special case for the method and its application. A attractive feature of our approach is that the sparsity pattern of each column of the preconditioner can be computed in parallel. Numerical comparisons are made with other incomplete Cholesky factorization preconditioners using problems from a range of practical applications. We demonstrate that the new preconditioner can outperform traditional level-based preconditioners and offer a parallel alternative to a serial limited-memory based approach.

We present an efficient block preconditioner for the two-dimensional biharmonic Dirichlet problem discretised by C1 bicubic Hermite finite elements. In this formulation each node in the mesh has four different degrees of freedom (DOFs). Grouping DOFs of the same type together leads to a natural blocking of the Galerkin coefficient matrix. Based on this block structure, we develop two preconditioners: a 2x2 block diagonal preconditioner (BD) and a block bordered diagonal (BBD) preconditioner. We prove mesh independent bounds for the spectra of the BD-preconditioned Galerkin matrix under certain conditions. The eigenvalue analysis is based on the fact that the proposed preconditioner, like the coefficient matrix itself, is symmetric positive definite and is assembled from element matrices. We demonstrate the effectiveness of an inexact version of the BBD preconditioner, which exhibits near optimal scaling in terms of computational cost with respect to the discrete problem size. Finally, we study robustness of this preconditioner with respect to element stretching, domain distortion and non-convex domains.

The tropical roots of $tp(x) = \max_{0\le j\le d}\|A_j\|x^j$ are points at which the maximum is attained at least twice. These roots, which can be computed in only $O(d)$ operations, can be good approximations to the moduli of the eigenvalues of the matrix polynomial $P(\lambda)=\sum_{j=0}^d \lambda^j A_j$, in particular when the norms of the matrices $A_j$ vary widely. Our aim is to investigate this observation and its applications. We start by providing annuli defined in terms of the tropical roots of $tp(x)$ that contain the eigenvalues of $P(\lambda)$. Our localization results yield conditions under which tropical roots offer order of magnitude approximations to the moduli of the eigenvalues of $P(\lambda)$. Our tropical localization of eigenvalues are less tight than eigenvalue localization results derived from a generalized matrix version of Pellet's theorem but they are easier to interpret. Tropical roots are already used to determine the starting points for matrix polynomial eigensolvers based on scalar polynomial root solvers such as the Ehrlich-Aberth method and our results further justify this choice. Our results provide the basis for analyzing the effect of Gaubert and Sharify's tropical scalings for $P(\lambda)$ on (a) the conditioning of linearizations of tropically scaled $P(\lambda)$ and (b) the backward stability of eigensolvers based on linearizations of tropically scaled $P(\lambda)$. We anticipate that the tropical roots of $tp(x)$, on which the tropical scalings are based, will help designing polynomial eigensolvers with better numerical properties than standard algorithms for polynomial eigenvalue problems such as that implemented in the MATLAB function \texttt{polyeig}.

We propose a novel approach to solve polynomial eigenvalue problems via linearization. The novelty lies in (a) our choice of linearization, which is constructed using input from tropical algebra and the notion of well-separated tropical roots, (b) an appropriate scaling applied to the linearization and (c) a modified stopping criterion for the QZ iterations that takes advantage of the properties of our scaled linearization. Numerical experiments show that our polynomial eigensolver computes all the finite and well-conditioned eigenvalues to high relative accuracy even when they are very different in magnitude.

We study the structured condition number of differentiable maps between smooth matrix manifolds, developing a theoretical framework that extends previous results for vector subspaces to any smooth manifold. We present algorithms to compute the structured condition number. As special cases of smooth manifolds, we analyze automorphism groups, and Lie and Jordan algebras associated with a scalar product. For such manifolds, we derive a lower bound on the structured condition number that is cheaper to compute than the structured condition number. We provide numerical comparisons between the structured and unstructured condition numbers for the principal matrix logarithm and principal matrix square root of matrices in automorphism groups as well as for the map between matrices in automorphism groups and their polar decomposition. We show that our lower bound can be used as a good estimate for the structured condition number when the matrix argument is well conditioned. We show that the structured and unstructured condition numbers can differ by many orders of magnitude, thus motivating the development of algorithms preserving structure.

A Hungarian scaling is a diagonal scaling of a matrix that is typically applied along with a permutation to a sparse linear system before calling a direct or iterative solver. A matrix that has been Hungarian scaled and reordered has all entries of modulus less than or equal to 1 and entries of modulus 1 on the diagonal. An important fact that has been overlooked by the previous research into Hungarian scaling of linear systems is that a given matrix typically has a range of possible Hungarian scalings and direct or iterative solvers may behave quite differently under each of these scalings. Since standard algorithms for computing Hungarian scalings return only one scaling, it is natural to ask whether a superior performing scaling can be obtained by searching within the set of all possible Hungarian scalings. To this end we propose a method for computing a Hungarian scaling that is optimal from the point of view of diagonal dominance. Our method uses max-balancing, which minimizes the largest off-diagonal entries in the matrix. Numerical experiments illustrate the increased diagonal dominance produced by max-balanced Hungarian scaling as well as the reduced need for row interchanges in Gaussian elimination with partial pivoting and the improved stability of LU factorizations without pivoting. We additionally find that applying the max-balancing scaling before computing incomplete LU preconditioners improves the convergence rate of certain iterative methods.

A square matrix can be reduced to simpler form via similarity transformations. Here ``simpler form'' may refer to diagonal (when possible), triangular (Schur), or Hessenberg form. Similar reductions exist for matrix pencils if we consider general equivalence transformations instead of similarity transformations. For both matrices and matrix pencils, well-established algorithms are available for each reduction, which are useful in various applications. For matrix polynomials, unimodular transformations can be used to achieve the reduced forms but we do not have a practical way to compute them. In this work we introduce a practical means to reduce a matrix polynomial with nonsingular leading coefficient to a simpler (diagonal, triangular, Hessenberg) form while preserving the degree and the eigenstructure. The key to our approach is to work with structure preserving similarity transformations applied to a linearization of the matrix polynomial instead of unimodular transformations applied directly to the matrix polynomial. As an applications, we illustrate how to use these reduced forms to solve parameterized linear systems.

Small- to medium-sized polynomial eigenvalue problems can be solved by linearizing the matrix polynomial and solving the resulting generalized eigenvalue problem using the QZ algorithm. The QZ algorithm, in turn, requires an initial reduction of a matrix pair to Hessenberg--triangular form. In this paper, we discuss the design and evaluation of high-performance parallel algorithms and software for Hessenberg--triangular reduction of a specific linearization of matrix polynomials of arbitrary degree. The proposed algorithm exploits the sparsity structure of the linearization to reduce the number of operations and improve the cache reuse compared to existing algorithms for unstructured inputs.Experiments on both a workstation and an HPC system demonstrate that our structure-exploiting parallel implementation can outperform both the general LAPACK routine DGGHRD and the prototype implementation DGGHR3 of a general blocked algorithm.

We present a new method for the a priori approximation of the orders of magnitude of the entries in the LU factors of a complex or real matrix $A$. This approximation can be used to quickly determine the positions of the largest entries in the LU factors of $A$ and these positions can then be used as the sparsity pattern for an incomplete LU factorization preconditioner. Our method uses max-plus algebra and is based solely on the moduli of the entries of $A$. We also present techniques for predicting which permutation matrices will be chosen by Gaussian elimination with partial pivoting. We exploit the strong connection between the field of Puiseux series and the max-plus semiring to prove properties of the max-plus LU factors. Experiments with a set of test matrices from the University of Florida sparse matrix collection show that our max-plus LU preconditioners outperform traditional level of fill methods and have similar performance to those preconditioners computed with more expensive threshold-based methods.

We present a new method for the a priori approximation of the orders of magnitude of the entries in the LU factors of a complex or real matrix $A$. This approximation can be used to quickly determine the positions of the largest entries in the LU factors of $A$ and these positions can then be used as the sparsity pattern for an incomplete LU factorization preconditioner. Our method uses max-plus algebra and is based solely on the moduli of the entries of $A$. We also present techniques for predicting which permutation matrices will be chosen by Gaussian elimination with partial pivoting. We exploit the strong connection between the field of Puiseux series and the max-plus semiring to prove properties of the max-plus LU factors. Experiments with a set of test matrices from the University of Florida sparse matrix collection show that our max-plus LU preconditioners outperform traditional level of fill methods and have similar performance to those preconditioners computed with more expensive threshold-based methods.

We present an efficient block preconditioner for the two-dimensional biharmonic Dirichlet problem discretised by C1 bicubic Hermite finite elements. In this formulation each node in the mesh has four different degrees of freedom (DOFs). Grouping DOFs of the same type together leads to a natural blocking of the Galerkin coefficient matrix. Based on this block structure, we develop two preconditioners: a 2x2 block diagonal preconditioner (BD) and a block bordered diagonal (BBD) preconditioner. We prove mesh independent bounds for the spectra of the BD-preconditioned Galerkin matrix under certain conditions. The eigenvalue analysis is based on the fact that the proposed preconditioner, like the coefficient matrix itself, is symmetric positive definite and is assembled from element matrices. We demonstrate the effectiveness of an inexact version of the BBD preconditioner, which exhibits near optimal scaling in terms of computational cost with respect to the discrete problem size. Finally, we study robustness of this preconditioner with respect to element stretching, domain distortion and non-convex domains.

The \textit{tropical scaling} algorithm has experimentally shown to generate accurate results in computing the eigenpairs of a matrix polynomial $P(\l)=\sum_{k=1}^\d \l^k A_k$. This algorithm scales $P(\l)$ by using certain quantities known as \textit{tropical roots}, then, for each tropical roots, it constructs a linearization of the tropically scaled polynomial, computes its eigenpairs and extracts eigenpairs of $P(\l)$ from those of the linearizations. In this work we analyse this algorithm in terms of backward error and conditioning. We show that when the tropical roots are well separated, for an eigenpair, the backward error of the scale matrix polynomial is bounded by the backward error of its linearization times the condition number of certain coefficients of $P(\l)$. These coefficients determine the abscissae of the nodes of the {\it Newton polygon} associated with the matrix polynomial. Similar results show that for an eigenvalue, the condition number of the linearization of a scaled matrix polynomial is bounded by the condition number of the scaled matrix polynomial times the condition number of these coefficients. Our results show that when these coefficients are well conditioned the eigenpairs can be computed without any numerical difficulty. These analysis is supported by the experiments which are provided at the end of the paper.