The Dynamic Mode Decomposition (DMD) is a method for computational analysis of nonlinear dynamical systems in data driven scenarios. Based on high fidelity numerical simulations or experimental data, the DMD can be used to reveal the latent structures in the dynamics or as a forecasting or a model order reduction tool. The theoretical underpinning of the DMD is the Koopman operator on a Hilbert space of observables of the dynamics under study. This paper describes a numerically robust and versatile variant of the DMD and its implementation using the state-of-the-art dense numerical linear algebra software package LAPACK . The features of the proposed software solution include residual bounds for the computed eigenpairs of the DMD matrix, eigenvectors refinements and computation of the eigenvectors of the Exact DMD, compressed DMD for efficient analysis of high dimensional problems that can be easily adapted for fast updates in a streaming DMD. Numerical analysis is the bedrock of numerical robustness and reliability of the software, that is tested following the highest standards and practices of LAPACK . Important numerical topics are discussed in detail and illustrated using numerous numerical examples.

The Dynamic Mode Decomposition (DMD) is a versatile and increasingly popular method for data driven analysis of dynamical systems that arise in a variety of applications in, e.g., computational fluid dynamics, robotics or machine learning. In the framework of numerical linear algebra, it is a data driven Rayleigh-Ritz procedure applied to a DMD matrix that is derived from the supplied data. In some applications, the physics of the underlying problem implies hermiticity of the DMD matrix, so the general DMD procedure is not computationally optimal. Furthermore, it does not guarantee important structural properties of the Hermitian eigenvalue problem and may return non-physical solutions. This paper proposes a software solution to the Hermitian (including the real symmetric) DMD matrices, accompanied with a numerical analysis that contains several fine and instructive numerical details. The eigenpairs are computed together with their residuals, and perturbation theory provides error bounds for the eigenvalues and eigenvectors. The software is developed and tested using the LAPACK package.

The quartic eigenvalue problem (λ 4 A +λ 3 B +λ 2 C +λ D + E ) x = 0 naturally arises in a plethora of applications, such as when solving the Orr–Sommerfeld equation in the stability analysis of the Poiseuille flow, in theoretical analysis and experimental design of locally resonant phononic plates, modeling a robot with electric motors in the joints, calibration of catadioptric vision system, or, for example, computation of the guided and leaky modes of a planar waveguide. This article proposes a new numerical method for the full solution (all eigenvalues and all left and right eigenvectors) that, starting with a suitable linearization, uses an initial, structure-preserving reduction designed to reveal and deflate a certain number of zero and infinite eigenvalues before the final linearization is forwarded to the QZ algorithm. The backward error in the reduction phase is bounded column wise in each coefficient matrix, which is advantageous if the coefficient matrices are graded. Numerical examples show that the proposed algorithm is capable of computing the eigenpairs with small residuals, and that it is competitive with the available state-of-the-art methods.

This article presents a new method for computing all eigenvalues and eigenvectors of quadratic matrix pencil Q (λ)=λ 2 M + λ C + K . It is an upgrade of the quadeig algorithm by Hammarlinget al., which attempts to reveal and remove by deflation a certain number of zero and infinite eigenvalues before QZ iterations. Proposed modifications of the quadeig framework are designed to enhance backward stability and to make the process of deflating infinite and zero eigenvalues more numerically robust. In particular, careful preprocessing allows scaling invariant/component-wise backward error and thus a better condition number. Further, using an upper triangular version of the Kronecker canonical form enables deflating additional infinite eigenvalues, in addition to those inferred from the rank of M . Theoretical analysis and empirical evidence from thorough testing of the software implementation confirm superior numerical performances of the proposed method.

A new software for computing the singular value decomposition (SVD) of real or complex matrices is proposed. The method implemented in the code xGESVDQ is essentially the QR SVD algorithm available as xGESVD in LAPACK. The novelty is an extra step, the QR factorization with column (or complete row and column) pivoting, also already available in LAPACK as xGEQP3. For experts in matrix computations, the combination of the QR factorization and an SVD computation routine is not new. However, what seems to be new and important for applications is that the resulting procedure is numerically superior to xGESVD and that it is capable of reaching the accuracy of the Jacobi SVD. Further, when combined with pivoted Cholesky factorization, xGESVDQ provides numerically accurate and fast solvers (designated as xPHEVC, xPSEVC) for the Hermitian positive definite eigenvalue problem. For instance, using accurately computed Cholesky factor, xPSEVC computes all eigenvalues of the 200 × 200 Hilbert matrix (whose spectral condition number is greater that 10 300 ) to nearly full machine precision. Furthermore, xGESVDQ can be used for accurate spectral decomposition of general (indefinite) Hermitian matrices.

This article proposes an efficient algorithm for reducing matrices to generalized Hessenberg form by unitary similarity, and recommends using it as a preprocessor in a variety of applications. To illustrate its usefulness, two cases from control theory are analyzed in detail: a solution procedure for a sequence of shifted linear systems with multiple right hand sides (e.g. evaluating the transfer function of a MIMO LTI dynamical system at many points) and computation of the staircase form. The proposed algorithm for the generalized Hessenberg reduction uses two levels of aggregation of Householder reflectors, thus allowing efficient BLAS 3-based computation. Another level of aggregation is introduced when solving many shifted systems by processing the shifts in batches. Numerical experiments confirm that the proposed methods have superior efficiency.

In this article, we propose a numerical algorithm for efficient and robust solution of a sequence of shifted Hessenberg linear systems. In particular, we show how the frequency response 𝒢(σ) = d - C ( A -σ I) -1 b in the single input case can be computed more efficiently than with other state-of-the-art methods. We also provide a backward stability analysis of the proposed algorithm.

This article reports an unexpected and rather erratic behavior of the LAPACK software implementation of the QR factorization with Businger-Golub column pivoting. It is shown that, due to finite precision arithmetic, the software implementation of the factorization can catastrophically fail to produce a properly structured triangular factor, thus leading to a potentially severe underestimate of a matrix's numerical rank. The 30-year old problem, dating back to LINPACK, has (undetectedly) badly affected many computational routines and software packages, as well as the study of rank-revealing QR factorizations. We combine computer experiments and numerical analysis to isolate, analyze, and fix the problem. Our modification of the current LAPACK xGEQP3 routine is already included in the LAPACK 3.1.0 release. The modified routine is numerically more robust and with a negligible overhead. We also provide a new, equally efficient, and provably numerically safe partial-column norm-updating strategy.