A general approach to multivariable recursive interpolation

We consider here the problem of constructing a general recursive algorithm to interpolate a given set of data with a rational function. While many algorithms of this kind already exist, they are either providing non-minimal degree solutions (like the Schur algorithm) or exhibit jumps in the degree of the interpolants (or of the partial realization, as the problem is called when the interpolation is at infinity, see Rissanen (SIAM J Control 9(3):420–430, 1971) and Gragg and Lindquist (in: Linear systems and control (special issue), linear algebra and its applications, vol 50. pp 277–319, 1983)). By imbedding the solution into a larger set of interpolants, we show that the increase in the degree of this representation is proportional to the increase in the length of the data. We provide an algorithm to interpolate multivariable tangential sets of data with arbitrary nodes, generalizing in a fundamental manner the results of Kuijper (Syst Control Lett 31:225–233, 1997). We use this new approach to discuss a special scalar case in detail. When the interpolation data are obtained from the Taylor-series expansion of a given function, then the Euclidean-type algorithm plays an important role.


Introduction
The concepts of controllability indexes and tools from module theory like polynomial basis have been widely used in system theory and they have been extensively applied to different representations of transfer functions (rational functions, quotient B Gy. Michaletzky michaletzky@caesar.elte.hu A. Gombani gombani@ieiit.cnr.it 1 IEIIT-CNR, Via Giovanni Gradenigo 6a, 35131 Padua, Italy 2 Eötvös Loránd University, Pázmány Péter sétány 1/C, Budapest 1117, Hungary of polynomials, etc., see Kailath [12]), and in Antoulas et al. [3] they have also been connected to an interpolation problem for given data by means of a rational function. In [18], the present authors developed an approach for the characterizations-under weaker conditions-of proper interpolants by means of rational functions using controllability indexes. We show here how this framework accommodates quite naturally a recursive algorithm to compute the interpolants.
The simplest (and well-known) instance of the problem is, given a sequence of complex values v 0 , . . . , v K , to find, for k such that 0 ≤ k ≤ K , all rational functions which interpolate those data at 0, that is, functions of the form f (s) = k i=0 s i v i +s k+1 g(s). This problem has a long history which (in its conceptual framework) goes back at least to Euler (see Wyman and Wyman [21] for a nice translation of Euler's paper and Meijering [17] for an exhaustive history of the problem) and, in its recursive form, has been widely studied especially in the case where there is one interpolation node at infinity (instead of 0: it is the well-known partial realization problem). In this context particularly relevant are the works of Ho and Kalman [11], of Rissanen [19] and of Gragg and Lindquist [10]. The common feature of all these algorithms is the intrinsic presence of sudden jumps in the degree of the interpolants as we add further data. A simple example of this behavior is a Fibonacci-like sequence which suddenly gets out of the recursion, like 1, 1, 2, 3, 5, 8, 13, 100: it can be seen that the degree of the recursion is 2 until it "hits" 100 and jumps to 4. This makes these algorithms quite unsuitable for parameterization and thus for applications in identification. Their reliance on the Euclidean algorithm also makes extensions to the multivariable case or to different interpolation nodes intrinsically very difficult. The generalization of the approach initiated by Loewner [15] in 1934 for scalar interpolation problem was first systematically applied to the multivariable situation in Anderson and Antoulas [1]. The so-called Loewner matrices (divided difference matrices, null-pole coupling matrices) play a key tool in Ball et al. [4], Mayo and Antoulas [16] and Lefteriu and Antoulas [14], mentioning only a few. While the multivariable multiple node problem has been studied by several authors (see, for example, Antoulas et al. [3], Ball et al. [4] and Gombani and Michaletzky [9] for Kimura-Georgiou like fixed degree interpolants), a recursive version of these characterizations are-in general-not discussed. A recursive algorithm, which only works for two-sided interpolation, is presented for the matrix case in Lefteriu and Antoulas [14]. We provide here a recursive scheme under weaker conditions than those required in Antoulas et al. [3] and Fuhrmann [8].
We show how considering the larger (quite canonical) family of solutions (actually an immersion into a larger space) considered in [18], provides an algorithm which exhibits an astonishing regularity and makes the extensions to different nodes and the multivariable case quite straightforward. This provides a substantial improvement on Kuijper [13], where a single pole at infinity is considered and a coprimeness assumption on the interpolating polynomials had to be made. This feature thus sparks renewed attention for the Behaviors approach to Systems Theory devised by Polderman and Willems [20], where lack of coprimeness of an AR representation was a fundamental and interesting feature of the models, but where the practical construction of recursive algorithms for such models relied on this very coprimeness assumption. Closer to our approach are the nice results in Boros et al. [5], where also polynomial families of solutions with the same regularity are obtained. Nevertheless, that approach does not guarantee uniqueness of the representation and its computation, in the case of confluent interpolation, is quite involved.
The fundamental idea is quite simple: most recursive techniques consider, for the scalar case, the smallest number for which some columns of data (in general arranged as a Hankel matrix) are linearly dependent. This means to consider enough elements so that we obtain a matrix with a one-dimensional kernel at each step. Adding a new row of data usually makes the matrix full rank and we thus have to add another column to this matrix to obtain a new kernel. The problems occur when the rank is not increased by the extra data for a while (the model already accommodates the new data) and then suddenly there is a jump in the rank. (We have to add several columns at once to have a nontrivial kernel.) We consider instead a slightly larger matrix, with a twodimensional kernel, for which we introduce an algorithm which always decreases by one the dimension of the kernel as we add data and thus has no jumps.
The paper is structured as follows: in Sect. 2 we introduce some general results for the tangential problem (see Antoulas et al. [3] or Ball et al. [4]) and discuss the existence of a polynomial solution. Since, as we said, the idea is quite simple, but the details for the multivariable case are rather intricate, in Sect. 3 we discuss the scalar rational interpolation problem and we lay the ground for characterizing all interpolants through the construction of a fat matrix (essentially a basis obtained from a nice selection with the addition of two extra columns) and connect it to a 2 × 2 polynomial matrix (see Antoulas et al. [3]) which we call fundamental solution: its main feature is that its rank is always 1 at the interpolation nodes and 2 everywhere else. (This solution is not unique, but all solutions are related by units.) As similar results are needed for the multivariable case, most proofs are deferred to that part. In Sect. 4 we exhibit a scalar interpolation algorithm in some detail, as we feel that it provides the essential ingredients for the multivariable case while providing a grasp of what the main idea is. In Sect. 5, we extend the construction of Sect. 3 to the case of tangential interpolating conditions. In Sect. 6 we show how to extend the recursive algorithm to the general tangential rational interpolation problem with arbitrary nodes. In Sect. 7 we apply our analysis to a scalar case having only one interpolation point at zero, and show how this special feature implies a very fine and detailed structure. In the "Appendix" we provide an example with the Fibonacci sequence.

Remark 1
The present approach seeks minimal degree interpolants and does not consider constrains on their norm and their analyticity, like in Nevanlinna-Pick or Caratheodory problems (see, for example, Dym [6] and Ball et al. [4]). It turns out, though that a similar, albeit non recursive, approach, can be used to tackle this kind of problems and that it yields a simpler solution to them. This will be discussed in a forthcoming paper.

A general tangential interpolation problem
Suppose that we are given a triplet (A, U , V ), where A, U , V are of size (K + 1) × (K +1), (K +1)×m and (K +1)× p, respectively, determining a so-called tangential interpolation problem. The ultimate goal is to characterize all matrix-valued rational functions Q(s) of size m × p possibly in the form Q(s) = β(s)α(s) −1 of a given McMillan degree (if any exists) such that is analytic at σ (A), the spectrum of A. Note that since the product Q(s)α(s) = β(s) is a polynomial a slightly different formulation of this problem is the following: Problem 1 For a given triplet (A, U , V ) find all the pairs of matrix polynomials (α, β) such that is a polynomial.
While the row dimensions of α and β are m and p, and their common column dimension could be an arbitrary natural number, it turns out that, to characterize all such pairs, we will need precisely α to be square and invertible.

Remark 2
Let us point out that, if β and α are right coprime and α is an invertible polynomial matrix, then Problem 1 is equivalent to the following: set Q(s) = β(s)α −1 (s). There exist a polynomial matrix p(s) such that Proof In fact, if N (s) in Eq. (2) above is a matrix polynomial then the right coprime property of β(s), α(s) implies that there exist matrix-valued polynomials φ(s), ψ(s) such that which is analytic on σ (A).
The converse statement is obvious.
As it is well known from the literature on interpolation (see, e.g., [3,4]) a coprime factorization of (s I − A) −1 [U , −V ] plays a crucial role. Since for the derivation of the present results-except in Proposition 1, for the time being the distinction between the tangential conditions U and the data V is irrelevant, we will set, from now on, W := [U , −V ] of dimension (K + 1) ×r and therefore we will consider the following polynomial coprime factorization with Φ, Γ right coprime. We will also refer to the interpolation problem as one defined by (A, W ) instead of (A, U , V ). Thus, clearly (3) can also be written as provided that Φ and Γ are right coprime and Γ (s) is invertible.
The following simple lemma-see Ex. 6.3-19 in Kailath [12]-provides a straightforward way to check whether, for a given factorization, this is the case: We will assume, from now on, that Φ and Γ are right coprime. There are a few properties of (3) which will be needed.

Remark 3
Since Φ, Γ are right coprime, (s I − A) and Γ (s) have the same-nonunity-invariant factors; especially, det Γ (s) = det(s I − A), and thus, the zeros of det Γ are contained in σ (A). Coprimeness also implies that if Γ (λ)ξ = 0 for some vector ξ ∈ C r , and λ ∈ C, then Φ(λ)ξ = 0; therefore, Φ(λ)ξ is in the kernel of (λI − A). With some abuse of notation, we say that a polynomial matrix γ generates solutions to Problem 1 defined by (A, W ) if there exists a polynomial matrix φ such that In other words, using γ (s) as a polynomial input to the transfer function (s I −A) −1 W the output-φ(s)-is polynomial, as well.
Although the next lemma is essentially the basic starting block of the derivation of all solutions to an interpolation problem in the form of linear fractional transformation (cf. [3,4,14]) for the readers convenience we include here a short proof of it.

Lemma 2
Assume that the pair (A, W ) is controllable and γ generates a solution to Problem 1 defined by (A, W ). Then, there exists a matrix polynomial π(s) such that γ (s) = Γ (s)π(s) (that is, Γ (s) generates the module of polynomials solving Problem 1).
The following proposition shows that under some mild conditions there exists a matrix polynomial solution to Problem 1.
Proof Consider the coprime factorization of (s (3) and define a partitioned form of Γ : corresponding to the partition W = [U , −V ]. First we show the controllability of the pair (A, U ) is equivalent to the left coprimeness of α U and α V . Assume that there exists a scalar s 0 ∈ C such that the matrix [α U (s 0 ), α V (s 0 )] is not of full row-rank. This implies that there exists a vector η such that In particular, ξ * U = 0. Thus, the pair (A, U ) is not a controllable pair.
Conversely, if the pair (A, U ) in not controllable, then there exists a left eigenvector ξ of A for which ξ * U = 0. Denoting the corresponding eigenvalue by s 0 we get that ξ * (s 0 I − A) = 0, implying that Since the vector ξ * [U , −V ] is nonzero due to the controllability of the pair (A, [U , −V ]) we have that ξ * V = 0. This implies that the polynomial matrices [α U and α V are not left coprime.
Next we show that the left coprimeness of α U and α V is equivalent to the existence of a matrix polynomial β 0 such that (I , β 0 ) provides a solution to Problem 1.
Since according to Lemma 2 any solution can be written in the form Γ π thus if (I , β 0 ) provides a solution then for some matrix polynomial π we have that The second entry shows that α U and α V are left coprime. Conversely, if α U , α V are left coprime, then there exist a polynomial pair π U , π V such that Introducing the notation We obtain that β 0 (s) I = Γ (s)π(s) .
Thus β 0 I provides a solution to Problem 1, concluding the proof of the proposition.

Remark 4
Let us point out the following immediate corollary of the previous Lemma. If the pair (A, U ) is controllable, then the rank drop of the polynomial matrix Γ (s) cannot be greater than the number of the columns of U at any s ∈ C.
The following proposition slightly generalizes the argument applied in the proof of Theorem 3.5 in [3].
Then all solutions β(s) α(s) can be written in the form for some matrix polynomial π g (s).
In particular, in this case where β 0 (s) is a particular polynomial solution (might be called Hermite interpolant), and (s I − A) −1 U Γ 1 (s) = Ψ 1 (s) is already a polynomial matrix.
Although all the ingredients of the recursive algorithm for constructing solutions of the tangential interpolation problem are quite simple, the details for the multivariate case are rather cumbersome, and therefore, we start our analysis with the scalar casepresenting the basic ideas in a simpler form and reverting to the general multivariate case in Sect. 5. But in order to avoid repetitions the detailed proofs will be given for the matrix-valued functions.

The template problem: scalar interpolation
We consider the situation when all the interpolation nodes are at the origin. Assume that we are given the data {v 0 , v 1 , . . . , v K } and want to characterize all rational functions Q(s) = β(s) α(s) of a given degree (if any exists) such that where Q 1 (s) is a rational function analytic in 0. We do not presently require that Q(s) is proper: a method to achieve this condition is thoroughly examined in [18]. Note that the product α(s)Q 1 (s) should obviously be a polynomial. Thus a slightly different formulation of the problem is the following.
Problem 2 Find all the pairs (α, β) such that is a polynomial.
(These two formulations coincide if α and β are coprime). A second well-known instance is, given points λ 0 , . . . λ K ∈ C distinct, and values v 0 , . . . v K , to find all rational functions Q(s) of a given degree (if any exists) such that Again, if we assume α and β coprime, we can rewrite the problem as Problem 3 Find all the coprime pairs (α, β) such that is a polynomial for all i : Both problems can be written in the same matricial form: in case of (9) define: where A is a (K + 1) × (K + 1) matrix and U , V are (K + 1)-dimensional vectors. For (11), set A := diag{λ 1 , . . . , λ K }, U := [1, 1, . . . , 1] T and V as above.
Both problems are thus immediately seen to be special cases of the Problem 4 Find all the pairs (α, β) such that G(s) We will thus focus on the general case (which includes the above examples) of a given matrix A of dimension (K +1)×(K +1) and vectors U , V of dimension K +1.
We will assume that the pair (A, U ) is controllable and, to construct recursively the interpolants, that the matrix A is lower triangular.
Since the recursive algorithm will be formulated using the coefficients of the pair α, β we need an equivalent formulation of Problem 4. We introduce the following notation: for any polynomial γ we denote by γ the column vector of the coefficients-starting with the constant term-of the polynomial γ . (In some cases we will have to increase the dimension of the corresponding vector, by adding extra zeros as entries. This is in coherence with considering higher order terms in the polynomial but with zero coefficients.) In case of the polynomi- We set U j to be the matrix of dimension (K + 1) × ( j + 1) of the form Similarly, we set V r to be the matrix of dimension (K + 1) × (r + 1) of the form Notice that U = U 0 and V = V 0 .
The proof is deferred to the general case result in Theorem 2 if each of its columns are solutions to the Problem 4 and for any other solution (β, α), there exist polynomials π U , π V such that To construct a coprime factorization (3) from the data, we denote by μ and ν the controllability indexes of U and V relatively to . . in this order. Then μ is the smallest number for which A μ U can be expressed as a linear combination of its preceding vectors. The ν is defined similarly.) In view of controllability, holds with α V ν = 1.

Lemma 4 Let the polynomials α U (s) and β U (s) be defined as
where the coefficients are those in (17); let α V (s), β V (s) be defined similarly from (18). Then is a fundamental solution to Problem 4.
Proof Let us point out that that for the polynomials constructed above the inequalities hold. Thus this restriction gives us a fundamental solution of special type, although it does not imply uniqueness. For a general fundamental solution these inequalities do not necessarily hold. From Lemma 3 we obtain that the coefficients , define polynomials providing solutions to Problem 4, that is there exists a matrix polynomial Ψ (s) of size (K + 1) × 2 such that Eq. (3) holds. Since we can assume that β U (s) and α V (s) are monic, the matrix of the coefficients of the highest column degrees of Γ (s) is upper triangular, with identity on the diagonal. Since μ + ν = K + 1, the degree of det Γ (s) is K + 1, that is it coincides with that of det(s I − A). Using (3), the controllability of the pair (A, [U , V ]) gives in view of Ex. 6.3-19 in Kailath [12] that Ψ (s) and Γ (s) are right coprime.
In order to ensure uniqueness, introduce the following notion:

Proposition 3 For any controllable set of data (A, U , V ) there exists a unique MF solution.
Proof the statement is essentially trivial because Lemmas 3 and 4 show that the MF solution is nothing else than a basis selection scheme. Another way to see this if, for example μ > ν, is that the vector A μ U can be uniquely expressed in terms of the columns of , since the matrix has full rank (similarly for A ν V ). For sake of completeness, we provide a proof in terms of matrices for the multivariable case in Proposition (4).

Scalar recursive interpolation
The interest of MF solutions is that its regularity and uniqueness allow for a straightforward recursion algorithm which, as we shall see, can be easily generalized to the multivariable case.
We consider the situation where we have a sequence of nested problems indexed by k ≤ K , where A k is lower triangular (k + 1) × (k + 1) and U k and V k are (k + 1)-dimensional vectors. We say that the problems are nested if, for each k, (20) and the pair (A k , U k ) is controllable for k ≤ K . Our goal is to present a recursive algorithm to compute a minimal degree solution for each k from the solution for k − 1. As in (3), for each k, we can consider right coprime polynomial matrices Γ k (s) and Ψ k (s) such that Partitioning As in (14) and (15), we can define, for each k, j, r integers, U j k and V r k as Let be the MF solution to the Problem 4 determined by the data The following remark contains the core idea of the recursion: has k + 1 rows (of data), k + 3 columns and thus, in view of controllability, it has a two-dimensional kernel. The recursion adds one row to this matrix: controllability will imply that its kernel is always one dimensional. Since it must be contained in the previous kernel (in view of the triangularity of A), it can be expressed as a linear combination of . This is the first step of the algorithm. To find the appropriate linear combinations the "error terms" should be considered arising from the solutions obtained in the k k step but applying those for the next interpolation data. Let us note that the construction of a recursive scheme for simultaneous left and right interpolation in [14] using Loewner matrix pencils is also based on the use of similar error terms.
The second step consists of adding a column to this extended matrix, so that its kernel has again dimension two and the property of being a minimal fundamental solution is preserved, as well.
In view of the above remark, us denote by [u and compute the error terms arising from using the linear combinations obtained in the k th step for the next interpolation values u k+1 , v k+1 , as well: Let us point out that if U k+1 = 0 then μ k+1 = μ k while V k+1 = 0 implies that ν k+1 = ν k . In view of controllability, this implies that either The above procedure allows to identify an element in the kernel of [U In view of controllability, at least one of the errors is not zero and thus at least one of Eqs. (27) and (28) is satisfied, identifying the kernel. In case both errors U k+1 and V k+1 are nonzero, we might have a choice between these representations of the kernel. In Theorem 1 we give explicitly which one should be selected, but in order to complete the recursive step, we will need to extend either U m k+1 or V p k+1 to make the kernel two dimensional and find an extra generating vector in this kernel. (We consider j, r generic in the next lemma.) So, let U j k+1 and V r k+1 be partitioned as (25). Then and thus A similar representation holds if we want to extend V r k+1 to V r +1 k+1 . We discuss now how to handle these extension and, crucially for defining the algorithm, which one to pick.
k+1 be partitioned as (25) and let A k+1 be partitioned as in (20). Then, for any vector γ of dimension m, it is for v r +1 k+1 in the partition equivalent to (25) of V r +1 k+1 . Proof Notice that, in view of (29), On the other hand, from (14), Subtracting (33) from (32) yields the result. The corresponding equality for v r +1 k+1 , V r +1 k+1 is proved similarly.
Notice that the expressions in (30) and (31) do not depend on the last data points u k+1 and v k+1 , respectively, a fact which is crucial to derive Eqs. (34) and (35) below.
We denote by We are now ready to construct an interpolant of (A k+1 , U k+1 , V k+1 ). from a min- Proof Since Γ k (s) is a minimal fundamental solution, the degree of α U k (s) is smaller than ν k and thus the vector representation of the first column in the identity (36) is , while the vector representation of the second column in the identity reads as which, in view of (27), is again easily seen to be in the kernel of the matrix Furthermore, it is a fundamental interpolant: in fact, if Γ (s) is a minimal fundamental interpolant, there exists a polynomial matrix P(s) such that Γ (s)P(s) = Γ k+1 (s). Since the determinants of Γ k+1 (s) (by construction) and Γ (s) (by definition) are both χ A k+1 (s) (χ A (s) being the characteristic polynomial of A), the matrix P(s) is a unit. A similar reasoning holds for (37), using (28) and (35).
So, the recursive step appears to be quite simple. The problem is that we obtain a fundamental solution which is not minimal. Since the proof of Lemma 6 makes use of (34) and (35), we need the degree of α U k (s) to be less than ν k or that of β V k (s) to be less than μ k . There is a small modification of the above procedure which ensures that this condition is satisfied, making the recursion complete. Notice that if both U k+1 and V k+1 are different from zero then the cases μ k ≤ ν k and μ k > ν k should be analyzed separately.
Theorem 1 Let (A k , U k , V k ) be a sequence of nested problems for k ≤ K and let be as in (22) and such that (19) are satisfied and let be as in (26).
In both cases, Γ k+1 (s) satisfies (19), that is it determines the minimal fundamental solution.
The proof will be given in Theorem 3, where the multivariable case is treated. Let us point out that Eq.

A general minimal degree solution for the tangential interpolation problem
Let us return to the interpolation Problem 1 on matrix valued functions.  pair (A, W ), that is, the smallest exponent κ i such that Notice that this implies that the vectors in ] are linearly independent. On the other hand, in view of controllability, they span C K +1 and thus r i=1 κ i = K + 1. We denote by γ i, j (s) the (i, j)-th entry of Γ (s) and set, for where γ i, j denotes the vector formed by the coefficients we add zero coefficients of the higher powers if the actual degree of γ i, j (s) is lower than κ i ). That is . . , r be scalar valued polynomials.
Then condition (5) is equivalent to with the corresponding in Eq. (5). Moreover, if γ 1, j , . . . , γ r , j have a common factor s − s 0 that is Proof A straightforward computation gives that: Thus, setting φ j (s) as in (42), we immediately get that (41) Conversely, if for some matrix polynomial Φ(s), Eq. (5) holds with Γ (s), then for each column φ j (s), for j = 1, . . . , r , we have, substituting in (43) The term on the left-hand side is constant, while the degree of the right-hand side-if it is not identically zero-is at least 1. Consequently, identity (41) should hold.
Thus we have that: Since (A − s 0 I ) is invertible, we reach the desired conclusion.
Immediate consequence of the calculation in the previous Theorem is the following corollary.

Corollary 1 Under the conditions of Theorem 2 the following representation holds
where Δ is any simple closed curve with counterclockwise orientation around the eigenvalues of A.
Proof Introducing the notation Equation (43) can be written in the form Multiplying both sides by (s I − A) −1 and integrating, we obtain: Using that 1 2πi Δ (s I − A) −1 ds = I we get (44).
As we have pointed out after Eq.

Remark 6
Let us observe that an immediate consequence of the degree constraints in the MF solution that it will be column reduced. This property plays a crucial role in the analysis of the minimal possible McMillan degree of the interpolants in [3].
The following proposition shows that the MF solution is essentially equivalent to a basis selection scheme corresponding to the controllability indexes of the pair (A, W ) (Cf. Section 6.7.2 in [12]). This statement follows from the more general analysis presented in [2] which is based on the so-called nice selections defined by Youngdiagrams, but we present here a short proof of it which fits more to our purposes.

Proposition 4 For any controllable set of data (A, W ) there exists a unique MF solution.
Proof By construction, the set constitutes a basis for C K +1 . This basis is uniquely determined by the order of the columns of W . Thus, for j = 1, . . . , r , there exist unique vectorsγ 1, j , . . . ,γ r , j of dimensions κ 1 , . . . , κ r , such that Set now, for j = 1, . . . , r , Then (47) can be rewritten as: Therefore, in view of Theorem (2), for each j = 1, . . . , r , the vector-polynomial γ j (s) corresponding to γ 1, j . . . , γ r , j as in (40) provides a solution to Problem 1.
Furthermore, the polynomials γ 1, j (s), γ 2, j (s), . . . , γ r , j (s) associated with the vectors γ j will have degree at most κ i − 1 if i = j and κ i if i = j. On the other hand, by construction, for j = 1, . . . , r only the components of (46) with exponent l ≤ κ j are needed, which means that deg γ i, j (s) ≤ κ j . Thus, the following inequalities hold: and γ j, j (s) will be monic. Thus conditions (45) are satisfied. Moreover, for each column j the degree of its j-th entry is exactly κ j and it is not greater than κ j if i < j and strictly less than j if i > j. Thus, setting Γ (s) := [γ 1 (s), γ 2 (s), . . . , γ r (s)], its column highest degree coefficient matrix is a triangular matrix with the diagonal equal to the identity. Therefore deg det Γ (s) = j κ j = K + 1. In particular, det Γ (s) is not identically zero. (Moreover, Γ is column reduced.) To prove uniqueness, observe that, ifΓ (s) is another solution to (4) which satisfies (45), then the differenceΓ (s) := Γ (s) −Γ (s) is such that, for each j = 1, . . . r " degγ j, j (s) < κ j . Thus, in view of (45), it is In particular On the other hand,Γ (s) still satisfies (4) and thus each of its columnsγ j (s) can be expressed in terms of the columns of the fundamental solution Γ (s): thus there exists a polynomialπ j (s) such thatγ j (s) = Γ (s)π j (s). Since Γ (s) is column reduced, it has the predictable degree property (see [12,, that is, whereπ l, j (s) is the l-th row ofπ j (s) for l = 1, . . . r . Let us introduce the notation τ 1 = max j κ j . From Eqs. (50) and (51) it is immediate thatπ i, j = 0, if κ i = τ 1 .
Consequently, the following representation holds: Let us observe that rearranging the columns (and the rows accordingly) of Γ according to the decreasing order of the controllability indexes the highest column degree coefficients matrix is transformed into an upper triangular matrix with diagonal entries equal to 1; thus, any of its principal submatrix is of full rank, and thus, the corresponding part of the matrix Γ has the predictable degree property. Reduce now the column vectors ofΓ keeping only those entriesγ i, j for which κ i < τ 1 . Now, let us introduce the notation: τ 2 = max {κ i |κ i < τ 1 }. Since the submatrix of Γ formed by the entries γ i,l for which κ i ≤ τ 2 and κ l ≤ τ 2 also has the degree predictable property we have that max i:κ i ≤τ 2 degγ i, j = max l:κ l ≤τ 2 ,π l, j =0 (deg γ l + degπ l, j ) .
The argument can be repeated for the matrix obtained fromΓ after eliminating those entriesγ i, j for which κ i < τ 2 , and so on, finally proving thatπ l, j (s) = 0 for l, j = 1, . . . , r , thusΓ (s) = 0, achieving the proof.

Remark 7
The uniqueness of the solution Γ (s) hinges on both of the conditions in (45) and thus it relies on the controllability indexes of (A, W ). If Γ (s) is another column-reduced matrix polynomial satisfying only the second condition in (45) then using Theorem 6.5-4. in [12] we see that Γ (s) can be obtained from Γ (s) multiplying it from the right by a unimodular matrix; consequently, Lemma 6.3-14. in [12] gives us that the column degrees of Γ (s) and Γ (s)-after an appropriate permutation of the columns-coincide.

General multivariable recursive interpolation
Similarly to the scalar case, given a controllable pair (A, W ), with A lower triangular of dimension (K +1)×(K +1) and W := [W 1 , W 2 , . . . , W r ] of dimension (K +1)×r , we can consider a sequence of reduced interpolation problems (A k , W k ) where A k is the submatrix of the first k + 1 rows and columns of A and W k = [W 1,k , . . . , W r ,k ] the truncation of W to its first k + 1 rows. According to Lemma 2 this means in particular to find right coprime polynomial matrices Φ k (s), Γ k (s) such that  (Let us point out that we assume (A k , W k ) controllable for each k. Now since A is lower triangular, this can be achieved if for example the first entry of a column of W and at least one element of each row of A before the diagonal are not zero. Another possibility is given when the diagonal elements of A are different and nonzero, and in each row of W there is a nonzero element.) Comparing the previous inequalities to the identities r j=1 κ j k = k + 1, r j=1 κ j k+1 = k + 2 we obtain that there is only one controllability index which changes during the step k → k + 1. To locate this let us introduce the following notation: For k, j fixed and i = 1, . . . , r let us denote by γ k i, j a column vector of dimension κ i k + 1 defined by the coefficients of the polynomial γ k i, j (s) (which is the (i, j) th entry of Γ k (s))-extended by extra zeros if necessary-and we set Notice that (see (48)). Furthermore, let w κ k k+1 := [w κ 1 k+1 , w κ 2 k+1 , . . . , w κ r k+1 ] denote the last row of the extended matrix ]. Let us observe that due to the assumption that A is lower triangular we have that and thus, similarly to the scalar case, ⎤ ⎥ ⎦ and eventually, for j = 1, . . . , r , Finally, set Let us observe that if j k+1 = 0 for some j, then κ j k+1 = κ j k . (The same linear combination which "worked" for the first k coordinates also gives zero at the last position.) Notice as well that, since W κ k k+1 has k + r + 1 columns and k + 2 rows, its kernel has dimension r − 1, in view of controllability. Thus k+1 = 0, since is of full rank. The following lemma shows how to locate which controllability index changes during the step k → k + 1.
we obtain in view of (56), (57) and (59). We claim that the controllability indexes of In view of the ordering, it is i > j and, therefore, the condition is satisfied also forγ i k , proving the claim.
Since, in view of controllability, the sum of the controllability indexes has to increase by 1, the only possibility is κ j k+1 = κ j k + 1. Thus, the controllability indexes of , . . . , κ r k . In the proof of the recursive scheme for the scalar case Lemma 6 played an important role which was based onto Lemma 5. The multivariate version of this latter one is formulated below.
j,k+1 be partitioned as (57) and let A k+1 be partitioned as in (55). Then, for any vector γ of dimension κ j k + 1, we have that The proof is similar to that of the scalar case in Lemma 5 and it will not be repeated here.
Again-using the observation that the last entry of the vector γ k i, j is zero for i = j, an immediate consequence of Lemma 8 is that, if γ k j ∈ ker W κ k k,i , then ⎡ where e j denotes the j th r -dimensional unit vector.
Proof Due to the definition of the index j to find the matrix polynomial Γ k+1 (s) satisfying (45) we have to consider special elements in the kernel of [W and, for i = j, and 2. Recursion step: Repeat while k < K (where K is the number of data): -Given W κ k k , γ k j for j = 1, . . . , r , compute W κ k k+1 using (57) and k+1 using (59). Set i to be the index corresponding to the smallest κ i k such that i k+1 in (59) is different from 0. Define, for j = 1, . . . , r , γ k+1 j using (66) and (67). Define W κ k+1 k+1 using (57) and (58). -Increase k.

Recursive interpolation of a given function-scalar case
In the remaining part of this paper we consider a very special situation, when all the interpolation nodes are at the origin. We are going to show that this problem exhibit several interesting properties not present in the general case.
We suppose that the data v 0 , . . . , v k we consider are derived from the power series expansion of a given rational function where we assume that α and β are coprime, α is monic and that α(0) = 0. We denote in this section by p and m the degrees of α and β, respectively. We would like to assess now how the controllability indexes μ k and ν k -giving the column degrees of minimal fundamental solutions-increase as the number of data v 0 , . . . , v k under consideration increases.
In this special scalar interpolation problem the matrices of the sequence of nested problems indexed by k, are given as where A k is lower triangular (k + 1) × (k + 1) and U k and V k are (k + 1)-dimensional vectors.
Due to the special structure of the matrix A the equations characterizing the polynomials β U , α U , β V , α V can be written in a simpler form. (Since for the time being we analyze the properties for a fixed k, we omit the index from the notation.) Thus the Similarly, the matrix Notice that U = U 0 and V = V 0 .
Denoting the controllability indexes of U and V by μ and ν, respectively, we see that Eq. (16) characterizing the coefficient vectors of the polynomials in the minimal fundamental solution can be written as finding the kernel of the matrix [U μ , −V ν ], which for μ ≥ ν can be written as follows: Let us point out that to determine μ the second block column should be written as a linear combination of the columns in the first and third block, while to find ν the last column should be expressed as a linear combination of all the previous ones. But the vectors in the first block are the unit vectors, thus for μ the part below the first horizontal line should be analyzed, for ν that one below the second horizontal line.
Using the special structure of (74) the behavior of the controllability indexes as k → ∞ is described by the following Lemma.

Lemma 9
Given the rational function f = β α define the A k , U k , V k matrices as in (71).
Denote by μ k , ν k the controllability indexes of (A k , Proof First let us assume that deg α ≥ deg β. That is, f is a proper rational function. In order to find the controllability index ν k the last column in (74) should be expressed as a linear combination of the "previous" ones. But the columns in the first block part of (74) contains unit vectors this means that we should look after linear combination between the vectors in the third part-below the second horizontal line.
Now, let us denote by p the degree of α. Then computing from the power series expansion of f = β α the coefficients of s k , k ≥ p + 1, in f (s)α(s) we obtain that Since α p = 0 the coefficient of the polynomial provides a linear combination we are looking for. Thus, we get if k is large enough then ν k ≤ p. Consequently, as k → ∞ the sequence ν k remains bounded, so μ k → ∞. At the same time, when we already have that μ k > ν k then columns in the third block in the matrix (74) are linearly independent, thus the matrix is of full column-rank. Consequently, in this case α V k = α and ν k = deg α. Next consider the case when p = deg α < deg β. In order to find the controllability index μ k the second column in (74) should be expressed as a linear combination of the "previous" ones. But the columns in the first block part of (74) contains unit vectors this means that we should look after linear combination expressing the vector in the second column using vectors from the third block column-but now below the first horizontal line.
Let us denote by m the degree of β. Expressing again the coefficients of s k in f (s)α(s) but now for k ≥ m we obtain that If ν k > μ k , then the columns of the matrix are linearly independent and thus in this case α U k = 1 β p α and μ k = deg β. This concludes the proof of the Lemma.
In the construction of the recursive scheme for interpolation the so-called error terms U , V played an important role.
Let U μ k k+1 , −V ν k k+1 be as in (25) and define the error terms arising from the linear combinations obtained in the k th step but including the next interpolation value v k+1 , as well In general, given two polynomials β k (s), α k (s) of degree μ k , ν k , respectively, and the corresponding vectors β k , α k , we can set k+1 := u The last row of the vector on the left hand side is k+1 ; thus, using that the last row of (s Proof From (78), we have: In the general recursion scheme presented in Theorem 1 the minimal fundamental solution Γ k+1 was defined using Γ k . In the present situation-as we are going to seeit will be natural to consider larger steps. To this end, given k 1 < k 2 and the associated minimal fundamental solutions Γ k 1 (s) and Γ k 2 (s), setΓ k 1 ,k 2 (s) := Γ −1 k 1 (s)Γ k 2 (s) (it clearly is a polynomial matrix). If k 2 = k 1 +1 then we write shortlyΓ k 2 (s) =Γ k 1 ,k 2 (s).
What is surprising is that this is true whenever one of the errors is 0, even as we go astray from the pattern 1, as the next result shows.

Theorem 4 Consider the recursive version of Problem 2 determined by the coefficients in the Taylor-series expansion of the rational function f (s) considered around the origin.
(i) Suppose that for some k we have that μ k = ν k + 1 and V k+1 = · · · = V k+l = 0 and V k+l+1 = 0. Let k := k + 2l: then μ k = ν k + 1 and k is the next smallest integer with this property. Furthermore, the matrixΓ k,k relating Γ k and Γ k from repeated iterations of (38) and (39) is lower triangular. That is, where l = ν k − ν k and p l is suitable polynomial of degree less than l. (ii) Similarly, suppose μ k = ν k and U k+1 = · · · = U k+l = 0 and U k+l+1 = 0. Let k := k + 2l: then μ k = ν k and k is the next smallest integer with this property. Moreover, the matrixΓ k,k relating Γ k and Γ k from repeated iterations of (38) and (39) is upper triangular. That is, where l = μ k − μ k , where p l is a suitable polynomial of degree no greater than l.
Proof Since μ k = ν k + 1 > ν k , but V k+1 = · · · = V k+l = 0, in accordance with Theorem 1, we need, in order to obtain Γ k,k , to use first (38) l times. It is immediate to verify that each of these factors is lower triangular. Thus also the matrix relating Γ k and Γ k+l is lower triangular. Furthermore, during these steps the second controllability index does not change, that is, ν k+l = ν k , while the first one in each step increases by one: μ k+l = μ k + l.
We claim now that, to obtain Γ k , we need l iterations of (39) which are also lower triangular. Since V k+l+1 = 0, for the step k + l + 1 we have to use (39). Now the column degrees of Γ k+l are ν k + l + 1, ν k . Therefore, deg β U k+l − deg β V k+l ≥ l + 1 and thus the coefficient β V k+l,μ k+l −1 = 0; consequently, also the last factor in (39) is lower triangular for this step. Thus, β V k+l+1 (s) = sβ V (s) k+l and α V k+l+1 (s) = sα V k+l (s) In view of Corollary 2, V k+l+2 = V k+l+1 . Thus it is different from zero. The same reasoning can be repeated whileΓ k+l+ j+1 (s) is lower triangular, that is as long as deg β U k+l+ j − deg β V k+l+ j > 1. This will happen, if μ k+l+ j > ν k+l+ j + 1. But during these steps the sequence of μ-s does not change. Thus μ k+l+ j = ν k + l + 1 , ν k+l+ j = ν k + j . Γ k+l+ j+1 (s) still will be lower triangular for j = 0, . . . , l − 1. In conclusion,Γ k,k (s) is lower triangular. Since it has determinant 2l, and each column has degree at least l, it can only have the form described in (79).
Let us point out that using again Corollary 2 for V k we obtain that V k +1 = 0 and thus μ k +1 = μ k and ν k +1 = ν k + 1. But now the factorΓ k +1 is not necessarily lower triangular.
The second statement can be proven similarly. But in this case the final application of Corollary 2 gives that μ k +1 = μ k + 1 and ν k +1 = ν k again pointing out that now the polynomial matrixΓ k +1 (s) is not necessarily upper triangular.
In conclusion, the possible paths the controllability indexes can exhibit are of the forms: The construction and the proof of Theorem 4 is based on Theorem 1 determining a recursion for the minimal fundamental solutions of Problem 2. The next theorem is based on a recursion for an appropriate sequence of fundamental (not necessarily minimal) solutions. This construction gives the possibility of explicitly determining the positions in the jumps of the degrees in the minimal solutions of the interpolation problem using the form of the given rational function f . The basic idea is that if for example β V , α V are already given and provide a solution for Problem 2 for v 0 , v 1 , . . . , v K , then to find the corresponding β U , α U which complete the pair to a minimal solution, a pair of polynomials should be determined in such a way that the determinant of the corresponding polynomial matrix formed from these polynomial vectors is equal to s K +1 .
In order to show this let us assume that f = β α , where deg α = p, deg β = m and α(0) = 0, and generate a sequence of polynomials using the following recursion: First let us denote by h 1 ≥ 0 the multiplicity of 0 as a possible zero of β and introduce the notation: β 1 (s)s h 1 = β(s). (That is β 1 (0) = 0.) Consider the solutions of the following equation where γ 1 and δ 1 are polynomials, deg γ 1 ≤ h 1 , while δ 1 (0) = 0. Notice that the coefficients of 1, s, . . . , s h 1 in s h 1 δ 1 are zero. Since according to the construction the constant term in β 1 is nonzero, this equation has a unique solution. Furthermore, Let us continue similarly-as in the Euclidean-algorithm-denoting by h 2 ≥ 1 the multiplicity of 0 as the zero of δ 1 and introducing the polynomial β 2 as solve the equation where γ 2 , δ 2 are polynomials, deg γ 2 ≤ h 2 and δ 2 (0) = 0. In general, if where δ j (0) = 0, then denote by h j+1 the multiplicity of 0 as the zero of δ j and set Then γ j+1 and δ j+1 are defined as solutions of where deg γ j+1 ≤ h j+1 and δ j+1 (0) = 0. We obtain that Consequently, the sequence deg β j , j ≥ 2 is strictly decreasing, so this algorithm terminates in finite steps. That is, for some r we have that Straightforward computation gives that . Now let us define a sequence a polynomial pairs as follows: for j = 1, . . . , r . We obtain that and In particular, Let us introduce the notation For the sake of simplicity we are going to assume that h 1 ≥ 1. In this case Let us point out that if h 1 = 0 then h 2 = min{h ≥ 1 | v h = 0}, thus the sequence of l 1 , l 2 , . . . essentially remains the same, except the 0 as a first element is added to it.
Theorem 5 Assume that f (s) = β(s) α(s) , where α(0) = 0 while β(0) = 0, with the Taylor-series expansion around the origin Apply the construction described above to obtain the polynomials N j , R j for j = 1, . . . , r . Then a fundamental solution for K = 2l j − 1 is given by Furthermore, both controllability indexes μ k , ν K are equal to l j , and coincide with the column degrees of this matrix.
Proof Let us note that the assumption β(0) = 0 implies that h 1 ≥ 1 and representation (82) implies that deg R j ≤ l j , deg N j ≤ l j First we are going to show by induction according to the number of products in the representation (82) that the polynomials N j , R j are coprime, and R j (0) = 0. Let us observe that since according to the construction β j−1 (0) = 0, β j (0) = 0 while δ j (0) = 0 we have from Eq. (81) that γ j (0) = 0. Thus s h j and γ j are coprime. Now let us introduce the notation for j = 2, . . . , r . Let us observe that N j (0) = 0. Then N j and R j can be expressed as where N (2) j and R (2) j are coprime using the induction hypothesis, while γ j (0) = 0 implying that N j and R j are coprime, as well. At the same time, since N (2) j (0) = 0 we obtain that R j (0) = 0 using the corresponding induction hypothesis on R (2) j . The representation (83) gives that Consequently, implying the N j R j is a solution of the Interpolation Problem 2 determined by the values v 0 = 0, v 1 , . . . , v l j +l j+1 −1 . In particular, the multiplicity of the origin as a zero of N j (s) is exactly h 1 .
Notice that the pair s h j N j−1 (s), s h j R j−1 (s) provides a solution for the Problem 2 determined by the values v 0 = 0, v 1 , . . . , v 2l j −1 .
Choosing K = 2 l j − 1 we obtain that both columns of the matrix where the matrix standing in the brackets is unimodular, thus its inverse is a polynomial matrix. Consequently, the matrix (84) is a fundamental solution, as well, concluding the proof of the first part. Furthermore, since the column degrees of the matrix (83) are no greater than l j while its determinant is s 2 l j , its highest column degree coefficient matrix should be of full rank and both column degrees should equal to l j . Using Lemma 6.3-14 in Kailath [12] we obtain that its column degrees and those of Γ K coincide. Thus l j = ν K and μ K = K + 1 − ν K = ν K .

Remark 8
In order to get the minimal fundamental solution Γ K from the matrix (84) the degrees of the (1, 2) and (2, 1) elements should be reduced and the (1, 1) and (2, 2) elements should be transformed to monic polynomials. Denoting by g j and f j the coefficients of s l j in R j (s) and N j (s), respectively, Eq. (86) implies that det − f j−1 f j −g j−1 g j = 1 . Consequently, where K = 2 l j − 1.
Now, combining Theorems 4 and 5 we can describe the paths of the controllability indexes generated by the minimal fundamental solutions to the restricted interpolation problems generated by the Taylor-series of a given function f . In fact, using the notations of Theorem 5 the following theorem holds. Theorem 6 Let introduce the notation k j = 2l j − 1. Then the set coincides with the set containing those values when for the recursive version of Problem 2 determined by the Taylor-series coefficients of the function f = β α the controllability indexes coincide. In particular, μ k j = ν k j = l j .
Furthermore, (i) if the coefficient g j of the term s l j in R j is zero, then μ k j +1 = μ k j , ν k j +1 = ν k j + 1 and the path of the pair of controllability indexes evolves according to part (i) of Theorem 4.
If moreover (a) h j+1 = 1 then in the next step both controllability indexes are equal: (This is the generic case.) (b) If g j = 0 and h j+1 > 1 then the path determined by the controllability indexes evolves according to part (ii) of Theorem 4.
Proof Theorem 5 implies that μ k j = ν k j = l j . Analyzing the path behavior of the controllability indexes we show that the next occasion when they are equal again is given by l j+1 , moreover at the same time we describe the form of the path between these two points. From Theorem 5 we have obtained that β U k j (s) = −g j s h j N j−1 (s) + g j−1 N j (s) α U k j (s) = −g j s h j R j−1 (s) + g j−1 R j (s) .
It was already pointed out that the pair N j , R j provide solutions to the Problem 2 exactly up to l j + l j+1 − 1. Thus, if the next interpolation value v l j +l j+1 is taken into consideration there will be a nonzero error term obtained via using this pair of polynomials. Let us denote this error by j . That is, Since h j+1 ≥ 1 and thus l j + l j+1 − 1 ≥ 2l j so the pair N j , R j still interpolates for the next value v 2l j , using Corollary 2 we obtain that U k j +1 = −g j−1 j−1 .

Similar computation gives that
V k j +1 = f j j−1 .
Consequently, U k j +1 = 0 if and only if g j = 0. Theorem 1 implies that in this case μ k j +1 = μ k j , and ν k j +1 = ν k j + 1. But if g j = 0 then the pair β U k j , α U k j are obtained from N j , R j via multiplication with a nonzero constant g j−1 , thus part ii) of Theorem 1 gives that in the present case β U k j +1 = β U k j , α U k j +1 = α U k j providing solutions to the interpolation problem determined by v 0 , v 1 , . . . , v k up to k ≤ 2 l j + h j+1 − 1. Now concluding the proof of the theorem.

Funding Open Access funding provided by Eötvös Loránd University
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Applying again Theorem 1 we get that Therefore, we cannot increase the index ν 5 . We thus have to increase μ 5 . While the Fibonacci sequence continues, this yields, for k ≥ 4, Γ k+1 (s) = Γ 4 (s) s k−3 0 0 1 . If for some k 0 the value v k 0 +1 is not a Fibonacci number (while all the previous are), then V k 0 +1 = 0 and we can increase ν k 0 by multiplying the second column of Γ k 0 (s) by s. A minimal interpolant will thus now have the degree of the first column of Γ k 0 +1 , which is k 0 .