RAM research project

RAM: Rational modelling: optimal conditioning and stable algorithms(01-01-2005 - 31-12-2008)

• [D1] The data

  How well the model should satisfy the conditions arising from the given data, depends on the context. When the data are less accurate, as is the case for measured data subject to noise, not even all data have to be used. In case of a continuous approximation problem, the correspondence of the model to the data is measured by means of a function norm. In many other circumstances, discrete data can be obtained very accurately, but sometimes at the expense of a very time-consuming computation process. A possible role of the model is then to give a sufficiently accurate but simplified representation of that numeric process. If the data points become gradually available, then they can be added to the modelling process in a recursive way. In all cases, the fact that the data points are fixed or that their location can be chosen during the computation, should be taken into account.

• [D2] The model

  Remains the choice of the model. The parameters of the model are computed from conditions which express the approximation of the data according to a chosen criterion. Polynomials and splines are rather popular for modelling and are easy to use. However, in the context of applications where a certain asymptotic behaviour has to be modelled, neither polynomials nor splines can be used since they tend to infinity. They can also not be used to capture drastic changes in function values. In these circumstances the use of rational models is recommended. When no specific behaviour is required, the rational model can reduce to a polynomial model. In a time-varying linear dynamical system, the coefficients of the model also vary with time. In this way a one-dimensional problem sometimes becomes higher dimensional.

Computational and numeric aspects

• [N1] The defining conditions

  Rational and polynomial models have in common the advantage that they can be computed by means of fast numeric techniques. These techniques make use of the structure, e.g. the displacement rank, of the linear system that expresses the modelling problem. Sometimes the computation of the model is also possible via a recursive technique. Then the complexity of the model and/or the size of the data set can be increased in a step by step manner.

  But for the advantages that rational models offer, there is a price to pay. The defining conditions are a priori non-linear and only give rise to a linear system after linearization. This gives rise to problems related to the existence and the convergence of the rational functions, etc.

• [N2] The chosen representation

  Several structured problems tend to be ill conditioned. This means that the computed model is very sensitive to small changes in the data. Small changes in the data can very often not be avoided, given the fact that the data can result from measurements or from a computationally intensive process which is subject to the inevitable accumulation of rounding errors.

  The conditioning of the problem is highly influenced by the way in which rational functions are represented. A rational function is in essence a pair of irreducible polynomials, each of which can be represented as a linear combination of either orthogonal polynomials, Newton-polynomials, the Bernstein basis and so on. Alternatives are orthogonal rational functions, continued fractions or recurrence relations, state-space formulas, a barycentric representation, etc. Because of the link with the conditioning, the choice of the representation is crucial. Yet determining an optimal representation is still an open problem, certainly in a multidimensional (several variables, several functions) context.

Problem formulation

Concisely formulated, the problem the teams plan to tackle in this project is this: how can a rational modelling problem be formulated such that it is optimally conditioned, especially in a multidimensional context? Given the very large range of applications, the answer to this problem is important.

The following example is prototypical for the problem formulated. The identification problem in the frequency domain is a rational approximation problem on the imaginary axis or the complex unit circle. The spectral transfer function is measured in a number of discrete points together with stochastic information of the measurement noise. The representation of the model as a couple of polynomials represented by monomials leads to unacceptable conditioning problems. The conditioning problem was largely solved by a quasi-optimal choice of the representation of the model in terms of orthogonal vector valued polynomials. The project promotors want to continue along this road, as explained in the project objectives. The team is especially interested in problems involving the modelling of multivariate data and/or matrix-valued data.

Participating teams

For this project, the participation of the following two teams is essential. They both have extensive experience in the field of modelling on one hand and in the use of rational functions on the other. In particular, the NALAG team of the K.U.Leuven (A. Bultheel, M. Van Barel) has specific experience in orthogonal rational functions, rational interpolation problems, multi-input multi-output systems, and the solution of structured linear systems. On the other hand, the CANT/ECT team from the University of Antwerp(A. Cuyt, B. Verdonk) is specialized in the solution of multidimensional rational modelling problems, in multivariate orthogonal polynomials an in scattered rational interpolation. Both teams were recently confronted with the problem formulated above and wish to join efforts to tackle this problem.

Objectives

• [O1] Objective 1: Analysis of the condition number

  First, a thorough analysis needs to be made of the causes of the ill conditioning. In this respect both the available data and the basis functions used, play an important role. A study, identifying an optimal orthogonal basis on one hand and the optimal placement of the data points on the other, is required.

  Optimal location of the data

  
  In the study of the optimal placement of the data points, the situation in which certain (directional) derivatives are given or can be obtained, must be considered. In a number of existing algorithms for (multidimensional) rational modelling, derivative information can not yet be handled.

  Furthermore, the effect of the starting values on the recursive computation of the rational model must be thorougly studied. Experiments have indicated that an algorithm which tries to compute subsequent data in (according to a certain criterion) optimally placed points, can be strongly influenced by the initial data chosen by the user. Especially in applications where the collection of data is a time-consuming process, it is important to be able to compute the desired model with a minimal number of data. The situation where a choice of initial data puts the modelling process on a wrong track needs to be avoided. In other cases, where the sampled data are fixed a priori, choosing data points freely is not an option.

  Optimal basis functions

  
  When the rational model is written down in terms of the classical monomial basis, in one or more dimensions, the problem is highly structured and in general ill conditioned. Switching to an orthogonal basis, or in several variables to the product of orthogonal basis functions or to spherical orthogonal polynomials, can lead to a substantial improvement in the conditioning. The problem is therefore "reduced" to the identification of an optimal basis, or in other words, the identification of the weight function which gives rise to the lowest condition number.

  In one variable, the basis of Bernstein polynomials is known to minimize the condition number for the solution of interpolation problems. In several variables the situation is not so straightforward since increasing the complexity of the model when a multidimensional Bernstein basis is used, increases the total degree, not the partial degree.

• [O2] Objective 2: Study of the close link with quadrature and cubature

  The whole problem is strongly related to the development of quadrature or cubature formulas of maximal order for the exact integration of rational functions. For numerical cubature or quadrature of integrands that are singular near the border of the region of integration, again it may be more appropriate to use rational approximants instead of polynomial approximants. In the polynomial case, Gauss quadrature formulas are known to have maximal degree of exactness. Similar techniques can be used to obtain quadrature formulas that have a maximal degree of exactness is certain spaces of rational functions. In the one-dimensional case this problem has been studied, but in the two or n-dimensional case, there are no rational analogs of Gaussian quadrature formulas known. In this case, we may fall back on the interpolatory type of cubature formulas, and then the problem is completely similar to the modelling problem that has just been sketched.

  The study of the link between rational modelling and integration formulas even has an impact on certain crucial questions in the domain of quantum computing! Going further into the matter would lead us too far.

• [O3] Objective 3: Stable algorithms

  
  

  For the basis

  
  Once the optimal basis has been identified, it still needs to be computed. For instance, for orthogonal rational functions this computation is usually based on a generalization of the classical 3-term recurrence relation for orthogonal polynomials. However, the conditioning problem of the matrix that we have mentioned above, is now hidden in the computation of the recurrence coefficients. There exist several algorithms to compute them theoretically, but these are numerically useless because of the ill conditioning. In the one-dimensional case, the teams have already developed a technique for the computation of rational basis functions orthogonal with respect to a discrete weight function. This technique is based on the solution of an inverse eigenvalue problem of a matrix which is the sum of a diagonal matrix and a semi-separable matrix.

  For the rational model

  
  Once the optimal basis has been indentified and computed, the rational model, which is the desired output of the modelling process, has to be computed as the solution of a structured system of equations that express the modelling conditions. When switching from the original formulation of the problem to the optimal basis, the structure may seemingly be destroyed, yet it is still latently available. The difficult part is to detect, and catch the structure in an efficient way and to process it in a numerically stable fashion. This structure is important since it gives rise to a reduction in the computational complexity. In turn, the use of an optimal basis makes sure that the problem is optimally conditioned.

  In addition, one can consider the application of appropriate preconditioners. The basis transformation which we have discussed so far corresponds to right preconditioning. But also left preconditioning can be considered to improve the conditioning of the system. While the preconditioning of linear systems has been thoroughly studied, the existing preconditioners are not applicable to the structure of the linear systems involved. Many of the currently available preconditioners are developed either for specific applications where certain assumptions can be made about the structure of the coefficient matrix (such as symmetry), or they are specifically aimed at certain matrix structures (such as Toeplitz with square sub-blocks). Neither of these are applicable to the specific structure of the systems with which the teams are confronted.