Large-scale computation of pseudospectra using ARPACK and eigs

ARPACK and its MATLAB counterpart, eigs, are software packages that calculate some eigenvalues of a large non-symmetric matrix by Arnoldi iteration with implicit restarts. We show that at a small additional cost, which diminishes relatively as the matrix dimension increases, good estimates of pseudospectra in addition to eigenvalues can be obtained as a by-product. Thus in large-scale eigenvalue calculations it is feasible to obtain routinely not just eigenvalue approximations, but also information as to whether or not the eigenvalues are likely to be physically significant. Examples are presented for matrices with dimension up to 200,000

Theoretical error estimates for computing the matrix logarithm by Pad\'e-type approximants

In this article, we focus on the error that is committed when computing the matrix logarithm using the Gauss--Legendre quadrature rules. These formulas can be interpreted as Pad\'e approximants of a suitable Gauss hypergeometric function. Empirical observation tells us that the convergence of these quadratures becomes slow when the matrix is not close to the identity matrix, thus suggesting the usage of an inverse scaling and squaring approach for obtaining a matrix with this property. The novelty of this work is the introduction of error estimates that can be used to select a priori both the number of Legendre points needed to obtain a given accuracy and the number of inverse scaling and squaring to be performed. We include some numerical experiments to show the reliability of the estimates introduced

A proposal for pseudospectra computation on graphic processor units

Introducción. El cálculo del pseudoespectro de matrices es requerido en muchas aplicaciones modeladas por ecuaciones diferenciales y discretizadas en tiempo y espacio. Este cálculo resulta ser muy costoso computacionalmente, sobre todo para matrices de gran magnitud, para las cuales se han implementado con éxito métodos altamente paralelizables ejecutados en máquinas de alto rendimiento. Objetivo. Se presenta un análisis exploratorio del cálculo del pseudoespectro y su posible implementación en una arquitectura híbrida CPU-GPU, en la cual el cómputo masivo y paralelo se realice en las unidades de procesamiento gráfico. Metodología. Para formular la propuesta se analizan algunas implementaciones de este cálculo que han resultado efectivas en máquinas de alto rendimientos, el uso de métodos basados en métodos de Krylov, y las capacidades de las unidades de procesamiento gráfico en el cómputo masivo. Resultados.  La propuesta resulta de interés debido a que la mayor parte de este cálculo puede realizarse en unidades de procesamiento gráfico, que en la actualidad son fáciles de adquirir a bajo costo, y están incluidas en la mayoría de las marcas y modelos presentes en el mercado. Conclusiones. En este documento se describe un esquema general para la paralelización del cálculo del pseudoespectro en una arquitectura híbrida CPU-GPU.Introduction. Computation of matrix pseudospectra is required in many applications when modeling by differential equations. This computation is really expensive, especially for large matrices, for which highly parallelizable algorithms have been successfully implemented on high performance computers. Objective. We present an exploratory analysis of the pseudospectra computation in a hybrid architecture CPU-GPU where the graphics processing unit performs the massive parallel computation. Methodology. A proposal is formulated after analyzing some parallel implementations on high performance computers, methods based on Krylov methods, and the capacities of the graphics processor units for massive computation in the large-scale setting. Results. The proposal is attractive since the graphics processing units currently can be found on a wide range of computers, or can be adapted to any computer at a very a low cost. Conclusions. In this document we describe a general scheme for the parallel computation of pseudospectra on a hybrid architecture CPU-GPU

Path-Following Method to Determine the Field of Values of a Matrix with High Accuracy

We describe a novel and efficient algorithm for calculating the field of values boundary, $\partial\textrm{W}(\cdot)$, of an arbitrary complex square matrix: the boundary is described by a system of ordinary differential equations which are solved using Runge--Kutta (Dormand--Prince) numerical integration to obtain control points with derivatives then finally Hermite interpolation is applied to produce a dense output. The algorithm computes $\partial\textrm{W}(\cdot)$ both efficiently and with low error. Formal error bounds are proven for specific classes of matrix. Furthermore, we summarise the existing state of the art and make comparisons with the new algorithm. Finally, numerical experiments are performed to quantify the cost-error trade-off between the new algorithm and existing algorithms

Numerical Computation of Resonances and Pseudospectra in Acoustic Scattering

Acoustic scattering is a well-known physical phenomenon which arises in a wide range of fields: when acoustic waves propagating in a medium impinge on a localised non-uniformity, such as a density fluctuation or an external obstacle, their trajectories are deviated and scattered waves are generated. A key role in scattering theory is played by resonances; these are particular scatterer-dependent non-physical ‘complex’ frequencies at which acoustic scattering exhibits exceptional behaviour. The study of acoustic resonances for a particular scatterer provides an insight in the behaviour that the acoustic scattering assumes at the near physical ‘real’ frequencies, and it is a fundamental step in many applications. Yet, the numerical computation of resonances and pseudospectra - a mathematical tool which can be used to study the influence of resonances on physical frequencies - remains very expensive. With the present Thesis we want to address this particular problem, by proposing numerical algorithms based on the Boundary Element Method (BEM) for computing resonances and pseudospectra and by analysing their efficiency and performance. Finally, we apply such algorithms to half a dozen of physically relevant scatterer, inspired from different fields where acoustic scattering plays a relevant role

Berechnung singulärer Punkte nichtlinearer Gleichungssysteme

Diese Arbeit beschäftigt sich mit der Berechnung singulärer Punkte nichtlinearer Gleichungssysteme F(x)=0. Dazu werden minimal erweiterte Systeme der Form F(x)+D*s=0, f(x)=0 betrachtet. Die allgemeine Vorgehensweise zur Berechnung singulärer Punkte mit solchen erweiterten Systemen wird geschlossen dargestellt. Dazu werden zuerst die (teilweise verallgemeinerten Ljapunov-Schmidt-)reduzierten Funktionen von Golubitsky und Schaeffer, Beyn, Jepson und Spence, Griewank und Reddien, Kunkel bzw. Govaerts verallgemeinert und zusammengefasst. Es wird die verallgemeinerte Kontaktäquivalenz all dieser verallgemeinerten reduzierten Funktionen und die Gleichheit der benötigten Regularitätsannahmen bewiesen. Für eine weitere, neu eingeführte reduzierte Funktion wird die in dieser Arbeit definierte Ableitungsäquivalenz zu den anderen reduzierten Funktionen gezeigt. Mit dieser neuen reduzierten Funktion wird eine Reihe singulärer Punkte klassifiziert. Aus dieser Klassifikation ergeben sich Funktionen f aus Ableitungen der neuen reduzierten Funktion. Mit den so eingeführten Funktionen f kann das zweistufiges Newtonverfahren nach Pönisch und Schwetlick effektiv angewendet werden. Alle benötigten Ableitungen werden mittels Automatischer Differentiation bestimmt. Die numerischen Ergebnisse für eine Reihe von Beispielen zeigen die Effizienz dieses Verfahrens. Beim Newtonverfahren werden lineare Gleichungssysteme mit geränderten Matrizen B gelöst. Es wird gezeigt, für welche Ränderungen die Konditionszahl von B minimal ist. Dies ist z.B. für gewisse Vielfache der Singulärvektoren zu den kleinsten Singulärwerten der Fall. Zur Bestimmung dieser Ränderungen wird die inverse Teilraumiteration für Singulärwerte in verschiedenen Algorithmen angewendet. Die Konvergenzeigenschaften werden untersucht. Für einen Algorithmus wird bewiesen, dass die Konditionszahlen der iterierten geränderten Matrizen monoton fallen. Die numerischen Experimente bestätigen diese Eigenschaften

Wave radiation in simple geophysical models

Wave radiation is an important process in many geophysical flows. In particular, it is by wave radiation that flows may adjust to a state for which the dynamics is slow. Such a state is described as “balanced”, meaning there is an approximate balance between the Coriolis force and horizontal pressure gradients, and between buoyancy and vertical pressure gradients. In this thesis, wave radiation processes relevant to these enormously complex flows are studied through the use of some highly simplified models, and a parallel aim is to develop accurate numerical techniques for doing so. This thesis is divided into three main parts. 1. We consider accurate numerical boundary conditions for various equations which support wave radiation to infinity. Particular attention is given to discretely non-reflecting boundary conditions, which are derived directly from a discretised scheme. Such a boundary condition is studied in the case of the 1-d Klein-Gordon equation. The limitations concerning the practical implementation of this scheme are explored and some possible improvements are suggested. A stability analysis is developed which yields a simple stability criterion that is useful when tuning the boundary condition. The practical use of higher-order boundary conditions for the 2-d shallow water equations is also explored; the accuracy of such a method is assessed when combined with a particular interior scheme, and an analysis based on matrix pseudospectra reveals something of the stability of such a method. 2. Large-scale atmospheric and oceanic flows are examples of systems with a wide timescale separation, determined by a small parameter. In addition they both undergo constant random forcing. The five component Lorenz-Krishnamurthy system is a system with a timescale separation controlled by a small parameter, and we employ it as a model of the forced ocean by further adding a random forcing of the slow variables, and introduce wave radiation to infinity by the addition of a dispersive PDE. The dynamics are reduced by deriving balance relations, and numerical experiments are used to assess the effects of energy radiation by fast waves. 3. We study quasimodes, which demonstrate the existence of associated Landau poles of a system. In this thesis, we consider a simple model of wave radiation that exhibits quasimodes, that allows us to derive some explicit analytical results, as opposed to physically realistic geophysical fluid systems for which such results are often unavailable, necessitating recourse to numerical techniques. The growth rates obtained for this system, which is an extension of one considered by Lamb, are confirmed using numerical experiments

Low-rank methods for parameter-dependent eigenvalue problems and matrix equations

The focus of this thesis is on developing efficient algorithms for two important problems arising in model reduction, estimation of the smallest eigenvalue for a parameter-dependent Hermitian matrix and solving large-scale linear matrix equations, by extracting and exploiting underlying low-rank properties. Availability of reliable and efficient algorithms for estimating the smallest eigenvalue of a parameter-dependent Hermitian matrix $A(\mu)$ for many parameter values $\mu$ is important in a variety of applications. Most notably, it plays a crucial role in \textit{a posteriori} estimation of reduced basis methods for parametrized partial differential equations. We propose a novel subspace approach, which builds upon the current state-of-the-art approach, the Successive Constraint Method (SCM), and improves it by additionally incorporating the sampled smallest eigenvectors and implicitly exploiting their smoothness properties. Like SCM, our approach also provides rigorous lower and upper bounds for the smallest eigenvalues on the parameter domain $D$. We present theoretical and experimental evidence to demonstrate that our approach represents a significant improvement over SCM in the sense that the bounds are often much tighter, at a negligible additional cost. We have successfully applied the approach to computation of the coercivity and the inf-sup constants, as well as computation of $\varepsilon$-pseudospectra. Solving an $m\times n$ linear matrix equation $A_1 X B_1^T + \cdots + A_K X B_K^T = C$ as an $m n \times m n$ linear system, typically limits the feasible values of $m,n$ to a few hundreds at most. We propose a new approach, which exploits the fact that the solution $X$ can often be well approximated by a low-rank matrix, and computes it by combining greedy low-rank techniques with Galerkin projection as well as preconditioned gradients. This can be implemented in a way where only linear systems of size $m \times m$ and $n \times n$ need to be solved. Moreover, these linear systems inherit the sparsity of the coefficient matrices, which allows to address linear matrix equations as large as $m = n = O(10^5)$. Numerical experiments demonstrate that the proposed methods perform well for generalized Lyapunov equations, as well as for the standard Lyapunov equations. Finally, we combine the ideas used for addressing matrix equations and parameter-dependent eigenvalue problems, and propose a low-rank reduced basis approach for solving parameter-dependent Lyapunov equations