Numerical Linear Algebra Techniques Research Articles

Control of port-Hamiltonian differential-algebraic systems and applications

We discuss the modelling framework of port-Hamiltonian descriptor systems and their use in numerical simulation and control. The structure is ideal for automated network-based modelling since it is invariant under power-conserving interconnection, congruence transformations and Galerkin projection. Moreover, stability and passivity properties are easily shown. Condensed forms under orthogonal transformations present easy analysis tools for existence, uniqueness, regularity and numerical methods to check these properties.After recalling the concepts for general linear and nonlinear descriptor systems, we demonstrate that many difficulties that arise in general descriptor systems can be easily overcome within the port-Hamiltonian framework. The properties of port-Hamiltonian descriptor systems are analysed, and time discretization and numerical linear algebra techniques are discussed. Structure-preserving regularization procedures for descriptor systems are presented to make them suitable for simulation and control. Model reduction techniques that preserve the structure and stabilization and optimal control techniques are discussed.The properties of port-Hamiltonian descriptor systems and their use in modelling simulation and control methods are illustrated with several examples from different physical domains. The survey concludes with open problems and research topics that deserve further attention.

Survey and Review

Gaussian (also known as normal) probability distributions in (\mathbb R^d\) play a central role in statistics and are important in many branches of pure and applied mathematics, statistical physics, and several other fields. In the univariate (d=1) case, the corresponding bell curve is known even to many people outside the sciences. At present, there is much interest in algorithms that generate random vectors (X) in (\mathbb R^d) distributed according to a target Gaussian density function. In the scalar (d=1) case, or when (d>1) but the scalar components of (X) are not correlated, there are several simple efficient algorithms, including the well-known Box--Muller method. In the general case, where (d>1) and there are correlations between the scalar components, the standard approach is to perform a linear change of variables (X = LY) so that the components (Y) are uncorrelated and may be sampled easily. Typically, the matrix (L) is obtained via a Cholesky factorization of the precision matrix of the target distribution, an approach which works well if (d) is not too large (for a laptop if, say, (d < 10^5)). However, many recent applications, including image analysis, spatial statistics, graphical structures, and others, operate with very large values of (d) and the Cholesky algorithm, with an (\mathcalO(d^3)) complexity and (\mathcalO(d^2)) memory requirements, may not be feasible. The Survey and Review paper in this issue, “High-Dimensional Gaussian Sampling: A Review and a Unifying Approach Based on a Stochastic Proximal Point Algorithm” by Maxime Vono, Nicolas Dobigeon, and Pierre Chainais, addresses the problem of obtaining Gaussian samples when the Cholesky factorization is not an option. The paper presents in a unified way and compares many algorithms suggested in different scientific communities. The algorithms may be grouped into two classes. In the first, numerical linear algebra techniques are used to reduce the computational complexity and/or the memory requirements. In the second, the samples are obtained via Markov chain Monte Carlo approaches; surprisingly, the resulting algorithms are very much related to classical iterative methods for the solution of linear systems, including Jacobi, Gauss--Seidel, and SOR. For this reason the paper will be relevant to readers interested in numerical linear algebra. For those whose work requires generating random samples, the paper includes a neat decision tree to choose, in a given application, among the many available algorithms. The authors have also made available software for all the methods considered in their survey.

Topical Issue Applied and Numerical Linear Algebra (2/2)

The present special issue of the GAMM Mitteilungen, which is the second of a two-part series, contains contributions on the topic of Applied and Numerical Linear Algebra, compiled by the GAMM Activity Group of the same name. The Activity Group has already contributed special issues to the GAMM Mitteilungen in 2004, 2006, and 2013. Because of the rapid development both in the theoretical foundations and the applicability of numerical linear algebra techniques throughout science and engineering, it is time again to survey the field and present the results to the readers of the GAMM Mitteilungen. We are happy that eight authors or teams of authors have accepted our invitation to report on recent research highlights in Applied Numerical Linear Algebra, and to point out the relevant literature as well as software. This work by Federico Poloni reviews a family of algorithms for Lyapunov- and Riccati-type equations which are all related by the idea of doubling. The algorithms are compared and their connections are highlighted. The paper also discusses open problems relating to their theory.

A Scalable Algorithm for Large-Scale Unsupervised Multi-View Partial Least Squares

We present an unsupervised multi-view partial least squares (PLS) by learning a common latent space from given multi-view data. Although PLS is a frequently used technique for analyzing relationships between two datasets, its extension to more than two views in unsupervised setting is seldom studied. In this article, we fill up the gap, and our model bears similarity to the extension of canonical correlation analysis (CCA) to more than two sets of variables and is built on the findings from analyzing PLS, CCA, and its variants. The resulting problem involves a set of orthogonality constraints on view-specific projection matrices, and is numerically challenging to existing methods that may have numerical instabilities and offer no orthogonality guarantee on view-specific projection matrices. To solve this problem, we propose a stable deflation algorithm that relies on proven numerical linear algebra techniques, can guarantee the orthogonality constraints, and simultaneously maximizes the covariance in the common space. We further adapt our algorithm to efficiently handle large-scale high-dimensional data. Extensive experiments have been conducted to evaluate the algorithm through performing two learning tasks, cross-modal retrieval, and multi-view feature extraction. The results demonstrate that the proposed algorithm outperforms the baselines and is scalable for large-scale high-dimensional datasets.

Finite convergence of an active signature method to local minima of piecewise linear functions

ABSTRACTWe previously derived first-order (KKT) and second-order (SOSC) optimality conditions for functions defined by evaluation programs involving smooth elementals and absolute values. For this class of problems we showed that the algorithm of successive abs-linear minimization with a proximal term (SALMIN) achieves a linear or even quadratic rate of convergence under suitable assumptions. An implementation of SALMIN called LiPsMin has been implemented and tested. It uses a bundle-type method direction finder for the inner loop, i.e. the minimization of the local piecewise linear model with a quadratic proximal term. As a consequence the inner solver can get stuck at stationary points of the piecewise linearization and the overall algorithm can only guarantee Clarke stationarity of its cluster points. Here we present a new method that computes a local minimizer of the proximal model objective. As a consequence all cluster points of the outer iteration must be first-order minimal, which is also known as criticality in nonsmooth optimization. The active signature strategy employed very much resembles the classical method for convex QOPs and utilizes the same numerical linear algebra techniques. The new custom-made algorithm provides opportunities for structure exploitation like warm starts in the context of the nonlinear, outer loop.

A stabilized normal form algorithm for generic systems of polynomial equations

We propose a numerical linear algebra based method to find the multiplication operators of the quotient ring C [ x ] ∕ I associated to a zero-dimensional ideal I generated by n C -polynomials in n variables. We assume that the polynomials are generic in the sense that the number of solutions in C n equals the Bézout number. The main contribution of this paper is an automated choice of basis for C [ x ] ∕ I , which is crucial for the feasibility of normal form methods in finite precision arithmetic. This choice is based on numerical linear algebra techniques and it depends on the given generators of I .

A comparative study on low-memory iterative solvers for FFT-based homogenization of periodic media

In this paper, we assess the performance of four iterative algorithms for solving non-symmetric rank-deficient linear systems arising in the FFT-based homogenization of heterogeneous materials defined by digital images. Our framework is based on the Fourier–Galerkin method with exact and approximate integrations that has recently been shown to generalize the Lippmann–Schwinger setting of the original work by Moulinec and Suquet from 1994. It follows from this variational format that the ensuing system of linear equations can be solved by general-purpose iterative algorithms for symmetric positive-definite systems, such as the Richardson, the Conjugate gradient, and the Chebyshev algorithms, that are compared here to the Eyre–Milton scheme — the most efficient specialized method currently available. Our numerical experiments, carried out for two-dimensional elliptic problems, reveal that the Conjugate gradient algorithm is the most efficient option, while the Eyre–Milton method performs comparably to the Chebyshev semi-iteration. The Richardson algorithm, equivalent to the still widely used original Moulinec–Suquet solver, exhibits the slowest convergence. Besides this, we hope that our study highlights the potential of the well-established techniques of numerical linear algebra to further increase the efficiency of FFT-based homogenization methods.

Fast low‐rank approximations of multidimensional integrals in ion‐atomic collisions modelling

SummaryAn efficient technique based on low‐rank separated approximations is proposed for computation of three‐dimensional integrals arising in the energy deposition model that describes ion‐atomic collisions. Direct tensor‐product quadrature requires grids of size 40003 that is unacceptable. Moreover, several of such integrals have to be computed simultaneously for different values of parameters. To reduce the complexity, we use the structure of the integrand and apply numerical linear algebra techniques for the construction of low‐rank approximation. The resulting algorithm is 103 faster than spectral quadratures in spherical coordinates used in the original DEPOSIT code. The approach can be generalized to other multidimensional problems in physics. Copyright © 2015 John Wiley &amp; Sons, Ltd.

Combining Algebraic and Numerical Techniques for Computing Matrix Determinant

Computing the sign of the determinant or the value of the determinant of an n × n matrix A is a classical well-know problem and it is a challenge for both numerical and algebraic methods. In this paper, we review, modify and combine various techniques of numerical linear algebra and rational algebraic computations (with no error) to achieve our main goal of decreasing the bit-precision for computing detA or its sign and enable us to obtain the solution with few arithmetic operations. In particular, we improved the precision bits of the p-adic lifting algorithm (H = 2h for a natural number h), which may exceed the computer precision β (see Section 5.2), to at most bits (see Section 6). The computational cost of the p-adic lifting can be performed in O(hn4). We reduced this cost to O(n3) by employing the faster p-adic lifting technique (see Section 5.3).

Pseudospectral methods for Boussinesq-type equations in an annular domain with applications to mid-sized lakes

We present a numerical solution procedure for a two-dimensional Boussinesq-type shallow water model on annular domains using pseudospectral discretization methods, including practical implementation details such as spectral filtering to prevent aliasing-driven instabilities and efficient numerical linear algebra techniques. The numerical model's potential for predicting and simulating wave motions in mid-sized lakes is illustrated with three test cases: (1) wave diffraction around an island and near-shore focusing; (2) the formation, propagation and destruction of wave trains and solitary-like waves in rotating basins; and (3) the influence of bottom bathymetry on wave formation and propagation from a Kelvin-seiche.

Computing polarizabilities without a Hamiltonian matrix

Using techniques of numerical linear algebra it is possible to calculate polarizabilities and hyperpolarizabilities without computing or storing a Hamiltonian matrix. This is done by exploiting the structure of the basis to evaluate matrix–vector products. The ideas are tested on a one dimensional Hubbard model.

Costate Computation by a Chebyshev Pseudospectral Method

AMONG the various pseudospectral (PS) methods for optimal control [1], only the Legendre PS method has been mathematically proven to guarantee the feasibility, consistency, and convergence of the approximations [2–5]. As exemplified by its experimental andflight applications in national programs [6–10], it is not surprising that the Legendre PS method has become the method of choice [11–19] in both industry and academia for solving optimal control problems. Efforts to improve the Legendre PS methods by using either other polynomials [20–22] or point distributions [23,24] have not yet resulted in any rigorous framework for convergence of these approximations [24,25]. Compared to Legendre PS methods, Chebyshev PS methods [21,22] for optimal control are somewhat more attractive for a number of reasons. When a function is approximated, it is well known that a Chebyshev expansion is very close to the best polynomial approximation in the infinity norm [26,27]. In addition, Chebyshev polynomials have an attractive computational advantage in terms of the computation of Chebyshev–Gauss–Lobatto (CGL) nodes. Unlike the Legendre–Gauss–Lobatto (LGL) nodes, CGL nodes can be evaluated in closed form [26]. Thus, a Chebyshev PS method offers the possibility of rapid computation because it does not require the use of advanced numerical linear algebra techniques that are necessary for the calculation of LGL nodes [21]. A similar numerical advantage applies to the computation of the derivative via a fast Chebyshev differentiation scheme that is similar to a fastFourier-transform (FFT) computation. In the same spirit, integration is also fast because of the connection between the Clenshaw–Curtis integration and the FFT [27]. Despite these attractive properties, Chebyshev PS methods have not advanced beyond the works of [21,22]. This is, in part, due to the absence of a covector mapping theorem that is crucial for the computation of the costates and other covectors. The computation of costates and other covectors is important in solving practical optimal control problems as it provides ameans for verification and validation of the computed solution [25]. Beyond verification and validation, information about covectors can also be used to facilitate the design of guidance and control algorithms [28]. In this Note, we fill the key gap of costate computation for Chebyshev PSmethods by furthering themethod of Fahroo andRoss [21]. We do this by combining some recent results from Clenshaw– Curtis integration [27], the unification principles proposed by Fahroo and Ross [1,29,30], and the new results of Gong et al. [23].

An Asymptotic Approach for Testing P0-Matrices

A direct approach to the P-matrix or P0-matrix problem is to evaluate all the principal minors of matrix A using standard numerical linear algebra techniques with O(2nn3) computational time complexity. The computational time complexity of the P-matrix problem has been reduced from O(2nn3) to O(2n) by applying recursively a criterion for P-matrices based on Schur complementation. But this algorithm can be not directly applied to test the P0-matrices because the Schur complementation can be not computed when some zero diagonal elements appear.This paper proposes an asymptotic approach for testing P0-matrices with O(2n) computational time complexity. Some numerical examples show that the proposed algorithm is effective for testing P0-matrices.

Algebraic Multigrid Solvers for Complex-Valued Matrices

In the mathematical modeling of real-life applications, systems of equations with complex coefficients often arise. While many techniques of numerical linear algebra, e.g., Krylov-subspace methods, extend directly to the case of complex-valued matrices, some of the most effective preconditioning techniques and linear solvers are limited to the real-valued case. Here, we consider the extension of the popular algebraic multigrid method to such complex-valued systems. The choices for this generalization are motivated by classical multigrid considerations, evaluated with the tools of local Fourier analysis, and verified on a selection of problems related to real-life applications.

Semidefinite Characterization and Computation of Zero-Dimensional Real Radical Ideals

For an ideal $I\subseteq\mathbb{R}[x]$ given by a set of generators, a new semidefinite characterization of its real radical $I(V_\mathbb{R}(I))$ is presented, provided it is zero-dimensional (even if $I$ is not). Moreover we propose an algorithm using numerical linear algebra and semidefinite optimization techniques, to compute all (finitely many) points of the real variety $V_\mathbb{R}(I)$ as well as a set of generators of the real radical ideal. The latter is obtained in the form of a border or Gr\"obner basis. The algorithm is based on moment relaxations and, in contrast to other existing methods, it exploits the real algebraic nature of the problem right from the beginning and avoids the computation of complex components.

Structured eigenvalue problems for rational gauss quadrature

The connection between Gauss quadrature rules and the algebraic eigenvalue problem for a Jacobi matrix was first exploited in the now classical paper by Golub and Welsch (Math. Comput. 23(106), 221–230, 1969). From then on many computational problems arising in the construction of (polynomial) Gauss quadrature formulas have been reduced to solving direct and inverse eigenvalue problems for symmetric tridiagonals. Over the last few years (rational) generalizations of the classical Gauss quadrature formulas have been studied, i.e., formulas integrating exactly in spaces of rational functions. This paper wants to illustrate that stable and efficient procedures based on structured numerical linear algebra techniques can also be devised for the solution of the eigenvalue problems arising in the field of rational Gauss quadrature.

Exponential data fitting applied to infiltration, hydrograph separation, and variogram fitting

Most lumped rainfall-runoff models separate the interflow and groundwater components from the measured runoff hydrograph in an attempt to model these as hydrologic reservoir units. Similarly, rainfall losses due to infiltration as well as other abstractions are separated from the measured rainfall hyetograph, which are then used as inputs to the various hydrologic reservoir units. This data pre-processing is necessary in order to use the linear unit hydrograph theory, as well as for maintaining a hydrologic budget between the surface and subsurface flow processes. Since infiltration determines the shape of the runoff hydrograph, it must be estimated as accurately as possible. When measured infiltration data is available, Horton’s exponential infiltration model is preferable due to its simplicity. However, estimating the parameters from Horton’s model constitutes a nonlinear least squares fitting problem. Hence, an iterative procedure that requires initialization is subject to convergence. In a similar context, the separation of direct runoff, interflow, and baseflow from the total hydrograph is typically done in an ad hoc manner. However, many practitioners use exponential models in a rather “layer peeling” fashion to perform this separation. In essence, this also constitutes an exponential data fitting problem. Likewise, certain variogram functions can be fitted using exponential data fitting techniques. In this paper we show that fitting a Hortonian model to experimental data, as well as performing hydrograph separation, and total hydrograph and variogram fitting can all be formulated as a system identification problem using Hankel-based realization algorithms. The main advantage is that the parameters can be estimated in a noniterative fashion, using robust numerical linear algebra techniques. As such, the system identification algorithms overcome the problem of convergence inherent in iterative techniques. In addition, the algorithms are robust to noise in the data since they optimally separate the signal and noise subspaces from the observed noisy data. The algorithms are tested with real data from field experiments performed in Surinam, as well as with real hydrograph data from a watershed in Louisiana. The system identification techniques presented herein can also be used with any other type of exponential data such as exponential decays from nuclear experiments, tracer studies, and compartmental analysis studies.

Computation of Periodic Solution Bifurcations in ODEs Using Bordered Systems

We consider numerical methods for the computation and continuation of the three generic secondary periodic solution bifurcations in autonomous ODEs, namely the fold, the period-doubling (or flip) bifurcation, and the torus (or Neimark--Sacker) bifurcation. In the fold and flip cases we append one scalar equation to the standard periodic BVP that defines the periodic solution; in the torus case four scalar equations are appended. Evaluation of these scalar equations and their derivatives requires the solution of linear BVPs, whose sparsity structure (after discretization) is identical to that of the linearization of the periodic BVP. Therefore the calculations can be done using existing numerical linear algebra techniques, such as those implemented in the software AUTO and COLSYS.

Applying numerical linear algebra techniques to analyzing algorithms in signal processing

Without theoretical evidence to support the use of an algorithm, one cannot be certain that if it solves one problem successfully it will necessarily solve all problems successfully. One can analyze an algorithm from many different settings, each having the possibility for different focal points for results and insights. In signal processing, an analysis into the stability of an algorithm can have as its focal point stability in the sense of Lyapunov or of numerical linear algebra (backward, forward), etc. Unfortunately, without a clear vocabulary to apply toward the focal point of interest, the interpretation is ambiguous. To avoid this problem, a clear distinction must be made that delineates each significant contribution to the process of solving the problem, from the definition of a problem to the derivation of the algorithm. It requires that perturbations assumed at each juncture be described in a clear language that offers no chance for confusion. We offer such a terminology and apply it to the fast transversal filter and its (existing) analyses, an algorithm that has long been known to diverge yet has nonetheless received much interest regarding the reasons behind its divergence as well as inquiries into how it might be modified to one that enjoys more dependable and robust behavior.

Numerical linear algebra in the integrity theory of the Global Positioning System

The Global Positioning System (GPS) is a satellite based navigation system. Since safety is the main concern for aircraft navigation, various means of monitoring the integrity (certainty of position) have been developed. This is an important area of research in the GPS community. In the following, it will be shown how some numerical linear algebra techniques can be applied to this interesting application. A typical model is presented. A uniform approach to derive the statistics for fault detection and isolation by orthogonal transformations is given. It is shown that the diagonal elements ω ii 2 of the orthogonal projection matrix onto the residual space are fundamental to the theory and understanding of integrity. ω ii can, for example, have a drastic effect on integrity when they are small. The sensitivity of related problems in this area are discussed.