A novel iterative method to approximate structured singular values

A novel method for approximating structured singular values (also known as mu-values) is proposed and investigated. These quantities constitute an important tool in the stability analysis of uncertain linear control systems as well as in structured eigenvalue perturbation theory. Our approach consists of an inner-outer iteration. In the outer iteration, a Newton method is used to adjust the perturbation level. The inner iteration solves a gradient system associated with an optimization problem on the manifold induced by the structure. Numerical results and comparison with the well-known Matlab function mussv, implemented in the Matlab Control Toolbox, illustrate the behavior of the method

Computing the structured pseudospectrum of a Toeplitz matrix and its extreme points

The computation of the structured pseudospectral abscissa and radius (with respect to the Frobenius norm) of a Toeplitz matrix is discussed and two algorithms based on a low rank property to construct extremal perturbations are presented. The algorithms are inspired by those considered in [SIAM J. Matrix Anal. Appl., 32 (2011), pp. 1166-1192] for the unstructured case, but their extension to structured pseudospectra and analysis presents several difficulties. Natural generalizations of the algorithms, allowing to draw significant sections of the structured pseudospectra in proximity of extremal points are also discussed. Since no algorithms are available in the literature to draw such structured pseudospectra, the approach we present seems promising to extend existing software tools (Eigtool, Seigtool) to structured pseudospectra representation for Toeplitz matrices. We discuss local convergence properties of the algorithms and show some applications to a few illustrative examples.Comment: 21 pages, 11 figure

Transient dynamics under structured perturbations: bridging unstructured and structured pseudospectra

The structured $\varepsilon$-stability radius is introduced as a quantity to assess the robustness of transient bounds of solutions to linear differential equations under structured perturbations of the matrix. This applies to general linear structures such as complex or real matrices with a given sparsity pattern or with restricted range and corange, or special classes such as Toeplitz matrices. The notion conceptually combines unstructured and structured pseudospectra in a joint pseudospectrum, allowing for the use of resolvent bounds as with unstructured pseudospectra and for structured perturbations as with structured pseudospectra. We propose and study an algorithm for computing the structured $\varepsilon$-stability radius. This algorithm solves eigenvalue optimization problems via suitably discretized rank-1 matrix differential equations that originate from a gradient system. The proposed algorithm has essentially the same computational cost as the known rank-1 algorithms for computing unstructured and structured stability radii. Numerical experiments illustrate the behavior of the algorithm

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

Dynamically Correct Formulations of the Linearised Navier-Stokes Equations

Motivated by the need to efficiently obtain low-order models of fluid flows around complex geometries for the purpose of feedback control system design, this paper considers the effect on system dynamics of basing plant models on different formulations of the linearised Navier-Stokes equations. We consider the dynamics of a single computational node formed by spatial discretisation of the governing equations in both primitive variables (momentum equation & continuity equation) and pressure Poisson equation (PPE) formulations. This reveals fundamental numerical differences at the nodal level, whose effects on the system dynamics at the full system level are exemplified by considering the corresponding formulations of a two-dimensional (2D) channel flow, subjected to a variety of different boundary conditions

Approximating the Real Structured Stability Radius with Frobenius Norm Bounded Perturbations

We propose a fast method to approximate the real stability radius of a linear dynamical system with output feedback, where the perturbations are restricted to be real valued and bounded with respect to the Frobenius norm. Our work builds on a number of scalable algorithms that have been proposed in recent years, ranging from methods that approximate the complex or real pseudospectral abscissa and radius of large sparse matrices (and generalizations of these methods for pseudospectra to spectral value sets) to algorithms for approximating the complex stability radius (the reciprocal of the $H_\infty$ norm). Although our algorithm is guaranteed to find only upper bounds to the real stability radius, it seems quite effective in practice. As far as we know, this is the first algorithm that addresses the Frobenius-norm version of this problem. Because the cost mainly consists of computing the eigenvalue with maximal real part for continuous-time systems (or modulus for discrete-time systems) of a sequence of matrices, our algorithm remains very efficient for large-scale systems provided that the system matrices are sparse

The Foundations of Infinite-Dimensional Spectral Computations

Spectral computations in infinite dimensions are ubiquitous in the sciences. However, their many applications and theoretical studies depend on computations which are infamously difficult. This thesis, therefore, addresses the broad question, “What is computationally possible within the field of spectral theory of separable Hilbert spaces?” The boundaries of what computers can achieve in computational spectral theory and mathematical physics are unknown, leaving many open questions that have been unsolved for decades. This thesis provides solutions to several such long-standing problems. To determine these boundaries, we use the Solvability Complexity Index (SCI) hierarchy, an idea which has its roots in Smale's comprehensive programme on the foundations of computational mathematics. The Smale programme led to a real-number counterpart of the Turing machine, yet left a substantial gap between theory and practice. The SCI hierarchy encompasses both these models and provides universal bounds on what is computationally possible. What makes spectral problems particularly delicate is that many of the problems can only be computed by using several limits, a phenomenon also shared in the foundations of polynomial root-finding as shown by McMullen. We develop and extend the SCI hierarchy to prove optimality of algorithms and construct a myriad of different methods for infinite-dimensional spectral problems, solving many computational spectral problems for the first time. For arguably almost any operator of applicable interest, we solve the long-standing computational spectral problem and construct algorithms that compute spectra with error control. This is done for partial differential operators with coefficients of locally bounded total variation and also for discrete infinite matrix operators. We also show how to compute spectral measures of normal operators (when the spectrum is a subset of a regular enough Jordan curve), including spectral measures of classes of self-adjoint operators with error control and the construction of high-order rational kernel methods. We classify the problems of computing measures, measure decompositions, types of spectra (pure point, absolutely continuous, singular continuous), functional calculus, and Radon--Nikodym derivatives in the SCI hierarchy. We construct algorithms for and classify; fractal dimensions of spectra, Lebesgue measures of spectra, spectral gaps, discrete spectra, eigenvalue multiplicities, capacity, different spectral radii and the problem of detecting algorithmic failure of previous methods (finite section method). The infinite-dimensional QR algorithm is also analysed, recovering extremal parts of spectra, corresponding eigenvectors, and invariant subspaces, with convergence rates and error control. Finally, we analyse pseudospectra of pseudoergodic operators (a generalisation of random operators) on vector-valued $l^p$ spaces. All of the algorithms developed in this thesis are sharp in the sense of the SCI hierarchy. In other words, we prove that they are optimal, realising the boundaries of what digital computers can achieve. They are also implementable and practical, and the majority are parallelisable. Extensive numerical examples are given throughout, demonstrating efficiency and tackling difficult problems taken from mathematics and also physical applications. In summary, this thesis allows scientists to rigorously and efficiently compute many spectral properties for the first time. The framework provided by this thesis also encompasses a vast number of areas in computational mathematics, including the classical problem of polynomial root-finding, as well as optimisation, neural networks, PDEs and computer-assisted proofs. This framework will be explored in the future work of the author within these settings

Structure-Preserving Model Reduction of Physical Network Systems

This paper considers physical network systems where the energy storage is naturally associated to the nodes of the graph, while the edges of the graph correspond to static couplings. The first sections deal with the linear case, covering examples such as mass-damper and hydraulic systems, which have a structure that is similar to symmetric consensus dynamics. The last section is concerned with a specific class of nonlinear physical network systems; namely detailed-balanced chemical reaction networks governed by mass action kinetics. In both cases, linear and nonlinear, the structure of the dynamics is similar, and is based on a weighted Laplacian matrix, together with an energy function capturing the energy storage at the nodes. We discuss two methods for structure-preserving model reduction. The first one is clustering; aggregating the nodes of the underlying graph to obtain a reduced graph. The second approach is based on neglecting the energy storage at some of the nodes, and subsequently eliminating those nodes (called Kron reduction).</p