Van Dooren Research Articles

Issues of terminology, gradient dynamics and the ease of sympatric speciation in Adaptive Dynamics

It has been plainly stated that Adaptive Dynamics is ‘not a scientific theory, but a mathematical framework’ (Metz, 2005); it is a ‘method but not a model’ (Kisdi & Gyllenberg, 2005). From a completely pragmatic point of view, there is, then, only one question that evolutionary biologists need to ask: Is Adaptive Dynamics a useful way of approaching problems? The 14 Commentaries present a variety of different views, indicating that there is no consensus in the answer to this question. Indeed, we view the previous pages, consisting of the Target Review and the 14 Commentaries, as a part of a much larger public discussion on Adaptive Dynamics, and we thank the other authors for their contributions to the discussion. We feel it necessary to acknowledge that the wording of our comment in Question 6 of our Target Review (Waxman & Gavrilets, 2005) concerning ‘hidden limitations and unconscious or implicit assumptions’ was not optimal. It was not our intention to imply that limitations and assumptions of Adaptive Dynamics are hidden or unconsciously made by practitioners of Adaptive Dynamics. Rather, the statement was intended for biologists who use Adaptive Dynamics, or who wish to understand the biological implications of its mathematical results, and may not be aware of the underlying limitations and assumptions. These limitations and assumptions are not, in our opinion, sufficiently emphasised in the literature and, equally importantly, the consequences of their violation are not made clear. The Commentary of Geritz & Gyllenberg (2005) makes it admirably clear what these limitations are but not the consequences of their violation. The Commentaries of Abrams (2005), Barton & Polechova (2005), Kokko (2005), Spencer & Feldman (2005), Butlin & Tregenza (2005), and Van Dooren (2005) provide additional discussion of the issue. We also feel it necessary to clarify our wording about the treatment of genetic drift in Adaptive Dynamics. In the second paragraph of Question 4.3, we said an implicit assumption of Adaptive Dynamics is that all beneficial mutations will always initially increase in frequency. However in the first paragraph, we noted that in the work of Dieckmann & Law (1996), fixation of new alleles occurs at a rate proportional to the probability of fixation of a mutant. Since the probability of fixation of all mutants is not unity, the two paragraphs are clearly inconsistent. Our reading of the Adaptive Dynamics literature, however, leads us to believe that it is common to make the assumption that beneficial mutations always initially increase in frequency (see Meszena, 2005). We are also of the opinion that the treatment of genetic drift, when included in Adaptive Dynamics, has not been particularly realistic. As an example, the probability of invasion of a deleterious mutation is invariably taken as zero in Adaptive Dynamics (Geritz & Gyllenberg, 2005; Meszena, 2005). However, if a mutation is deleterious, but only very slightly so, and the population is finite, then such near neutral mutants can invade and become fixed in the population (Kimura, 1957). At the end of our Review, we stated our earnest and legitimate wish to see more details of numerical simulations, but this was somehow interpreted, in some Commentaries, as an attempt to denigrate Adaptive Dynamics, because of the lack of analytical results, or to attach specific significance to simulations in Adaptive Dynamics. This was not our intention. Although we appreciate all of the Commentaries made, we will not answer them directly. Rather we discuss three controversial points raised in more than one Commentary. A separate response is, however, also necessary for a number of peripheral points raised by Dieckmann & Doebeli (2005). Our response to these points is contained in Appendix B.

THE STABLE EMBEDDING PROBLEM

We compute a stable polynomial matrix embedding a stabilizable one. The algorithm resembles the one described by Beelen and Van Dooren for the unimodular embedding problem. We also desribe the numerical problems associated with this kind of algorithms, and point out why the stable embedding algorithm has better prospects than its unimodular counterpart.

Solving Polynomials with Small Leading Coefficients

We explore the computation of roots of polynomials via eigenvalue problems. In particular, we look at the case when the leading coefficient is relatively very small. We argue that the companion matrix algorithm (used, for instance, by the Matlab {\tt roots} function) is inaccurate in this case. The accuracy problem is addressed by using matrix pencils instead. This improvement can be predicted from the backward error bound of Edelman and Murakami (for companion matrices) versus the bound of Van Dooren and Dewilde (for pencils). We then show how to extend the accurate algorithm to Bézier polynomials and present computational experiments.

Improper Applications of Proportional Reasoning

One of the major goals of elementary and middle-grades mathematics education is for students to obtain a deep understanding of the proportional model in a variety of forms and applications. However, the reinforcement of proportionality at numerous occasions in school mathematics, along with the teaching of some standardized methods for solving proportionality problems, appear to lead to a resistant tendency in some students and adults to see and apply proportions everywhere. This same application occurs in situations where another method of solution is appropriate. Along with mastering the proportional scheme, its misuse seems to appear, as well. This overgeneralization of proportion has many faces: It has been found at different age levels and in a variety of mathematical domains, such as elementary arithmetic (Cramer, Post, and Currier 1993), algebra (Matz 1982), geometry (De Bock, Verschaffel, and Janssens 1998, 2002) and probability (Van Dooren et al. 2002).

Computation of canonical matrices for chains and cycles of linear mappings

Van Dooren [Linear Algebra Appl. 27 (1979) 103] constructed an algorithm for the computation of all irregular summands in Kronecker’s canonical form of a matrix pencil. The algorithm is numerically stable since it uses only unitary transformations. We construct a unitary algorithm for computation of the canonical form of the matrices of a chain of linear mappings V 1—V 2—⋯—V t andextend Van Dooren’s algorithm to the matrices of a cycle of linear mappings whereall V i are complex vector spaces and each line denotes → or ←.

The fixed poles of the disturbance decoupling problem and almost stability subspace $\mathcal V^\star_{b,g}({\rm ker}(C))$

This paper is concerned with the fixed poles of the disturbance decoupling problem for linear time-invariant systems by state feedback and measurement feedback. Algebraic characterizations for these fixed poles are given using the reducing subspace technique in numerical linear algebra. These algebraic characterizations lead that the fixed poles can be directly computed by numerically reliable algorithms. As an additional application of the reducing subspace technique, we also develop a numerically stable method for computing the almost stability subspace $\mathcal V^\star_{b,g}({\rm ker}(C))$ , which is an extension of the method given by Van Dooren for computing the invariant subspaces in geometric control theory.

Problems and Techniques

The Rayleigh quotient iteration (RQI) that we know and love is used to calculate an eigenvalue-eigenvector pair of a symmetric matrix. There are numerous situations when multiple eigenvalues and the subspace spanned by the corresponding eigenvectors are required. Our first article, by Absil, Mahony, Sepulchre, and Van Dooren, describes a generalization of RQI to address this problem. While other algorithms for this problem exist, they either converge at a slow linear \textit{O}(\textit{n}) rate in the dimension \textit{n} of the matrix or have a high \textit{O}(\textit{n}$^{\mbox{\font size{6}{6}\selectfont 3}}$) cost per iteration. The present generalized RQI procedure has cubic convergence and an \textit{O}(\textit{np}$^{\mbox{\font size{6}{6}\selectfont 2}}$) unit cost for invariant subspaces of dimension \textit{p}. The key involves the solution of a Sylvester system working with a condensed form of the original matrix and the theory relies onthe structure of Grassmann manifolds. The current situation in Afghanistan would seem to make the article by Hassan Sedaghat on the dynamics of a discrete-time model of combat particularly relevant. Sedaghat examines the behavior of a "ground version" of a combat model proposed by Epstein by neglecting air support. He emphasizes that realistic combat models used for "war game" simulations would include this as well as deterministic and stochastic effects due to, e.g., weather, terrain, psychology, and logistics. Nevertheless, the simplified deterministic ground model, consisting of three nonlinear difference equations, has a rich structure that may form the core of a more complex model. The key ingredient of the ground model is a withdrawal strategy that introduces a jump discontinuity into the model. Sedaghat's asymptotic and transient analyses of the ground model reveal that steady, oscillatory, and chaotic behaviors are all possible. The latter effect sounds realistic to me! Please enjoy these very interesting articles.

Parameter depending state space descriptions of index-2-matrix polynomials

To a quadratic matrix polynomial P with coefficients in R n×n , which originated from an electrical network and depending on a parameter vector q∈ R ϱ , a matrix A( q) and a parameter set Ω are assigned such that for all q∈Ω the eigenvalues of A( q) coincide with the zeros of det P(.; q). To find A( q) and the corresponding parameter set Ω two algorithms are proposed where the first one is similar to the Algorithm 3.6 of Van Dooren [Linear Algebra Appl. 27 (1979) 103]. The reason to compute A( q) parameter depending is given by the desire to apply the matrix perturbation theory of Stewart and Sun [Matrix Perturbation Theory (Academic Press, 1990)] to study the influence of q on the zeros of the determinant of P. The assumption that P is derived from the Laplace transform of a DAE system describing an electrical network implies that its coefficients admit representations containing sums of the form ∑ v k q k w k T, where v k and w k are the unit vectors or the nonvanishing difference of two such vectors. This circumstance is decisive for the efficiency of our two algorithms.

Introduction to the Special Section on Sparse and Structured Matrices and Their Applications

The field of sparse matrices is a broad and important area of the computational sciences that includes structured matrices and those with seemingly little or no structure. The relevance of the field is highlighted by the wide range of application areas that require the exploitation of matrix sparsity and structure in order to achieve a solution given real-world constraints on computing resources and/or time. Applications in which sparse matrices appear include structural analysis, computational fluid dynamics, economic modeling, financial analysis, numerical optimization, statistical modeling, power network analysis, electromagnetics, meteorology, medical imaging, data mining, and many more. A number of significant advancements in sparse matrix computations have been made in recent years. These advances have led to new challenges as multidisciplinary problems are now ambitiously posed, and they symbolize the growth in the area and demonstrate the dependence of the future of many fields on sparse matrix methods. Moreover, they give rise to a host of new problems that have yet to be addressed by the sparse matrix community. The Second SIAM Conference on Sparse Matrices, which was held in Coeur d'Alene, Idaho, October 9--11, 1996, was organized specifically to address some of the challenging issues in sparse matrix computations. This special issue of SIAM Journal on Matrix Analysis and Applications consists of some of the outstanding papers that were presented at the Sparse Matrix conference. Because of an error in scheduling, the paper titled ``Using Generalized Cayley Transformations within an Inexact Rational Krylov Sequence Method" by Richard Lehoucq, which should have appeared in this section, was published in SIAM Journal on Matrix Analysis and Applications, volume 20, number 1. We are grateful to Paul van Dooren, the Editor-in-Chief of SIAM Journal on Matrix Analysis and Applications, and the SIAM office for agreeing to publish a special section on sparse matrices. We would like to thank Roland Freund, Anne Greenbaum, Joseph Liu, and Zdenek Strakos for serving on the Guest Editorial Board; their effort and cooperation in handling the papers were much appreciated. Finally, we would like to express thanks to all the authors who have submitted papers to the special section; it would not have become a reality without their support.

Successes and failures of current treatment of heart failure

Successes and failures of current treatment of heart failure

Globally convergent algorithms for robust pole assignment by state feedback

It is observed that an algorithm proposed in 1985 by Kautsky, Nichols, and Van Dooren (KNV) amounts to maximizing, at each iteration, the determinant of the candidate closed-loop eigenvector matrix X with respect to one of its columns (with unit-length constraint), subject to the constraint that it remains an achievable closed-loop eigenvector matrix. This interpretation is used to prove convergence of the KNV algorithm. It is then shown that a more efficient algorithm is obtained if det (X) is concurrently maximized with respect to two columns of X, and such a scheme is easily extended to the case where the eigenvalues to be assigned include complex conjugate pairs. Variations exploiting the availability of multiple processors are suggested. Convergence properties of the proposed algorithms are established. Their superiority is demonstrated numerically.

Polynomial Roots from Companion Matrix Eigenvalues

In classical linear algebra, the eigenvalues of a matrix are sometimes defined as the roots of the characteristic polynomial. An algorithm to compute the roots of a polynomial by computing the eigenvalues of the corresponding companion matrix turns the tables on the usual definition. We derive a first-order error analysis of this algorithm that sheds light on both the underlying geometry of the problem as well as the numerical error of the algorithm. Our error analysis expands on work by Van Dooren and Dewilde in that it states that the algorithm is backward normwise stable in a sense that must be defined carefully. Regarding the stronger concept of a small componentwise backward error, our analysis predicts a small such error on a test suite of eight random polynomials suggested by Toh and Trefethen. However, we construct examples for which a small componentwise relative backward error is neither predicted nor obtained in practice. We extend our results to polynomial matrices, where the result is essentially the same, but the geometry becomes more complicated.

Polynomial roots from companion matrix eigenvalues

In classical linear algebra, the eigenvalues of a matrix are sometimes defined as the roots of the characteristic polynomial. An algorithm to compute the roots of a polynomial by computing the eigenvalues of the corresponding companion matrix turns the tables on the usual definition. We derive a first-order error analysis of this algorithm that sheds light on both the underlying geometry of the problem as well as the numerical error of the algorithm. Our error analysis expands on work by Van Dooren and Dewilde in that it states that the algorithm is backward normwise stable in a sense that must be defined carefully. Regarding the stronger concept of a small componentwise backward error, our analysis predicts a small such error on a test suite of eight random polynomials suggested by Toh and Trefethen. However, we construct examples for which a small componentwise relative backward error is neither predicted nor obtained in practice. We extend our results to polynomial matrices, where the result is essentially the same, but the geometry becomes more complicated.

A Jacobi-Type Systolic Algorithm for Riccati and Lyapunov Equations

A systolic algorithm is devised for solving a class of Riccati and Lyapunov equations, which makes use of a factored version of the matrix sign recursions due to Charlier and Van Dooren ( MCSS 2 (1989), 109-136). The original algorithm is worked into an alternative Jacobi-type algorithm, which is readily implemented on a systolic array. Compared to the array of Charlier and Van Dooren, the present algorithm/array is conceptually simpler and, furthermore, roughly three times more efficient.

Computation of system zeros with balancing

The numerical computation of the transmission zeros of a linear, time-invariant system requires the solution of a generalized eigenvalue problem ( GEP) Mx= λNx with N=block diag ( I,0). The Emami-Naeini Van Dooren (ENVD) algorithm for the computation of system zeros avoids the computation of spurious zeros by isolating the finite zeros of a system prior to numerical solution the corresponding GEP. Unfortunately, the ENVD algorithm can be sensitive to the dynamic range of system coefficients. We present a numerical procedure for reduction of the dynamic range of coefficients in a zero-computation problem that preserves the underlying structure N=block diag ( I,0) of the zero-computation GEP.

The Laurent expansion of pencils that are singular at the origin

This paper gives two new proofs of a theorem of Langenhop on the Laurent expansion of a matrix pencil that is singular at the origin. The second, based on a decomposition of Van Dooren, leads to a computational algorithm.

Disturbance decoupling with stability and impulse-free response in singular systems

The problem of disturbance decoupling with stability and impulse-free response in square singular systems E x ̇ = Ax + Bu + Gq with output y = Dx , by applying a state feedback u = Fx , is considered. By employing Van Dooren's recursive algorithm for the computation of the subspace supremal ( A, B)-invariant subspace in Kernel D and using Kalman's decomposition theorem, the paper's approach presents a complete solution to the underlying problem. Regarding state space systems, it is shown that the solution to the underlying problem is obtained almost directly from Van Dooren's recursive algorithm.

Rational matrix GCDs and the design of squaring-down compensators-a state-space theory

A state-space construction for rational matrix greatest common divisors (GCDs) of rational transfer matrices is given. It is shown how the GCD results can be used to solve the problem of designing stable minimum-phase squaring-down compensators for multi-variable plants. One application is a direct state-space construction for such compensators and a state-space solution to 'fat-plant' H-infinity control problems. The results make use of the concepts of strongly observable systems and maximally observable systems and build upon the concepts introduced in the state-space GCD extraction results for polynomial matrices of L.M. Silverman and P. Van Dooren (1979).< <ETX xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">&gt;</ETX>

Feedback System Design and Structural Properties: The Role and Use of Rational Matrix GCD's

A state-space construction for rational matrix greatest common divisors (GCO's) is given. It is shown how the GCO results can be used to solve the problem of designing stable minimum-phase squaring-down compensators for multivariable plants, leading to a direct state-space construction for such compensators and a state-space solution to "fat-plant- H-infinity control problems. The theory also provides a state-space construction for minimal bases for rational vector spac(ts in the sense of Forney. The results make use of the concepts of strongly observable and maximally unobservable systems initiated in the respective works of Silverman and of Wonham and build upon the concepts introduced in the state-space GCO extraction results for polynomial matrices of Silverman and Van Dooren.

Cholesky Factor Updating Techniques for Rank 2 Matrix Modifications

Gill, Golub, Murray, and Saunders have described five methods by which the Cholesky factors of a positive definite matrix may be updated when the matrix is subjected to a symmetric rank 1 modification. For a negative rank 1 update, a modification of one of their methods was given by Lawson and Hanson and analyzed by Bojanczyk, Brent, van Dooren, and de Hoog. In many minimization algorithms, symmetric rank 2 modifications are found.This paper shows how each of the rank 1 methods gives rise to a single-application rank 2 method. For some of the methods, this involves a new Householder transformation technique designed to eliminate elements of two vectors at once using a rank 1 correction of the identity matrix.The authors’ experiments on scalar, vector, and shared-memory multiple-instructions multiple-data machines show that it is more economical to perform rank 2 updates rather than two rank 1 updates. In their comparison, the authors do not consider pipelining two applications of the rank 1 algorithms, wh...